text:
# FrequencyDomainAnalysis.jl
# Tools for spectral density estimation and analysis of phase relationships
# between sets of signals.

# Copyright (C) 2013   Simon Kornblith

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.

# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

# import Base.accumulate

# export PowerSpectrum, PowerSpectrumVariance, CrossSpectrum, Coherence, Coherency, PLV, PPC, PLI,
#        PLI2Unbiased, WPLI, WPLI2Debiased, ShiftPredictor, allpairs

# Get all pairs of channels
function allpairs(n)
    pairs = Array{Int}(2, binomial(n, 2))
    k = 0
    for i = 1:n-1, j = i+1:n
        k += 1
        pairs[1, k] = i
        pairs[2, k] = j
    end
    pairs
end

# Statistics computed on transformed data
abstract type TransformStatistic{T<:Real} end
abstract type PairwiseTransformStatistic{T<:Real} <: TransformStatistic{T} end

# Generate a new PairwiseTransformStatistic, including constructors
macro pairwisestat(name, xtype)
    esc(quote
        type $name{T<:Real} <: PairwiseTransformStatistic{T}
            pairs::Array{Int, 2}
            x::$xtype
            nepochs::Vector{Int}
            tmp::BitVector
            $name() = new()
            $name(pairs::Array{Int, 2}) = new(pairs)
        end
        $name() = $name{Float64}()
    end)
end

# Most PairwiseTransformStatistics will initialize their fields the same way
function init{T}(s::PairwiseTransformStatistic{T}, nout, nchannels, ntapers)
    if !isdefined(s, :pairs); s.pairs = allpairs(nchannels); end
    # s.x = zeros(eltype(fieldtype(s, :x)), nout, size(s.pairs, 2))
    s.x = zeros(Float64, nout, size(s.pairs, 2))

    s.nepochs = zeros(Int, size(s.pairs, 2))
    s.tmp = falses(nchannels)
end

# Create accumulate function that loops over pairs of channels, performing some transform for each
macro accumulatebypair(stat, arr, freqindex, pairindex, ch1ft, ch2ft, code)
    quote
        # Split out the two FFTs to allow efficient computation in shift predictor case
        function $(esc(:accumulatepairs))(s::$stat, fftout1, fftout2, itaper)
            $arr = s.x
            tmp = s.tmp
            pairs = s.pairs
            nepochs = s.nepochs
            @inbounds begin
                # Find channels with NaNs
                for i = 1:size(fftout1, 2)
                    good = true
                    for j = 1:size(fftout1, 1)
                        good &= !isnan(real(fftout1[j, i])) & !isnan(real(fftout2[j, i]))
                    end
                    tmp[i] = good
                end

                for $pairindex = 1:size(pairs, 2)
                    ch1 = pairs[1, $pairindex]
                    ch2 = pairs[2, $pairindex]

                    # Skip channels with NaNs
                    if !tmp[ch1] || !tmp[ch2] continue end
                    nepochs[$pairindex] += 1

                    for $freqindex = 1:size(fftout1, 1)
                        $ch1ft = fftout1[$freqindex, ch1]
                        $ch2ft = fftout2[$freqindex, ch2]
                        $code
                    end
                end
            end
        end
    end
end

accum(s::PairwiseTransformStatistic, fftout, itaper) =
    accumulatepairs(s, fftout, fftout, itaper)

#
# Power spectrum
#
type PowerSpectrum{T<:Real} <: TransformStatistic{T}
    x::Array{T, 2}
    nepochs::Vector{Int}
    PowerSpectrum() = new()
end
PowerSpectrum() = PowerSpectrum{Float64}()
function init{T}(s::PowerSpectrum{T}, nout, nchannels, ntapers)
    s.x = zeros(T, nout, nchannels)
    s.nepochs = zeros(Int, nchannels)
end
function accum(s::PowerSpectrum, fftout, itaper)
    A = s.x
    @inbounds for i = 1:size(fftout, 2)
        # Skip channels with NaNs
        good = true
        for j = 1:size(fftout, 1)
            good &= !isnan(real(fftout[j, i]))
        end
        if !good continue end

        s.nepochs[i] += 1
        for j = 1:size(fftout, 1)
            A[j, i] += abs2(fftout[j, i])
        end
    end
end
finish(s::PowerSpectrum) = scale!(s.x, 1./s.nepochs)

#
# Variance of power spectrum across trials
#
type PowerSpectrumVariance{T<:Real} <: TransformStatistic{T}
    x::Array{T, 3}
    ntrials::Vector{Int}
    ntapers::Int
    PowerSpectrumVariance() = new()
end
PowerSpectrumVariance() = PowerSpectrum{Float64}()
function init{T}(s::PowerSpectrumVariance{T}, nout, nchannels, ntapers)
    s.x = zeros(T, 3, nout, nchannels)
    s.ntrials = zeros(Int, nchannels)
    s.ntapers = ntapers
end
function accum(s::PowerSpectrumVariance, fftout, itaper)
    @inbounds begin
        A = s.x
        for i = 1:size(fftout, 2), j = 1:size(fftout, 1)
            A[1, j, i] += abs2(fftout[j, i])
        end
        if itaper == s.ntapers
            for i = 1:size(A, 3)
                # Skip channels with NaNs
                good = true
                for j = 1:size(A, 2)
                    good &= !isnan(real(A[1, j, i]))
                end
                if !good continue end

                # http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
                n = (s.ntrials[i] += 1) # n = n + 1
                for j = 1:size(A, 2)
                    x = A[1, j, i]
                    mean = A[2, j, i]
                    delta = x - mean
                    mean = mean + delta/n
                    A[3, j, i] += delta*(x - mean) # M2 = M2 + delta*(x - mean)
                    A[2, j, i] = mean
                end
            end
        end
    end
end
finish(s::PowerSpectrumVariance) = scale!(squeeze(s.x[3, :, :], 1), 1./(s.ntrials - 1))

#
# Cross spectrum
#
@pairwisestat CrossSpectrum Matrix{Complex{T}}
@accumulatebypair CrossSpectrum A j i x y begin
    A[j, i] += conj(x)*y
end
finish(s::CrossSpectrum) = scale!(s.x, 1./s.nepochs)

#
# Coherency and coherence
#
for sym in (:Coherency, :Coherence)
    @eval begin
        type $sym{T<:Real} <: PairwiseTransformStatistic{T}
            pairs::Array{Int, 2}
            psd::PowerSpectrum{T}
            xspec::CrossSpectrum{T}
            $sym() = new()
            $sym(pairs::Array{Int, 2}) = new(pairs)
        end
        $sym() = $sym{Float64}()
    end
end
function init{T}(s::Union{Coherency{T}, Coherence{T}}, nout, nchannels, ntapers)
    if !isdefined(s, :pairs); s.pairs = allpairs(nchannels); end
    s.psd = PowerSpectrum{T}()
    s.xspec = CrossSpectrum{T}(s.pairs)
    init(s.psd, nout, nchannels, ntapers)
    init(s.xspec, nout, nchannels, ntapers)
end
function accumulatepairs(s::Union{Coherency, Coherence}, fftout1, fftout2, itaper)
    accum(s.psd, fftout1, itaper)
    accumulatepairs(s.xspec, fftout1, fftout2, itaper)
end
function finish(s::Coherency)
    psd = finish(s.psd)
    xspec = finish(s.xspec)
    pairs = s.pairs
    for i = 1:size(pairs, 2)
        ch1 = pairs[1, i]
        ch2 = pairs[2, i]
        for j = 1:size(xspec, 1)
            xspec[j, i] = xspec[j, i]/sqrt(psd[j, ch1]*psd[j, ch2])
        end
    end
    xspec
end
function finish{T}(s::Coherence{T})
    psd = finish(s.psd)
    xspec = finish(s.xspec)
    out = zeros(T, size(xspec, 1), size(xspec, 2))
    pairs = s.pairs
    for i = 1:size(pairs, 2)
        ch1 = pairs[1, i]
        ch2 = pairs[2, i]
        for j = 1:size(xspec, 1)
            out[j, i] = abs(xspec[j, i])/sqrt(psd[j, ch1]*psd[j, ch2])
        end
    end
    out
end

#
# Phase locking value and pairwise phase consistency
#
# For PLV, see Lachaux, J.-P., Rodriguez, E., Martinerie, J., & Varela,
# F. J. (1999). Measuring phase synchrony in brain signals. Human Brain
# Mapping, 8(4), 194–208.
# doi:10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C
#
# For PPC, see Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P.,
# & Pennartz, C. M. A. (2010). The pairwise phase consistency: A
# bias-free measure of rhythmic neuronal synchronization. NeuroImage,
# 51(1), 112–122. doi:10.1016/j.neuroimage.2010.01.073
@pairwisestat PLV Matrix{Complex{T}}
@pairwisestat PPC Matrix{Complex{T}}
@accumulatebypair Union{PLV, PPC} A j i x y begin
    # Add phase difference between pair as a unit vector
    z = conj(x)*y
    A[j, i] += z/abs(z)
    # Faster, but less precise
    #A[j, i] += z*(1/sqrt(abs2(real(z))+abs2(imag(z))))
end
finish(s::PLV) = scale!(abs(s.x), 1./s.nepochs)
function finish{T}(s::PPC{T})
    out = zeros(T, size(s.x, 1), size(s.x, 2))
    for i = 1:size(s.x, 2)
        nepochs = s.nepochs[i]
        for j = 1:size(s.x, 1)
            # This is equivalent to the formulation in Vinck et al. (2010), since
            # 2*sum(unique pairs) = sum(trials)^2-n.
            out[j, i] = (abs2(s.x[j, i])-nepochs)/(nepochs*(nepochs-1))
        end
    end
    out
end

#
# Phase lag index and unbiased PLI^2
#
# For PLI, see Stam, C. J., Nolte, G., & Daffertshofer, A. (2007).
# Phase lag index: Assessment of functional connectivity from multi
# channel EEG and MEG with diminished bias from common sources.
# Human Brain Mapping, 28(11), 1178–1193. doi:10.1002/hbm.20346
#
# For unbiased PLI^2, see Vinck, M., Oostenveld, R., van
# Wingerden, M., Battaglia, F., & Pennartz, C. M. A. (2011). An
# improved index of phase-synchronization for electrophysiological data
# in the presence of volume-conduction, noise and sample-size bias.
# NeuroImage, 55(4), 1548–1565. doi:10.1016/j.neuroimage.2011.01.055
@pairwisestat PLI Matrix{Int}
@pairwisestat PLI2Unbiased Matrix{Int}
@accumulatebypair Union{PLI, PLI2Unbiased} A j i x y begin
    z = imag(conj(x)*y)
    if z != 0
        A[j, i] += 2*(z > 0)-1
    end
end
function finish{T}(s::PLI{T})
    out = zeros(T, size(s.x, 1), size(s.x, 2))
    for i = 1:size(s.x, 2)
        nepochs = s.nepochs[i]
        for j = 1:size(s.x, 1)
            out[j. i] = abs(s.x[j, i])/nepochs
        end
    end
    out
end
function finish{T}(s::PLI2Unbiased{T})
    out = zeros(T, size(s.x, 1), size(s.x, 2))
    for i = 1:size(s.x, 2)
        nepochs = s.nepochs[i]
        for j = 1:size(s.x, 1)
            out[j, i] = (nepochs * abs2(s.x[j, i]/nepochs) - 1)/(nepochs - 1)
        end
    end
    out
end

#
# Weighted phase lag index
#
# See Vinck et al. (2011) as above.
@pairwisestat WPLI Array{T,3}
# We need 2 fields per freq/channel in s.x
function init{T}(s::WPLI{T}, nout, nchannels, ntapers)
    if !isdefined(s, :pairs); s.pairs = allpairs(nchannels); end
    s.x = zeros(eltype(fieldtype(s, :x)), 2, nout, size(s.pairs, 2))
    s.nepochs = zeros(Int, size(s.pairs, 2))
    s.tmp = falses(nchannels)
end
@accumulatebypair WPLI A j i x y begin
    z = imag(conj(x)*y)
    A[1, j, i] += z
    A[2, j, i] += abs(z)
end
function finish{T}(s::WPLI{T}, nepochs)
    out = zeros(T, size(s.x, 2), size(s.x, 3))
    for i = 1:size(out, 2), j = 1:size(out, 1)
        out[j, i] = abs(s.x[1, j, i])/s.x[2, j, i]
    end
    out
end

#
# Debiased (i.e. still somewhat biased) WPLI^2
#
# See Vinck et al. (2011) as above.
@pairwisestat WPLI2Debiased Array{T,3}
# We need 3 fields per freq/channel in s.x
function init{T}(s::WPLI2Debiased{T}, nout, nchannels, ntapers)
    if !isdefined(s, :pairs); s.pairs = allpairs(nchannels); end
    s.x = zeros(eltype(fieldtype(s, :x)), 3, nout, size(s.pairs, 2))
    s.nepochs = zeros(Int, size(s.pairs, 2))
    s.tmp = falses(nchannels)
