text:
#=
Optimized search functions for UTF16Str

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

function _srch_cp(::Fwd, ::CodeUnitMulti, str::Str{UTF16CSE}, cp::AbsChar, pos, len)
    @preserve str begin
        beg = pointer(str)
        (ch = codepoint(cp)) <= 0x0ffff && return _srch_codeunit(Fwd(), beg, ch%UInt16, pos, len)
        pnt = bytoff(beg, pos)
        fin = bytoff(beg, len)
        lead, trail = get_utf16(ch)
        while pnt < fin
            # This checks for the trailing surrogate, since they are much less like to match
            # (there are only a few leading surrogates for the assigned areas in the non-BMP
            # planes)
            get_codeunit(pnt) == trail && get_codeunit(pnt - 2) == lead &&
                return chrdiff(pnt, beg)
            pnt += sizeof(UInt16)
        end
        0
    end
end

function _srch_cp(::Rev, ::CodeUnitMulti, str::Str{UTF16CSE}, cp::AbsChar, pos, len)
    @preserve str begin
        beg = pointer(str)
        (ch = codepoint(cp)) <= 0x0ffff && _srch_codeunit(Rev(), beg, ch%UInt16, pos)
        pnt = bytoff(beg, pos)
        lead, trail = get_utf16(ch)
        beg += sizeof(UInt16) # Need to have room for lead
        while pnt > beg
            get_codeunit(pnt -= 2) == trail && get_codeunit(pnt -= 2) == lead &&
                return chrdiff(pnt, beg) + 1
        end
        0
    end
end


text:
module StatsBase

import Base: length, size, isempty, values, sum, show, maximum, minimum, extrema
import Base.Cartesian: @nloops, @nref, @nextract
using Base: @irrational, @propagate_inbounds
using DataAPI
import DataAPI: describe
import DataStructures: heapify!, heappop!, percolate_down!
using SortingAlgorithms
using Missings

using Statistics
using LinearAlgebra
using Random
using Printf
using SparseArrays
import Random: rand, rand!
import LinearAlgebra: BlasReal, BlasFloat
import Statistics: mean, mean!, var, varm, varm!, std, stdm, cov, covm,
                   cor, corm, cov2cor!, unscaled_covzm, quantile, sqrt!,
                   median, middle
import StatsAPI: pairwise, pairwise!

    ## tackle compatibility issues

export

    ## weights
    AbstractWeights,    # abstract type to represent any weight vector
    Weights,            # to represent a generic weight vector
    AnalyticWeights,    # to represent an analytic/precision/reliability weight vector
    FrequencyWeights,   # to representing a frequency/case/repeat weight vector
    ProbabilityWeights, # to representing a probability/sampling weight vector
    UnitWeights,        # to representing a uniform weight vector
    weights,            # construct a generic Weights vector
    aweights,           # construct an AnalyticWeights vector
    fweights,           # construct a FrequencyWeights vector
    pweights,           # construct a ProbabilityWeights vector
    eweights,           # construct an exponential Weights vector
    uweights,           # construct an UnitWeights vector
    wsum,               # weighted sum with vector as second argument
    wsum!,              # weighted sum across dimensions with provided storage

    ## moments
    skewness,       # (standardized) skewness
    kurtosis,       # (excessive) kurtosis
    moment,         # central moment of given order
    mean_and_var,   # (mean, var)
    mean_and_std,   # (mean, std)
    mean_and_cov,   # (mean, cov)

    ## scalarstats
    geomean,     # geometric mean
    harmmean,    # harmonic mean
    genmean,     # generalized/power mean
    middle,      # the mean of two real numbers
    mode,        # find a mode from data (the first one)
    modes,       # find all modes from data

    zscore,      # compute Z-scores
    zscore!,     # compute Z-scores inplace or to a pre-allocated array

    percentile,  # quantile using percentage (instead of fraction) as argument
    nquantile,   # quantiles at [0:n]/n

    span,        # The range minimum(x):maximum(x)
    variation,   # ratio of standard deviation to mean
    sem,         # standard error of the mean, i.e. sqrt(var / n)
    mad,         # median absolute deviation
    iqr,         # interquatile range

    genvar,      # generalized variance
    totalvar,    # total variation

    entropy,        # the entropy of a probability vector
    renyientropy,   # the Rényi (generalised) entropy of a probability vector
    crossentropy,   # cross entropy between two probability vectors
    kldivergence,   # K-L divergence between two probability vectors

    summarystats,   # summary statistics
    describe,       # print the summary statistics

    # deviation
    counteq,        # count the number of equal pairs
    countne,        # count the number of non-equal pairs
    sqL2dist,       # squared L2 distance between two arrays
    L2dist,         # L2 distance between two arrays
    L1dist,         # L1 distance between two arrays
    Linfdist,       # L-inf distance between two arrays
    gkldiv,         # (Generalized) Kullback-Leibler divergence between two vectors
    meanad,         # mean absolute deviation
    maxad,          # maximum absolute deviation
    msd,            # mean squared deviation
    rmsd,           # root mean squared deviation
    psnr,           # peak signal-to-noise ratio (in dB)

    # cov
    scattermat,     # scatter matrix (i.e. unnormalized covariance)
    cov2cor,        # converts a covariance matrix to a correlation matrix
    cor2cov,        # converts a correlation matrix to a covariance matrix
    CovarianceEstimator, # abstract type for covariance estimators
    SimpleCovariance,    # simple covariance estimator

    ## counts
    addcounts!,     # add counts to an accumulating array or map
    counts,         # count integer values in given arrays
    proportions,    # proportions of integer values in given arrays
                    # (normalized version of counts)
    countmap,       # count distinct values and return a map
    proportionmap,  # proportions of distinct values returned as a map

    ## ranking
    ordinalrank,    # ordinal ranking ("1234" ranking)
    competerank,    # competition ranking ("1 2 2 4" ranking)
    denserank,      # dense ranking ("1 2 2 3" ranking)
    tiedrank,       # tied ranking ("1 2.5 2.5 4" ranking)

    ## rankcorr
    corspearman,       # spearman's rank correlation
    corkendall,        # kendall's rank correlation

    ## partialcor
    partialcor,        # partial correlation

    ## signalcorr
    autocov!, autocov,      # auto covariance
    autocor!, autocor,      # auto correlation
    crosscov!, crosscov,    # cross covariance
    crosscor!, crosscor,    # cross correlation
    pacf!, pacf,            # partial auto-correlation

    ## sampling
    samplepair,     # draw a pair of distinct elements
    sample,         # sampling from a population
    sample!,        # sampling from a population, with pre-allocated output
    wsample,        # sampling from a population with weights
    wsample!,       # weighted sampling, with pre-allocated output

    ## empirical
    ecdf,           # empirical cumulative distribution function
    ECDF,           # type for empirical cumulative distribution function

    AbstractHistogram,
    Histogram,
    midpoints,
    # histrange,

    ## robust
    trim,           # trimmed set
    trim!,          # trimmed set
    winsor,         # Winsorized set
    winsor!,        # Winsorized set
    trimvar,        # variance of the mean of a trimmed set

    ## misc
    rle,            # run-length encoding
    inverse_rle,    # inverse run-length encoding
    indexmap,       # construct a map from element to index
    levelsmap,      # construct a map from n unique elements to [1, ..., n]
    indicatormat,   # construct indicator matrix
    pairwise,       # pairwise application of functions
    pairwise!,      # pairwise! application of functions

    # statistical models
    CoefTable,
    StatisticalModel,
    RegressionModel,

    adjr2,
    adjr²,
    aic,
    aicc,
    bic,
    coef,
    coefnames,
    coeftable,
    confint,
    cooksdistance,
    crossmodelmatrix,
    deviance,
    dof,
    dof_residual,
    fit,
    fit!,
    fitted,
    informationmatrix,
    isfitted,
    islinear,
    leverage,
    loglikelihood,
    meanresponse,
    modelmatrix,
    mss,
    response,
    responsename,
    nobs,
    nulldeviance,
    nullloglikelihood,
    rss,
    score,
    stderror,
    vcov,
    predict,
    predict!,
    residuals,
    r2,
    r²,

    ConvergenceException,

    # data standardization
    standardize,
    AbstractDataTransform, # the type to represent a abstract data transformation
    ZScoreTransform,       # the type to represent a z-score data transformation
    UnitRangeTransform     # the type to represent a 0-1 data transformation

# source files

include("common.jl")
include("weights.jl")
include("moments.jl")
include("scalarstats.jl")
include("robust.jl")
include("deviation.jl")
include("cov.jl")
include("counts.jl")
include("ranking.jl")
include("toeplitzsolvers.jl")
include("rankcorr.jl")
include("signalcorr.jl")
include("partialcor.jl")
include("empirical.jl")
include("hist.jl")
include("pairwise.jl")
include("misc.jl")

include("sampling.jl")
include("statmodels.jl")

include("transformations.jl")

include("deprecates.jl")

end # module

ext: jl


text:
import Distributions
import LinearAlgebra

const VISUAL_RANDOM_STRUCTURE = true

@testset "$(format("Random Structure: FixedRankMatrix[Constructor]"))" begin
    # test data type
    @test isa(FixedRankMatrix{Distributions.Poisson},DataType) == true
    @test isa(FixedRankMatrix{Distributions.Poisson}(),DataType) == false
    # test default constructor
    @test FixedRankMatrix{Distributions.Poisson}(Distributions.Poisson();rank=4).rank==4
    # test shorthanded constructor
    @test FixedRankMatrix(Distributions.Poisson(5),rank=2).rank == 2
    @test typeof(FixedRankMatrix(Distributions.Poisson(5),rank=2).dist) == Distributions.Poisson{Float64}
    # test mixed distribution

end

@testset "$(format("Random Structure: FixedRankMatrix[overload: Base.rand]"))" begin
    let
        tc1 = rand(FixedRankMatrix(Distributions.Gaussian(0,1),rank =3),
                   10,10)
        @test_logs (:warn,"Rank is not specified. Using formula: rank= ⌊0.3 * (row ∧ col)⌋") rand(FixedRankMatrix(Distributions.Gaussian(0,1)), 100,30)
        tc2 = rand(FixedRankMatrix(Distributions.Gaussian(0,1)),
                   100,30)
        @test check(:rank,tc1,3)
        @test check(:rank,tc2,9)

        tc3 = rand([(FixedRankMatrix(Distributions.Gaussian(0,1),  rank=2),10,3),
                    (FixedRankMatrix(Distributions.Poisson(5),     rank=2),10,3),
                    (FixedRankMatrix(Distributions.Bernoulli(0.5), rank=2),10,4)])
        @test LinearAlgebra.rank(tc3) == 6
    end
end

@testset "$(format("Random Structure: FixedRankMatrix[overload: Concepts.provide]"))" begin
    # test provide interface
    let
        tc = provide(FixedRankMatrix(Distributions.Gaussian(0,1),rank=5),
                     row = 10,col=10)
        @test check(:dimension,tc,(10,10))
        @test check(:rank,tc,5)
    end
end

@testset "$(format("Random Structure: GaussianMatrix"))" begin
    let
        tc = GaussianMatrix(20,20,rank=10,μ=0.0,σ=1.0)
        @test check(:dimension,tc,(20,20))
        @test check(:rank,tc,10)
        # test optional rank parameter
        tc1 = GaussianMatrix(20,20;μ=0,σ=1)
        @test check(:dimension,tc1,(20,20))
        @test check(:rank,tc1,20)
    end
end

@testset "$(format("Random Structure: PoissonMatrix"))" begin
    let
        tc = PoissonMatrix(20,20,rank=10,λ=3)
        @test check(:dimension,tc,(20,20))
        @test check(:rank,tc,10)
        # test optional rank parameter
        tc1 = PoissonMatrix(20,20;λ=3)
        @test check(:dimension,tc1,(20,20))
        @test check(:rank,tc1,20)
    end
end

@testset "$(format("Random Structure: BernoulliMatrix"))" begin
    let
        tc = BernoulliMatrix(20,20,rank=10,p=0.5)
        @test check(:dimension,tc,(20,20))
        @test check(:rank,tc,10)
        # test optional rank parameter
        tc1 = BernoulliMatrix(20,20;p=0.5)
        @test check(:dimension,tc1,(20,20))
        @test check(:rank,tc1,20)
    end
end


text:
using TokenizeMeta
using Test