end
@accumulatebypair WPLI2Debiased A j i x y begin
    z = imag(conj(x)*y)
    A[1, j, i] += z
    A[2, j, i] += abs(z)
    A[3, j, i] += abs2(z)
end
function finish{T}(s::WPLI2Debiased{T})
    out = zeros(T, size(s.x, 2), size(s.x, 3))
    for i = 1:size(out, 2), j = 1:size(out, 1)
        imcsd = s.x[1, j, i]
        absimcsd = s.x[2, j, i]
        sqimcsd = s.x[3, j, i]
        out[j, i] = (abs2(imcsd) - sqimcsd)/(abs2(absimcsd) - sqimcsd)
    end
    out
end

#
# Shift predictors
#
type ShiftPredictor{T<:Real,S<:PairwiseTransformStatistic} <: TransformStatistic{T}
    stat::S
    first::Array{Complex{T}, 3}
    previous::Array{Complex{T}, 3}
    isfirst::Bool
end
ShiftPredictor{T}(stat::PairwiseTransformStatistic{T}) =
    ShiftPredictor{T,typeof(stat)}(stat, Array{Complex{T}}(0, 0, 0), Array{Complex{T}}(0, 0, 0), true)

function init{T}(s::ShiftPredictor{T}, nout, nchannels, ntapers)
    s.first = Array{Complex{T}}(nout, nchannels, ntapers)
    s.previous = Array{Complex{T}}(nout, nchannels, ntapers)
    s.isfirst = true
    init(s.stat, nout, nchannels, ntapers)
end
function accum(s::ShiftPredictor, fftout, itaper)
    first = pointer_to_array(pointer(s.first, size(s.first, 1)*size(s.first, 2)*(itaper - 1)+1),
                             (size(s.first, 1), size(s.first, 2)), false)
    previous = pointer_to_array(pointer(s.previous, size(s.previous, 1)*size(s.previous, 2)*(itaper - 1)+1),
                                (size(s.previous, 1), size(s.previous, 2)), false)
    if s.isfirst
        copy!(first, fftout)
        s.isfirst = itaper != size(s.first, 3)
    else
        accumulatepairs(s.stat, previous, fftout, itaper)
    end
    copy!(previous, fftout)
end
function finish(s::ShiftPredictor)
    if !s.isfirst
        for itaper = 1:size(s.first, 3)
            first = pointer_to_array(pointer(s.first, size(s.first, 1)*size(s.first, 2)*(itaper - 1)+1),
                                     (size(s.first, 1), size(s.first, 2)), false)
            previous = pointer_to_array(pointer(s.previous, size(s.previous, 1)*size(s.previous, 2)*(itaper - 1)+1),
                                        (size(s.previous, 1), size(s.previous, 2)), false)
            accumulatepairs(s.stat, previous, first, itaper)
        end
    end
    finish(s.stat)
end


text:
# tests information stored as choice point info

#
# sequence choice points
#

@generator CPRepsLiteralMinMaxGen begin # prefix CP (for Choice Point info) to avoid type name clashes
    start() = reps(a,2,4)
    a() = gensym()
end

@generator CPRepsNonLiteralMinGen begin
    start() = reps(a,1+1,4)
    a() = gensym()
end

@generator CPRepsNonLiteralMaxGen begin
    start() = reps(a,2,2+2)
    a() = gensym()
end

@generator CPRepsNonLiteralMinMaxGen begin
    start() = reps(a,1+1,2+2)
    a() = gensym()
end

@generator CPRepsNoMaxGen begin
    start() = reps(a,2)
    a() = gensym()
end

@generator CPMultGen begin
    start() = mult(a)
    a() = gensym()
end

@generator CPPlusGen begin
    start() = plus(a)
    a() = gensym()
end

@generator CPRuleGen begin
    start() = x()
    x() = 'a'
    x() = 'b'
    x() = 'c'
end

@generator CPChooseIntLiteralMinMaxGen begin
    start() = choose(Int,2,4)
end

@generator CPChooseIntNonLiteralMinGen begin
    start() = choose(Int,1+1,4)
end

@generator CPChooseIntNonLiteralMaxGen begin
    start() = choose(Int,2,2+2)
end

@generator CPChooseIntNonLiteralMinMaxGen begin
    start() = choose(Int,1+1,2+2)
end

@generator CPChooseIntNoMaxGen begin
    start() = choose(Int,2)
end

@generator CPChooseIntNoMinMaxGen begin
    start() = choose(Int)
end

@generator CPChooseFloat64LiteralMinMaxGen begin
    start() = choose(Float64,2.1,4.2)
end

@generator CPChooseBoolGen begin
    start() = choose(Bool)
end

@generator CPMultipleCPGen begin
    start() = mult(a)
    a() = choose(Int,2,4)
    a() = 'a'
end

@generator CPBoolGen begin
    start() = choose(Bool)
end

@generator CPIntGen begin
    start() = choose(Int,5,9)
end

@generator CPMainGen(boolGen, intGen) begin
    start() = map(i->choose(boolGen), 1:choose(intGen))
end

@testset "choice point info" begin

	@testset "sequence choice point info" begin

		@testset "reps with literal min and max" begin
		    gn = CPRepsLiteralMinMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test cp1[:min] == 2
		    @test cp1[:max] == 4
		end

		@testset "reps with non-literal min" begin
		    gn = CPRepsNonLiteralMinGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test !haskey(cp1,:min)
		    @test cp1[:max] == 4
		end

		@testset "reps with non-literal max" begin
		    gn = CPRepsNonLiteralMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test cp1[:min] == 2
		    @test !haskey(cp1,:max)
		end

		@testset "reps with non-literal min and max" begin
		    gn = CPRepsNonLiteralMinMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test !haskey(cp1,:min)
		    @test !haskey(cp1,:max)
		end

		@testset "reps with no max" begin
		    gn = CPRepsNoMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test cp1[:min] == 2
		    @test cp1[:max] == typemax(Int)
		end

		@testset "mult" begin
		    gn = CPMultGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test cp1[:min] == 0
		    @test cp1[:max] == typemax(Int)
		end

		@testset "plus" begin
		    gn = CPPlusGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :sequence
		    @test cp1[:min] == 1
		    @test cp1[:max] == typemax(Int)
		end

	end

	@testset "rule choice point info" begin

		@testset "rule choice" begin
		    gn = CPRuleGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :rule
		    @test cp1[:min] == 1
		    @test cp1[:max] == 3
		    @test cp1[:rulename] == :x
		end

	end

	@testset "value choice point info" begin

		@testset "choose Int with literal min and max" begin
		    gn = CPChooseIntLiteralMinMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Int
		    @test cp1[:min] == 2
		    @test cp1[:max] == 4
		end

		@testset "choose Int with non literal min" begin
		    gn = CPChooseIntNonLiteralMinGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Int
		    @test !haskey(cp1,:min)
		    @test cp1[:max] == 4
		end

		@testset "choose Int with non literal max" begin
		    gn = CPChooseIntNonLiteralMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Int
		    @test cp1[:min] == 2
		    @test !haskey(cp1,:max)
		end

		@testset "choose Int with non literal min and max" begin
		    gn = CPChooseIntNonLiteralMinMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Int
		    @test !haskey(cp1,:min)
		    @test !haskey(cp1,:max)
		end

		@testset "choose Int with no max" begin
		    gn = CPChooseIntNoMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Int
		    @test cp1[:min] == 2
		    @test cp1[:max] == typemax(Int)
		end

		@testset "choose Int with no min nor max" begin
		    gn = CPChooseIntNoMinMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Int
		    @test cp1[:min] == typemin(Int)
		    @test cp1[:max] == typemax(Int)
		end

		@testset "choose Float64 with literal min and max" begin
		    gn = CPChooseFloat64LiteralMinMaxGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Float64
		    @test cp1[:min] == 2.1
		    @test cp1[:max] == 4.2
		end

		@testset "choose Bool" begin
		    gn = CPChooseBoolGen()
		    cpi = gn.choicepointinfo
		    cpids = collect(keys(cpi))
		    @test length(cpi) == 1
		    cp1 = cpi[cpids[1]]
		    @test cp1[:type] == :value
		    @test cp1[:datatype] == Bool
		    @test cp1[:min] == false
		    @test cp1[:max] == true
		end

	end

	@testset "multiple choice points" begin

		@testset "generator sequence, rule, and value choice points" begin
		    gn = CPMultipleCPGen()
		    cpi = gn.choicepointinfo
		    @test length(cpi) == 3
		    cptypes = [cpdict[:type] for cpdict in values(cpi)]
		    @test :sequence in cptypes
		    @test :rule in cptypes
		    @test :value in cptypes
		end

	end

	@testset "choicepointinfo function" begin

		@testset "returns info for the generator itself when no subgenerators" begin
		    gn = CPMultipleCPGen()
		    cpi = gn.choicepointinfo
		    @test choicepointinfo(gn) == cpi
		end

		@testset "returns info for the generator itself only even where there are subgenerators" begin
		    bgn = CPBoolGen()
		    ign = CPIntGen()
		    gn = CPMainGen(bgn,ign)
		    bcpi = bgn.choicepointinfo
		    icpi = ign.choicepointinfo
		    cpi = gn.choicepointinfo
		    @test choicepointinfo(gn) == cpi
		end

	end

end

ext: jl


text:
"""
  MshFileReader
function for reading gmsh's mesh format files http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
    Input:
        - mshFilename: the name of the .msh file to be read
    Output:
        - nodesMat: matrix with 4 columns: [x y z physicalTag]
        - conecMat: matrix with 5 columns: [ n1 n2 n3 n4 physicalTag ]
            for elements with less than four nodes 0 is used as node index.
        - physicalNames: cell with strings of physical names.

    Assumptions:
        - physical names are saved as strings
        - maximum of one physical tag per entity
        - maximum number of nodes per element: 4 (linear tetrahedron)
"""
function MshFileReader( mshFilename; verbosityBoolean::Bool = true )

    mshFileLines = open( readlines, mshFilename, "r")

    # ----------------------------------------------------
    # verifies if the msh version is 4.1
    currLine = 2
    cmp( mshFileLines[currLine][1:3], "4.1" ) != 0 && error("version of file error")

    # ----------------------------------------------------
    # reads physical names
    currLine += 2
    if cmp( mshFileLines[ currLine ][2:14], "PhysicalNames" ) == 0
        verbosityBoolean && println("Physical names found. Reading ... ")

        currLine += 1
        numPhysicalNames = parse( Int32, mshFileLines[ currLine ] )

        numsPhysProp   = zeros( Int32, numPhysicalNames )
        physicalNames  = Vector{String}(undef, numPhysicalNames )

        for i in 1:numPhysicalNames
            currLine += 1
            _, aux, aux2 = split( mshFileLines[ currLine ] )
            numsPhysProp[i]  = parse( Int32, aux ) # convert string to integer
            physicalNames[i] = aux2[2:end-1]       # remove quotes from start and end

            verbosityBoolean && println("  property ", i, ": ", physicalNames[i] )
        end

    end

    # ----------------------------------------------------
    currLine += 2
    if cmp( mshFileLines[ currLine ][2:7], "Entiti" ) == 0
        verbosityBoolean && println("Entities found. Reading ... ")

        currLine += 1
        entNumsPerDim = parse.( Int32, split( mshFileLines[ currLine ]) )

        #    matsPhysicalPropsPerEntity = cell(4) ;
        vecsPhysicalPropsPerEntity = Vector{Matrix{Int32}}(undef,4)
        [ vecsPhysicalPropsPerEntity[i] = zeros(entNumsPerDim[i],2) for i in 1:4 ]

        for indDim in 1:4
             colNumTags  = 1+3+3*(indDim>1)+1 # this is the column number of the entry with the number of physical properties (or tags)
             colTags     = 1+3+3*(indDim>1)+2 # this is the column number of the entry with the first physical property

            if entNumsPerDim[indDim] > 0
                for i in 1:entNumsPerDim[indDim]
                    currLine += 1
                    aux = parse.( Int32, split( mshFileLines[currLine] ) )
                    if aux[colNumTags] > 0 # if there are physicalTags
                        vecsPhysicalPropsPerEntity[indDim][i,:] = [ aux[1] aux[colTags] ] # assigns [element_tag physical_tag ]
                    else
                        vecsPhysicalPropsPerEntity[indDim][i,:] = [ aux[1] 0            ]
                    end
                end
            end
        end
        verbosityBoolean && print( "\n vec phys: ", vecsPhysicalPropsPerEntity,"\n\n")
    end # if - entities block

    # ----------------------------------------------------
    currLine += 2
    if cmp( mshFileLines[ currLine ][2:5], "Node" ) == 0
        currLine += 1
        numEntBlocks, numNodes = parse.( Int32, split( mshFileLines[ currLine ] ) )[1:2]

        nodesCoordMat = zeros( numNodes, 3 )

        # vector used to store the number of physical property of each node
        #   positive is assigned, zero if not.
        nodePhysNums  = zeros(Int32, numNodes )

        for block in 1:numEntBlocks
            currLine += 1
            aux = parse.( Int32, split( mshFileLines[ currLine ] ) )
            entityDimen   = aux[1]
            entityTag     = aux[2]
            numNodesInEnt = aux[4]
            if numNodesInEnt > 0 # if there are nodes in the block
                currLine += 1
                nodesTags = parse.( Int32, split( mshFileLines[ currLine ] ) )

                for indnode in nodesTags
                    currLine += 1
                    nodesCoordMat[ indnode, :] = parse.( Float64, split( mshFileLines[ currLine ] ) )
                end

                # only saves physical tags for nodes defined as nodes (no inheritance)
                if entityDimen == 0 # node created manualy
                    dimenEntTags = vecsPhysicalPropsPerEntity[entityDimen+1][:,1]
                    indEnt = findall( x->x==entityTag, dimenEntTags )
                    nodePhysNums[ nodesTags] = vecsPhysicalPropsPerEntity[entityDimen+1][indEnt ,2]
                end

            end
        end
    end # if - nodes block

    # ----------------------------------------------------
    currLine += 2
    if cmp( mshFileLines[ currLine ][2:5], "Elem" ) == 0
        verbosityBoolean && println("Elements found. Reading ... ")
        currLine += 1
        numEntBlocks, numElems, _, _ = parse.( Int32, split( mshFileLines[ currLine ] ) )

        # vector used to store the number of physical property of each element
        #   positive is assigned, zero if not.
        elemPhysNums = zeros(Int32, numElems )

        connectivity = Vector{Vector{Int64}}(undef, numElems)

        for block in 1:numEntBlocks
            currLine += 1
            dimOfBlock, entityTag, _, numElemInBlock = parse.( Int32, split( mshFileLines[ currLine ] ) )