@testset "TokenizeMeta.jl" begin

    t1 = collect(tokenize("Int64"))[1]
    @test Tokenize.Tokens.tevalfast(t1) == Int64
    @test Tokenize.Tokens.tgetfield(t1) == Int64
    @test Tokenize.Tokens.teval(t1) == Int64
    @test Tokenize.Tokens.ttypeof(t1) == DataType
    @test Tokenize.Tokens.tisa(t1, DataType) == true

    t2 = collect(tokenize("[5]"))
    @test Tokenize.Tokens.tevalfast(t2) == [5]
    @test Tokenize.Tokens.teval(t2) == [5]
    @test Tokenize.Tokens.ttypeof(t2) == Array{Int64,1}
    @test Tokenize.Tokens.tisa(t2, Array{Int64}) == true

end


text:
# Core functions for evaluating PDDL formulae and state transitions

"""
    satisfy(formulae, state, domain=nothing; mode=:any)

Returns whether `formulae` can be satisfied by facts in the given `state`,
and any axioms / derived predicates in `domain`. If `formulae` contains
unbound variables, `mode=:any` returns the first substitution found,
while `mode=:all` returns all substitutions.
"""
function satisfy(formulae::Vector{<:Term}, state::State,
                 domain::Union{Domain,Nothing}=nothing; mode::Symbol=:any)
    # Do quick check as to whether formulae are in the set of facts
    function in_facts(f::Term)
        if !isempty(state.fluents)
            f = eval_term(f, Subst(), state.fluents) end
        if f in state.facts || f in state.types
            return true end
        if !isnothing(domain) && f in get_const_facts(domain)
            return true end
        if f.name in Julog.comp_ops || f.name in keys(state.fluents)
            return eval_term(f, Subst(), state.fluents).name == true end
        return false
    end
    if all(f -> f.name == :not ? !in_facts(f.args[1]) : in_facts(f), formulae)
        return true, [Subst()] end
    # Initialize Julog knowledge base
    clauses = isnothing(domain) ? Clause[] : get_clauses(domain)
    clauses = Clause[clauses; collect(state.types); collect(state.facts)]
    # Pass in fluents and function definitions as a dictionary of functions
    funcs = isnothing(domain) ?
        state.fluents : merge(state.fluents, domain.funcdefs)
    return resolve(formulae, clauses; funcs=funcs, mode=mode)
end

satisfy(formula::Term, state::State, domain::Union{Domain,Nothing}=nothing;
        options...) = satisfy(Term[formula], state, domain; options...)

"""
    evaluate(formula, state, domain=nothing; as_const=true)

Evaluates `formula` as fully as possible with respect to the fluents defined
in `state`, along with any axioms / derived predicates defined in `domain`.
Returns a Julog `Const` if `as_const=true`, otherwise return an unwrapped
Julia value.
"""
function evaluate(formula::Term, state::State,
                  domain::Union{Domain,Nothing}=nothing; as_const::Bool=true)
    # Evaluate formula as fully as possible
    funcs = isnothing(domain) ?
        state.fluents : merge(state.fluents, domain.funcdefs)
    val = eval_term(formula, Subst(), funcs)
    # Return if formula evaluates to a Const (unwrapping if as_const=false)
    if isa(val, Const) && !isa(val.name, Symbol)
        return as_const ? val : val.name end
    # If val is not a Const, check if holds true in the state
    sat, _ = satisfy(Term[val], state, domain)
    return as_const ? Const(sat) : sat # Wrap in Const if as_const=true
end

"""
    find_matches(formula, state, domain=nothing)

Returns a list of all matching substitutions of `formula` with respect to
a given `state` and `domain`.
"""
function find_matches(formula::Term, state::State,
                      domain::Union{Domain,Nothing}=nothing)
    if formula.name in keys(state.fluents)
        clauses = Vector{Clause}(get_fluents(state))
        _, subst = resolve(formula, clauses; mode=:all)
    else
        clauses = isnothing(domain) ? Clause[] : get_clauses(domain)
        clauses = Clause[clauses; collect(state.types); collect(state.facts)]
        funcs = state.fluents
        _, subst = resolve(formula, clauses; funcs=funcs, mode=:all)
    end
    matches = Term[substitute(formula, s) for s in subst]
    return matches
end

"Construct initial state from problem definition."
function init_state(problem::Problem)
    types = Term[@julog($ty(:o)) for (o, ty) in problem.objtypes]
    state = State(problem.init, types)
    return state
end

"Construct goal state from problem definition."
function goal_state(problem::Problem)
    types = Term[@julog($ty(:o)) for (o, ty) in problem.objtypes]
    state = State(flatten_conjs(problem.goal), types)
    return state
end

"Construct initial state from problem definition."
initialize(problem::Problem) = init_state(problem)

"""
    transition(domain, state, action::Term; kwargs...)
    transition(domain, state, actions::Set{<:Term}; kwargs...)

Returns the successor to `state` in the given `domain` after applying a single
`action` or a set of `actions` in parallel, along with any events triggered
by the effects of those actions. Keyword arguments specify whether to `check`
if action preconditions hold, and the `fail_mode` (`:error` or `:no_op`)
if they do not.
"""
function transition(domain::Domain, state::State, action::Term;
                    check::Bool=true, fail_mode::Symbol=:error)
    state = execute(action, state, domain; check=check, fail_mode=fail_mode)
    if length(domain.events) > 0
        state = trigger(domain.events, state, domain)
    end
    return state
end

function transition(domain::Domain, state::State, actions::Set{<:Term};
                    check::Bool=true, fail_mode::Symbol=:error)
    # Execute all actions in parallel
    state = execpar(actions, state, domain, check=check, fail_mode=fail_mode)
    if length(domain.events) > 0
        state = trigger(domain.events, state, domain)
    end
    return state
end

"""
    simulate(domain, state, actions; kwargs...)

Returns the state trajectory that results from applying a sequence of `actions`
to an initial `state` in a given `domain`. Keyword arguments specify whether
to `check` if action preconditions hold, the `fail_mode` (`:error` or `:no_op`)
if they do not, and a `callback` function to apply after each step.
"""
function simulate(domain::Domain, state::State, actions::Vector{<:Term};
                  check::Bool=true, fail_mode::Symbol=:error,
                  callback::Function=(d,s,a)->nothing)
    trajectory = State[state]
    callback(domain, state, Const(:start))
    for act in actions
        state = transition(domain, state, act; check=check, fail_mode=fail_mode)
        push!(trajectory, state)
        callback(domain, state, act)
    end
    return trajectory
end

function simulate(domain::Domain, state::State, actions::Vector{Set{<:Term}};
                  check::Bool=true, fail_mode::Symbol=:error,
                  callback::Function=(d,s,a)->nothing)
    trajectory = State[state]
    callback(domain, state, Set([Const(:start)]))
    for acts in actions
        state = transition(domain, state, acts;
                           check=check, fail_mode=fail_mode)
        push!(trajectory, state)
        callback(domain, state, acts)
    end
    return trajectory
end


text:
"""
This module is used in DeConv, ConvMix and Coupling
"""
module Conv
using FFTW
using DSP
using LinearAlgebra
using LinearMaps

using DSP: nextfastfft
using Misfits

struct D end
struct G end
struct S end

export P_conv
"""
Model : d = convolution(g, s)
d, g and s can have arbitary +ve and -ve lags
"""
mutable struct P_conv{T<:Real, Nd, Ng, Ns}
	dsize::Vector{Int}
	gsize::Vector{Int}
	ssize::Vector{Int}
	"length after zero padding"
	np2::Int # nfastfft length
	d::Array{T,Nd}
	g::Array{T,Ng}
	s::Array{T,Ns}
	dpad::Array{T,Nd}
	gpad::Array{T,Ng}
	spad::Array{T,Ns}
	dfreq::Array{Complex{T},Nd}
	gfreq::Array{Complex{T},Ng}
	sfreq::Array{Complex{T},Ns}
	dfftp::FFTW.rFFTWPlan
	difftp::FFTW.Plan
	gfftp::FFTW.rFFTWPlan
	gifftp::FFTW.Plan
	sfftp::FFTW.rFFTWPlan
	sifftp::FFTW.Plan
	"+ve and -ve lags of g"
	glags::Vector{Int}
	"+ve and -ve lags of s"
	slags::Vector{Int}
	"+ve and -ve lags of d"
	dlags::Vector{Int}
end

function Base.fill!(pa::P_conv, k)
	for iff in [:d, :g, :s]
		m=getfield(pa,iff)
		for i in eachindex(m)
			m[i]=eltype(pa.d)(k)
		end
	end
end

function P_conv(T=Float64;
	       dsize::Vector{Int64}=[1],
	       ssize::Vector{Int64}=[1],
	       gsize::Vector{Int64}=[1],
	       d=zeros(T,dsize...),
	       s=zeros(T,ssize...),
	       g=zeros(T,gsize...),
	       slags=nothing,
	       dlags=nothing,
	       glags=nothing,
	       np2=nextfastfft(ssize[1]+gsize[1]-1), # fft dimension for plan, such that circular convolution is same as linear convolution
	       fftwflag=FFTW.ESTIMATE
	       )
	@assert size(s) == tuple(ssize...)
	@assert size(g) == tuple(gsize...)
	@assert size(d) == tuple(dsize...)

	Nd=length(dsize)
	Ns=length(ssize)
	Ng=length(gsize)
	nts=ssize[1]
	ntg=gsize[1]
	ntd=dsize[1]

	# default lags
	if(slags===nothing)
		# equal +ve and -ve lags for s
		nwplags=div(nts-1,2)
		nwnlags=nts-1-nwplags
		slags=[nwplags, nwnlags]
	end
	if(dlags===nothing)
		# no negative lags
		nsplags=ntd-1
		nsnlags=ntd-1-nsplags
		dlags=[nsplags, nsnlags]
	end
	if(glags===nothing)
		# no negative lags
		nrplags=ntg-1
		nrnlags=ntg-1-nrplags
		glags=[nrplags, nrnlags]
	end

	# check if nlags are consistent with the first dimension of inputs
	(sum(glags)+1 ≠ size(g,1)) && error("glags")
	(sum(dlags)+1 ≠ size(d,1)) && error("dlags")
	(sum(slags)+1 ≠ size(s,1)) && error("slags")

	FFTW.set_num_threads(Sys.CPU_THREADS)
	nrfft=div(np2,2)+1

	dfftp=plan_rfft(zeros(T, np2, dsize[2:end]...),[1], flags=fftwflag,  timelimit=Inf)
	difftp=plan_irfft(complex.(zeros(T, nrfft, dsize[2:end]...)),np2,[1], flags=fftwflag,  timelimit=Inf)

	if(gsize[2:end] ≠  dsize[2:end])
		gfftp=plan_rfft(zeros(T, np2, gsize[2:end]...),[1], flags=fftwflag,  timelimit=Inf)
		gifftp=plan_irfft(complex.(zeros(T, nrfft, gsize[2:end]...)),np2,[1], flags=fftwflag,  timelimit=Inf)
	else
		gfftp=dfftp
		gifftp=difftp
	end

	if(ssize[2:end] ≠ dsize[2:end])
		sfftp=plan_rfft(zeros(T, np2, ssize[2:end]...),[1], flags=fftwflag,  timelimit=Inf)
		sifftp=plan_irfft(complex.(zeros(T, nrfft, ssize[2:end]...)),np2,[1], flags=fftwflag,  timelimit=Inf)
	else
		sfftp=dfftp
		sifftp=difftp

	end

	# preallocate freq domain vectors after rfft
	dfreq=complex.(zeros(T,nrfft,dsize[2:end]...))
	gfreq=complex.(zeros(T,nrfft,gsize[2:end]...))
	sfreq=complex.(zeros(T,nrfft,ssize[2:end]...))

	# preallocate padded arrays
	dpad=(zeros(T,np2,dsize[2:end]...))
	gpad=(zeros(T,np2,gsize[2:end]...))
	spad=(zeros(T,np2,ssize[2:end]...))

	pa=P_conv(dsize, gsize, ssize, np2, d, g, s,
	  dpad, gpad, spad,
	  dfreq, gfreq, sfreq,
	  dfftp, difftp,
	  gfftp, gifftp,
	  sfftp, sifftp,
	  glags, slags, dlags)
end

"""
Convolution that allocates `P_conv` internally.
Don't use inside loops, use `mod!` instead.
"""
function conv!(
	   d::AbstractArray{T,Nd},
	   g::AbstractArray{T,Ng},
	   s::AbstractArray{T,Ns}, attrib::Union{D,G,S};
	       slags=nothing,
	       dlags=nothing,
	       glags=nothing,
	   ) where {T,Nd,Ng,Ns}
	dsize=[size(d)...]
	gsize=[size(g)...]
	ssize=[size(s)...]

	# allocation of freq matrices
	pa=P_conv(dsize=dsize, ssize=ssize, gsize=gsize, g=g, s=s, d=d, slags=slags, dlags=dlags,glags=glags)

	# using pa, return d, g, s according to attrib
	mod!(pa, attrib)
end

function initialize_d!(pa::P_conv, d=pa.d)
	T=eltype(pa.d)
	fill!(pa.dfreq,complex(T(0)))
	fill!(pa.dpad,T(0))
	pad!(d, pa.dpad, pa.dlags[1], pa.dlags[2], pa.np2)
end

function initialize_g!(pa::P_conv, g=pa.g)
	T=eltype(pa.g)
	fill!(pa.gfreq, complex(T(0)))
	fill!(pa.gpad,T(0))
	pad!(g, pa.gpad, pa.glags[1], pa.glags[2], pa.np2)
end

function initialize_s!(pa::P_conv, s=pa.s)
	T=eltype(pa.s)
	fill!(pa.sfreq, complex(T(0)))
	fill!(pa.spad, T(0))
	pad!(s, pa.spad, pa.slags[1], pa.slags[2], pa.np2)
end

# initialize freq vectors
function initialize_all!(pa::P_conv, d=pa.d, g=pa.g, s=pa.s)
	initialize_d!(pa, d);	initialize_g!(pa, g);    initialize_s!(pa, s)
end

"""
Convolution modelling with no allocations at all.
By default, the fields `g`, `d` and `s` in pa are modified accordingly.
Otherwise use keyword arguments to input them.
"""
function mod!(pa::P_conv{T,N,N,N}, ::D; g=pa.g, d=pa.d, s=pa.s) where {N,T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	for i in eachindex(pa.dfreq)
		@inbounds pa.dfreq[i] = pa.gfreq[i] * pa.sfreq[i]
	end
	mul!(pa.dpad, pa.difftp, pa.dfreq)
	truncate!(d, pa.dpad, pa.dlags[1], pa.dlags[2], pa.np2)
	return pa
end

function mod!(pa::P_conv{T,N,N,N}, ::G; g=pa.g, d=pa.d, s=pa.s) where {N,T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.sfreq)
	for i in eachindex(pa.dfreq)
		@inbounds pa.gfreq[i] = pa.dfreq[i] * pa.sfreq[i]
	end
	mul!(pa.gpad, pa.gifftp, pa.gfreq)
	truncate!(g, pa.gpad, pa.glags[1], pa.glags[2], pa.np2)
	return pa
end

function mod!(pa::P_conv{T,N,N,N}, ::S; g=pa.g, d=pa.d, s=pa.s) where {N,T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.gfreq)
	for i in eachindex(pa.dfreq)
		@inbounds pa.sfreq[i] = pa.dfreq[i] * pa.gfreq[i]
	end
	mul!(pa.spad, pa.sifftp, pa.sfreq)
	truncate!(s, pa.spad, pa.slags[1], pa.slags[2], pa.np2)
	return pa
end

# dsum=true
function mod!(pa::P_conv{T,1,2,2}, ::D; g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	for i in CartesianIndices(size(pa.gfreq))
		@inbounds pa.dfreq[i[1]] += pa.gfreq[i] * pa.sfreq[i]
	end
	mul!(pa.dpad, pa.difftp, pa.dfreq)
	truncate!(d, pa.dpad, pa.dlags[1], pa.dlags[2], pa.np2)
	return pa
end

function mod!(pa::P_conv{T,1,2,2}, ::G; g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.sfreq)
	for i in CartesianIndices(size(pa.gfreq))
		@inbounds pa.gfreq[i] = pa.dfreq[i[1]] * pa.sfreq[i]
	end
	mul!(pa.gpad, pa.gifftp, pa.gfreq)
	truncate!(g, pa.gpad, pa.glags[1], pa.glags[2], pa.np2)
	return pa
end

function mod!(pa::P_conv{T,1,2,2}, ::S; g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.gfreq)
	for i in CartesianIndices(size(pa.gfreq))
		@inbounds pa.sfreq[i] = pa.dfreq[i[1]] * pa.gfreq[i]
	end
	mul!(pa.spad, pa.sifftp, pa.sfreq)
	truncate!(s, pa.spad, pa.slags[1], pa.slags[2], pa.np2)
	return pa
end

# gsum=true
function mod!(pa::P_conv{T,2,1,2}, ::D;  g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	for i in CartesianIndices(size(pa.sfreq))
		@inbounds pa.dfreq[i] = pa.gfreq[i[1]] * pa.sfreq[i]
	end
	mul!(pa.dpad, pa.difftp, pa.dfreq)
	truncate!(d, pa.dpad, pa.dlags[1], pa.dlags[2], pa.np2)
	return pa
end
function mod!(pa::P_conv{T,2,1,2}, ::G;  g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.sfreq)
	for i in CartesianIndices(size(pa.sfreq))
		@inbounds pa.gfreq[i[1]] += pa.dfreq[i] * pa.sfreq[i]
	end
	mul!(pa.gpad, pa.gifftp, pa.gfreq)
	truncate!(g, pa.gpad, pa.glags[1], pa.glags[2], pa.np2)
	return pa
end

function mod!(pa::P_conv{T,2,1,2}, ::S;  g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.gfreq)
	for i in CartesianIndices(size(pa.sfreq))
		@inbounds pa.sfreq[i] = pa.dfreq[i] * pa.gfreq[i[1]]
	end
	mul!(pa.spad, pa.sifftp, pa.sfreq)
	truncate!(s, pa.spad, pa.slags[1], pa.slags[2], pa.np2)
	return pa
end

# ssum=true
function mod!(pa::P_conv{T,2,2,1}, ::D; g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	for i in CartesianIndices(size(pa.gfreq))
		@inbounds pa.dfreq[i] = pa.gfreq[i] * pa.sfreq[i[1]]
	end
	mul!(pa.dpad, pa.difftp, pa.dfreq)
	truncate!(d, pa.dpad, pa.dlags[1], pa.dlags[2], pa.np2)
	return pa
end
function mod!(pa::P_conv{T,2,2,1}, ::G; g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.sfreq, pa.sfftp, pa.spad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.sfreq)
	for i in CartesianIndices(size(pa.gfreq))
		@inbounds pa.gfreq[i] = pa.dfreq[i] * pa.sfreq[i[1]]
	end
	mul!(pa.gpad, pa.gifftp, pa.gfreq)
	truncate!(g, pa.gpad, pa.glags[1], pa.glags[2], pa.np2)
	return pa
end
function mod!(pa::P_conv{T,2,2,1}, ::S; g=pa.g, d=pa.d, s=pa.s) where {T<:Real}
	initialize_all!(pa, d, g, s)
	mul!(pa.gfreq, pa.gfftp, pa.gpad)
	mul!(pa.dfreq, pa.dfftp, pa.dpad)
	conj!(pa.gfreq)
	for i in CartesianIndices(size(pa.gfreq))
		@inbounds pa.sfreq[i[1]] += pa.dfreq[i] * pa.gfreq[i]
	end
	mul!(pa.spad, pa.sifftp, pa.sfreq)
	truncate!(s, pa.spad, pa.slags[1], pa.slags[2], pa.np2)
	return pa
end

include("operators.jl")
include("Pad.jl")
include("Xcorr.jl")
include("Misfits.jl")
include("getprop.jl")

end


text:
using SoleBase
using Test

@testset "SoleBase.jl" begin
    # Write your tests here.
end


text:
module StructuredLaTeXTables
using DataFrames

function tabular_header(data_labels, meta_labels, levels, swap_headers)
    number_meta_cols = length(meta_labels)
    number_data_cols = length(data_labels)
    number_levels = length(levels)

    if swap_headers
        # Use levels as the main headers.
        main_labels = levels
        sub_labels = data_labels
        number_main_cols = number_levels
        number_sub_cols = number_data_cols
    else
        # Use levels as the subheaders.
        main_labels = data_labels
        sub_labels = levels
        number_main_cols = number_data_cols
        number_sub_cols = number_levels
    end

    # Begin a new tabular environment with the given column specification.
    spec1 = repeat("l", number_meta_cols)
    spec2 = repeat(repeat("c", number_sub_cols), number_main_cols)
    column_specification = string("{", spec1, spec2, "}")
    output_string = "\\begin{tabular}$(column_specification)"

    # Build the table header with the given nested structure.
    meta_header = repeat("& ", number_meta_cols-1)
    data_header = ""
    horizontal_rules = repeat("& ", number_meta_cols-1)
    start_idx = 1 + number_meta_cols
    stop_idx = start_idx + number_sub_cols - 1
    for label in main_labels
        data_header *= "& \\multicolumn{$(number_sub_cols)}{c}{$(label)} "
        horizontal_rules *= "\\cmidrule(lr){$(start_idx)-$(stop_idx)}"
        start_idx += number_sub_cols
        stop_idx += number_sub_cols
    end
    header = strip(string(meta_header, " ", data_header))
    horizontal_rules = strip(horizontal_rules)

    # Build the table subheader.
    meta_subheader = join(meta_labels, " & ")
    data_subheader = repeat(string(" & ", join(sub_labels, " & ")), number_main_cols)
    subheader = strip(string(meta_subheader, " ", data_subheader))

    # Put everything together and add a horizontal line to mark the end of the header.
    output_string *= "\n\\toprule"
    output_string *= "\n$(header) \\\\"
    output_string *= "\n$(horizontal_rules)"
    output_string *= "\n$(subheader) \\\\"
    output_string *= "\n\\midrule"

    return output_string
end

function tabular_body(df, data_cols, meta_cols, swap_headers)
    output_string = ""

    # Group DataFrame by the selected metadata and iterate over its SubDataFrames.
    for sub in groupby(df, meta_cols)
        # Extract metadata for a particular record.
        row1_meta, row1_data = join(sub[1, meta_cols], " & "), ""

        if swap_headers
            # Build the record by the *rows* of the SubDataFrame.
            for row in eachrow(sub)
                data1_formatted = map(col -> string(row[col]), data_cols)
                row1_data = string(row1_data, " & ", join(data1_formatted, " & "))
            end
        else
            # Build the record by the *columns* of the SubDataFrame.
            for col in data_cols
                data_raw = sub[!, col]
                data1_formatted = map(x -> string(x), data_raw)
                row1_data = string(row1_data, " & ", join(data1_formatted, " & "))
            end
        end
        row1_meta = strip(row1_meta)
        row1_data = strip(row1_data)
        output_string *= "\n$(row1_meta) $(row1_data) \\\\"
    end

    return output_string
end

tabular_footer() = string("\n", raw"\bottomrule", "\n", raw"\end{tabular}")

"""
`make_latex_tabular(df::DataFrame; [data], [metadata], [emphasis], [swap_headers=false])`

Represent structured data from a `DataFrame` as a string containing a `tabular` LaTeX environment with nested columns.

This assumes `data` columns have associated `metadata` columns.

# Options

- `data`: Column-label pairs indicating what data goes in the main body of the table. The labels are used to construct nested columns in the main body.
- `metadata`: Column-label pairs indicating what metadata goes on the left of the table.
- `emphasis`: A column, which should appear in the `metadata`, to emphasize in the main body. The levels for `df[!,emphasis]` are used to construct nested columns in the main body.
- `swap_headers`: Indicate whether to nest `emphasis` levels under `data` columns (`false`) or to nest `data` columns under `emphasis` levels (`true`).
"""
function make_latex_tabular(df::DataFrame;
    data::Vector{Pair{Symbol,String}}=Pair{Symbol,String}[],
    metadata::Vector{Pair{Symbol,String}}=Pair{Symbol,String}[],
    emphasis::Symbol=:none,
    swap_headers::Bool=false,
    )
    data_cols, data_labels = map(first, data), map(last, data)
    meta_cols, meta_labels = map(first, metadata), map(last, metadata)
    levels = unique(df[!, emphasis])

    idx = findfirst(isequal(emphasis), meta_cols)
    deleteat!(meta_cols, idx)
    deleteat!(meta_labels, idx)

    number_meta_cols = length(meta_labels)
    number_levels = length(levels)
    blanks = ["" for _ in 1:number_levels]

    output_string = tabular_header(data_labels, meta_labels, levels, swap_headers)
    output_string *= tabular_body(df, data_cols, meta_cols, swap_headers)
    output_string *= tabular_footer()

    return output_string
end

export make_latex_tabular

end # module


text:
using CodeTransformation
using Test

@testset "CodeTransformation.jl" begin
    @test CodeTransformation.getmodule(typeof(sin)) === Base
    @test CodeTransformation.getmodule(sin) === Base

    let
        # Test example from doctring to addmethod!
        g(x) = x + 13
        ci = code_lowered(g)[1]
        function f end
        addmethod!(f, (Any,), ci)
        @test f(1) === 14

        # Alternative syntax
        function f2 end
        @test CodeTransformation.makesig(f2, (Any,)) === Tuple{typeof(f2), Any}
        addmethod!(Tuple{typeof(f2), Any}, ci)
        @test f2(1) === 14

    end

    let
        # Test example from doctring to codetransform!
        g(x) = x + 13
        function e end
        codetransform!(g => e) do ci
            for ex in ci.code
                if ex isa Expr
                    map!(x -> x === 13 ? 7 : x, ex.args, ex.args)
                end
            end
            ci
        end
        @test e(1) === 8
        @test g(1) === 14
    end

    let
        a = Vector{T} where T
        b = CodeTransformation.typevars(a)
        @test b isa Tuple
        @test length(b) == 1
        @test b[1] isa TypeVar
        @test b[1].name == :T
        @test b[1].lb === Union{}
        @test b[1].ub === Any
    end

end


text:
using StatsPlots
using Clustering
using Distances
using Distributions
using DataFrames

import Base: ∈
import Base: ∩

include("../src/utilities/utils.jl")

# a = rand(Uniform(0,.8), 2, 20)
# b = rand(Uniform(1.2,2), 2, 20)
# c = a + cat(1.2ones(1,20),zeros(1,20), dims=1)
# X = cat(a,b,c, dims=2)

# scatt = scatter(X[1,:], X[2,:], label=false)

# D_a = pairwise(Euclidean(), X, X)
# SL = hclust(D_a, linkage=:single)
# CL = hclust(D_a, linkage=:complete)
# AV = hclust(D_a, linkage=:average)
# plot(SL)
# plot(CL)
# plot(AV)

# propertynames(Hclust)

# propertynames(SL)
# SL.merges
# SL.heights
# SL.order
# SL.linkage
# SL.height
# SL.labels
# SL.merge
# SL.method

# plot(SL.heights)
# SL.heights[56]
# SL.heights[57]
# SL.heights[58]
# SL.heights[59]

# kmeans_ = kmeans(reshape(SL.heights, (1,59)), 2)
# kmeans_ = kmeans(reshape(CL.heights, (1,59)), 2)
# propertynames(kmeans_)
# @show kmeans_

# cutree(SL, k=2)

# D_SL = pairwise(Euclidean(), SL.heights)
# dbscan(reshape(SL.heights, (1,59)), 0.8maximum(SL.heights))
# dbscan(reshape(SL.heights, (1,59)), 0.05maximum(SL.heights))
# dbscan(reshape(SL.heights, (1,59)), 0.2maximum(SL.heights))

# dbscan_ = dbscan(reshape(CL.heights, (1,59)), 0.1maximum(CL.heights))
# dbscan_[1]
# dbscan_[2]
# dbscan(reshape(CL.heights, (1,59)), 0.08maximum(CL.heights))
# dbscan(reshape(CL.heights, (1,59)), 0.2maximum(CL.heights))

# cutree(CL, h=CL.heights[58])
# x = Interval
# x = Interval(1,2)

# struct OrderedPair
#     x::Real
#     y::Real
#     OrderedPair(x,y) = x > y ? error("out of order") : new(x,y)
# end
# OrderedPair(3,1)

# I = Interval(0.2, 4.0)
# -1 ∈ I
# 1 ∈ I
# J = Interval(-1., 3.6)

# # function intersection(I::Interval, J::Interval)
# #     min = maximum([I.a, J.a])
# #     max = minimum([I.b, J.b])

# #     return Interval(min, max)
# # end

# ∩(I,J)

begin
using StatsPlots
using Clustering
using Distances
using Distributions
using DataFrames

import Base: ∈
import Base: ∩

include("../src/utilities/utils.jl")
Data = noisy_circle_data(500, 0.03)

function proj_1d(X)
    min = minimum(X[1, :])
    output = X[1, :] .- min

    return output
end

struct Interval{T<:Real}
    a::T
    b::T
    # Interval(a,b) = a > b ? error("out of order") : new(a,b)
end

mutable struct MapperInformations
    data_number
    original_coordinate
    original_subset
    clustered_subset
    filtered_coordinate
    interval_subset
end

function ∈(x::T, interval::Interval) where T<:Real
    if interval.a ≤ x ≤ interval.b
        return true
    else
        return false
    end
end

function ∩(I::Interval, J::Interval)
    min = maximum([I.a, J.a])
    max = minimum([I.b, J.b])

    return Interval(min, max)
end

function filtered_coordinate(x, p)
    x.filtered_coordinate = p
end

#Y_ = MapperInformations(:None, :None, :None, :None, :None)
end

function mapper(X, filter, interval, overlab)

    original_X = []
    for i in 1:size(X,2)
        push!(original_X, X[:,i])
    end
    # findfirst(==(Info[500].original_coordinate), original_X)

    Info = [MapperInformations(i, X[:, i], Int64[], Tuple{Int64, Int64}[], :None, :None) for i ∈ 1:length(X[1,:])]

    filtered_X = filter(X)
    filtered_coordinate.(Info, filtered_X)
    Info

    min_X = minimum(filtered_X)
    max_X = maximum(filtered_X)
    filter_range = Interval(min_X, max_X)

    length_ = (max_X - min_X)/interval
    ϵ = overlab*length_

    Intervals = [∩(Interval(min_X +(i-1)*length_ - ϵ, min_X +i*length_ + ϵ), filter_range) for i ∈ 1:interval]

    X_subset = interval_test(Info, Intervals)

    Info

    X_subset_clustering = domain_clustering(Info, X_subset)

    # function myshowall(io, x, limit = false)
    #     println(io, summary(x), ":")
    #     Base.print_matrix(IOContext(io, :limit => limit), x)
    # end
    # myshowall(Base.stdout, Info, false)

    # X_subset_clustering[1]
    # X_subset_clustering[2]
    # X_subset_clustering[3]
    # X_subset_clustering[4]
    # X_subset_clustering[5]

    link = []

    for i ∈ 1:length(X_subset)-1
        for j ∈ 1:length(X_subset_clustering[i])
            for k ∈ 1:length(X_subset_clustering[i+1])
                if length(∩(X_subset_clustering[i][j], X_subset_clustering[i+1][k])) > 0
                    push!(link, [[i,j], [i+1,k]])
                end
            end
        end
    end

    center_ = coor_mean(X_subset_clustering, Info)

    center_number = 0
    for i ∈ 1:length(X_subset)
        for j ∈ 1:length(X_subset_clustering[i])
            center_number +=1
        end
    end

    center_matrix = zeros(2, center_number)
    k_temp = 1
    for i ∈ 1:length(X_subset)
        for j ∈ 1:length(X_subset_clustering[i])
            center_matrix[:,k_temp] = center_[i][j]
            k_temp += 1
        end
    end

    zcolor_temp = filter(cat(Data, center_matrix, dims=2))

    zcolor = [[] for i ∈ 1:length(X_subset)]
    k_temp = 1
    for i ∈ 1:length(X_subset)
        for j ∈ 1:length(X_subset_clustering[i])
            push!(zcolor[i], zcolor_temp[size(Data)[2]+k_temp])
            k_temp +=1
        end
    end

    return link, X_subset_clustering, Info, center_, zcolor
end

function mapper_plot(link, X_subset_clustering, center_, zcolor, xlim, ylim)

    linked_center = zeros(length(link), 4)
    for i ∈ 1:length(link)
        linked_center[i, 1:2] = center_[link[i][1][1]][link[i][1][2]]
        linked_center[i, 3:4] = center_[link[i][2][1]][link[i][2][2]]
    end
    linked_center

    p = plot([linked_center[1,1], linked_center[1,3]], [linked_center[1,2], linked_center[1,4]], color=:black, xlim=xlim, ylim=ylim, label="")
    for i ∈ 2:length(link)
        plot!([linked_center[i,1], linked_center[i,3]], [linked_center[i,2], linked_center[i,4]], color=:black, label="")
    end

    for i ∈ 1:length(X_subset_clustering)
        for j ∈ 1:length(X_subset_clustering[i])
            scatter!([center_[i][j][1]], [center_[i][j][2]], markersize=25*(452.106/181)*sqrt(length(X_subset_clustering[i][j])/500), zcolor=zcolor[i][j], label="", color=:jet, clim=(0,Inf))
            length_temp = length(X_subset_clustering[i][j])
            annotate!([center_[i][j][1]], [center_[i][j][2]], "$length_temp")
        end
    end

    return p
end

d1 = scatter(Data[1,:], Data[2,:], label="")
d2 = scatter(Data[1,:], Data[2,:], zcolor=proj_1d(Data), color=:jet, label="")

link, X_subset_clustering, Info, center_, zcolor = mapper(Data, proj_1d, 5, 0.2)
p = mapper_plot(link, X_subset_clustering, center_, zcolor, (-1.5, 1.5), (-1.5,1.5))

plot(d1, d2, p, size=(1400,800))

# plot(d, p, size=(1000,400))
# include("../src/utilities/tDATAsets.jl")

# Data = rand(Sphere(2), 500)
# Data = cat(rand(Infty(), 500), rand(Sphere(2), 500).+1, dims=2)

# link, X_subset_clustering, Info = mapper(Data, proj_1d, 5, 0.3)
# p = mapper_plot(link, X_subset_clustering, Info, (-1.5, 1.5), (-1.5,1.5))

# plot(d, p, size=(1000,400))

# d = scatter(Data[1,:], Data[2,:])
# link, X_subset_clustering, Info = mapper(Data, proj_1d, 10, 0.2)
# p = mapper_plot(link, X_subset_clustering, Info)

# plot(d, p, size=(1000,400))


text:
module GenRedner

using GenPyTorch; TorchGenerativeFunction, TorchArg
using PyCall

# coordinate frame: https://learnopencv.com/wp-content/uploads/2020/02/world-camera-image-coordinates.png

function get_depth_renderer(;num_samples=4, print_timing=false)

py"""
import pyredner
import torch
import torch.nn as nn

pyredner.set_print_timing($(print_timing))

def construct_intrinsic_mat(intrinsics_vector, width, height):
    # see https://github.com/BachiLi/redner/blob/96d3e27cdd39129b1dc00e0271d3e27c0c5cec8f/pyredner/camera.py#L39-L50
    # for the specification of the intrinsics matrix that this function returns
    fx = intrinsics_vector[0:1]
    fy = intrinsics_vector[1:2]
    cx = intrinsics_vector[2:3]
    cy = intrinsics_vector[3:4]
    a = fx * 2 / width
    b = -fy * 2 / width
    x0 = (2 * cx / width) - 1
    y0 = -(2 * cy / width) + (height / width)
    zero = torch.zeros_like(a)
    one = torch.ones_like(a)
    skew = zero
    row0 = torch.cat([a, skew, x0])
    row1 = torch.cat([ zero,   b, y0])
    row2 = torch.cat([ zero,    zero,  one])
    intrinsic_mat = torch.stack([row0, row1, row2]).contiguous()
    return intrinsic_mat

class RednerDepthRenderer(torch.nn.Module):

    def __init__(self, num_samples):
        super(RednerDepthRenderer, self).__init__()
        self.num_samples = num_samples

    def forward(self, vertices, indices, intrinsics, dims):

        # vertices is num_vertices x 3 tensor of floats
        # indices is num_triangles x 3 tensor of ints, and starts from one

        # shift indices down by one, so they start from zero
        indices = indices - 1

        (width, height) = dims
        intrinsic_mat = construct_intrinsic_mat(intrinsics, width, height)

        # the mesh vertices must already be in the camera's 3D coordinate frame
        cam_to_world = torch.eye(4)

        camera = pyredner.Camera(
            resolution=(height, width),
            cam_to_world=cam_to_world,
            clip_near=0.0001,
            intrinsic_mat=intrinsic_mat)

        objects = [pyredner.Object(vertices.contiguous(), indices.contiguous(), pyredner.Material())]

        scene = pyredner.Scene(camera=camera, objects=objects)

        # return a height x width x 3 tensor of floats
        # the channels give x, y, and z coordinates in the 3D camera frame, respectively
        # so the third channel points[:,:,2] is the depth
        points = pyredner.render_g_buffer(scene, [pyredner.channels.position], self.num_samples)
        assert tuple(points.size()) == (height, width, 3)
        return points
"""

    return TorchGenerativeFunction(
    py"RednerDepthRenderer($(num_samples))", [
        TorchArg(true, py"torch.float"), # vertices
        TorchArg(false, py"torch.int"), # indices
        TorchArg(true, py"torch.float"), # intrinsics
        TorchArg(false, py"torch.int"), # dims
    ], 1)

end

export get_depth_renderer

end # module


text:
m = model(data...) # instantiate model on the data
q = vi(m, vi_alg)  # perform VI on `m` using the VI method `vi_alg`, which returns a `VariationalPosterior`

using Random
using Turing
using Turing: Variational

Random.seed!(42);

# generate data
x = randn(2000);

@model model(x) = begin
    s ~ InverseGamma(2, 3)
    m ~ Normal(0.0, sqrt(s))
    for i = 1:length(x)
        x[i] ~ Normal(m, sqrt(s))
    end
end;

# Instantiate model
m = model(x);

samples_nuts = sample(m, NUTS(200, 0.65), 10000);

@doc(Variational.vi)

@doc(Variational.meanfield)

@doc(Variational.ADVI)

# ADVI
advi = ADVI(10, 1000)
q = vi(m, advi);

q isa MultivariateDistribution

rand(q)

logpdf(q, rand(q))

var(x), mean(x)

(mean(rand(q, 1000); dims = 2)..., )

samples = rand(q, 10000);

# setup for plotting
using Plots, LaTeXStrings, StatsPlots

p1 = histogram(samples[1, :], bins=100, normed=true, alpha=0.2, color = :blue, label = "")
density!(samples[1, :], label = "s (ADVI)", color = :blue, linewidth = 2)
density!(samples_nuts, :s; label = "s (NUTS)", color = :green, linewidth = 2)
vline!([var(x)], label = "s (data)", color = :black)
vline!([mean(samples[1, :])], color = :blue, label ="")

p2 = histogram(samples[2, :], bins=100, normed=true, alpha=0.2, color = :blue, label = "")
density!(samples[2, :], label = "m (ADVI)", color = :blue, linewidth = 2)
density!(samples_nuts, :m; label = "m (NUTS)", color = :green, linewidth = 2)
vline!([mean(x)], color = :black, label = "m (data)")
vline!([mean(samples[2, :])], color = :blue, label="")

plot(p1, p2, layout=(2, 1), size=(900, 500))

# used to compute closed form expression of posterior
using ConjugatePriors

# closed form computation
# notation mapping has been verified by explicitly computing expressions
# in "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy
μ₀ = 0.0 # => μ
κ₀ = 1.0 # => ν, which scales the precision of the Normal
α₀ = 2.0 # => "shape"
β₀ = 3.0 # => "rate", which is 1 / θ, where θ is "scale"

# prior
pri = NormalGamma(μ₀, κ₀, α₀, β₀)

# posterior
post = posterior(pri, Normal, x)

# marginal distribution of τ = 1 / σ²
# Eq. (90) in "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy
# `scale(post)` = θ
p_τ = Gamma(post.shape, scale(post))
p_σ²_pdf = z -> pdf(p_τ, 1 / z) # τ => 1 / σ²

# marginal of μ
# Eq. (91) in "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy
p_μ = TDist(2 * post.shape)

μₙ = post.mu    # μ → μ
κₙ = post.nu    # κ → ν
αₙ = post.shape # α → shape
βₙ = post.rate  # β → rate

# numerically more stable but doesn't seem to have effect; issue is probably internal to
# `pdf` which needs to compute ≈ Γ(1000)
p_μ_pdf = z -> exp(logpdf(p_μ, (z - μₙ) * exp(- 0.5 * log(βₙ) + 0.5 * log(αₙ) + 0.5 * log(κₙ))))

# posterior plots
p1 = plot();
histogram!(samples[1, :], bins=100, normed=true, alpha=0.2, color = :blue, label = "")
density!(samples[1, :], label = "s (ADVI)", color = :blue)
density!(samples_nuts, :s; label = "s (NUTS)", color = :green)
vline!([mean(samples[1, :])], linewidth = 1.5, color = :blue, label ="")

# normalize using Riemann approx. because of (almost certainly) numerical issues
Δ = 0.001
r = 0.75:0.001:1.50
norm_const = sum(p_σ²_pdf.(r) .* Δ)
plot!(r, p_σ²_pdf, label = "s (posterior)", color = :red);
vline!([var(x)], label = "s (data)", linewidth = 1.5, color = :black, alpha = 0.7);
xlims!(0.75, 1.35);

p2 = plot();
histogram!(samples[2, :], bins=100, normed=true, alpha=0.2, color = :blue, label = "")
density!(samples[2, :], label = "m (ADVI)", color = :blue)
density!(samples_nuts, :m; label = "m (NUTS)", color = :green)
vline!([mean(samples[2, :])], linewidth = 1.5, color = :blue, label="")

# normalize using Riemann approx. because of (almost certainly) numerical issues
Δ = 0.0001
r = -0.1 + mean(x):Δ:0.1 + mean(x)
norm_const = sum(p_μ_pdf.(r) .* Δ)
plot!(r, z -> p_μ_pdf(z) / norm_const, label = "m (posterior)", color = :red);
vline!([mean(x)], label = "m (data)", linewidth = 1.5, color = :black, alpha = 0.7);

xlims!(-0.25, 0.25);

p = plot(p1, p2; layout=(2, 1), size=(900, 500))

Random.seed!(1);

# Import RDatasets.
using RDatasets

# Hide the progress prompt while sampling.
Turing.setprogress!(false);

# Import the "Default" dataset.
data = RDatasets.dataset("datasets", "mtcars");

# Show the first six rows of the dataset.
first(data, 6)

# Function to split samples.
function split_data(df, at = 0.70)
    r = size(df,1)
    index = Int(round(r * at))
    train = df[1:index, :]
    test  = df[(index+1):end, :]
    return train, test
end

# A handy helper function to rescale our dataset.
function standardize(x)
    return (x .- mean(x, dims=1)) ./ std(x, dims=1), x
end

# Another helper function to unstandardize our datasets.
function unstandardize(x, orig)
    return (x .+ mean(orig, dims=1)) .* std(orig, dims=1)
end

# Remove the model column.
select!(data, Not(:Model))

# Standardize our dataset.
(std_data, data_arr) = standardize(Matrix(data))

# Split our dataset 70%/30% into training/test sets.
train, test = split_data(std_data, 0.7)

# Save dataframe versions of our dataset.
train_cut = DataFrame(train, names(data))
test_cut = DataFrame(test, names(data))

# Create our labels. These are the values we are trying to predict.
train_label = train_cut[:, :MPG]
test_label = test_cut[:, :MPG]

# Get the list of columns to keep.
remove_names = filter(x->!in(x, [:MPG, :Model]), names(data))

# Filter the test and train sets.
train = Matrix(train_cut[:,remove_names]);
test = Matrix(test_cut[:,remove_names]);

# Bayesian linear regression.
@model linear_regression(x, y, n_obs, n_vars, ::Type{T}=Vector{Float64}) where {T} = begin
    # Set variance prior.
    σ₂ ~ truncated(Normal(0,100), 0, Inf)

    # Set intercept prior.
    intercept ~ Normal(0, 3)

    # Set the priors on our coefficients.
    coefficients ~ MvNormal(zeros(n_vars), 10 * ones(n_vars))

    # Calculate all the mu terms.
    mu = intercept .+ x * coefficients
    y ~ MvNormal(mu, σ₂)
end;

n_obs, n_vars = size(train)
m = linear_regression(train, train_label, n_obs, n_vars);

q0 = Variational.meanfield(m)
typeof(q0)

advi = ADVI(10, 10_000)

using Flux, Turing
using Turing.Variational

vi(m, advi; optimizer = Flux.ADAM())

opt = Variational.DecayedADAGrad(1e-2, 1.1, 0.9)

q = vi(m, advi, q0; optimizer = opt)
typeof(q)

z = rand(q, 10_000);

avg = vec(mean(z; dims = 2))

_, sym2range = bijector(m, Val(true));
sym2range

avg[union(sym2range[:σ₂]...)]

avg[union(sym2range[:intercept]...)]

avg[union(sym2range[:coefficients]...)]

function plot_variational_marginals(z, sym2range)
    ps = []

    for (i, sym) in enumerate(keys(sym2range))
        indices = union(sym2range[sym]...)  # <= array of ranges
        if sum(length.(indices)) > 1
            offset = 1
            for r in indices
                for j in r
                    p = density(z[j, :], title = "$(sym)[$offset]", titlefontsize = 10, label = "")
                    push!(ps, p)

                    offset += 1
                end
            end
        else
            p = density(z[first(indices), :], title = "$(sym)", titlefontsize = 10, label = "")
            push!(ps, p)
        end
    end

    return plot(ps..., layout = (length(ps), 1), size = (500, 1500))
end

plot_variational_marginals(z, sym2range)

chain = sample(m, NUTS(0.65), 10_000);

plot(chain)

vi_mean = vec(mean(z; dims = 2))[[union(sym2range[:coefficients]...)..., union(sym2range[:intercept]...)..., union(sym2range[:σ₂]...)...]]

mean(chain).nt.mean

sum(abs2, mean(chain).nt.mean .- vi_mean)

# Import the GLM package.
using GLM

# Perform multivariate OLS.
ols = lm(@formula(MPG ~ Cyl + Disp + HP + DRat + WT + QSec + VS + AM + Gear + Carb), train_cut)

# Store our predictions in the original dataframe.
train_cut.OLSPrediction = unstandardize(GLM.predict(ols), data.MPG);
test_cut.OLSPrediction = unstandardize(GLM.predict(ols, test_cut), data.MPG);

# Make a prediction given an input vector.
function prediction_chain(chain, x)
    p = get_params(chain)
    α = mean(p.intercept)
    β = collect(mean.(p.coefficients))
    return  α .+ x * β
end

# Make a prediction using samples from the variational posterior given an input vector.
function prediction(samples::AbstractVector, sym2ranges, x)
    α = mean(samples[union(sym2ranges[:intercept]...)])
    β = vec(mean(samples[union(sym2ranges[:coefficients]...)]; dims = 2))
    return  α .+ x * β
end

function prediction(samples::AbstractMatrix, sym2ranges, x)
    α = mean(samples[union(sym2ranges[:intercept]...), :])
    β = vec(mean(samples[union(sym2ranges[:coefficients]...), :]; dims = 2))
    return  α .+ x * β
end

# Unstandardize the dependent variable.
train_cut.MPG = unstandardize(train_cut.MPG, data.MPG);
test_cut.MPG = unstandardize(test_cut.MPG, data.MPG);

# Show the first side rows of the modified dataframe.
first(test_cut, 6)

z = rand(q, 10_000);

# Calculate the predictions for the training and testing sets using the samples `z` from variational posterior
train_cut.VIPredictions = unstandardize(prediction(z, sym2range, train), data.MPG);
test_cut.VIPredictions = unstandardize(prediction(z, sym2range, test), data.MPG);

train_cut.BayesPredictions = unstandardize(prediction_chain(chain, train), data.MPG);
test_cut.BayesPredictions = unstandardize(prediction_chain(chain, test), data.MPG);

vi_loss1 = mean((train_cut.VIPredictions - train_cut.MPG).^2)
bayes_loss1 = mean((train_cut.BayesPredictions - train_cut.MPG).^2)
ols_loss1 = mean((train_cut.OLSPrediction - train_cut.MPG).^2)

vi_loss2 = mean((test_cut.VIPredictions - test_cut.MPG).^2)
bayes_loss2 = mean((test_cut.BayesPredictions - test_cut.MPG).^2)
ols_loss2 = mean((test_cut.OLSPrediction - test_cut.MPG).^2)

println("Training set:
    VI loss: $vi_loss1
    Bayes loss: $bayes_loss1
    OLS loss: $ols_loss1
Test set:
    VI loss: $vi_loss2
    Bayes loss: $bayes_loss2
    OLS loss: $ols_loss2")

z = rand(q, 1000);
preds = hcat([unstandardize(prediction(z[:, i], sym2range, test), data.MPG) for i = 1:size(z, 2)]...);

scatter(1:size(test, 1), mean(preds; dims = 2), yerr=std(preds; dims = 2), label="prediction (mean ± std)", size = (900, 500), markersize = 8)
scatter!(1:size(test, 1), unstandardize(test_label, data.MPG), label="true")
xaxis!(1:size(test, 1))
ylims!(95, 140)
title!("Mean-field ADVI (Normal)")

preds = hcat([unstandardize(prediction_chain(chain[i], test), data.MPG) for i = 1:5:size(chain, 1)]...);

scatter(1:size(test, 1), mean(preds; dims = 2), yerr=std(preds; dims = 2), label="prediction (mean ± std)", size = (900, 500), markersize = 8)
scatter!(1:size(test, 1), unstandardize(test_label, data.MPG), label="true")
xaxis!(1:size(test, 1))
ylims!(95, 140)
title!("MCMC (NUTS)")

using Bijectors

using Bijectors: Scale, Shift

d = length(q)
base_dist = Turing.DistributionsAD.TuringDiagMvNormal(zeros(d), ones(d))

to_constrained = inv(bijector(m));

function getq(θ)
    d = length(θ) ÷ 2
    A = @inbounds θ[1:d]
    b = @inbounds θ[d + 1: 2 * d]

    b = to_constrained ∘ Shift(b; dim = Val(1)) ∘ Scale(exp.(A); dim = Val(1))

    return transformed(base_dist, b)
end

q_mf_normal = vi(m, advi, getq, randn(2 * d));

p1 = plot_variational_marginals(rand(q_mf_normal, 10_000), sym2range) # MvDiagNormal + Affine transformation + to_constrained
p2 = plot_variational_marginals(rand(q, 10_000), sym2range)  # Turing.meanfield(m)

plot(p1, p2, layout = (1, 2), size = (800, 2000))

using LinearAlgebra

# Using `ComponentArrays.jl` together with `UnPack.jl` makes our lives much easier.
using ComponentArrays, UnPack

proto_arr = ComponentArray(
    L = zeros(d, d),
    b = zeros(d)
)
proto_axes = proto_arr |> getaxes
num_params = length(proto_arr)

function getq(θ)
    L, b = begin
        @unpack L, b = ComponentArray(θ, proto_axes)
        LowerTriangular(L), b
    end
    # For this to represent a covariance matrix we need to ensure that the diagonal is positive.
    # We can enforce this by zeroing out the diagonal and then adding back the diagonal exponentiated.
    D = Diagonal(diag(L))
    A = L - D + exp(D) # exp for Diagonal is the same as exponentiating only the diagonal entries

    b = to_constrained ∘ Shift(b; dim = Val(1)) ∘ Scale(A; dim = Val(1))

    return transformed(base_dist, b)
end

advi = ADVI(10, 20_000)

q_full_normal = vi(m, advi, getq, randn(num_params); optimizer = Variational.DecayedADAGrad(1e-2));

A = q_full_normal.transform.ts[1].a

heatmap(cov(A * A'))

zs = rand(q_full_normal, 10_000);

p1 = plot_variational_marginals(rand(q_mf_normal, 10_000), sym2range)
p2 = plot_variational_marginals(rand(q_full_normal, 10_000), sym2range)

plot(p1, p2, layout = (1, 2), size = (800, 2000))

# Unfortunately, it seems like this has quite a high variance which is likely to be due to numerical instability,
# so we consider a larger number of samples. If we get a couple of outliers due to numerical issues,
# these kind affect the mean prediction greatly.
z = rand(q_full_normal, 10_000);

train_cut.VIFullPredictions = unstandardize(prediction(z, sym2range, train), data.MPG);
test_cut.VIFullPredictions = unstandardize(prediction(z, sym2range, test), data.MPG);

vi_loss1 = mean((train_cut.VIPredictions - train_cut.MPG).^2)
vifull_loss1 = mean((train_cut.VIFullPredictions - train_cut.MPG).^2)
bayes_loss1 = mean((train_cut.BayesPredictions - train_cut.MPG).^2)
ols_loss1 = mean((train_cut.OLSPrediction - train_cut.MPG).^2)

vi_loss2 = mean((test_cut.VIPredictions - test_cut.MPG).^2)
vifull_loss2 = mean((test_cut.VIFullPredictions - test_cut.MPG).^2)
bayes_loss2 = mean((test_cut.BayesPredictions - test_cut.MPG).^2)
ols_loss2 = mean((test_cut.OLSPrediction - test_cut.MPG).^2)

println("Training set:
    VI loss: $vi_loss1
    Bayes loss: $bayes_loss1
    OLS loss: $ols_loss1
Test set:
    VI loss: $vi_loss2
    Bayes loss: $bayes_loss2
    OLS loss: $ols_loss2")

z = rand(q_mf_normal, 1000);
preds = hcat([unstandardize(prediction(z[:, i], sym2range, test), data.MPG) for i = 1:size(z, 2)]...);

p1 = scatter(1:size(test, 1), mean(preds; dims = 2), yerr=std(preds; dims = 2), label="prediction (mean ± std)", size = (900, 500), markersize = 8)
scatter!(1:size(test, 1), unstandardize(test_label, data.MPG), label="true")
xaxis!(1:size(test, 1))
ylims!(95, 140)
title!("Mean-field ADVI (Normal)")

z = rand(q_full_normal, 1000);
preds = hcat([unstandardize(prediction(z[:, i], sym2range, test), data.MPG) for i = 1:size(z, 2)]...);

p2 = scatter(1:size(test, 1), mean(preds; dims = 2), yerr=std(preds; dims = 2), label="prediction (mean ± std)", size = (900, 500), markersize = 8)
scatter!(1:size(test, 1), unstandardize(test_label, data.MPG), label="true")
xaxis!(1:size(test, 1))
ylims!(95, 140)
title!("Full ADVI (Normal)")

preds = hcat([unstandardize(prediction_chain(chain[i], test), data.MPG) for i = 1:5:size(chain, 1)]...);

p3 = scatter(1:size(test, 1), mean(preds; dims = 2), yerr=std(preds; dims = 2), label="prediction (mean ± std)", size = (900, 500), markersize = 8)
scatter!(1:size(test, 1), unstandardize(test_label, data.MPG), label="true")
xaxis!(1:size(test, 1))
ylims!(95, 140)
title!("MCMC (NUTS)")

plot(p1, p2, p3, layout = (1, 3), size = (900, 250), label="")

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: leachim/TuringTutorials

max_issues_repo_name: leachim/TuringTutorials

max_forks_repo_name: leachim/TuringTutorials


text:
"""
    BayesianSVM()

The [Bayesian SVM](https://arxiv.org/abs/1707.05532) is a Bayesian interpretation of the classical SVM.
```math
p(y|f) \\propto \\exp(2 \\max(1-yf, 0))
````

---

For the analytic version of the likelihood, it is augmented via:

```math
p(y|f, ω) = \\frac{1}{\\sqrt(2\\pi\\omega) \\exp(-\\frac{(1+\\omega-yf)^2}{2\\omega}))
```

where ``ω ∼ 𝟙[0,∞)`` has an improper prior (his posterior is however has a valid distribution, a Generalized Inverse Gaussian). For reference [see this paper](http://ecmlpkdd2017.ijs.si/papers/paperID502.pdf)
"""
struct BayesianSVM{T<:Real,A<:AbstractVector{T}} <: ClassificationLikelihood{T}
    ω::A
    θ::A
    function BayesianSVM{T}() where {T<:Real}
        return new{T,Vector{T}}()
    end
    function BayesianSVM{T}(ω::A, θ::A) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(ω, θ)
    end
end

function BayesianSVM()
    return BayesianSVM{Float64}()
end

implemented(::BayesianSVM, ::AnalyticVI) = true

function init_likelihood(
    ::BayesianSVM{T}, ::AbstractInference{T}, ::Int, nSamplesUsed::Int
) where {T}
    return BayesianSVM{T}(rand(T, nSamplesUsed), zeros(T, nSamplesUsed))
end
function (::BayesianSVM)(y::Real, f::Real)
    return svmlikelihood(y * f)
end

function Base.show(io::IO, ::BayesianSVM{T}) where {T}
    return print(io, "Bayesian SVM")
end

"""Return likelihood equivalent to SVM hinge loss"""
function svmlikelihood(f::Real)
    pos = svmpseudolikelihood(f)
    return pos ./ (pos .+ svmpseudolikelihood(-f))
end

"""Return the pseudo likelihood of the SVM hinge loss"""
function svmpseudolikelihood(f::Real)
    return exp(-2.0 * max.(1.0 - f, 0))
end

function compute_proba(
    ::BayesianSVM{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    N = length(μ)
    pred = zeros(T, N)
    sig_pred = zeros(T, N)
    for i in 1:N
        x = pred_nodes .* sqrt(max(σ²[i], zero(T))) .+ μ[i]
        pred[i] = dot(pred_weights, svmlikelihood.(x))
        sig_pred[i] = max(dot(pred_weights, svmlikelihood.(x) .^ 2) - pred[i]^2, zero(T))
    end
    return pred, sig_pred
end

## Updates

function local_updates!(
    l::BayesianSVM{T}, y::AbstractVector, μ::AbstractVector, diagΣ::AbstractVector
) where {T}
    @. l.ω = abs2(one(T) - y * μ) + diagΣ
    @. l.θ = inv(sqrt(l.ω))
end

@inline function ∇E_μ(l::BayesianSVM{T}, ::AOptimizer, y::AbstractVector) where {T}
    return (y .* (l.θ .+ one(T)),)
end

@inline ∇E_Σ(l::BayesianSVM{T}, ::AOptimizer, ::AbstractVector) where {T} = (0.5 .* l.θ,)

## Lower bounds

function expec_loglikelihood(
    l::BayesianSVM{T},
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
) where {T}
    tot = -(0.5 * length(y) * logtwo)
    tot += dot(μ, y)
    tot += -0.5 * dot(l.θ, diag_cov) + dot(l.θ, abs2.(one(T) .- y .* μ))
    return tot
end

AugmentedKL(l::BayesianSVM, ::AbstractVector) = Zygote.@ignore(GIGEntropy(l))

function GIGEntropy(l::BayesianSVM)
    return 0.5 * sum(log.(l.ω)) + sum(log.(2.0 * besselk.(0.5, sqrt.(l.ω)))) -
           0.5 * sum(sqrt.(l.ω))
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

include("tree_mesh.jl")
include("structured_mesh.jl")
include("surface_interpolant.jl")
include("unstructured_mesh.jl")
include("face_interpolant.jl")
include("transfinite_mappings_3d.jl")
include("p4est_mesh.jl")
include("mesh_io.jl")
include("dgmulti_meshes.jl")

end # @muladd


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.plaindigraph" begin
    g = plaindigraph(5, [(1, 2), (3, 4), (4, 3), (4, 5)])
    @test neighbors(g, 4)[4] == 5
    @test outneighbors(g, 1)[1] == 2
    @test inneighbors(g, 2)[1] == 1
    @test getedge(g, 2) == (3, 4)
    @test getedge(g, 3) == (4, 3)
    @test !hasedge(g, 2, 1)
    @test_throws ErrorException nodeattr(g, 1)
    @test_throws ErrorException edgeattr(g, 1)
    @test_throws ErrorException edgeattr(g, 3, 4)

    @test issetequal(adjacencies(g, 4), [3, 5])
    @test issetequal(successors(g, 3), [4])
    @test issetequal(predecessors(g, 5), [4])
    @test issetequal(incidences(g, 4), [2, 3, 4])
    @test issetequal(outincidences(g, 4), [3, 4])
    @test issetequal(in_incidences(g, 4), [2])
    @test issetequal(nodeset(g), 1:5)
    @test issetequal(edgeset(g), 1:4)
    @test_throws ErrorException nodeattrs(g)
    @test_throws ErrorException edgeattrs(g)

    @test nodecount(g) == 5
    @test edgecount(g) == 4
    @test neighborcount(g, 2) == 1
    @test degree(g, 4) == 3
    @test outdegree(g, 4) == 2
    @test indegree(g, 4) == 1

    addnode!(g)
    @test degree(g, 6) == 0
    addedge!(g, 4, 6)
    @test outdegree(g, 4) == 3
    @test_throws ErrorException nodeattrtype(g)
    @test_throws ErrorException edgeattrtype(g)

    emp = plaindigraph()
    @test nodecount(emp) == 0

    # TODO build from attributed graph
    # TODO build from subgraph
    # TODO unlink
end # graph.plaindigraph


text:
include("activate.jl")

import Random
using SparseArrays, LinearAlgebra
using Test
using Plots

const nfeat = 100 * 5
const nobs = 500

# rank of the real data
const r = 30
const Xreal = Matrix{Float64}(undef, nobs, nfeat)

const X_gen_cols = randn(nfeat, r)
const X_gen_rows = randn(r, nobs)
const svals = 100 * rand(r)
for i in 1:nobs
    for j in 1:nfeat
        Xreal[i, j] = sum(X_gen_cols[j, k] * X_gen_rows[k, i] * svals[k] for k in 1:r)
    end
end

nucnorm(Xmat) = sum(abs(σi) for σi in LinearAlgebra.svdvals(Xmat))

@test rank(Xreal) == r

# 0.2 of entries missing
const missing_entries = unique!([(rand(1:nobs), rand(1:nfeat)) for _ in 1:10000])
const present_entries = [(i, j) for i in 1:nobs, j in 1:nfeat if (i, j) ∉ missing_entries]

f(X) = 0.5 * sum((X[i, j] - Xreal[i, j])^2 for (i, j) in present_entries)

function grad!(storage, X)
    storage .= 0
    for (i, j) in present_entries
        storage[i, j] = X[i, j] - Xreal[i, j]
    end
    return nothing
end

const lmo = FrankWolfe.NuclearNormLMO(275_000.0)
const x0 = FrankWolfe.compute_extreme_point(lmo, zero(Xreal))

FrankWolfe.benchmark_oracles(f, grad!, () -> randn(size(Xreal)), lmo; k=100)

# gradient descent
gradient = similar(x0)
xgd = Matrix(x0)
for _ in 1:5000
    @info f(xgd)
    grad!(gradient, xgd)
    xgd .-= 0.01 * gradient
    if norm(gradient) ≤ sqrt(eps())
        break
    end
end

grad!(gradient, x0)
v0 = FrankWolfe.compute_extreme_point(lmo, gradient)
@test dot(v0 - x0, gradient) < 0

const k = 500

xfin, vmin, _, _, traj_data = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e7,
    max_iteration=k,
    print_iter=k / 10,
    trajectory=true,
    verbose=true,
    linesearch_tol=1e-7,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,
    gradient=spzeros(size(x0)...),
)

xfinAFW, vmin, _, _, traj_data = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e7,
    max_iteration=k,
    print_iter=k / 10,
    trajectory=true,
    verbose=true,
    linesearch_tol=1e-7,
    lazy=true,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,#,
)

xfinBCG, vmin, _, _, traj_data = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0;
    epsilon=1e7,
    max_iteration=k,
    print_iter=k / 10,
    trajectory=true,
    verbose=true,
    linesearch_tol=1e-7,
    line_search=FrankWolfe.Adaptive(),
    emphasis=FrankWolfe.memory,#,
)

plot(svdvals(xfin), label="FW", width=3, yaxis=:log)
plot!(svdvals(xfinAFW), label="AFW", width=3, yaxis=:log)
plot!(svdvals(xfinBCG), label="BCG", width=3, yaxis=:log)
plot!(svdvals(xgd), label="Gradient descent", width=3, yaxis=:log)
plot!(svdvals(Xreal), label="Real matrix", linestyle=:dash, width=3, color=:black)
title!("Singular values")


text:
using OVERT
using OVERTVerify
include("../src/overt_parser_minimal.jl")
include("../src/overt_to_mip.jl")

@testset "overt_2_mip_test1" begin
    oA_test = OverApproximation()
    oA_test.approx_eq = [:(x1 == 1-x), :(x2 == 2.0*max(0, x) + 3.0*max(0, x1))]
    oA_test.ranges = Dict(:x => [-1., 2.], :x1 => [-1., 2.], :x2 => [2., 6.])
    overt_mip_model = OvertMIP(oA_test)
    JuMP.optimize!(overt_mip_model.model)
    for (x, x_var) in pairs(overt_mip_model.vars_dict)
       l, u = oA_test.ranges[x]
       @test l <= value(x_var)
       @test u >= value(x_var)
    end

    @objective(overt_mip_model.model, Min, overt_mip_model.vars_dict[:x2])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 2.

    @objective(overt_mip_model.model, Max, overt_mip_model.vars_dict[:x2])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 6.

    @objective(overt_mip_model.model, Min, overt_mip_model.vars_dict[:x2] - overt_mip_model.vars_dict[:x])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 1.

    @objective(overt_mip_model.model, Max, overt_mip_model.vars_dict[:x2] - overt_mip_model.vars_dict[:x])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 7.
end

@testset "overt_2_mip_test2" begin
    oA_test = OverApproximation()
    oA_test.approx_eq = [:(x1 == 1-x), :(x2 == 1.0*max(0, min(x,x1)) + 1.0*max(0, min(x,x1)))]
    oA_test.ranges = Dict(:x => [-1., 2.], :x1 => [-1., 2.], :x2 => [0., 0.5])
    overt_mip_model = OvertMIP(oA_test)
    JuMP.optimize!(overt_mip_model.model)
    for (x, x_var) in pairs(overt_mip_model.vars_dict)
       l, u = oA_test.ranges[x]
       @test l <= value(x_var)
       @test u >= value(x_var)
    end

    @objective(overt_mip_model.model, Min, overt_mip_model.vars_dict[:x2])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 0.

    @objective(overt_mip_model.model, Max, overt_mip_model.vars_dict[:x2])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 0.5

    @objective(overt_mip_model.model, Min, overt_mip_model.vars_dict[:x2] - overt_mip_model.vars_dict[:x])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ -2.

    @objective(overt_mip_model.model, Max, overt_mip_model.vars_dict[:x2] - overt_mip_model.vars_dict[:x])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 1.
end

@testset "overt_2_mip_test3" begin
    oA_test = OverApproximation()
    oA_test.approx_eq = [:(x1 == 2x), :(x2 == 1-x), :(x3 == 2.0*max(0, x1) + 1.0*max(0, x2) + -3.0*max(0, min(x1,x2)))]
    oA_test.ranges = Dict(:x => [-1., 2.], :x1 => [-2., 4.], :x2 => [-1., 2.], :x3 => [0., 8.])
    overt_mip_model = OvertMIP(oA_test)
    JuMP.optimize!(overt_mip_model.model)
    for (x, x_var) in pairs(overt_mip_model.vars_dict)
       l, u = oA_test.ranges[x]
       @test l <= value(x_var)
       @test u >= value(x_var)
    end

    @objective(overt_mip_model.model, Min, overt_mip_model.vars_dict[:x3])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 0.

    @objective(overt_mip_model.model, Max, overt_mip_model.vars_dict[:x3])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈  8.

    @objective(overt_mip_model.model, Min, overt_mip_model.vars_dict[:x3] - overt_mip_model.vars_dict[:x])
    JuMP.optimize!(overt_mip_model.model)
    @test 3.0*objective_value(overt_mip_model.model) ≈ -1.

    @objective(overt_mip_model.model, Max, overt_mip_model.vars_dict[:x3] - overt_mip_model.vars_dict[:x])
    JuMP.optimize!(overt_mip_model.model)
    @test objective_value(overt_mip_model.model) ≈ 6.
end


text:
# Matrix setup
#==========================================================================================#

"""
	influence_coefficient(r, horseshoe, normal, V_hat, symmetry)

Compute the influence coefficient of the velocity of a Horseshoe with trailing lines in a given direction ``\\hat V`` at a point ``r`` projected to a normal vector.
"""
influence_coefficient(horseshoe :: Horseshoe, r, normal, V_hat, finite_core = false) = dot(velocity(r, horseshoe, 1., V_hat, finite_core), normal)

"""
	influence_matrix(horseshoes, collocation_points, normals, V_hat)

Assemble the Aerodynamic Influence Coefficient (AIC) matrix given horseshoes, collocation points, associated normal vectors, and a unit vector representing the freestream.
"""
influence_matrix(horseshoes :: Vector{<: Horseshoe}, collocation_points, normals, V_hat, finite_core = false) = [ influence_coefficient(horsie_j, r_i, n_i, V_hat, finite_core) for (r_i, n_i) ∈ zip(collocation_points, normals), horsie_j ∈ horseshoes ]

"""
	boundary_condition(velocities, normals)

Compute the projection of a velocity vector with respect to normal vectors of panels. Corresponds to construction of the boundary condition for the RHS of the AIC system.
"""
boundary_condition(velocities, normals) = dot.(velocities, normals)

## Pre-allocated versions
#====================================================#

function matrix_assembly!(AIC, RHS, horseshoes :: Vector{<: Horseshoe}, collocation_points, normals, V, Ω)
    for i ∈ 1:length(collocation_points)
        for j ∈ 1:length(horseshoes)
            AIC[i,j] = dot(total_horseshoe_velocity(r1(collocation_points[i], horseshoes[j]), r2(collocation_points[i], horseshoes[j]), 1., -normalize(V)), normals[i])
        end
        RHS[i] = dot(V + Ω × collocation_points[i], normals[i])
    end
    nothing
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
module POMDPSimulators

using POMDPs
using Random
using ProgressMeter
using DataFrames
using POMDPPolicies
using POMDPModelTools
using BeliefUpdaters
using Distributed

import POMDPs: simulate, discount

export RolloutSimulator
include("rollout.jl")

export
    SimHistory,
    POMDPHistory,
    MDPHistory,
    AbstractPOMDPHistory,
    AbstractMDPHistory,
    HistoryIterator,
    eachstep,
    state_hist,
    action_hist,
    observation_hist,
    belief_hist,
    reward_hist,
    info_hist,
    ainfo_hist,
    uinfo_hist,
    exception,
    backtrace,
    undiscounted_reward,
    discounted_reward,
    n_steps,
    step_tuple
include("history.jl")

export sim
include("sim.jl")

export HistoryRecorder
include("history_recorder.jl")

export
    StepSimulator,
    stepthrough
include("stepthrough.jl")

export
    Sim,
    run,
    run_parallel,
    problem
include("parallel.jl")

end # module


text:
using SamplingDataFrame
using Test
using DataFrames

@testset "SamplingDataFrame.jl" begin
    # Write your tests here.
    #データセット作成
    x=1:10
    y1=2x
    y2=map(x->2x-1,x)
    label1=["a","a","a","b","b","b","b","c","c","c"]
    label2=["1","2","3","1","2","3","1","2","3","1"]
    df1=DataFrame(x=x,y1=y1,y2=y2,label1=label1,label2=label2)

    #行番号
    trainingdata,testdata=TestDataSetsNumbering(df1,PERCENTAGE=0.6)
    @test testdata==setdiff(1:size(df1,1),trainingdata)

    # データセット作成
    train,test=CreateTestDataSets(df1,PERCENTAGE=0.6)
    trainingdata=train[:,:x]
    testdata=setdiff(1:size(df1,1),trainingdata)
    @test test ==df1[testdata,:]
end


text:
module ClassPolynomials

import Nemo, GZip

DIR = joinpath(dirname(pathof(ClassPolynomials)), "..", "deps")

export ClassicalModularPolynomial, AtkinModularPolynomial, EtaModularPolynomial, WeberModularPolynomial

function ClassicalModularPolynomial(ell)
    parent, (X, Y) = Nemo.PolynomialRing(Nemo.ZZ, ["X", "Y"])
    ClassicalModularPolynomial(ell, X, Y)
end

function ClassicalModularPolynomial(ell, X, Y)
    _parse(joinpath("ModPolCls", "$ell.pol.gz")) do deg1, deg2, c
        d1, d2 = parse(Int, deg1), parse(Int, deg2)
        mon = X^d1 * Y^d2
        if d1 != d2
            mon += X^d2 * Y^d1
        end
        parse(BigInt, c) * mon
    end
end

function EtaModularPolynomial(ell)
    parent, (X, Y) = Nemo.PolynomialRing(Nemo.ZZ, ["X", "Y"])
    EtaModularPolynomial(ell, X, Y)
end

function EtaModularPolynomial(ell, X, Y)
    _parse(joinpath("ModPolEta", "$ell.pol.gz")) do deg1, deg2, c
        parse(BigInt, c) * X^parse(Int, deg1) * Y^parse(Int, deg2)
    end
end

function AtkinModularPolynomial(ell)
    parent, (X, Y) = Nemo.PolynomialRing(Nemo.ZZ, ["X", "Y"])
    AtkinModularPolynomial(ell, X, Y)
end

function AtkinModularPolynomial(ell, X, Y)
    _parse(joinpath("ModPolAtk", "$ell.pol.gz")) do deg1, deg2, c
        parse(BigInt, c) * X^parse(Int, deg1) * Y^parse(Int, deg2)
    end
end

function WeberModularPolynomial(ell)
    parent, (X, Y) = Nemo.PolynomialRing(Nemo.ZZ, ["X", "Y"])
    WeberModularPolynomial(ell, X, Y)
end

function WeberModularPolynomial(ell, X, Y)
    _parse(joinpath("ModPolWeb", "$ell.pol.gz")) do deg1, deg2, c
        d1, d2 = parse(Int, deg1), parse(Int, deg2)
        mon = X^d1 * Y^d2
        if d1 != d2
            mon += X^d2 * Y^d1
        end
        parse(BigInt, c) * mon
    end
end

function _parse(monomial, file)
    data = GZip.open(joinpath(DIR, file))
    res = 0
    for line in eachline(data)
        res += monomial(split(line)...)
    end
    close(data)
    return res
end

end # module


text:
#= ForwardDiff is removed as a dependency in this version.
import MechanicalSketch
import MechanicalSketch: PALETTE, O, HE, WI, EM, FS, PT, ∙, settext
import MechanicalSketch: @import_expand, set_scale_sketch
import MechanicalSketch: generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: clamped_velocity_matrix
import MechanicalSketch: place_image
import MechanicalSketch: LicSceneOpts
import MechanicalSketch: ColorSchemes, color_with_lumin
import MechanicalSketch: BinLegend, setfont, FS, fontsize, background, sethue, dimension_aligned, line
import MechanicalSketch: Movie, Scene, animate, color_matrix_current, draw_legend, @layer
import MechanicalSketch: fontsize, fontface, setfont, color_with_lumin
import MechanicalSketch: convolution_matrix, ColorSchemes
import ColorSchemes:     isoluminant_cgo_80_c38_n256

let

if !@isdefined m²
    @import_expand(~m, s, °)
end
include("test_functions_33.jl")

# Scaling and placement
physwidth = 10.0m
physheight = 4.0m
framedheight = physheight * 1.2
totheight = framedheight * 3
set_scale_sketch(totheight, HE)
totwidth = totheight * WI / HE
Δx = totwidth / 4
Δy = framedheight
Oadj = O + (0, EM /3)
cent1 = Oadj + (-Δx, Δy)
cent2 = Oadj + ( Δx, Δy)
cent3 = Oadj + (-Δx, 0.0m)
cent4 = Oadj + ( Δx, 0.0m)
cent5 = Oadj + (-Δx, -Δy)

# Reused velocity field from earlier tests
max_velocity = 0.5m/s
velocity_matrix = clamped_velocity_matrix(ϕ_33; physwidth, physheight, cutoff = max_velocity);
speedmatrix = hypot.(velocity_matrix)
c1 = convolution_matrix(velocity_matrix)
lum = 50
cols = isoluminant_cgo_80_c38_n256.colors .|> co -> color_with_lumin(co, lum)
legen = BinLegend(maxlegend = max_velocity, binwidth = 0.1m/s, colorscheme = cols, name = :Speed)

function scenepic(scene, framenumber)
    colmat = color_matrix_current(scene, framenumber)
    place_image(scene.opts.O, colmat)
end

function sc1(scene, framenumber)
    ulp, lrp = scenepic(scene, framenumber)
    str = "    1: Luminosity variation = 1 / 6"
    settext(str, ulp, markup = true)

end
function sc2(scene, framenumber)
    ulp, _ = scenepic(scene, framenumber)
    str = "    2: Luminosity variation = 2 / 6"
    settext(str, ulp, markup = true)
end

function sc3(scene, framenumber)
    ulp, _ = scenepic(scene, framenumber)
    str = "    3: Luminosity variation = 3 / 6"
    settext(str, ulp, markup = true)
end
function sc4(scene, framenumber)
    ulp, _ = scenepic(scene, framenumber)
    str = "    4: Luminosity variation = 4 / 6"
    settext(str, ulp, markup = true)
end

function sc5(scene, framenumber)
    ulp, lrp = scenepic(scene, framenumber)
    str = "    5: Luminosity variation = 5 / 6"
    settext(str, ulp, markup = true)
    # Centre of this scene's rectangle
    O = scene.opts.O
    # Top left of legend
    legendpos = O + ((lrp[1] - ulp[1]) / 2 + EM, (ulp[2] - lrp[2]) / 2)
    @layer begin
        sethue("white")
        # Font for the 'toy' text interface
        fontface("Calibri")
        # 1 pt font size = 12/72 inch - by the book.
        # Letter spacing works differently here than in Word, so we adjust a little.
        fontsize(FS)
        draw_legend(legendpos, scene.opts.legend)
    end
end

function backdrop(scene, framenumber)
    setfont("Calibri", FS)
    fontsize(FS)
    O = scene.opts.O
    background(color_with_lumin(PALETTE[6], 30))
    sethue(PALETTE[3])
    Δt = scene.opts.cycle_duration
    frames_per_cycle = Δt ∙ scene.opts.framerate
    # Where we are in the repeating cycle, [0, 1 - frame_duration]
    normalized_time = framenumber / frames_per_cycle
    t = Δt ∙ normalized_time

    dimension_aligned(O + (-2.0m, physheight / 2), O + (-1.0m, physheight / 2),
        fromextension = (0, EM), toextension = (0, EM), offset = -EM)
    dimension_aligned(O + (-3.0m, physheight / 2), O + (-2.0m, physheight / 2),
        fromextension = (0, EM), toextension = (0, EM), offset = -EM)
    dimension_aligned(O + (-4.0m, physheight / 2), O + (-3.0m, physheight / 2),
        fromextension = (0, EM), toextension = (0, EM), offset = -EM)
    for i = 0:4
        x = t ∙ 0.5m/s + i ∙ 1.0m - 4.0m
        line(O + (x, 1.1 * physheight / 2), O + (x, physheight / 2), :stroke)
    end
end
movie = Movie(WI, HE, "Five flow fields with varying luminosity variation", 15:15)
Δt = 2.0s
framerate = 30/s
endframe = Δt ∙ framerate - 1
o1 = LicSceneOpts(Δt,  cent1; phase_magnitude_matrix = c1, data = speedmatrix, legend = legen, luminosity_variation = 1 / 6)
o2 = LicSceneOpts(Δt,  cent2; phase_magnitude_matrix = c1, data = speedmatrix, legend = legen, luminosity_variation = 2 / 6)
o3 = LicSceneOpts(Δt,  cent3; phase_magnitude_matrix = c1, data = speedmatrix, legend = legen, luminosity_variation = 3 / 6)
o4 = LicSceneOpts(Δt,  cent4; phase_magnitude_matrix = c1, data = speedmatrix, legend = legen, luminosity_variation = 4 / 6)
o5 = LicSceneOpts(Δt,  cent5; phase_magnitude_matrix = c1, data = speedmatrix, legend = legen, luminosity_variation = 5 / 6)

scenes = [  Scene(movie, backdrop, 0:endframe, optarg = LicSceneOpts(Δt, O)),
            Scene(movie, sc1, 0:endframe, optarg = o1),
            Scene(movie, sc2, 0:endframe, optarg = o2),
            Scene(movie, sc3, 0:endframe, optarg = o3),
            Scene(movie, sc4, 0:endframe, optarg = o4),
            Scene(movie, sc5, 0:endframe, optarg = o5)];

animate(movie, scenes;
    creategif = true,
    pathname = joinpath(@__DIR__, "test_34.png"),
    framerate = framerate * s)

end # Let
=#


text:
# here we compare various deconvolution options in terms of image quality
using IndexFunArrays, LinearAlgebra, Random, Noise, TestImages, FourierTools
using DeconvOptim, FFTW, Optim, LineSearches
using View5D, Plots

obj = 100f0 .* Float32.(testimage("resolution_test_512"))
# simulate a simple PSF
sz = size(obj);
R_max = sz[1] ./ 12.0;
psf = generate_psf(sz, R_max);

# simulate a perfect image
conv_img = DeconvOptim.conv(obj, psf);

# set a fixed point for the measured data quality
max_photons = 1000

Random.seed!(42)
measured = poisson(conv_img, max_photons);

opt_options = nothing

iterations = 100
function get_data(summary)
    return (summary["best_ncc_img"], summary["best_nvar_img"])
end

function show_ncc!(summary, title="", dpi=300)
        nccs = summary["nccs"]
        plt = plot!(nccs, label=title*" NCC")
        col = plt[1][end].plotattributes[:markercolor]
        vline!([summary["best_ncc_idx"]], line=:dash, color=col, label=title*"_best NCC", dpi=dpi)
        println("best ncc index is $(summary["best_ncc_idx"])")
        println("best ncc is $(summary["best_ncc"])")
        println("best ncc loss is $(summary["losses"][summary["best_ncc_idx"]])")
        println("final ncc is $(summary["nccs"][end])")
        println("final loss is $(summary["losses"][end])")
        xlabel!("iteration")
        ylabel!("normalized cross correlation")
end

function show_nvar!(summary, title="")
    nvars = summary["nvars"]
    nvars_norm = nvars ./ nvars[1]
    plt = plot!(nvars_norm, label=title*" NVAR")
    col = plt[1][end].plotattributes[:markercolor]
    vline!([summary["best_nvar_idx"]], line=:dash, color=col, label=title*"_best NVAR")
    println("best nvar index is $(summary["best_nvar_idx"])")
    println("best nvar is $(summary["best_nvar"])")
    println("best nvar loss is $(summary["losses"][summary["best_nvar_idx"]])")
    println("final nvar is $(summary["nvars"][end])")
    println("final loss is $(summary["losses"][end])")
    xlabel!("iteration")
    ylabel!("normalized variance")
end

function show_loss!(summary, addCurves="", lowest_loss=minimum(summary["losses"]))
        losses = summary["losses"]
        log_losses = log.(losses .- lowest_loss .+ 1)
        rel_losses = log_losses # .- maximum(log_losses)
        plot!(rel_losses[1:end-1], label=addCurves)
        xlabel!("iteration")
        ylabel!("log loss")
end

opt_options, noreg_summary = DeconvOptim.options_trace_deconv(obj, iterations, Non_negative());

res_noreg = deconvolution(measured, psf;
        regularizer=nothing, mapping=Non_negative(),
        opt_options=opt_options, debug_f=nothing)

opt_grad, grad_summary = DeconvOptim.options_trace_deconv(obj, iterations, Non_negative(), );
res_grad = deconvolution(measured, psf;
        regularizer=nothing,  mapping=Non_negative(),
        opt=GradientDescent(), opt_options=opt_grad, debug_f=nothing)

opt_gr, gr_summary = DeconvOptim.options_trace_deconv(obj, iterations, Non_negative());
res_gr = deconvolution(measured, psf;
        regularizer=DeconvOptim.GR(), λ=1e-3,
        mapping=Non_negative(), opt_options=opt_gr, debug_f=nothing)

opt_tv, tv_summary = DeconvOptim.options_trace_deconv(obj, iterations, Non_negative())
res_tv = deconvolution(measured, psf;
        regularizer=DeconvOptim.TV(), λ=1e-3,
        mapping=Non_negative(), opt_options=opt_tv, debug_f=nothing)

plot()
title!("Regularization")
do_display=false
show_ncc!(noreg_summary, "NoReg")
show_nvar!(noreg_summary, "NoReg")
show_ncc!(gr_summary, "GR")
show_nvar!(gr_summary, "GR")
show_ncc!(tv_summary, "TV")
show_nvar!(tv_summary, "TV")

plot()
title!("Optimization")
show_loss!(noreg_summary, "LBFGS", minimum(noreg_summary["losses"]))
show_loss!(grad_summary, "SteepestDecent", minimum(noreg_summary["losses"]))
show_loss!(tv_summary, "TV", minimum(noreg_summary["losses"]))
show_loss!(gr_summary, "GR", minimum(noreg_summary["losses"]))

best_ncc_img, best_nvar_img = get_data(tv_summary)
@vt obj
@vt measured
@vt best_ncc_img
@vt best_nvar_img
@vt res_noreg


text:
@testset "products" begin
  res = products()

  @test isa(res, DataFrame)
end

@testset "product" begin
  res = product("btc-usd")

  @test isa(res, DataFrame)
end


text:
using BenchmarkTools
using Printf

# Get the number of iterations from the command line
n, = size(ARGS)
if (n < 1)
    println("Usage: test_markov_chain.jl N")
    println("       ---> Please specify the number of iterations.")
    exit()
end

N = parse(Int, ARGS[1])

"""
    function markov_optimized(x,N)

Operate the Markov chain on two identical inputs x a total of N times.
"""
function markov_optimized(x, N)
    x1, x2 = x
    f(x1, x2) = exp(sin(x1*5) - x1^2 - x2^2)
    p = f(x1, x2)
    for n = 1:N
        y1 = x1 + 0.01*randn(x1)
        y2 = x2 + 0.01*randn(x2)
        q = f(y1, y2)
        if p*rand() < q
            x1 = y1
            x2 = y2
            p = q
        end
    end
end

"""
    function markov(x, y, N)

Operate the Markov chain on two inputs x, y a total of N times.
"""
function markov(x, y, N)
    f(x, y) = exp(sin(x*5) - x^2 - y^2)
    p = f(x, y)
    for n = 1:N
        x2 = x + 0.01*randn()
        y2 = y + 0.01*randn()
        p2 = f(x2, y2)
        if rand() < p2/p
            x = x2
            y = y2
            p = p2
        end
    end
    x, y
end

# trigger JIT
markov(0.0, 0.0, 10)

println("--------------------------")
println(@sprintf "Markov Chain calculations: %d" N)
println("--------------------------")

@btime markov(0.0, 0.0, N)

#println("--------------------------")
#println(@sprintf "Markov Chain calculations - Opt: %d" N)
#println("--------------------------")
#
@btime markov_optimized([0.0 0.0],N)