print("\n block:",block)

            elemInds = []
            if numElemInBlock > 0 # if there are elements in the block
                for i in 1:numElemInBlock

                    print("\n element:",i)

                    currLine += 1
                    aux = parse.( Int32, split( mshFileLines[ currLine ] ) )
                    print(aux)
                    connectivity[ aux[1] ] = aux[2:end]
                    push!(elemInds, aux[1] )
                end
            end
print("\n\nelem inds", elemInds)
            dimenEntTags = vecsPhysicalPropsPerEntity[ dimOfBlock+1 ][:,1]
            indEnt = findall( x->x==entityTag, dimenEntTags )
            elemPhysNums[ elemInds ] .= vecsPhysicalPropsPerEntity[ dimOfBlock+1 ][indEnt ,2]

        end
    end # if - elements block

    return nodesCoordMat, connectivity, physicalNames, elemPhysNums
end # function


text:
# function that returns the joint tree for an n-rotor vtol vehicle
# include("simulate.jl")
function gen_nRotor(n::Integer)
    # Inertial and mass properties from:
    # Hossain, M. Raju, and Nicholas Krouglicof. "Multi-body dynamics modeling & control of quadrotor helicopter using bond graph." Proceedings of the International Conference on Bond Graph Modeling and Simulation. 2016.

    mb = 0.30155; mp = 0.011721;

    Ib = zeros(3,3);
    Ib[1,1] = 6.4e-4; Ib[2,2] = 6.4e-4; Ib[3,3] = 1.22e-3;

    Ip = zeros(3,3);
    Ip[1,1] = 6.4061e-7; Ip[2,2] = 6.2732e-6; Ip[3,3] = 3.008e-5;

    InFrame = InertialFrameAsRB()
    QuadFrame = RigidBody(mb,Ib,2)
    x0Quad = [zeros(3);[1;zeros(3)];zeros(3);zeros(4)];
    initialiseRigidBody!(QuadFrame,x0Quad);

    Props = Vector{RigidBody}(undef,n);
    j = Vector{Joint}(undef, n+1);
    j[1] = Joint(InFrame,QuadFrame,zeros(3),zeros(3));

    th = 2*pi/n;
    l = 1.0;
    axis = [0.0 0.0 1.0][:];
    for i=1:n
        Props[i] = RigidBody(mp, Ip, 2+i);
        ang =(i-1)*th; # angle
        r0Prop_i = [l*cos(ang); l*sin(ang); 0.0];
        β0Prop_i = abs.([cos(ang/2); sin(ang/2).*[0.0;0.0;1.0]]);
        x0Prop_i = [r0Prop_i; β0Prop_i; zeros(7)];
        initialiseRigidBody!(Props[i],x0Prop_i);

        j[i+1] = Joint(QuadFrame, Props[i], r0Prop_i, zeros(3), type = "Revolute", axis = axis);
    end
    return j
end
## simulate quadcopter
# clearconsole();
# qRotor_j = gen_nRotor(4);
#
# tEnd = 0.1;
# tSpan = 0.01;
# g =([0.0,0.0,-9.81]); # Gravity Vector.
# tSim, solFinal = simulate(tEnd, tSpan, qRotor_j, g=g);
#
# solQuad = solFinal[1];
# solProp1 = solFinal[2];
# solProp2 = solFinal[3];
# solProp3 = solFinal[4];
# solProp4 = solFinal[5];

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
using Plots, BenchmarkTools
plotly()

Model = fooSdErgm1
T = 300
N=10
GBA = 1
GW = 2#N#N#1# N#N # 20
groupsInds = [Utilities.distributeAinVecN(Array{Int64}(1:GW),N), Utilities.distributeAinVecN(Array{Int64}(1:GBA),GW)]

UMdgp,dgpdegs = linSpacedPar(fooSdErgm1,N; NgroupsW = GW)
Bdgp = 0.9*ones(GBA)
Adgp = 0.1*ones(GBA)
Wdgp = UMdgp.*(1-Bdgp[groupsInds[2]])
scoreScal = "FISHER-EWMA"#""#
Mod = SdErgm1(ones(2,2),[Wdgp , Bdgp, Adgp],groupsInds,scoreScal)
Mod,Y_T,Fitness_T = sampl(Mod,T)
degs_T = squeeze(sum(Y_T,2),2)
#time = @elapsed  out = estimate(Mod)
#out
#println((T,N,GW ,round(time),out[1][2],out[1][3]))

Fitness_T = Float64.(Fitness_T)
corrMat=  squeeze(StatsBase.crosscor( Fitness_T,Fitness_T,[0]),1)
# ppar = plot(Fitness_T,title = "$(round.(squeeze(mean(Fitness_T,1),1),1))  $(round.(UMdgp,1))")
# pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
# plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))

#
p1 = contour(corrMat,fill=true,yflip = true, title = "Scaling = " *scoreScal )
    p2 = histogram(corrMat[.!Bool.(eye(corrMat))] )
    plot(p1,p2,layout= (2,1),size=(600,600),legend = :none)
    corrMat=  squeeze(StatsBase.crosscor( degs_T,degs_T,[0]),1)
    p1 = contour(corrMat,fill=true,yflip = true, title = " degs Corr with Scaling = " *scoreScal )
    p2 = histogram(corrMat[.!Bool.(eye(corrMat))] )
    plot(p1,p2,layout= (2,1),size=(600,600),legend = :none)    #

# ppar = plot(Fitness_T,title = "$(round.(squeeze(mean(Fitness_T,1),1),1))  $(round.(UMdgp,1))")
# pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
#     plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))
#
#     vresPar = [Mod.Par[1];Mod.Par[2];Mod.Par[3]]
#     Fitness_T_filt, tmp = score_driven_filter_or_dgp( Mod, vresPar)
#     ppar_filt = plot(Fitness_T_filt,title = "$(round.(squeeze(mean(Fitness_T_filt,1),1),1))  $(round.(UMdgp,1))")
#     plot(ppar,ppar_filt ,layout = (2,1),legend=:none,size = (1200,600))

## Test Single snapshot estimates
#@time Fitness_T_estSnap = estimateSnapSeq(Mod)
    # UMsingSnap = squeeze(mean(Fitness_T_estSnap,1),1)
    # UMdgp
    # ppar = plot(Fitness_T_estSnap,title = "$(round.(squeeze(mean(Fitness_T_estSnap,1),1),1))  $(round(UMdgp,1))")
    # pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
    # plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))

## Test Target Estimate

estOutSS = estimateSnapSeq(Mod)

#println(Mod.obsT[110,:])

estOutTarg = estimateTargeting(Mod)
Mod_estTarg = SdErgm1(Mod.obsT,estOutTarg[1],groupsInd,scoreScal)

vresParTarg = [Mod_estTarg.Par[1];Mod_estTarg.Par[2];Mod_estTarg.Par[3]]
Fitness_T_estTarg, tmp = score_driven_filter_or_dgp( Mod, vresParTarg)

ppar = plot(Fitness_T_estTarg,title = "$(round.(squeeze(mean(Fitness_T_estTarg,1),1),1))  $(round.(UMdgp,1))")
pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))

## Test estimate

estOut = estimate(Mod)
Mod_est = SdErgm1(Mod.obsT,estOut[1],groupsInd,scoreScal)

vresPar = [Mod_est.Par[1];Mod_est.Par[2];Mod_est.Par[3]]
Fitness_T_est, tmp = score_driven_filter_or_dgp( Mod, vresPar)

ppar = plot(Fitness_T_est,title = "$(round.(squeeze(mean(Fitness_T_est,1),1),1))  $(round.(UMdgp,1))")
pdegs = plot(degs_T,title = "$(dgpdegs)  $(round.(squeeze(mean(degs_T,1),1),1))")
plot(ppar,pdegs,layout = (2,1),legend=:none,size = (1200,600))

#

##
using JLD2

N_est = 300
Tvals =  [1000]
N_ind = length(Tvals)[1]

n  =1
N = 10 ;
GBAdgp =1; GWdgp = N  ;
GBAest =1; GWest = GWdgp; useStartVal = false

useStartVal? startValFlag = "T" : startValFlag = "F";# false #false #
# save
save_fold = "./data/estimatesTest/"
file_name = "EstTest_$(N)_$(GWdgp)_$(GBAdgp)_$(Tvals)_$(N_est)_$(useStartVal)_Miss_$(GWest)_$(GBAest)_rescaling.jld"
save_path = save_fold*file_name#
estDataT = @load(save_path,groupsIndsEst,groupsIndsDgp,allObs,
        simPar_1,simPar_2,simPar_3,estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)

vestPar = [estPar_1[1,n,:]; estPar_2[n];estPar_3[n]]
Fitness_Tfil ,like= score_driven_filter_or_dgp( Mod,vestPar)

    pparfil = plot(Fitness_Tfil)
    pparfil = plot(Fitness_Tfil)
    plot(ppar,pparfil,layout = (2,1),legend=:none,size = (1200,600))
##
vResGasPar = [Mod.Par[1];Mod.Par[2];Mod.Par[3];]
Fitness_Tfil ,like= score_driven_filter_or_dgp( Mod,vResGasPar)

    pparfil = plot(Fitness_Tfil)
    pparfil = plot(Fitness_Tfil)
    plot(ppar,pparfil,layout = (2,1),legend=:none,size = (1200,600))

##

#


text:
import Base.+, Base.*, Base.length, Base.product
using LinearAlgebra

# Abstract base class and definitions
abstract type Term end

function get_terms(term::Term)
    return tuple(term)
end

function get_real_coefficients(term::Term)
    return (Array{Float64}(undef, 0), Array{Float64}(undef, 0))
end

function get_complex_coefficients(term::Term)
    return (Array{Float64}(undef, 0), Array{Float64}(undef, 0), Array{Float64}(undef, 0), Array{Float64}(undef, 0))
end

function get_all_coefficients(term::Term)
    rc = get_real_coefficients(term)
    cc = get_complex_coefficients(term)
    return (rc[1], rc[2], cc[1], cc[2], cc[3], cc[4])
end

function get_value(term::Term, tau)
    coeffs = get_all_coefficients(term)
    t = abs.(tau)
    k = zeros(size(tau))
    for i in 1:length(coeffs[1])
        k .+= coeffs[1][i] .* exp.(-coeffs[2][i] .* t)
    end
    for i in 1:length(coeffs[3])
        k .+= (coeffs[3][i] .* cos.(coeffs[6][i] .* t) .+ coeffs[4][i] .* sin.(coeffs[6][i] .* t)) .* exp.(-coeffs[5][i] .* t)
    end
    return k
end

function get_psd(term::Term, omega)
    coeffs = get_all_coefficients(term)
    w2 = omega.^2
    w4 = w2.^2
    p = zeros(w2)
    for i in 1:length(coeffs[1])
        a = coeffs[1][i]
        c = coeffs[2][i]
        p = p + a*c ./ (c*c + w2)
    end
    for i in 1:length(coeffs[3])
        a = coeffs[3][i]
        b = coeffs[4][i]
        c = coeffs[5][i]
        d = coeffs[6][i]
        w02 = c*c+d*d
        p = p + ((a*c+b*d)*w02+(a*c-b*d).*w2) ./ (w4 + 2.0*(c*c-d*d).*w2+w02*w02)
    end
    return sqrt(2.0 / pi) .* p
end

function length(term::Term)
    return 0
end

function get_parameter_vector(term::Term)
    return Array{Float64}(undef, 0)
end

function set_parameter_vector!(term::Term, vector::Array)
end

# A sum of terms
mutable struct TermSum <: Term
    terms
end

function get_terms(term::TermSum)
    return term.terms
end

function get_all_coefficients(term_sum::TermSum)
    a_real = Array{Float64}(undef, 0)
    c_real = Array{Float64}(undef, 0)
    a_complex = Array{Float64}(undef, 0)
    b_complex = Array{Float64}(undef, 0)
    c_complex = Array{Float64}(undef, 0)
    d_complex = Array{Float64}(undef, 0)
    for term in term_sum.terms
        coeffs = get_all_coefficients(term)
#        a_real = cat(1, a_real, coeffs[1])
#        c_real = cat(1, c_real, coeffs[2])
#        a_complex = cat(1, a_complex, coeffs[3])
#        b_complex = cat(1, b_complex, coeffs[4])
#        c_complex = cat(1, c_complex, coeffs[5])
#        d_complex = cat(1, d_complex, coeffs[6])
        a_real = vcat(a_real, coeffs[1])
        c_real = vcat(c_real, coeffs[2])
        a_complex = vcat(a_complex, coeffs[3])
        b_complex = vcat(b_complex, coeffs[4])
        c_complex = vcat(c_complex, coeffs[5])
        d_complex = vcat(d_complex, coeffs[6])
    end
    return (a_real, c_real, a_complex, b_complex, c_complex, d_complex)
end

function +(t1::Term, t2::Term)
    return TermSum((get_terms(t1)..., get_terms(t2)...))
end

function length(term_sum::TermSum)
    return +(map(length, term_sum.terms)...)
end

function get_parameter_vector(term_sum::TermSum)
    return cat(1, map(get_parameter_vector, term_sum.terms)...)
end

function set_parameter_vector!(term_sum::TermSum, vector::Array)
    index = 1
    for term in term_sum.terms
        len = length(term)
        set_parameter_vector!(term, vector[index:index+len-1])
        index = index + len
    end
end

# A product of two terms
mutable struct TermProduct <: Term
    term1
    term2
end

function chain(a...)
    return (el for it in a for el in it)
end

function get_all_coefficients(term_sum::TermProduct)
    c1 = get_all_coefficients(term_sum.term1)
    nr1, nc1 = length(c1[1]), length(c1[3])
    c2 = get_all_coefficients(term_sum.term2)
    nr2, nc2 = length(c2[1]), length(c2[3])

    # First compute real terms
    nr = nr1 * nr2
    ar = Array{Float64}(nr)
    cr = Array{Float64}(nr)
    gen = product(zip(c1[1], c1[2]), zip(c2[1], c2[2]))
    for (i, ((aj, cj), (ak, ck))) in enumerate(gen)
        ar[i] = aj * ak
        cr[i] = cj + ck
    end

    # Then the complex terms
    nc = nr1 * nc2 + nc1 * nr2 + 2 * nc1 * nc2
    ac = Array{Float64}(nc)
    bc = Array{Float64}(nc)
    cc = Array{Float64}(nc)
    dc = Array{Float64}(nc)

    # real * complex
    gen = product(zip(c1[1], c1[2]), zip(c2[3:end]...))
    gen = chain(gen, product(zip(c2[1], c2[2]), zip(c1[3:end]...)))
    for (i, ((aj, cj), (ak, bk, ck, dk))) in enumerate(gen)
        ac[i] = aj * ak
        bc[i] = aj * bk
        cc[i] = cj + ck
        dc[i] = dk
    end

    # complex * complex
    i0 = nr1 * nc2 + nc1 * nr2 + 1
    gen = product(zip(c1[3:end]...), zip(c2[3:end]...))
    for (i, ((aj, bj, cj, dj), (ak, bk, ck, dk))) in enumerate(gen)
        ac[i0 + 2*(i-1)] = 0.5 * (aj * ak + bj * bk)
        bc[i0 + 2*(i-1)] = 0.5 * (bj * ak - aj * bk)
        cc[i0 + 2*(i-1)] = cj + ck
        dc[i0 + 2*(i-1)] = dj - dk

        ac[i0 + 2*(i-1) + 1] = 0.5 * (aj * ak - bj * bk)
        bc[i0 + 2*(i-1) + 1] = 0.5 * (bj * ak + aj * bk)
        cc[i0 + 2*(i-1) + 1] = cj + ck
        dc[i0 + 2*(i-1) + 1] = dj + dk
    end

    return ar, cr, ac, bc, cc, dc
end

function *(t1::Term, t2::Term)
    return TermProduct(t1, t2)
end

function length(term_prod::TermProduct)
    return length(term_prod.term1) + length(term_prod.term2)
end

function get_parameter_vector(term_prod::TermProduct)
    return cat(1, get_parameter_vector(term_prod.term1), get_parameter_vector(term_prod.term2))
end

function set_parameter_vector!(term_prod::TermProduct, vector::Array)
    l = length(term_prod.term1)
    set_parameter_vector!(term_prod.term1, vector[1:l])
    set_parameter_vector!(term_prod.term2, vector[l+1:end])
end

# A real term where b=0 and d=0
mutable struct RealTerm <: Term
    log_a::Float64
    log_c::Float64
end

function get_real_coefficients(term::RealTerm)
    return [exp(term.log_a)], [exp(term.log_c)]
end

function length(term::RealTerm)
    return 2
end

function get_parameter_vector(term::RealTerm)
    return [term.log_a, term.log_c]
end

function set_parameter_vector!(term::RealTerm, vector::Array)
    term.log_a = vector[1]
    term.log_c = vector[2]
end

# General celerite term
mutable struct ComplexTerm <: Term
    log_a::Float64
    log_b::Float64
    log_c::Float64
    log_d::Float64
end

function get_complex_coefficients(term::ComplexTerm)
    return [exp(term.log_a)], [exp(term.log_b)], [exp(term.log_c)], [exp(term.log_d)]
end

function length(term::ComplexTerm)
    return 4
end

function get_parameter_vector(term::ComplexTerm)
    return [term.log_a, term.log_b, term.log_c, term.log_d]
end

function set_parameter_vector!(term::ComplexTerm, vector::Array)
    term.log_a = vector[1]
    term.log_b = vector[2]
    term.log_c = vector[3]
    term.log_d = vector[4]
end

# General celerite term with negative parameters
mutable struct ComplexTermNeg <: Term
    a::Float64
    b::Float64
    c::Float64
    d::Float64
end

function get_complex_coefficients(term::ComplexTermNeg)
    return [term.a], [term.b], [term.c], [term.d]
end

function length(term::ComplexTermNeg)
    return 4
end

function get_parameter_vector(term::ComplexTermNeg)
    return [term.a, term.b, term.c, term.d]
end

function set_parameter_vector!(term::ComplexTermNeg, vector::Array)
    term.a = vector[1]
    term.b = vector[2]
    term.c = vector[3]
    term.d = vector[4]
end

# A SHO term
mutable struct SHOTerm <: Term
    log_S0
    log_Q
    log_omega0
end

function get_real_coefficients(term::SHOTerm)
    Q = exp(term.log_Q)
    if Q >= 0.5
        return Array{Float64}(undef, 0), Array{Float64}(undef, 0)
    end
    S0 = exp(term.log_S0)
    w0 = exp(term.log_omega0)
    f = sqrt(1.0 - 4.0 * Q^2)
    return (
        0.5*S0*w0*Q*[1.0+1.0/f, 1.0-1.0/f],
        0.5*w0/Q*[1.0-f, 1.0+f]
    )
end

function get_complex_coefficients(term::SHOTerm)
    Q = exp(term.log_Q)
    if Q < 0.5
        return Array{Float64}(undef, 0), Array{Float64}(undef, 0), Array{Float64}(undef, 0), Array{Float64}(undef, 0)
    end
    S0 = exp(term.log_S0)
    w0 = exp(term.log_omega0)
    f = sqrt(4.0 * Q^2 - 1.0)
    return (
        [S0 * w0 * Q],
        [S0 * w0 * Q / f],
        [0.5 * w0 / Q],
        [0.5 * w0 / Q * f],
    )
end

function length(term::SHOTerm)
    return 3
end

function get_parameter_vector(term::SHOTerm)
    return [term.log_S0, term.log_Q, term.log_omega0]
end

function set_parameter_vector!(term::SHOTerm, vector::Array)
    term.log_S0 = vector[1]
    term.log_Q = vector[2]
    term.log_omega0 = vector[3]
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the /(x, y) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

@testset "/(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "/(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "/(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::Blade, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::Blade, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "/(B::Blade, C::Scalar)" begin
    @test_skip 1
end

@testset "/(B::Blade, C::One)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))
    C = One()

    @test B / C === B
end

@testset "/(B::Blade, C::Zero)" begin
    @test_skip 1
end

@testset "/(B::Blade, C::Real)" begin
    @test_skip 1
end

@testset "/(B::Blade, C::Vector)" begin
    @test_skip 1
end

# ------ B::Pseudoscalar

@testset "/(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::Pseudoscalar, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::Pseudoscalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_dim = 10

    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value_1)

    # --- Tests

    # B != C, C.value != Inf
    test_value_2 = abs(test_value_1) + rand() + 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == test_value_1 / test_value_2

    # B == C, C.value != Inf
    C = Pseudoscalar(test_dim, test_value_1)

    B_slash_C = B / C
    @test B_slash_C isa One

    # B != C, C.value == Inf
    test_value_2 = Inf
    C = Pseudoscalar(test_dim, test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == C, C.value == Inf
    test_value = Inf
    B = Pseudoscalar(test_dim, test_value)
    C = Pseudoscalar(test_dim, test_value)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end

@testset "/(B::Pseudoscalar, C::Scalar)" begin
    # --- Preparations

    test_dim = 10
    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value_1)

    # --- Tests

    # B != Inf, C != Inf
    test_value_2 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test B_slash_C.value == test_value_1 / test_value_2

    # B != Inf, C == Inf
    test_value_2 = Inf
    C = Scalar(test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == Inf, C == Inf
    test_value = Inf
    B = Pseudoscalar(test_dim, test_value)
    C = Scalar(test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test isnan(B_slash_C.value)
end

@testset "/(B::Pseudoscalar, C::One)" begin
    test_dim = 10
    test_value = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value)

    C = One()

    @test B / C === B
end

@testset "/(B::Pseudoscalar, C::Zero)" begin
    test_dim = 10
    test_value = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value)

    C = Zero()

    expected_result = sign(B) > 0 ? Inf : -Inf

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test B_slash_C.value == expected_result
end

@testset "/(B::Pseudoscalar, C::Real)" begin
    # --- Preparations

    test_dim = 10
    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    B = Pseudoscalar(test_dim, test_value_1)

    # --- Tests

    # B.value != Inf, C != Inf
    test_value_2 = get_random_value(1)  # add 1 to keep value away from 0
    C = test_value_2

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test B_slash_C.value == test_value_1 / test_value_2

    # B.value != Inf, C == Inf
    C = Inf

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B.value == Inf, C == Inf
    B = Pseudoscalar(test_dim, Inf)
    C = Inf

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test isnan(B_slash_C.value)
end

@testset "/(B::Pseudoscalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::Scalar

@testset "/(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::Scalar, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::Scalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # B != Inf, C.value != Inf
    test_value_2 = get_random_value(1)  # add 1 to keep value away from 0

    for test_dim in 5:8
        C = Pseudoscalar(test_dim, test_value_2)

        expected_result = mod(test_dim, 4) < 2 ?
            Pseudoscalar(test_dim, test_value_1 / test_value_2) :
            Pseudoscalar(test_dim, -test_value_1 / test_value_2)

        B_slash_C = B / C
        @test B_slash_C isa Pseudoscalar
        @test B_slash_C.dim == expected_result.dim
        @test B_slash_C.value == expected_result.value
    end

    # B != Inf, C.value == Inf
    test_dim = 10
    test_value_2 = Inf
    C = Pseudoscalar(test_dim, test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == Inf, C.value == Inf
    test_dim = 10
    test_value = Inf
    B = Scalar(test_value)
    C = Pseudoscalar(test_dim, test_value)

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test isnan(B_slash_C.value)
end

@testset "/(B::Scalar, C::Scalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # B != C, B != Inf, C != Inf
    test_value_2 = abs(test_value_1) + rand() + 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == test_value_1 / test_value_2

    # B == C, B != Inf
    C = Scalar(test_value_1)

    B_slash_C = B / C
    @test B_slash_C isa One

    # B != Inf, C == Inf
    C = Scalar(Inf)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == Inf, C == Inf
    test_value = Inf
    B = Scalar(test_value)
    C = Scalar(test_value)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end

@testset "/(B::Scalar, C::One)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = One()

    B_slash_C = B / C
    @test B_slash_C === B
end

@testset "/(B::Scalar, C::Zero)" begin
    # --- Preparations

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1

    C = Zero()

    # --- Tests

    # B > 0
    B = Scalar(abs(test_value))
    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == Inf

    # B < 0
    B = Scalar(-abs(test_value))
    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == -Inf
end

@testset "/(B::Scalar, C::Real)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    # --- Tests

    # B != C, B != Inf, C != Inf, C != 0
    test_value_2 = abs(test_value_1) + rand() + 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == test_value_1 / test_value_2

    # B == C, B != Inf, C != 0
    C = test_value_1

    B_slash_C = B / C
    @test B_slash_C isa One

    # B != Inf, C == Inf
    C = Inf

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B > 0, B != Inf, C == 0
    B = Scalar(abs(test_value_1))
    C = 0

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == Inf

    # B < 0, B != -Inf, C == 0
    B = Scalar(-abs(test_value_1))
    C = 0

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == -Inf

    # B == Inf, C == Inf
    test_value = Inf
    B = Scalar(test_value)
    C = test_value

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end

@testset "/(B::Scalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::One

@testset "/(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::One, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::One, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0

    B = One()

    # --- Tests

    # C.value != Inf
    for test_dim in 5:8
        C = Pseudoscalar(test_dim, test_value_1)

        expected_result = mod(test_dim, 4) < 2 ?
            Pseudoscalar(test_dim, 1/ test_value_1) :
            Pseudoscalar(test_dim, -1/ test_value_1)

        B_slash_C = B / C
        @test B_slash_C isa Pseudoscalar
        @test B_slash_C.dim == expected_result.dim
        @test B_slash_C.value == expected_result.value
    end

    # C.value == Inf
    test_dim = 10
    test_value_2 = Inf
    C = Pseudoscalar(test_dim, test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Zero
end

@testset "/(B::One, C::Scalar)" begin
    # C != Inf
    B = One()

    test_value = get_random_value(2)  # add 1 to keep value away from 0 and 1
    C = Scalar(test_value)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == 1 / test_value

    # C == Inf
    C = Scalar(Inf)

    B_slash_C = B / C
    @test B_slash_C isa Zero
end

@testset "/(B::One, C::One)" begin
    B = One()
    C = One()
    @test isone(B / C)
end

@testset "/(B::One, C::Zero)" begin
    B = One()
    C = Zero()

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C == Inf
end

@testset "/(B::One, C::Real)" begin
    B = One()

    # B != C, C != Inf, C != 0
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = test_value

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == 1 / test_value

    # B == C
    C = 1

    B_slash_C = B / C
    @test B_slash_C isa One

    # C == Inf
    C = Inf

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # C == 0
    C = 0

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == Inf
end

@testset "/(B::One, C::Vector)" begin
    @test_skip 1
end

# ------ B::Zero

@testset "/(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::Zero, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::Zero, C::Pseudoscalar)" begin
    B = Zero()

    test_dim = 10
    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    C = Pseudoscalar(test_dim, test_value_1)

    @test iszero(B / C)
end

@testset "/(B::Zero, C::Scalar)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = Zero()

    @test iszero(C / B)
end

@testset "/(B::Zero, C::One)" begin
    B = Zero()
    C = One()
    @test iszero(B / C)
end

@testset "/(B::Zero, C::Zero)" begin
    B = Zero()
    C = Zero()

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end

@testset "/(B::Zero, C::Real)" begin
    # --- Preparations

    B = Zero()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1

    # --- Tests

    # C != 0
    C = test_value

    @test iszero(B / C)

    # C == 0
    C = 0

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end

@testset "/(B::Zero, C::Vector)" begin
    @test_skip 1
end

# ------ B::Real

@testset "/(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::Real, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::Real, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = test_value_1

    # --- Tests

    # B != Inf, B != 0, C.value != Inf
    test_value_2 = get_random_value(1)  # add 1 to keep value away from 0

    for test_dim in 5:8
        C = Pseudoscalar(test_dim, test_value_2)

        expected_result = mod(test_dim, 4) < 2 ?
            Pseudoscalar(test_dim, test_value_1 / test_value_2) :
            Pseudoscalar(test_dim, -test_value_1 / test_value_2)

        B_slash_C = B / C
        @test B_slash_C isa Pseudoscalar
        @test B_slash_C.dim == expected_result.dim
        @test B_slash_C.value == expected_result.value
    end

    # B != Inf, B != 0, C.value == Inf
    test_dim = 10
    test_value_2 = Inf
    C = Pseudoscalar(test_dim, test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == 0
    B = 0

    test_dim = 10
    C = Pseudoscalar(test_dim, test_value_1)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == Inf, C.value == Inf
    test_value = Inf
    test_dim = 10
    B = test_value
    C = Pseudoscalar(test_dim, test_value)

    B_slash_C = B / C
    @test B_slash_C isa Pseudoscalar
    @test B_slash_C.dim == test_dim
    @test isnan(B_slash_C.value)
end

@testset "/(B::Real, C::Scalar)" begin
    # --- Preparations

    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = test_value_1

    # --- Tests

    # B != C, B != Inf, C != Inf
    test_value_2 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == test_value_1 / test_value_2

    # B == C, B != Inf
    C = Scalar(test_value_1)

    B_slash_C = B / C
    @test B_slash_C isa One

    # B != Inf, C == Inf
    test_value_2 = Inf
    C = Scalar(test_value_2)

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == 0
    B = 0

    B_slash_C = B / C
    @test B_slash_C isa Zero

    # B == Inf, C == Inf
    test_value = Inf
    B = test_value
    C = Scalar(test_value)

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end

@testset "/(B::Real, C::One)" begin
    # --- Preparations

    C = One()

    # --- Tests

    # B != C, B != 0
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = test_value

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C.value == test_value

    # B == C
    B = 1

    B_slash_C = B / C
    @test B_slash_C isa One

    # B == 0
    B = 0

    B_slash_C = B / C
    @test B_slash_C isa Zero
end

@testset "/(B::Real, C::Zero)" begin
    # --- Preparations

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1

    C = Zero()

    # --- Tests

    # B > 0
    B = abs(test_value)
    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C == Inf

    # B = 0
    B = 0
    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)

    # B < 0
    B = -abs(test_value)
    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test B_slash_C == -Inf
end

# ------ B::Vector

@testset "/(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "/(B::Vector, C::Blade)" begin
    @test_skip 1
end

@testset "/(B::Vector, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "/(B::Vector, C::Scalar)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    # B is not zero
    B = rand(5)

    B_slash_C = B / C
    @test B_slash_C isa Blade
    @test B_slash_C == Blade(B) / test_value

    # B is zero
    B = [0; 0; 0; 0; 0]
    @test iszero(B / C);
end

@testset "/(B::Vector, C::One)" begin
    C = One()

    # B is not zero
    B = rand(5)

    B_slash_C = B / C
    @test B_slash_C isa Blade
    @test B_slash_C == Blade(B)

    # B is zero
    B = [0; 0; 0; 0; 0]
    @test iszero(B / C)
end

@testset "/(B::Vector, C::Zero)" begin
    C = Zero()

    # B is not zero
    B = rand(5)

    B_slash_C = B / C
    @test B_slash_C isa Blade
    @test B_slash_C == Blade(B) / 0

    # B is zero
    B = [0; 0; 0; 0; 0]

    B_slash_C = B / C
    @test B_slash_C isa Scalar
    @test isnan(B_slash_C.value)
end


text:
"""
    bump()

Lorem ipsum dolor sit amet.
"""
function bump(args::Dict)
  bump()
end

function bump()

  package_name = rsplit(pwd(), "/"; limit=2)[2]

  package_dir = rsplit(dirname(@__FILE__), "/"; limit=9)[1]

  target_name = replace(package_name, ".jl", "")

  target_path = "$package_dir/$target_name"

  if target_path != pwd()
    cd("..") do
      cp(package_name, target_path, remove_destination=true)
    end
  end

  return target_name

end


text:
using SymbolicNumericIntegration
using Symbolics

using SymbolicUtils
using SymbolicUtils.Rewriters

using Test
using PyCall
import SymPy

include("axiom.jl")

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

@variables x β

"""
    a list of basic standard integral tests
    based on http://integral-table.com/ with modifications
"""
basic_integrals = [
    # Basic Forms
    1,
    x^2,
    4x^3,
    # Integrals of Rational Functions
    1 / x,
    1 / (2x + 5),
    1 / (x + 1)^2,
    (x + 3)^3,
    x * (x - 2)^4,
    1 / (1 + x^2),
    1 / (9 + x^2),
    x / (4 + x^2),
    x^2 / (16 + x^2),
    x^3 / (1 + x^2),
    1 / (x^2 - 5x + 6),
    1 / (x^2 + x + 1),
    x / (x + 4)^2,
    x / (x^2 + x + 1),
    # Integrals with Roots
    sqrt(x - 2),
    1 / sqrt(x - 1),
    1 / sqrt(x + 1),
    1 / sqrt(4 - x),
    x * sqrt(x - 3),
    sqrt(2x + 5),
    (3x - 1)^1.5,
    x / sqrt(x - 1),
    x / sqrt(x + 1),
    sqrt(x / (4 - x)),
    sqrt(x / (4 + x)),
    x * sqrt(2x + 3),
    sqrt(x * (x + 2)),
    sqrt(x^3 * (x + 3)),
    sqrt(x^2 + 4),
    sqrt(x^2 - 4),
    sqrt(4 - x^2),
    x * sqrt(x^2 + 9),
    x * sqrt(x^2 - 9),
    1 / sqrt(x^2 + 4),
    1 / sqrt(x^2 - 4),
    1 / sqrt(4 - x^2),
    x / sqrt(x^2 + 4),
    x / sqrt(x^2 - 4),
    x / sqrt(4 - x^2),
    x^2 / sqrt(x^2 + 4),
    x^2 / sqrt(x^2 - 4),
    sqrt(x^2 - 5x + 6),
    x * sqrt(x^2 - 5x + 6),
    1 / sqrt(x^2 - 5x + 6),
    1 / (4 + x^2)^1.5,
    # Integrals with Logarithms
    log(x),
    x * log(x),
    x^2 * log(x),
    log(2x) / x,
    log(x) / x^2,
    log(2x + 1),
    log(x^2 + 4),
    log(x^2 - 4),
    log(x^2 - 5x + 6),
    x * log(x + 2),
    x * log(9 - 4x^2),
    log(x)^2,
    log(x)^3,
    x * log(x)^2,
    x^2 * log(x)^2,
    # Integrals with Exponentials
    exp(x),
    sqrt(x) * exp(x),
    x * exp(x),
    x * exp(3x),
    x^2 * exp(x),
    x^2 * exp(5x),
    x^3 * exp(x),
    x^3 * exp(2x),
    exp(x^2),
    x * exp(x^2),
    # Integrals with Trigonometric Functions
    sin(4x),
    sin(x)^2,
    sin(x)^3,
    cos(3x),
    cos(x)^2,
    cos(2x)^3,
    sin(x) * cos(x),
    sin(3x) * cos(5x),
    sin(x)^2 * cos(x),
    sin(3x)^2 * cos(x),
    sin(x) * cos(x)^2,
    sin(x) * cos(5x)^2,
    sin(x)^2 * cos(x),
    sin(x)^2 * cos(x)^2,
    sin(4x)^2 * cos(4x^2),
    tan(x),
    tan(7x),
    tan(x)^2,
    tan(x)^3,
    sec(x),
    sec(x) * tan(x),
    sec(x)^2 * tan(x),
    csc(x),
    sec(x) * csc(x),
    # Products of Trigonometric Functions and Monomials
    x * cos(x),
    x * cos(3x),
    x^2 * cos(x),
    x^2 * cos(5x),
    x * sin(x),
    x * sin(3x),
    x^2 * sin(x),
    x^2 * sin(5x),
    x * cos(x)^2,
    x * sin(x)^2,
    x * tan(x)^2,
    x * sec(x)^2,
    x^3 * sin(x),
    x^4 * cos(2x),
    sin(x)^2 * cos(x)^3,
    # Products of Trigonometric Functions and Exponentials
    exp(x) * sin(x),
    exp(3x) * sin(2x),
    exp(x) * cos(x),
    exp(2x) * cos(7x),
    x * exp(x) * sin(x),
    x * exp(x) * cos(x),
    # Integrals of Hyperbolic Functions
    cosh(x),
    exp(x) * cosh(x),
    sinh(3x),
    exp(2x) * sinh(3x),
    tanh(x),
    exp(x) * tanh(x),
    cos(x) * cosh(x),
    cos(x) * sinh(x),
    sin(x) * cosh(x),
    sin(x) * sinh(x),
    sinh(x) * cosh(x),
    sinh(3x) * cosh(5x),
    # Misc
    exp(x) / (1 + exp(x)),
    cos(exp(x)) * sin(exp(x)) * exp(x),
    cos(exp(x))^2 * sin(exp(x)) * exp(x),
    1 / (x * log(x)),
    (log(x - 1) + (x - 1)^-1) * log(x),
    1 / (exp(x) - 1),
    1 / (exp(x) + 5),
    sqrt(x) * log(x),
    log(log(x)) / x,
    x^3 * exp(x^2),
    sin(log(x)),
    x * cos(x) * exp(x),
    log(x - 1)^2,
    1 / (exp(2x) - 1),
    exp(x) / (exp(2x) - 1),
    x / (exp(2x) - 1),
    # derivative-divide examples (Lamangna 7.10.2)
    exp(x) * exp(exp(x)),
    exp(sqrt(x)) / sqrt(x),
    log(log(x)) / (x * log(x)),
    log(cos(x)) * tan(x),
    # rothstein-Trager examples (Lamangna 7.10.9)
    1 / (x^3 - x),
    1 / (x^3 + 1),
    1 / (x^2 - 8),
    (x + 1) / (x^2 + 1),
    x / (x^4 - 4),
    x^3 / (x^4 + 1),
    1 / (x^4 + 1),
    # bypass = true
    β,      # turn of bypass = true
    exp(x) / x - exp(x) / x^2,
    cos(x) / x - sin(x) / x^2,
    1 / log(x) - 1 / log(x)^2,
]

function test_integrals(; kw...)
    args = Dict(kw)
    misses = []
    k = 1

    for eq in basic_integrals
        if isequal(eq, β)
            printstyled("**** bypass on ****\n"; color=:red)
            bypass = true
            args[:bypass] = true
        else
            printstyled(k, ": "; color=:blue)
            k += 1
            printstyled(eq, " =>\n"; color=:green)
            solved, unsolved = SymbolicNumericIntegration.integrate(eq; args...)
            printstyled('\t', solved; color=:cyan)
            if isequal(unsolved, 0)
                println()
            else
                printstyled(" + ∫ ", unsolved, '\n'; color=:red)
                push!(misses, eq)
            end
        end
    end

    n = length(misses)
    if n > 0
        println("**** missess (n=$n) *****")
    end
    for eq in misses
        printstyled(eq, '\n'; color=:red)
    end
    return n
end

@testset "integral" begin
    n = test_integrals(; symbolic=false, verbose=false, homotopy=true, num_steps=2, num_trials=10)
    @test n > 0
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

import Base.Threads: @threads

function _find_violations(
    model::JuMP.Model;
    max_per_line::Int,
    max_per_period::Int,
)
    instance = model[:instance]
    net_injection = model[:net_injection]
    overflow = model[:overflow]
    length(instance.buses) > 1 || return []
    violations = []
    @info "Verifying transmission limits..."
    time_screening = @elapsed begin
        non_slack_buses = [b for b in instance.buses if b.offset > 0]
        net_injection_values = [
            value(net_injection[b.name, t]) for b in non_slack_buses,
            t in 1:instance.time
        ]
        overflow_values = [
            value(overflow[lm.name, t]) for lm in instance.lines,
            t in 1:instance.time
        ]
        violations = UnitCommitment._find_violations(
            instance = instance,
            net_injections = net_injection_values,
            overflow = overflow_values,
            isf = model[:isf],
            lodf = model[:lodf],
            max_per_line = max_per_line,
            max_per_period = max_per_period,
        )
    end
    @info @sprintf(
        "Verified transmission limits in %.2f seconds",
        time_screening
    )
    return violations
end

"""
    function _find_violations(
        instance::UnitCommitmentInstance,
        net_injections::Array{Float64, 2};
        isf::Array{Float64,2},
        lodf::Array{Float64,2},
        max_per_line::Int,
        max_per_period::Int,
    )::Array{_Violation, 1}

Find transmission constraint violations (both pre-contingency, as well as
post-contingency).

The argument `net_injection` should be a (B-1) x T matrix, where B is the
number of buses and T is the number of time periods. The arguments `isf` and
`lodf` can be computed using UnitCommitment.injection_shift_factors and
UnitCommitment.line_outage_factors. The argument `overflow` specifies how much
flow above the transmission limits (in MW) is allowed. It should be an L x T
matrix, where L is the number of transmission lines.
"""
function _find_violations(;
    instance::UnitCommitmentInstance,
    net_injections::Array{Float64,2},
    overflow::Array{Float64,2},
    isf::Array{Float64,2},
    lodf::Array{Float64,2},
    max_per_line::Int,
    max_per_period::Int,
)::Array{_Violation,1}
    B = length(instance.buses) - 1
    L = length(instance.lines)
    T = instance.time
    K = nthreads()

    size(net_injections) == (B, T) || error("net_injections has incorrect size")
    size(isf) == (L, B) || error("isf has incorrect size")
    size(lodf) == (L, L) || error("lodf has incorrect size")

    filters = Dict(
        t => _ViolationFilter(
            max_total = max_per_period,
            max_per_line = max_per_line,
        ) for t in 1:T
    )

    pre_flow::Array{Float64} = zeros(L, K)           # pre_flow[lm, thread]
    post_flow::Array{Float64} = zeros(L, L, K)       # post_flow[lm, lc, thread]
    pre_v::Array{Float64} = zeros(L, K)              # pre_v[lm, thread]
    post_v::Array{Float64} = zeros(L, L, K)          # post_v[lm, lc, thread]

    normal_limits::Array{Float64,2} = [
        l.normal_flow_limit[t] + overflow[l.offset, t] for
        l in instance.lines, t in 1:T
    ]

    emergency_limits::Array{Float64,2} = [
        l.emergency_flow_limit[t] + overflow[l.offset, t] for
        l in instance.lines, t in 1:T
    ]

    is_vulnerable::Array{Bool} = zeros(Bool, L)
    for c in instance.contingencies
        is_vulnerable[c.lines[1].offset] = true
    end

    @threads for t in 1:T
        k = threadid()

        # Pre-contingency flows
        pre_flow[:, k] = isf * net_injections[:, t]

        # Post-contingency flows
        for lc in 1:L, lm in 1:L
            post_flow[lm, lc, k] =
                pre_flow[lm, k] + pre_flow[lc, k] * lodf[lm, lc]
        end

        # Pre-contingency violations
        for lm in 1:L
            pre_v[lm, k] = max(
                0.0,
                pre_flow[lm, k] - normal_limits[lm, t],
                -pre_flow[lm, k] - normal_limits[lm, t],
            )
        end

        # Post-contingency violations
        for lc in 1:L, lm in 1:L
            post_v[lm, lc, k] = max(
                0.0,
                post_flow[lm, lc, k] - emergency_limits[lm, t],
                -post_flow[lm, lc, k] - emergency_limits[lm, t],
            )
        end

        # Offer pre-contingency violations
        for lm in 1:L
            if pre_v[lm, k] > 1e-5
                _offer(
                    filters[t],
                    _Violation(
                        time = t,
                        monitored_line = instance.lines[lm],
                        outage_line = nothing,
                        amount = pre_v[lm, k],
                    ),
                )
            end
        end

        # Offer post-contingency violations
        for lm in 1:L, lc in 1:L
            if post_v[lm, lc, k] > 1e-5 && is_vulnerable[lc]
                _offer(
                    filters[t],
                    _Violation(
                        time = t,
                        monitored_line = instance.lines[lm],
                        outage_line = instance.lines[lc],
                        amount = post_v[lm, lc, k],
                    ),
                )
            end
        end
    end

    violations = _Violation[]
    for t in 1:instance.time
        append!(violations, _query(filters[t]))
    end

    return violations
end


text:
abstract type AbstractUQ end

"""
    struct UQ1D{
        A<:Integer,
        B<:AbstractString,
        C<:AbstractString,
        D<:AbstractOrthoPoly,
        E<:Union{AbstractVector,Tuple},
        F<:AbstractVector,
        G<:AbstractVector,
        H<:AbstractMatrix,
        I<:AbstractArray,
        J<:AbstractVector,
    } <: AbstractUQ
        nr::A
        nm::A
        nq::A
        method::B
        optype::C
        op::D
        p::E
        phiRan::F
        t1Product::G
        t2Product::H
        t3Product::I
        pce::J
        pceSample::J
    end

Struct of UQ setup

"""
struct UQ1D{
    A<:Integer,
    B<:AbstractString,
    C<:AbstractString,
    D<:AbstractOrthoPoly,
    E<:Union{AbstractVector,Tuple},
    F<:AbstractMatrix,
    G<:AbstractVector,
    H<:AbstractMatrix,
    I<:AbstractArray,
    J<:AbstractVector,
} <: AbstractUQ
    nr::A
    nm::A
    nq::A
    method::B
    optype::C
    op::D
    p::E
    phiRan::F
    t1Product::G
    t2Product::H
    t3Product::I
    pce::J
    pceSample::J
end

function UQ1D(
    NR::Int,
    NREC::Int,
    P1::AbstractFloat,
    P2::AbstractFloat,
    TYPE = "uniform"::AbstractString,
    METHOD = "collocation"::AbstractString,
)
    method = METHOD
    nr = NR
    nRec = NREC
    optype = TYPE
    nm = nr # 1D

    if TYPE == "gauss"
        op = GaussOrthoPoly(nr, Nrec = nRec, addQuadrature = true)
    elseif TYPE == "uniform"
        # uniform ∈ [0,1]
        # op = Uniform01OrthoPoly(nr, Nrec=nRec, addQuadrature=true)

        # uniform ∈ [-1, 1]
        op = Uniform_11OrthoPoly(nr, Nrec = nRec, addQuadrature = true)
        #supp = (-1.0, 1.0)
        #uni_meas = Measure("uni_meas", x -> 0.5, supp, true, Dict())
        #op = OrthoPoly("uni_op", nr, uni_meas; Nrec = nRec)
    elseif TYPE == "custom"
        supp = (P1, P2)
        uni_meas = Measure("uni_meas", x -> 1 / (P2 - P1), supp, true, Dict())
        op = OrthoPoly("uni_op", nr, uni_meas; Nrec = nRec)
    else
        @warn "polynomial chaos unavailable"
    end

    nq = op.quad.Nquad

    p1 = P1
    p2 = P2
    p = [p1, p2]

    phiRan = evaluate(collect(0:op.deg), op.quad.nodes, op)

    t1 = PolyChaos.Tensor(1, op) # < \phi_i >
    t2 = PolyChaos.Tensor(2, op) # < \phi_i \phi_j >
    t3 = PolyChaos.Tensor(3, op) # < \phi_i \phi_j \phi_k >

    t1Product = OffsetArray{Float64}(undef, 0:nr)
    t2Product = OffsetArray{Float64}(undef, 0:nr, 0:nr)
    t3Product = OffsetArray{Float64}(undef, 0:nr, 0:nr, 0:nr)
    for i = 0:nr
        t1Product[i] = t1.get([i])
    end
    for i = 0:nr
        for j = 0:nr
            t2Product[i, j] = t2.get([i, j])
        end
    end
    for i = 0:nr
        for j = 0:nr
            for k = 0:nr
                t3Product[i, j, k] = t3.get([i, j, k])
            end
        end
    end

    if TYPE == "gauss"
        pce = [convert2affinePCE(p1, p2, op); zeros(nr - 1)]
    elseif TYPE == "uniform"
        # pce = [ convert2affinePCE(p1, p2, op); zeros(nr-1) ] # uniform ∈ [0, 1]
        pce = [[0.5 * (p1 + p2), 0.5 * (p2 - p1)]; zeros(nr - 1)] # uniform ∈ [-1, 1]
    else
        pce = [[0.5 * (p1 + p2), 0.5 * (p2 - p1)]; zeros(nr - 1)] # uniform ∈ [-1, 1]
    end

    # pceSample = [1.0] # test
    # pceSample= samplePCE(2000, pce, op) # Monte-Carlo
    pceSample = evaluatePCE(pce, op.quad.nodes, op) # collocation

    return UQ1D(
        nr,
        nm,
        nq,
        method,
        optype,
        op,
        p,
        phiRan,
        t1Product,
        t2Product,
        t3Product,
        pce,
        pceSample,
    )
end

"""
    struct UQ2D{
        A<:Integer,
        B<:AbstractString,
        C<:Union{AbstractVector,Tuple},
        D<:AbstractVector,
        E<:AbstractArray{<:AbstractFloat,2},
        F<:AbstractVector,
        G<:AbstractMatrix,
        H<:AbstractArray,
        I<:AbstractMatrix,
        J<:AbstractVector,
    } <: AbstractUQ
        nr::A
        nm::A
        nq::A
        method::B
        optype::C
        op::MultiOrthoPoly
        p::D
        phiRan::E
        t1Product::F
        t2Product::G
        t3Product::H
        points::I
        weights::J
    end

Struct of UQ setup

"""
struct UQ2D{
    A<:Integer,
    B<:AbstractString,
    C<:Union{AbstractVector,Tuple},
    D<:AbstractVector,
    E<:AbstractArray{<:AbstractFloat,2},
    F<:AbstractVector,
    G<:AbstractMatrix,
    H<:AbstractArray,
    I<:AbstractMatrix,
    J<:AbstractVector,
} <: AbstractUQ
    nr::A
    nm::A
    nq::A
    method::B
    optype::C
    op::MultiOrthoPoly
    p::D
    phiRan::E
    t1Product::F
    t2Product::G
    t3Product::H
    points::I
    weights::J
end

function UQ2D(
    NR::Integer,
    NREC::Integer,
    P::AbstractVector,
    TYPE = ["uniform", "uniform"],
    METHOD = "collocation",
)

    ops = map(TYPE) do x
        if x == "uniform"
            return Uniform_11OrthoPoly(NR, Nrec = NREC, addQuadrature = true)
        elseif x == "gauss"
            return GaussOrthoPoly(NR, Nrec = NREC, addQuadrature = true)
        else
            throw("No default polynomials available")
        end
    end

    phi = MultiOrthoPoly(ops, 4)
    nm = size(phi.ind, 1) - 1

    t1 = PolyChaos.Tensor(1, phi)
    t2 = PolyChaos.Tensor(2, phi)
    t3 = PolyChaos.Tensor(3, phi)

    t1Product = OffsetArray{Float64}(undef, 0:nm)
    t2Product = OffsetArray{Float64}(undef, 0:nm, 0:nm)
    t3Product = OffsetArray{Float64}(undef, 0:nm, 0:nm, 0:nm)
    for i = 0:nm
        t1Product[i] = t1.get([i])
    end
    for i = 0:nm
        for j = 0:nm
            t2Product[i, j] = t2.get([i, j])
        end
    end
    for i = 0:nm
        for j = 0:nm
            for k = 0:nm
                t3Product[i, j, k] = t3.get([i, j, k])
            end
        end
    end

    p = [[P[1], P[2]], [P[3], P[4]]]

    nq = ops[1].quad.Nquad * ops[2].quad.Nquad
    weights = zeros(nq)
    points = zeros(nq, 2)
    for i = 1:ops[1].quad.Nquad, j = 1:ops[2].quad.Nquad
        idx = ops[1].quad.Nquad * (j - 1) + i

        points[idx, 1] = ops[1].quad.nodes[i]
        points[idx, 2] = ops[2].quad.nodes[j]
        weights[idx] = ops[1].quad.weights[i] * ops[2].quad.weights[j]
    end

    phiRan = evaluate(phi.ind, points, phi) |> permutedims

    return UQ2D(
        NR,
        nm,
        nq,
        METHOD,
        TYPE,
        phi,
        p,
        phiRan,
        t1Product,
        t2Product,
        t3Product,
        points,
        weights,
    )

end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Random
using ElasticArrays, ArraysOfArrays

@testset "array_shape" begin
    @test_throws ArgumentError @inferred(ValueShapes._valshapeoftype(Vector{Int}))

    @test @inferred(valshape(rand(3))) == ArrayShape{Float64,1}((3,))
    @test @inferred(valshape(rand(3, 4, 5))) == ArrayShape{Float64,3}((3, 4, 5))

    @test @inferred(ValueShapes.default_unshaped_eltype(ArrayShape{Complex,3}((3, 4, 5)))) == Float64
    @test @inferred(ValueShapes.default_unshaped_eltype(ArrayShape{Complex{Float32},3}((3, 4, 5)))) == Float32

    @test @inferred(ValueShapes.shaped_type(ArrayShape{Real}(2, 3, 4))) == Array{Float64, 3}
    @test @inferred(ValueShapes.shaped_type(ArrayShape{Complex}(2))) == Array{Complex{Float64}, 1}
    @test @inferred(ValueShapes.shaped_type(ArrayShape{Complex{Real}}(2), Float32)) == Array{Complex{Float32}, 1}
    @test @inferred(ValueShapes.shaped_type(ArrayShape{Complex{Int16}}(2, 3))) == Array{Complex{Int16}, 2}

    @test @inferred(totalndof(ArrayShape{Float64,1}((3,)))) == 3
    @test @inferred(totalndof(ArrayShape{Complex,3}((3, 4, 5)))) == 120

    @test size(@inferred(Vector{Float64}(undef, ArrayShape{Complex}((2, 1, 3))))) == (2 * 2*1*3,)
    @test size(flatview(@inferred(VectorOfSimilarVectors{Float32}(ArrayShape{Complex}((2, 1, 3)))))) == (2 * 2*1*3, 0)

    shape = ArrayShape{Real}(2,3)
    A = @inferred(shape(undef))
    @test typeof(A) == Array{Float64,2}
    @test size(A) == (2, 3)
    @test @inferred(length(shape)) == shape.dims[1]*shape.dims[2]
    @test @inferred(length(shape)) == prod(shape.dims)
    @test @inferred(valshape(A)) == ArrayShape{Float64}(2,3)
    @test @inferred(shape([1, 2, 3, 4, 5, 6])) == [1 3 5; 2 4 6]
    @test_throws ArgumentError shape([1, 2, 3, 4, 5, 6, 7])

    @test @inferred(ArrayShape{Float64}(4,5) <= ArrayShape{Real}(4,5)) == true
    @test @inferred(ArrayShape{Real}(4,5) >= ArrayShape{Float64}(4,5)) == true
    @test (ArrayShape{Float64}(4,5) <= ArrayShape{Integer}(4,5)) == false
    @test (ArrayShape{Float64}(2,2) <= ArrayShape{Float64}(3,3)) == false

    @test eltype(@inferred(Vector{Float32}(undef, shape))) == Float32
    @test eltype(eltype(@inferred(VectorOfSimilarVectors{Float32}(shape)))) == Float32

    @test valshape(shape.(push!(@inferred(VectorOfSimilarVectors{Float64}(shape)), @inferred(Vector{Float64}(undef, shape))))[1]) == valshape(shape(undef))

    let A = [2.2, 4.4, 3.3]
        @test @inferred(unshaped(A, ArrayShape{Real}(3))) === A
        @test_throws ArgumentError unshaped(A, ArrayShape{Real}(2))
        @test_throws ArgumentError unshaped(A, ArrayShape{Integer}(3))
    end
    @test @inferred(unshaped([1 2; 3 4; 5 6], ArrayShape{Real}(3,2))) == [1, 3, 5, 2, 4, 6]
    let shape = ArrayShape{Real}(3,2), UA = [1, 3, 5, 2, 4, 6]
        @test @inferred(unshaped(shape(UA), shape)) == UA
    end

    let
        A = collect(1:8)
        ac = ValueAccessor(ArrayShape{Real}(3), 2)

        @test @inferred(getindex(A, ac)) == [3, 4, 5]
        @test @inferred(view(A, ac)) == [3, 4, 5]
        @test @inferred(setindex!(A, [7, 2, 6], ac)) === A
        @test A[ac] == [7, 2, 6]
    end

    let
        A = collect(1:6*6)
        ac = ValueAccessor(ArrayShape{Real}(2,3), 17)

        @test @inferred(getindex(A, ac)) == [18 20 22; 19 21 23]
        @test @inferred(view(A, ac)) == [18 20 22; 19 21 23]
        @test @inferred(setindex!(A, [6 4 3; 2 1 5], ac)) === A
        @test A[ac] == [6 4 3; 2 1 5]
    end

    let
        A = collect(reshape(1:6*6, 6, 6))
        ac1 = ValueAccessor(ArrayShape{Real}(2), 2)
        ac2 = ValueAccessor(ArrayShape{Real}(3), 3)

        @test @inferred(getindex(A, ac1, ac2)) == [21 27 33; 22 28 34]
        @test @inferred(view(A, ac1, ac2)) == [21 27 33; 22 28 34]
        @test @inferred(setindex!(A, [6 4 3; 2 1 5], ac1, ac2)) === A
        @test A[ac1, ac2] == [6 4 3; 2 1 5]
    end

    let A = collect(reshape(1:5*5*5, 5, 5, 5))
        ac1 = ValueAccessor(ArrayShape{Real}(size(A)[1]), 0)
        ac2 = ValueAccessor(ArrayShape{Real}(1), 0)
        ac3 = ValueAccessor(ArrayShape{Real}(1), 1)
        ac4 = ValueAccessor(ArrayShape{Real}(1), 2)

        @test @inferred(getindex(A, ac1, ac1, ac1)) == A
        @test @inferred(getindex(A, ac2, ac2, ac2))[1] == A[1]
        @test @inferred(getindex(A, ac4, ac4, ac4))[1] == 63

        @test @inferred(view(A, ac1, ac1, ac1)) == A
        @test @inferred(view(A, ac2, ac2, ac2))[1] == A[1]
        @test @inferred(view(A, ac4, ac4, ac4))[1] == 63

        first_layer = @inferred(getindex(A, ac1, ac1, ac2))
        setindex!(A, first_layer, ac1, ac1, ac3)
        @test @inferred(A[:,:,1]) == @inferred(A[:,:,2])
    end

    let d1 = [11, 12, 13, 14], d2 = [21, 22]
        d = vcat(d1, d2)
        reshaped = shape(d)
        @test reshaped == reshape(d, (2,3))
    end
end

@testset "broadcasting and copy" begin
    data1d = [rand(4), rand(4), rand(4), rand(4)]
    data2d = [[rand(4,)] [rand(4,)] [rand(4,)]]

    VoV = VectorOfVectors(data1d)

    shape1 = ArrayShape{Float64, 1}((4,))
    shape2 = ArrayShape{Float64}(1,4)
    shape3 = ArrayShape{Float64}(2,2)

    shape1_VoV_bcast = broadcast(shape1, VoV)
    shape2_VoV_bcast = broadcast(shape2, VoV)
    shape3_VoV_bcast = broadcast(shape3, VoV)

    shape1_bcast = broadcast(shape1, data1d)
    shape1_data_bcast = broadcast(shape1, data1d)
    shape1_data_dcast = shape1.(data1d)
    shape2_data_bcast = broadcast(shape2, data1d)
    shape3_data_bcast = broadcast(shape3, data1d)

    @test shape1_data_dcast == shape1_data_bcast
    @test isapprox(shape1_VoV_bcast, shape1_data_bcast)
    @test isapprox(shape2_VoV_bcast, shape2_data_bcast)
    @test isapprox(shape3_VoV_bcast, shape3_data_bcast)

    AoSV = ArrayOfSimilarVectors{Float64}(data1d)
    AoSA = ArrayOfSimilarArrays{Float64}(data2d)

    shaped_AoSV_bcast = broadcast(shape1, AoSV)
    shaped_AoSV_dcast = shape1.(AoSV)

    shaped_AoSA_bcast = broadcast(shape1, AoSA)
    shaped_AoSA_dcast = shape1.(AoSA)

    @test shaped_AoSV_bcast == AoSV
    @test shaped_AoSV_dcast == shaped_AoSV_bcast
    @test shape1.(VoV) == VoV

    shape1_bcast = broadcast(shape1, data2d)
    @test shape1_bcast == shape1.(data2d)
    for i in 1:length(data2d)
        @test shape1_bcast[i] == data2d[i]
    end
    unshaped1 = unshaped.(shaped_AoSA_bcast)
    @test unshaped1 == data2d
end


text:
"""
Non-blocking poll for events. Should return an element of type `Nothing` if no events were caught.
"""
poll_for_event(wm::AbstractWindowManager) = not_implemented_for(wm)

"""
Wait until the next event occurs.
"""
wait_for_event(wm::AbstractWindowManager) = not_implemented_for(wm)


text:
# Handling regression files

"""
    $(SIGNATURES)

Retrieve a regression coefficient of the form `:afqt3`.
Ignores default dummies (j < 2).
Optional: allow for missing regressors.
Assumes that regressors have been renamed.

## Example
```
get_regr_coef(rt, :afqt, 2) # -> coefficient value for `:afqt2`
get_regr_coef(rt, :parental, 1) == 0.0
```
"""
function get_regr_coef(rt :: RegressionTable, colCat :: Symbol, j :: Integer;
    errorIfMissing :: Bool = true)

    if j < 2
        rCoef = 0.0;
    else
        rName = regressor_name(colCat, j);
        if has_regressor(rt, rName)
            rCoef = get_coefficient(rt, rName);
        elseif errorIfMissing
            error("Cannot retrieve $rName from $rt");
        else
            rCoef = 0.0;
        end
    end
    @assert isa(rCoef, Float64)
    return rCoef
end

function get_intercept(rt :: RegressionTable)
    return get_coefficient(rt, :cons)
end

## ----------  Regressor names

"""
	$(SIGNATURES)

Column headers for RegressionTable objects. Regardless of names in regression files (which have a tendency to change over time), this is what will be returned.

## Example
```
regressor_name(:gpa, 2) == :afqt2
```
"""
function regressor_name(colCat :: Symbol, j :: Integer)
    return Symbol(regressor_string(colCat, j))
end

"""
	$(SIGNATURES)

Retrieve regressor names as Symbols for indices `1 : jMax`.
"""
regressor_names(varName :: Symbol, jMax :: Integer) =
    [regressor_name(varName, j)  for j = 1 : jMax];

regressor_names(varName :: Symbol, jV :: AbstractVector) =
    [regressor_name(varName, j)  for j ∈ jV];

"""
	$(SIGNATURES)

Make name of a dummy regressor, such as `:school2`.

## Example
```
    regressor_string(:gpa, 2) == "afqt2"
```
"""
function regressor_string(varName :: Symbol, idx :: Integer)
    if varName ∈ (:gpa, :afqt)
        regName = "afqt$idx";
    elseif varName == :parental
        regName = "parental$idx"
    elseif varName ∈ (:quality, :school)
        regName = "$varName$idx";
    elseif varName == :lastColl
        regName = "last_type$idx";
    elseif varName == :const
        @assert idx <= 0  "Constant does not have index"
        regName = "cons"
    else
        error("Invalid varName: $varName")
    end
    return regName
end

regressor_strings(varName :: Symbol, idx) =
    [regressor_string(varName, idx1)  for idx1 in idx];

"""
	$(SIGNATURES)

Name of intercept regressor in raw data files.
"""
intercept_name() = :cons;

# All regressors for given fields. As returned to the outside, not as loaded from disk.
group_regressors(ds :: DataSettings, group :: Symbol) =
    regressor_names(group,  1 : n_groups(ds, group));

gpa_regressors(ds :: DataSettings) =
    group_regressors(ds, :gpa);
parental_regressors(ds :: DataSettings) =
    group_regressors(ds, :parental);
quality_regressors(ds :: DataSettings) =
    group_regressors(ds, :quality);
school_regressors(ds :: DataSettings) =
    group_regressors(ds, :school);

## ----------  Names in original files
# collected here because they tend to change over time.

"""
    $(SIGNATURES)

Column headers for regression files. These are the headers that are in the raw files. What is returned is determined by `regressor_name`.
Different from cross-tabs.
"""
function regr_col_header(colCat :: Symbol, j :: Integer)
    if colCat ∈ (:gpa, :afqt)
        hd = "afqt$j";
    elseif colCat == :parental
        hd = "inc_quartile$j";
    elseif colCat == :school
        hd = "school$j";
    elseif colCat == :lastColl
        hd = "last_type$j";
    elseif colCat == :exper
        if j == 1
            hd = "exp";
        else
            hd = "exp$j";
        end
    else
        error("Invalid: $colCat")
    end
    return Symbol(hd)
end

regr_col_headers(colCat :: Symbol, jMax :: Integer) =
    [regr_col_header(colCat, j)  for j = 1 : jMax];

# Valid column headers (without indices) in original files.
# Some files have :afqt; others have :quartile.
valid_col_headers() =
    [:afqt, :quartile, :inc_quartile, :quality, :school, :last_type, :exp, :cons];

# Check that regressor names in a loaded file are valid.
# Before renaming.
function validate_regressor_names(nameV :: AbstractVector{Symbol})
    isValid = true;
    for rName ∈ nameV
        isValid = isValid  &&  validate_regressor_name(rName);
    end
    return isValid
end

function validate_regressor_name(rName :: Symbol)
    isValid = true;
    group, idx = break_regressor_name(rName);
    if !(Symbol(group) ∈ valid_col_headers())
        isValid = false;
        @warn "Invalid column header: $rName"
    end
    return isValid
end

# Break a regressor name of the form `:afqt3` into "afqt" and 3.
# The index is `nothing` if it is missing
function break_regressor_name(rName)
    rStr = string(rName);
    m = match(r"([A-Za-z_]*)([0-9]*)",  rStr);
    group = m.captures[1];
    if isempty(m.captures[2])
        idx = nothing;
    else
        idx = parse(Int, m.captures[2]);
    end
    return group, idx
end

# Rename regressors so that output names match loaded names
# Does nothing if regressor does not exist or does not need to be renamed.
function rename_regressors(rt :: RegressionTable, colCat :: Symbol)
    oldName = regr_col_header(colCat, 2);
    newName = regressor_name(colCat, 2);
    if (oldName != newName)  &&  has_regressor(rt, oldName)
        for j = 1 : 10
            oldName = regr_col_header(colCat, j);
            newName = regressor_name(colCat, j);
            if has_regressor(rt, oldName)
                EconometricsLH.rename_regressor(rt, oldName, newName);
            end
        end
    end
    return nothing
end

# Ensure that all regressors match output regressor names
function rename_regressors(rt :: RegressionTable)
    for colCat ∈ (:afqt, :parental, :school)
        rename_regressors(rt, colCat);
    end
end

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


text:
@doc doc"""
`Q` state Potts model with energy $E = -\sum_{i,j} \delta_{\sigma_i, \sigma_j}$,
where $\sigma_i$ takes an integer value from $1$ to $Q$ and $\delta$ is a Kronecker's delta.
Order parameter (total magnetization) is defined as
\begin{equation}
    M = \frac{Q-1}{Q}N_1 - \frac{1}{Q}(N-N_1),
\end{equation}
where $N$ is the number of sites and $N_1$ is the number of $\sigma=1$ spins.
"""
mutable struct Potts <: Model
    lat :: Lattice
    Q :: Int
    spins :: Matrix{Int}
    rng :: Random.MersenneTwister

    function Potts(lat::Lattice, Q::Integer, rng::Random.AbstractRNG)
        if Q < 2
            error("Q should be 2 or more: Q = $Q")
        end
        spins = rand(rng, 1:Q, 1, numsites(lat))
        return new(lat, Q, spins, rng)
    end
end

Potts(lat::Lattice, Q::Integer) = Potts(lat, Q, Random.seed!(Random.MersenneTwister(0)))
Potts(lat::Lattice, Q::Integer, seed) = Potts(lat, Q, Random.seed!(Random.MersenneTwister(0), seed...))

@doc doc"""
    Potts(param)

Generates `Potts` using `param["Lattice"]`, `param["Q"]`, and `param["Seed"]` (if defined).
Each spin will be initialized randomly and independently.
"""
function Potts(param::Parameter)
    lat = generatelattice(param)
    Q = param["Q"]
    if "Seed" in keys(param)
        return Potts(lat, Q, param["Seed"])
    else
        return Potts(lat, Q)
    end
end

include("update.jl")
include("estimator.jl")
include("postproc.jl")


text: hull2array(H) = [collect(H[i][1] for i in 1:length(H)) collect(H[i][2] for i in 1:length(H))]

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
module MultivariateStats
    using LinearAlgebra
    using StatsBase: SimpleCovariance, CovarianceEstimator
    import Statistics: mean, var, cov, covm
    import Base: length, size, show, dump
    import StatsBase: fit, predict, ConvergenceException
    import SparseArrays
    import LinearAlgebra: eigvals

    export

    ## common
    evaluate,           # evaluate discriminant function values (imported from Base)
    predict,            # use a model to predict responses (imported from StatsBase)
    fit,                # fit a model to data (imported from StatsBase)
    centralize,         # subtract a mean vector from each column
    decentralize,       # add a mean vector to each column
    indim,              # the input dimension of a model
    outdim,             # the output dimension of a model
    projection,         # the projection matrix
    reconstruct,        # reconstruct the input (approximately) given the output
    transform,          # apply a model to transform a vector or a matrix

    # lreg
    llsq,               # Linear Least Square regression
    ridge,              # Ridge regression

    # whiten
    Whitening,          # Type: Whitening transformation

    invsqrtm,           # Compute inverse of matrix square root, i.e. inv(sqrtm(A))
    invsqrtm!,          # Compute inverse of matrix square root inplace
    cov_whitening,      # Compute a whitening transform based on covariance
    cov_whitening!,     # Compute a whitening transform based on covariance (input will be overwritten)
    invsqrtm,           # Compute C^{-1/2}, i.e. inv(sqrtm(C))

    ## pca
    PCA,                # Type: Principal Component Analysis model

    pcacov,             # PCA based on covariance
    pcasvd,             # PCA based on singular value decomposition of input data
    principalratio,     # the ratio of variances preserved in the principal subspace
    principalvar,       # the variance along a specific principal direction
    principalvars,      # the variances along all principal directions

    tprincipalvar,      # total principal variance, i.e. sum(principalvars(M))
    tresidualvar,       # total residual variance
    tvar,               # total variance

    ## ppca
    PPCA,               # Type: the Probabilistic PCA model

    ppcaml,             # Maximum likelihood probabilistic PCA
    ppcaem,             # EM algorithm for probabilistic PCA
    bayespca,           # Bayesian PCA
    loadings,           # factor loadings matrix

    ## kpca
    KernelPCA,          # Type: the Kernel PCA model

    ## cca
    CCA,                # Type: Correlation Component Analysis model

    ccacov,             # CCA based on covariances
    ccasvd,             # CCA based on singular value decomposition of input data

    xindim,             # input dimension of X
    yindim,             # input dimension of Y
    xmean,              # sample mean of X
    ymean,              # sample mean of Y
    xprojection,        # projection matrix for X
    yprojection,        # projection matrix for Y
    xtransform,         # transform for X
    ytransform,         # transform for Y
    correlations,       # correlations of all projected directions

    ## cmds
    MDS,
    classical_mds,      # perform classical MDS over a given distance matrix
    eigvals,            # eignenvalues of the transformation
    stress,             # stress evaluation

    gram2dmat, gram2dmat!,  # Gram matrix => Distance matrix
    dmat2gram, dmat2gram!,  # Distance matrix => Gram matrix

    ## lda
    Discriminant,           # Abstract Type: for all discriminant functionals
    LinearDiscriminant,     # Type: Linear Discriminant functional
    MulticlassLDAStats,     # Type: Statistics required for training multi-class LDA
    MulticlassLDA,          # Type: Multi-class LDA model
    SubspaceLDA,            # Type: LDA model for high-dimensional spaces

    ldacov,                 # Linear discriminant analysis based on covariances

    classweights,           # class-specific weights
    classmeans,             # class-specific means
    withclass_scatter,      # with-class scatter matrix
    betweenclass_scatter,   # between-class scatter matrix
    multiclass_lda_stats,   # compute statistics for multiclass LDA training
    multiclass_lda,         # train multi-class LDA based on statistics
    mclda_solve,            # solve multi-class LDA projection given scatter matrices
    mclda_solve!,           # solve multi-class LDA projection (inputs are overriden)

    ## ica
    ICA,                    # Type: the Fast ICA model

    icagfun,                # a function to get a ICA approx neg-entropy functor
    fastica!,               # core algorithm function for the Fast ICA

    ## fa
    FactorAnalysis,         # Type: the Factor Analysis model

    faem,                   # Maximum likelihood probabilistic PCA
    facm                    # EM algorithm for probabilistic PCA

    ## source files
    include("common.jl")
    include("lreg.jl")
    include("whiten.jl")
    include("pca.jl")
    include("ppca.jl")
    include("kpca.jl")
    include("cca.jl")
    include("cmds.jl")
    include("lda.jl")
    include("ica.jl")
    include("fa.jl")

end # module