text:
"""
`prior(m::AbstractModel{T})`

Calculates log joint prior density of m.parameters.
"""
function prior{T<:AbstractFloat}(m::AbstractModel{T})
    free_params = Base.filter(θ -> !θ.fixed, m.parameters)
    logpdfs = map(logpdf, free_params)
    return sum(logpdfs)
end

"""
```
posterior{T<:AbstractFloat}(m::AbstractModel{T}, data::Matrix{T};
                            mh::Bool = false, catch_errors::Bool = false)
```

Calculates and returns the log of the posterior distribution for `m.parameters`:

```
log posterior  = log likelihood + log prior + const
log Pr(Θ|data) = log Pr(data|Θ) + log Pr(Θ) + const
```

### Arguments

- `m`: the model object
- `data`: matrix of data for observables

### Optional Arguments

- `mh`: Whether `metropolis_hastings` is the caller. If `mh = true`,
  `catch_errors` is set to `true` (see below)
- `catch_errors`: Whether or not to catch errors of type `GensysError` and
  `ParamBoundsError`
"""
function posterior{T<:AbstractFloat}(m::AbstractModel{T},
                                     data::Matrix{T};
                                     mh::Bool = false,
                                     catch_errors::Bool = false)
    catch_errors = catch_errors || mh
    post = likelihood(m, data; mh = mh, catch_errors = catch_errors) + prior(m)
    return post
end

"""
```
posterior!{T<:AbstractFloat}(m::AbstractModel{T}, parameters::Vector{T}, data::Matrix{T};
                             mh::Bool = false, catch_errors::Bool = false)
```

Evaluates the log posterior density at `parameters`.

### Arguments

- `m`: The model object
- `parameters`: New values for the model parameters
- `data`: Matrix of input data for observables

### Optional Arguments

- `mh`: Whether `metropolis_hastings` is the caller. If `mh = true`,
  `catch_errors` is set to `true` (see below)
- `catch_errors`: Whether or not to catch errors of type `GensysError` and
  `ParamBoundsError`
"""
function posterior!{T<:AbstractFloat}(m::AbstractModel{T},
                                      parameters::Vector{T},
                                      data::Matrix{T};
                                      mh::Bool = false,
                                      catch_errors::Bool = false)
    catch_errors = catch_errors || mh
    if mh
        try
            update!(m, parameters)
        catch err
            if isa(err, ParamBoundsError)
                return -Inf
            else
                throw(err)
            end
        end
    else
        update!(m, parameters)
    end
    return posterior(m, data; mh=mh, catch_errors=catch_errors)

end

"""
```
likelihood{T<:AbstractFloat}(m::AbstractModel, data::Matrix{T};
                             mh::Bool = false, catch_errors::Bool = false)
```

Evaluate the DSGE likelihood function. Can handle two-part estimation where the observed
sample contains both a normal stretch of time (in which interest rates are positive) and
a stretch of time in which interest rates reach the zero lower bound. If there is a
zero-lower-bound period, then we filter over the 2 periods separately. Otherwise, we
filter over the main sample all at once.

### Arguments

- `m`: The model object
- `data`: matrix of data for observables

### Optional Arguments

- `mh`: Whether `metropolis_hastings` is the caller. If `mh = true`,
  `catch_errors` is set to `true` (see below)
- `catch_errors`: Whether or not to catch errors of type `GensysError`
"""
function likelihood{T<:AbstractFloat}(m::AbstractModel,
                                      data::Matrix{T};
                                      mh::Bool = false,
                                      catch_errors::Bool = false)
    catch_errors = catch_errors || mh

    # During Metropolis-Hastings, return -∞ if any parameters are not within their bounds
    if mh
        for θ in m.parameters
            (left, right) = θ.valuebounds
            if !θ.fixed && !(left <= θ.value <= right)
                return -Inf
            end
        end
    end

    # Compute state-space system
    system = try
        compute_system(m)
    catch err
        if catch_errors && isa(err, GensysError)
            return -Inf
        else
            rethrow(err)
        end
    end

    # Return total log-likelihood, excluding the presample
    try
        kal = filter(m, data, system; outputs = [:loglh], include_presample = false)
        return kal[:total_loglh]
    catch err
        if catch_errors && isa(err, DomainError)
            warn("Log of incremental likelihood is negative; returning -Inf")
            return -Inf
        else
            rethrow(err)
        end
    end
end


text:
function TabularClassificationSingle(
        blocks::Tuple{<:TableRow, <:Label},
        data)
    tabledata, targetdata = data
    tabledata isa TableDataset || error("`data` needs to be a tuple of a `TableDataset` and targets")

    return SupervisedTask(
        blocks,
        (
            setup(TabularPreprocessing, blocks[1], tabledata),
            OneHot()
        )
    )
end

"""
    TabularClassificationSingle(blocks, data)

Learning task for single-label tabular classification. Continuous columns are
normalized and missing values are filled, categorical columns are label encoded
taking into account any missing values which might be present. The target value
is predicted from `classes`. `blocks` should be an input and target block
`(TableRow(...), Label(...))`.

    TabularClassificationSingle(classes, tabledata [; catcols, contcols])

Construct learning task with `classes` to classify into and a `TableDataset`
`tabledata`. The column names can be passed in or guessed from the data.
"""
function TabularClassificationSingle(
        classes::AbstractVector,
        tabledata::TableDataset;
        catcols = nothing,
        contcols = nothing)

    blocks = (
        setup(TableRow, tabledata; catcols = catcols, contcols = contcols),
        Label(classes)
    )
    return TabularClassificationSingle(blocks, (tabledata, nothing))
end

# ## Tests

@testset "TabularClassificationSingle [task]" begin
    df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"], C = ["P", "F", "P", "F"])
    td = TableDataset(df)

    task = TabularClassificationSingle(["P", "F"], td; catcols=(:B,), contcols=(:A,), )
    testencoding(getencodings(task), getblocks(task).sample)
    FastAI.checktask_core(task)
    @test_nowarn tasklossfn(task)
    @test_nowarn taskmodel(task)

    @testset "`encodeinput`" begin
        row = mockblock(getblocks(task)[1])

        xtrain = encodeinput(task, Training(), row)
        @test length(xtrain[1]) == length(getblocks(task).input.catcols)
        @test length(xtrain[2]) == length(getblocks(task).input.contcols)

        @test eltype(xtrain[1]) <: Number
    end

    @testset "`encodetarget`" begin
        category = "P"
        y = encodetarget(task, Training(), category)
        @test y ≈ [1, 0]
    end

    @test_nowarn TabularClassificationSingle(["P", "F"], td)
    @test TabularClassificationSingle(["P", "F"], td).blocks[1].catcols == (:B, :C)

    FastAI.test_task_show(task, ShowText(Base.DevNull()))
end

ext: jl


text:
import Base.+, Base.*

# Basic idea: Data structure to store a differential 1-form e.g. a contact form
struct OneForm
	# mapping from R^3 to R^3 that represents a contact from. Input are three arguments and a single 3-component vector should be returned
	rep
end

# the names should tell you everything
dx = OneForm((x,y,z) -> [1 0 0])
dy = OneForm((x,y,z) -> [0 1 0])
dz = OneForm((x,y,z) -> [0 0 1])

# addition of two contact forms
function (+)(a::OneForm, b::OneForm)
	return OneForm((x, y, z) -> a.rep(x,y,z) + b.rep(x,y,z))
end

# addition of a real number and a contact form
function (+)(a::Real, b::OneForm)
	return OneForm((x, y, z) -> a + b.rep(x,y,z))
end

# addition of a real number and a contact form
function (+)(a::OneForm, b::Real)
	return OneForm((x, y, z) -> a.rep(x,y,z) + b)
end

# addition of a real matrix and a contact form
function (+)(a::Matrix, b::OneForm)
	return OneForm((x, y, z) -> a + b.rep(x,y,z))
end

# addition of a real matrix and a contact form
function (+)(a::OneForm, b::Matrix)
	return OneForm((x, y, z) -> a.rep(x,y,z) + b)
end

# multiplication of a R^3 -> R function and a contact form
function (*)(f::Function, a::OneForm)
	return OneForm((x,y,z) -> f(x,y,z) * a.rep(x,y,z))
end

# multiplication of a real number and a contact form
function (*)(a::Real, b::OneForm)
	return OneForm((x,y,z) -> a * b.rep(x,y,z))
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

using Plots
using ForwardDiff
using Combinatorics

# # Generate 1 soliton

# +
v=1.
c=3.
θ = - 2c/√v

function u(x,t)
    v/2 * (sech(√v / 2 * (x-v*t-θ)))^2
end

# +
using Plots

anim = @animate for t in 0.:0.1:10
    x=collect(-10:0.01:10)
    y=u.(x,t)
    plot(x,y,ylim=[0,1.])
end
# -

gif(anim)

# # Generate 2 soliton

# +
p = [1.,2.]
θ = [0, 1]
function τ(x,t)
    η=[
        p[1]*x - p[1]^3 * t - θ[1],
        p[2]*x - p[2]^3 * t - θ[2]
    ]
    A12 = ((p[1]-p[2])/(p[1]+p[2]))^2
    return 1 + exp(η[1])+exp(η[2]) + A12*exp(η[1]+η[2])
end

logτ(x,t)=log(τ(x,t))
logτ(xt) = logτ(xt[1], xt[2])
hessian = xt -> ForwardDiff.hessian(logτ,xt)
u(x,t)=hessian([x,t])[1]
anim = @animate for t in 0:0.05:10
    x=-10:0.1:10
    y=u.(x,t)
    plot(x,y,ylim=[0,3.])
end

gif(anim)
# -

# # n-soliton n =3

# +
θ = [0.0,5.0,-15.0]
p = [1.5,2.0,2.5]
N = length(p)

A = zeros(N,N)
for (i, j) in combinations(1:N,2)
    A[i,j] = ((p[i]-p[j])/(p[i]+p[j]))^2 |> log
end

function τ(x,t)
    η = @. p*x - p^3*t - θ
    ret = 0.
    for μ in Iterators.product(repeat([[0,1]],N)...)
        s1 = sum(μ .* η)
        s2 = 0.
        for (i,j) in combinations(1:N,2)
            s2= μ[i]*μ[j]*A[i,j]
        end
        ret += exp(s1 + s2)
    end
    ret
end

logτ(x,t)=log(τ(x,t))
logτ(xt) = logτ(xt[1], xt[2])
hessian = xt -> ForwardDiff.hessian(logτ,xt)
u(x,t)=hessian([x,t])[1]
anim = @animate for t in 0:0.05:10
    x=-10:0.1:10
    y=u.(x,t)
    plot(x,y,ylim=[0,3.],label="3-soliton")
end

gif(anim)
# -

# # Reference
#
# - 箱玉系の数理


text:
############# Problem Set 1 #############
############# Giuseppe Grasso, Econ PhD student at LISER/University of Luxembourg (https://liser.elsevierpure.com/en/persons/giuseppe-grasso)
println("Hello, this is my first Julia script!")

############# 1. Initializing variables and practice with basic matrix operations

#1.a
using Random
Random.seed!(1234)
using Distributions
A=rand(Uniform(-5,10),10,7)
B=rand(Normal(-2,15),10,7)
using LinearAlgebra
C=hcat(A[1:5,1:5],B[1:5,6:7])
#currently giving up on D

#1.b
length(A)

#1.c skipped

#1.d
reshape(B,10,7)==
# too late to start this assignment. will work on it later at some point


text:
using Test, DiscreteExteriorCalculus
const DEC = DiscreteExteriorCalculus

@testset "parametrize_subspace" begin
    M, v = [[0,1] [-1,1] [0,1]], [-.5,1]
    a, b = DEC.parametrize_subspace(M, v)
    @test M * a ≈ v
    @test M * (a + b * [1]) ≈ v
    @test M * (a + b * [-2]) ≈ v

    M, v = [[.5, 1] [-.5, 1] [-.5, 1]], [0.0, 1]
    a, b = DEC.parametrize_subspace(M, v)
    @test M * a ≈ v
    @test M * (a + b * [1]) ≈ v
    @test M * (a + b * [-2]) ≈ v

    M, v = [[.5, 1] [-.5, 1]], [0.0, 1]
    a, b = DEC.parametrize_subspace(M, v)
    @test size(b) == (2,0)
    @test M * a ≈ v
end

@testset "sign_of_permutation" begin
    @test DEC.sign_of_permutation([1],[1]) == 1
    @test DEC.sign_of_permutation([1,2],[2,1]) == -1
    @test DEC.sign_of_permutation([1,2],[1,2]) == 1
end


text:
module LinearWrapper

using JavaCall

const AbstractFieldMatrix = @jimport org.hipparchus.linear.AbstractFieldMatrix
const AbstractRealMatrix = @jimport org.hipparchus.linear.AbstractRealMatrix
const AnyMatrix = @jimport org.hipparchus.linear.AnyMatrix
const Array2DRowFieldMatrix = @jimport org.hipparchus.linear.Array2DRowFieldMatrix
const Array2DRowRealMatrix = @jimport org.hipparchus.linear.Array2DRowRealMatrix
const ArrayFieldVector = @jimport org.hipparchus.linear.ArrayFieldVector
const ArrayRealVector = @jimport org.hipparchus.linear.ArrayRealVector
const BiDiagonalTransformer = @jimport org.hipparchus.linear.BiDiagonalTransformer
const BlockFieldMatrix = @jimport org.hipparchus.linear.BlockFieldMatrix
const BlockRealMatrix = @jimport org.hipparchus.linear.BlockRealMatrix
const CholeskyDecomposer = @jimport org.hipparchus.linear.CholeskyDecomposer
const CholeskyDecomposition = @jimport org.hipparchus.linear.CholeskyDecomposition
const Class = @jimport java.lang.Class
const ComplexEigenDecomposition = @jimport org.hipparchus.linear.ComplexEigenDecomposition
const ConjugateGradient = @jimport org.hipparchus.linear.ConjugateGradient
const DecompositionSolver = @jimport org.hipparchus.linear.DecompositionSolver
const DefaultFieldMatrixChangingVisitor = @jimport org.hipparchus.linear.DefaultFieldMatrixChangingVisitor
const DefaultFieldMatrixPreservingVisitor = @jimport org.hipparchus.linear.DefaultFieldMatrixPreservingVisitor
const DefaultIterativeLinearSolverEvent = @jimport org.hipparchus.linear.DefaultIterativeLinearSolverEvent
const DefaultRealMatrixChangingVisitor = @jimport org.hipparchus.linear.DefaultRealMatrixChangingVisitor
const DefaultRealMatrixPreservingVisitor = @jimport org.hipparchus.linear.DefaultRealMatrixPreservingVisitor
const DiagonalMatrix = @jimport org.hipparchus.linear.DiagonalMatrix
const Double = @jimport java.lang.Double
const EigenDecomposition = @jimport org.hipparchus.linear.EigenDecomposition
const EventObject = @jimport java.util.EventObject
const Field = @jimport org.hipparchus.Field
const FieldDecompositionSolver = @jimport org.hipparchus.linear.FieldDecompositionSolver
const FieldElement = @jimport org.hipparchus.FieldElement
const FieldLUDecomposition = @jimport org.hipparchus.linear.FieldLUDecomposition
const FieldMatrix = @jimport org.hipparchus.linear.FieldMatrix
const FieldMatrixChangingVisitor = @jimport org.hipparchus.linear.FieldMatrixChangingVisitor
const FieldMatrixPreservingVisitor = @jimport org.hipparchus.linear.FieldMatrixPreservingVisitor
const FieldPosition = @jimport java.text.FieldPosition
const FieldQRDecomposition = @jimport org.hipparchus.linear.FieldQRDecomposition
const FieldVector = @jimport org.hipparchus.linear.FieldVector
const FieldVectorChangingVisitor = @jimport org.hipparchus.linear.FieldVectorChangingVisitor
const FieldVectorPreservingVisitor = @jimport org.hipparchus.linear.FieldVectorPreservingVisitor
const Function = @jimport java.util.function.Function
const HessenbergTransformer = @jimport org.hipparchus.linear.HessenbergTransformer
const IterationEvent = @jimport org.hipparchus.util.IterationEvent
const IterationManager = @jimport org.hipparchus.util.IterationManager
const IterativeLinearSolver = @jimport org.hipparchus.linear.IterativeLinearSolver
const IterativeLinearSolverEvent = @jimport org.hipparchus.linear.IterativeLinearSolverEvent
const JComplex = @jimport org.hipparchus.complex.Complex
const JIterator = @jimport java.util.Iterator
const JString = @jimport java.lang.String
const JacobiPreconditioner = @jimport org.hipparchus.linear.JacobiPreconditioner
const LUDecomposer = @jimport org.hipparchus.linear.LUDecomposer
const LUDecomposition = @jimport org.hipparchus.linear.LUDecomposition
const Locale = @jimport java.util.Locale
const MatrixDecomposer = @jimport org.hipparchus.linear.MatrixDecomposer
const MatrixUtils = @jimport org.hipparchus.linear.MatrixUtils
const NumberFormat = @jimport java.text.NumberFormat
const Object = @jimport java.lang.Object
const ObjectInputStream = @jimport java.io.ObjectInputStream
const ObjectOutputStream = @jimport java.io.ObjectOutputStream
const OpenMapRealMatrix = @jimport org.hipparchus.linear.OpenMapRealMatrix
const OpenMapRealVector = @jimport org.hipparchus.linear.OpenMapRealVector
const OrderedComplexEigenDecomposition = @jimport org.hipparchus.linear.OrderedComplexEigenDecomposition
const OrderedEigenDecomposition = @jimport org.hipparchus.linear.OrderedEigenDecomposition
const ParsePosition = @jimport java.text.ParsePosition
const PreconditionedIterativeLinearSolver = @jimport org.hipparchus.linear.PreconditionedIterativeLinearSolver
const QRDecomposer = @jimport org.hipparchus.linear.QRDecomposer
const QRDecomposition = @jimport org.hipparchus.linear.QRDecomposition
const RRQRDecomposition = @jimport org.hipparchus.linear.RRQRDecomposition
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RealLinearOperator = @jimport org.hipparchus.linear.RealLinearOperator
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const RealMatrixChangingVisitor = @jimport org.hipparchus.linear.RealMatrixChangingVisitor
const RealMatrixFormat = @jimport org.hipparchus.linear.RealMatrixFormat
const RealMatrixPreservingVisitor = @jimport org.hipparchus.linear.RealMatrixPreservingVisitor
const RealVector = @jimport org.hipparchus.linear.RealVector
const RealVectorChangingVisitor = @jimport org.hipparchus.linear.RealVectorChangingVisitor
const RealVectorFormat = @jimport org.hipparchus.linear.RealVectorFormat
const RealVectorPreservingVisitor = @jimport org.hipparchus.linear.RealVectorPreservingVisitor
const RectangularCholeskyDecomposition = @jimport org.hipparchus.linear.RectangularCholeskyDecomposition
const RiccatiEquationSolver = @jimport org.hipparchus.linear.RiccatiEquationSolver
const RiccatiEquationSolverImpl = @jimport org.hipparchus.linear.RiccatiEquationSolverImpl
const SchurTransformer = @jimport org.hipparchus.linear.SchurTransformer
const SingularValueDecomposer = @jimport org.hipparchus.linear.SingularValueDecomposer
const SingularValueDecomposition = @jimport org.hipparchus.linear.SingularValueDecomposition
const SparseFieldMatrix = @jimport org.hipparchus.linear.SparseFieldMatrix
const SparseFieldVector = @jimport org.hipparchus.linear.SparseFieldVector
const SparseRealVector = @jimport org.hipparchus.linear.SparseRealVector
const StringBuffer = @jimport java.lang.StringBuffer
const SymmLQ = @jimport org.hipparchus.linear.SymmLQ
const TriDiagonalTransformer = @jimport org.hipparchus.linear.TriDiagonalTransformer
const UnivariateFunction = @jimport org.hipparchus.analysis.UnivariateFunction

export AbstractFieldMatrix
export AbstractRealMatrix
export AnyMatrix
export Array2DRowFieldMatrix
export Array2DRowRealMatrix
export ArrayFieldVector
export ArrayRealVector
export BiDiagonalTransformer
export BlockFieldMatrix
export BlockRealMatrix
export CholeskyDecomposer
export CholeskyDecomposition
export ComplexEigenDecomposition
export ConjugateGradient
export DecompositionSolver
export DefaultFieldMatrixChangingVisitor
export DefaultFieldMatrixPreservingVisitor
export DefaultIterativeLinearSolverEvent
export DefaultRealMatrixChangingVisitor
export DefaultRealMatrixPreservingVisitor
export DiagonalMatrix
export EigenDecomposition
export FieldDecompositionSolver
export FieldLUDecomposition
export FieldMatrix
export FieldMatrixChangingVisitor
export FieldMatrixPreservingVisitor
export FieldQRDecomposition
export FieldVector
export FieldVectorChangingVisitor
export FieldVectorPreservingVisitor
export HessenbergTransformer
export IterativeLinearSolver
export IterativeLinearSolverEvent
export JacobiPreconditioner
export LUDecomposer
export LUDecomposition
export MatrixDecomposer
export MatrixUtils
export OpenMapRealMatrix
export OpenMapRealVector
export OrderedComplexEigenDecomposition
export OrderedEigenDecomposition
export PreconditionedIterativeLinearSolver
export QRDecomposer
export QRDecomposition
export RRQRDecomposition
export RealLinearOperator
export RealMatrix
export RealMatrixChangingVisitor
export RealMatrixFormat
export RealMatrixPreservingVisitor
export RealVector
export RealVectorChangingVisitor
export RealVectorFormat
export RealVectorPreservingVisitor
export RectangularCholeskyDecomposition
export RiccatiEquationSolver
export RiccatiEquationSolverImpl
export SchurTransformer
export SingularValueDecomposer
export SingularValueDecomposition
export SparseFieldMatrix
export SparseFieldVector
export SparseRealMatrix
export SparseRealVector
export SymmLQ
export TriDiagonalTransformer
export _end
export add
export add_to_entry
export append
export big_fraction_matrix_to_real_matrix
export block_inverse
export check_addition_compatible
export check_column_index
export check_matrix_index
export check_multiplication_compatible
export check_row_index
export check_same_column_dimension
export check_same_row_dimension
export check_sub_matrix_index
export check_subtraction_compatible
export check_symmetric
export combine
export combine_to_self
export copy
export copy_sub_matrix
export cosine
export create
export create_blocks_layout
export create_column_field_matrix
export create_column_real_matrix
export create_field_diagonal_matrix
export create_field_identity_matrix
export create_field_matrix
export create_field_vector
export create_matrix
export create_real_diagonal_matrix
export create_real_identity_matrix
export create_real_matrix
export create_real_vector
export create_row_field_matrix
export create_row_real_matrix
export decompose
export deserialize_real_matrix
export deserialize_real_vector
export dot_product
export ebe_divide
export ebe_multiply
export equals
export format
export fraction_matrix_to_real_matrix
export get_available_locales
export get_b
export get_class
export get_column
export get_column_dimension
export get_column_matrix
export get_column_separator
export get_column_vector
export get_condition_number
export get_covariance
export get_d
export get_data
export get_data_ref
export get_determinant
export get_dimension
export get_distance
export get_eigenvalues
export get_eigenvector
export get_entry
export get_field
export get_format
export get_frobenius_norm
export get_h
export get_imag_eigenvalue
export get_imag_eigenvalues
export get_inverse
export get_inverse_condition_number
export get_iteration_manager
export get_iterations
export get_k
export get_l
export get_l1_distance
export get_l1_norm
export get_l_inf_distance
export get_l_inf_norm
export get_lt
export get_max_index
export get_max_value
export get_min_index
export get_min_value
export get_norm
export get_norm1
export get_norm_infty
export get_norm_of_residual
export get_p
export get_pivot
export get_prefix
export get_pt
export get_q
export get_qt
export get_r
export get_rank
export get_real_eigenvalue
export get_real_eigenvalues
export get_real_matrix_format
export get_real_vector_format
export get_residual
export get_right_hand_side_vector
export get_root_matrix
export get_row
export get_row_dimension
export get_row_matrix
export get_row_prefix
export get_row_separator
export get_row_suffix
export get_row_vector
export get_s
export get_separator
export get_singular_values
export get_solution
export get_solver
export get_source
export get_sparsity
export get_square_root
export get_sub_matrix
export get_sub_vector
export get_suffix
export get_t
export get_trace
export get_u
export get_ut
export get_v
export get_vt
export has_complex_eigenvalues
export hash_code
export inverse
export is_infinite
export is_na_n
export is_non_singular
export is_singular
export is_square
export is_symmetric
export is_transposable
export iterator
export kronecker_product
export map
export map_add
export map_add_to_self
export map_divide
export map_divide_to_self
export map_inv
export map_inv_to_self
export map_multiply
export map_multiply_to_self
export map_subtract
export map_subtract_to_self
export map_to_self
export matrix_exponential
export multiply
export multiply_entry
export multiply_transposed
export notify
export notify_all
export operate
export operate_transpose
export outer_product
export parse
export power
export pre_multiply
export projection
export provides_residual
export scalar_add
export scalar_multiply
export serialize_real_matrix
export serialize_real_vector
export set
export set_column
export set_column_matrix
export set_column_vector
export set_entry
export set_row
export set_row_matrix
export set_row_vector
export set_sub_matrix
export set_sub_vector
export should_check
export solve
export solve_in_place
export solve_lower_triangular_system
export solve_upper_triangular_system
export sparse_iterator
export sqrt
export stack
export start
export subtract
export to_array
export to_blocks_layout
export to_string
export transpose
export transpose_multiply
export unit_vector
export unitize
export unmodifiable_real_vector
export unstack_square
export visit
export wait
export walk_in_column_order
export walk_in_default_order
export walk_in_optimized_order
export walk_in_row_order

include("abstract_field_matrix.jl")
include("abstract_real_matrix.jl")
include("any_matrix.jl")
include("array2_d_row_field_matrix.jl")
include("array2_d_row_real_matrix.jl")
include("array_field_vector.jl")
include("array_real_vector.jl")
include("bi_diagonal_transformer.jl")
include("block_field_matrix.jl")
include("block_real_matrix.jl")
include("cholesky_decomposer.jl")
include("cholesky_decomposition.jl")
include("complex_eigen_decomposition.jl")
include("conjugate_gradient.jl")
include("decomposition_solver.jl")
include("default_field_matrix_changing_visitor.jl")
include("default_field_matrix_preserving_visitor.jl")
include("default_iterative_linear_solver_event.jl")
include("default_real_matrix_changing_visitor.jl")
include("default_real_matrix_preserving_visitor.jl")
include("diagonal_matrix.jl")
include("eigen_decomposition.jl")
include("field_decomposition_solver.jl")
include("field_lu_decomposition.jl")
include("field_matrix.jl")
include("field_matrix_changing_visitor.jl")
include("field_matrix_preserving_visitor.jl")
include("field_qr_decomposition.jl")
include("field_vector.jl")
include("field_vector_changing_visitor.jl")
include("field_vector_preserving_visitor.jl")
include("hessenberg_transformer.jl")
include("iterative_linear_solver.jl")
include("iterative_linear_solver_event.jl")
include("jacobi_preconditioner.jl")
include("lu_decomposer.jl")
include("lu_decomposition.jl")
include("matrix_decomposer.jl")
include("matrix_utils.jl")
include("open_map_real_matrix.jl")
include("open_map_real_vector.jl")
include("ordered_complex_eigen_decomposition.jl")
include("ordered_eigen_decomposition.jl")
include("preconditioned_iterative_linear_solver.jl")
include("qr_decomposer.jl")
include("qr_decomposition.jl")
include("real_linear_operator.jl")
include("real_matrix.jl")
include("real_matrix_changing_visitor.jl")
include("real_matrix_format.jl")
include("real_matrix_preserving_visitor.jl")
include("real_vector.jl")
include("real_vector_changing_visitor.jl")
include("real_vector_format.jl")
include("real_vector_preserving_visitor.jl")
include("rectangular_cholesky_decomposition.jl")
include("riccati_equation_solver.jl")
include("riccati_equation_solver_impl.jl")
include("rrqr_decomposition.jl")
include("schur_transformer.jl")
include("singular_value_decomposer.jl")
include("singular_value_decomposition.jl")
include("sparse_field_matrix.jl")
include("sparse_field_vector.jl")
include("sparse_real_matrix.jl")
include("sparse_real_vector.jl")
include("symm_l_q.jl")
include("tri_diagonal_transformer.jl")

end


text:
#unum-unum.jl

#contains information about the unum type and helper functions directly related to constructor.

#the unum type is an abstract type.  We'll be overloading the call function later
#so we can do "pseudo-constructions" on this type.
doc"""
`Unum{ESS,FSS}` creates a Unum with esizesize ESS and fsizesize FSS.

NB:  Internally this may cast to a different Unum type (UnumLarge or UnumSmall)
for performance purposes.  The `Unum{ESS,FSS}(...)` constructor is always safe
and purposes where computational performance is an issue, use of the internal
types is recommended.
"""
abstract Unum{ESS, FSS} <: Real
export Unum

#general parameter checking for unums.
function __general_unum_check(ESS, FSS, exponent::UInt64, esize::UInt16, fsize::UInt16)
  (ESS > 6) && throw(ArgumentError("ESS == $ESS > 6 disallowed in current implementation"))
  (FSS > 11) && throw(ArgumentError("FSS == $FSS > 11 disallowed in current implementation"))
  _mfs = max_fsize(FSS)
  (fsize > _mfs) && throw(ArgumentError("fsize == $fsize > $_mfs maximum for FSS == $FSS."))
  _mes = max_esize(ESS)
  (esize > _mes) && throw(ArgumentError("esize == $esize > $_mes maximum for ESS == $ESS."))
  _mbe = max_biased_exponent(esize)
  (exponent > _mbe) && throw(ArgumentError("exponent == $exponent > $_mbe maximum for esize == $esize."))
  nothing
end

doc"""
`Unums.UnumSmall{ESS,FSS}` is the internal type for Unums with FSS < 7.  These
numbers require a single unsigned 64-bit integer to store their fractions.
"""
@dev_check type UnumSmall{ESS, FSS} <: Unum{ESS, FSS}
  exponent::UInt64
  fraction::UInt64
  flags::UInt16
  esize::UInt16
  fsize::UInt16
end

#parameter checking.  The call to this check is autogenerated by the @dev_check macro
function __check_UnumSmall(ESS, FSS, exponent::UInt64, fraction::UInt64, flags::UInt16, esize::UInt16, fsize::UInt16)
  (FSS > 6) && throw(ArgumentError("UnumSmall internal class is inappropriate for FSS == $FSS > 6."))
  __general_unum_check(ESS, FSS, exponent, esize, fsize)
  nothing
end

doc"""
`Unums.UnumLarge{ESS,FSS}` is the internal type for Unums with FSS > 6.  These
numbers require an array of 64-bit integers to store their fractions.
"""
@dev_check type UnumLarge{ESS, FSS} <: Unum{ESS, FSS}
  exponent::UInt64
  fraction::ArrayNum{FSS}
  flags::UInt16
  esize::UInt16
  fsize::UInt16
end

#parameter checking.  The call to this check is autogenerated by the @dev_check macro
function __check_UnumLarge{FSS}(_ESS, _FSS, exponent::UInt64, fraction::ArrayNum{FSS}, flags::UInt16, esize::UInt16, fsize::UInt16)
  _FSS != FSS && throw(ArgumentError("ArrayNum FSS must match Unum FSS"))
  (FSS < 7) && throw(ArgumentError("UnumLarge internal class is inappropriate for FSS == $FSS < 7."))
  __general_unum_check(_ESS, _FSS, exponent, esize, fsize)
  nothing
end

#universal copy constructors
UnumSmall{ESS,FSS}(x::UnumSmall{ESS,FSS}) = UnumSmall{ESS,FSS}(x.exponent, x.fraction, x.flags, x.esize, x.fsize)
UnumLarge{ESS,FSS}(x::UnumLarge{ESS,FSS}) = UnumLarge{ESS,FSS}(x.exponent, copy(x.fraction), x.flags, x.esize, x.fsize)
#overide call to provide copy constructors that look like Unum{ESS,FSS} instead of UnumSmall or UnumLarge
@universal function (::Type{Unum{ESS,FSS}})(x::Unum)
  (FSS < 7) ? (UnumSmall(x)) : (UnumLarge(x))
end
@universal function Base.copy(x::Unum)
  (FSS < 7) ? (UnumSmall(x)) : (UnumLarge(x))
end

#override call to allow direct instantiation using the Unum{ESS,FSS} pseudo-constructor.
#The Unum{ESS,FSS} constructor also does trim operation that sets the ubit correctly.
@generated function (::Type{Unum{ESS,FSS}}){ESS,FSS}(exponent::UInt64, fraction, flags::UInt16, esize::UInt16, fsize::UInt16)
  arrayconversion = :(fraction)
  if (fraction == UInt64)
    utype = UnumSmall
  elseif (fraction == Array{UInt64, 1})
    utype = UnumLarge
    arrayconversion = :(ArrayNum{FSS}(fraction))
  elseif (fraction == ArrayNum{FSS})
    utype = UnumLarge
  else
    throw(ArgumentError("fraction must be UInt64, ArrayNum, or Array{UInt64,1} "))
  end

  quote
     #mask out values outside of the flag range.
     flags &= UNUM_FLAG_MASK
     x = $utype{ESS,FSS}(exponent, $arrayconversion, flags, esize, fsize)
     trim_and_set_ubit!(x, fsize)
     x
  end
end

doc"""
  `unum(Unum{ESS,FSS}, exponent, fraction, esize, fsize, flags)` is a guaranteed-safe
  constructor for unums.  Use this if you don't trust your incoming data.  You can
  also use the easy constructor `unum(Unum{ESS,FSS}, exponent, fraction, flags)`
  passing a signed integer exponent, which automatically creates the size parameters
  in a safe fashion.  `unum(x::Unum{ESS,FSS})`, unlike the equivalent `Unum{ESS,FSS}`
  constructor forms, is a safe copy constructor.
"""
function unum{ESS,FSS}(::Type{Unum{ESS,FSS}}, exponent::UInt64, fraction, flags::UInt16, esize::UInt16, fsize::UInt16)
  if !options[:devmode]
    __general_unum_check(ESS, FSS, exponent, esize, fsize)
    (FSS > 6) && __check_ArrayNum(FSS, fraction)
  end
  Unum{ESS,FSS}(exponent, fraction, flags, esize, fsize)
end
function unum{ESS,FSS}(::Type{Unum{ESS,FSS}}, exponent::Int64, fraction, flags::UInt16)
  (esize, exp) = encode_exp(exponent)
  fsize = __minimum_data_width(fraction)
  unum(Unum{ESS,FSS}, exp, fraction, flags, esize, fsize)
end
@universal function unum(x::Unum)
  __general_unum_check(ESS, FSS, x.exponent, x.esize, x.fsize)
  (FSS > 6) && __check_ArrayNum(x.fraction)
  (FSS < 7) ? UnumSmall{ESS,FSS}(x.exponent, x.fraction, x.flags, x.esize, x.fsize) : UnumLarge{ESS,FSS}(x.exponent, x.fraction, x.flags, x.esize, x.fsize)
end

doc"""
  `Unums.buildunum(::Type{[UnumLarge|UnumSmall]}, exponent::Int64, fraction, flags, fsize)` builds a
  unum out of the respective components.  `buildunum` is intended for internal
  use and instead of throwing errors when the values are outside of the standard
  ranges, it will create constants mmr/sss.
"""
@universal function buildunum(T::Type{Unum}, true_exponent::Int64, fraction, flags::UInt16, fsize::UInt16)

  #determine properties of the destination type.
  min_exp_subnormal = min_exponent(ESS, FSS)
  min_exp_normal = min_exponent(ESS)
  max_exp = max_exponent(ESS)
  mfsize = max_fsize(FSS)

  #eject easy candidates.
  (true_exponent < min_exp_subnormal) && return sss(U, flags & UNUM_SIGN_MASK)
  (true_exponent > max_exp) && return mmr(U, flags & UNUM_SIGN_MASK)

  #check if we need to make this a subnormal number.
  if (true_exponent < min_exp_normal)
    #figure the needed difference.
    shiftvalue = to16(min_exp_normal - true_exponent)
    #load the values into a holding unum.
    result = U(z64, fraction, flags, max_esize(ESS), min(fsize, mfsize))
    rsh_and_set_ubit!(result, shiftvalue)
    frac_set_bit!(result, shiftvalue)
  else
    shiftvalue = z16
    (esize, exponent) = encode_exp(true_exponent)
    result = U(exponent, fraction, flags, esize, min(fsize, mfsize))
  end

  trim_and_set_ubit!(result, min(fsize + shiftvalue, mfsize))
  exact_trim!(result)
  result
end

export unum

#a blank constructor uses the default environment setting
function (::Type{Unum})(args...)
  environment()(args...)
end

#masks for the unum flags variable.
const UNUM_SIGN_MASK = UInt16(0x0002)
const UNUM_UBIT_MASK = UInt16(0x0001)
const UNUM_FLAG_MASK = UInt16(0x0003)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using ..CartesianDomains: RectangularDomain

""" Like a linear index... """
function linear_coordinate(elem, domain, cpu_x, cpu_y, cpu_z)

    corner_x = cpu_x[1, elem.index]
    corner_y = cpu_y[1, elem.index]
    corner_z = cpu_z[1, elem.index]

    Δx = corner_x - domain.x[1]
    Δy = corner_y - domain.y[1]
    Δz = corner_z - domain.z[1]

    coordinate = (
        Δz / domain.L.z * domain.Ne.z * domain.Ne.y * domain.L.x +
        Δy / domain.L.y * domain.Ne.y * domain.L.x +
        Δx
    )

    return coordinate
end

"""
    SpectralElementField(domain::RectangularDomain, grid, realdata::AbstractArray)

Returns a `SpectralElementField` whose `elements` provide a Cartesian `view`
into `realdata`, assuming that `realdata` lives on `domain::RectangularDomain`.
"""
function SpectralElementField(
    domain::RectangularDomain{FT},
    grid::DiscontinuousSpectralElementGrid,
    realdata::AbstractArray,
    data::Union{AbstractArray, Nothing} = nothing,
    realelems::Union{UnitRange, Nothing} = nothing,
) where {FT}

    # Build element list
    Te = prod(domain.Ne) # total number of elements

    element_list = [RectangularElement(domain, grid, realdata, i) for i in 1:Te]

    # Transfer coordinate data to CPU for element sorting
    volume_geometry = grid.vgeo

    x = convert(Array, view(volume_geometry, :, grid.x1id, :))
    y = convert(Array, view(volume_geometry, :, grid.x2id, :))
    z = convert(Array, view(volume_geometry, :, grid.x3id, :))

    # Sort elements by linear coordinate for x, y, z structuring
    sort!(element_list, by = elem -> linear_coordinate(elem, domain, x, y, z))

    # Reshape and permute dims to obtain array where i, j, k correspond to x, y, z
    Ne = domain.Ne
    element_array = reshape(element_list, Ne.x, Ne.y, Ne.z)

    return SpectralElementField(
        element_array,
        domain,
        grid,
        realdata,
        data,
        realelems,
    )
end


text:
# Functions related to the counting bound.
# This version uses analytic approximations from Wikipedia,
# and works with log-transformed numbers.
# (All logs are natural, unless otherwise noted).

module LogCountingBound

export logCountingBound, logApproxNumMaximalCliques,
  approxLogNChooseK, approxBoundAlmostOneExp, logDiff
	# ??? export fewer of these?

"""
Approximate log of binomial(n, k), when n >> k.
From Wikipedia.
"""
function approxLogNChooseK1(n, k)
  k * log(n/k - 0.5) + k - 0.5 * log(2 * pi * k)
end

"""
Better approximation of log(binomial(n, k)) when n >> k.
Also from Wikipedia.
"""
function approxLogNChooseK(n, k)
	((n + 0.5) * log((n + 0.5) / (n + 0.5 - k))
		+ k * log((n + 0.5 - k) / k)
		- 0.5 * log(2 * pi * k))
end

"""
Upper-bound on

log( (1 - 1/a)^b )

given the log of a and b. This is more accurate when
a is large, and b is small.
  logA, logB: log of a and b, respectively
  Returns: log of the above expression
"""
function approxBoundAlmostOneExp(logA, logB)
  - exp( logB - logA )
end

"""
Given log(a) and b, computes log( a +/- b ).
	logA: a log-transformed number
	b: a (non-log-transformed) number, smaller than e^a
	plus: if true, add b; if false, subtract b
	Returns: log(a +/- b), or NaN if something goes negative.
(This is computed using techniques adapted from those used
to deal with log-likelihoods e.g. in NLP).
"""
function logPlusMinus(logA, b, plus)
	if b <= 0
		return NaN
	end
	# compute b, as a fraction of a
	bScaled = exp( log(b) - logA )
	# add or subtract that from a
	diff = if plus; 1 + bScaled; else 1 - bScaled; end
	if diff <= 0
		return NaN
	end
	# return log of that, back on original scale
	log( diff ) + logA
end

"""
Counting bound, based on Shannon's argument.
  m: number of edges
  logW: log of number of 'wires' -- that is,
		log(log2(number of functions)), where
		'log2(number of functions)' is the number of bits
		needed to specify a function
  Returns: log of average number of NAND gates (with unbounded
		fan-in) required to compute any of those functions.
    (This may not be an integer).
"""
function logCountingBound(m, logW)
  m = BigFloat(m)
  logW = BigFloat(logW)
  b = m - 0.5
	# this presumably should be this:
  # the "-1" here is because this is the average, not the max.
  # bound = sqrt(2*w + b*b) - (b - 1)
	# but we use log-transformed numbers in several places
	x = logPlusMinus(log(2) + logW, b*b, true) / 2
	bound = logPlusMinus(x, b - 1, false)
	bound
end

"""
  Log of the approximate number of maximal hypercliques of some size.
  Note that k < r < n .
  Also, the precision of what's returned can be set by setprecision().
  k: number of vertices per hyperedge
  r: number of vertices in the clique
  n: vertices in the larger graph
  Returns: (natural log of) the number of maximal hypercliques. This
    is approximate, because it's the expected number, and it uses an
    approximation for (1 - 1/a)^b. (Again, presumably it's more
    accurate for larger numbers).
  We assume that binomial(r, k-1) and binomial(r, k) are
		computable exactly, but that n may be much larger, and so
		'n choose r' or 'n choose k' require approximation.
"""
function logApproxNumMaximalCliques(k, r, n)
  k = BigInt(k)
  r = BigInt(r)
  n = BigInt(n)

  # log of probability that one of those is
	# not covered by a larger clique
  logA = log(BigInt(2)) * binomial(r, k-1)
	logP = approxBoundAlmostOneExp(logA, log(n-r))

	logNumMaximalCliques = (
    # number of possible cliques of this size
		approxLogNChooseK(n, r)
    # P(clique is present ...)
		+ (- log(BigInt(2)) * binomial(r, k))
    # P( ... and not covered by a larger clique)
    + logP)

  logNumMaximalCliques
end

end

max_stars_repo_name: joshtburdick/misc

max_issues_repo_name: joshtburdick/misc

max_forks_repo_name: joshtburdick/misc


text:
using FFTW, StatsFuns, Statistics, SpecialFunctions

# Some functions for continous wavelet transfrom and cross-spectra
#
#
# This file is based on a straightforward Julia translation of the
# code for 1d Wavelet transform from Torrence and Compo
#
# There is no intention to have a fancy, performant, or elegant "translation"
# it just should be able to "do stuff"
# The copyright remains by Torrence and Compo (see below)
# This code is provided without any warranties
# by Bernd Blasius (bernd.blasius@gmail.com)
#
#  HERE GOES THE ORIGINAL COPYRIGHT
#----------------------------------------------------------------------------
#   Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
#
#   This software may be used, copied, or redistributed as long as it is not
#   sold and this copyright notice is reproduced on each copy made. This
#   routine is provided as is without any express or implied warranties
#   whatsoever.
#
# Notice: Please acknowledge the use of the above software in any publications:
#    ``Wavelet software was provided by C. Torrence and G. Compo,
#      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
#
# Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
#            Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
#
# Please send a copy of such publications to either C. Torrence or G. Compo:
#  Dr. Christopher Torrence               Dr. Gilbert P. Compo
#  Research Systems, Inc.                 Climate Diagnostics Center
#  4990 Pearl East Circle                 325 Broadway R/CDC1
#  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
#  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
#----------------------------------------------------------------------------

#WAVE_SIGNIF  Significance testing for the 1D Wavelet transform WAVELET
#
#   [signif,fft_theor] = ...
#      wave_signif(y,dt,scale,sigtest,lag1,siglvl,dof,mother,param)
#
# INPUTS:
#    y = the time series, or, the VARIANCE of the time series.
#        (If this is a single number, it is assumed to be the variance...)
#    dt = amount of time between each Y value, i.e. the sampling time.
#    scale = the vector of scale indices, from previous call to WAVELET.
#
# OUTPUTS:
#    signif = significance levels as a function of SCALE
#    fft_theor = output theoretical red-noise spectrum as fn of PERIOD
#
# OPTIONAL INPUTS:
# Note * setting any of the following to -1 will cause the default value to be used
#
#    sigtest = 0, 1, or 2.    If omitted, then assume 0.
#
#         If 0 (the default), then just do a regular chi-square test,
#             i.e. Eqn (18) from Torrence & Compo.
#         If 1, then do a "time-average" test, i.e. Eqn (23).
#             In this case, dof should be set to NA, the number
#             of local wavelet spectra that were averaged together.
#             For the Global Wavelet Spectrum, this would be NA=N,
#             where N is the number of points in your time series.
#         If 2, then do a "scale-average" test, i.e. Eqns (25)-(28).
#             In this case, dof should be set to a
#             two-element vector [S1,S2], which gives the scale
#             range that was averaged together.
#             e.g. if one scale-averaged scales between 2 and 8,
#             then dof=[2,8].
#
#    lag1 = LAG 1 Autocorrelation, used for signif levels. Default is 0.0
#
#    siglvl = significance level to use. Default is 0.95
#
#    dof = degrees-of-freedom for signif test.
#         IF sigtest=0, then (automatically) dof = 2 (or 1 for mother="DOG")
#         IF sigtest=1, then dof = na, the number of times averaged together.
#         IF sigtest=2, then dof = [s1,s2], the range of scales averaged.
#
#      Note: IF sigtest=1, then dof can be a vector (same length as scale's),
#            in which case na is assumed to vary with scale.
#            This allows one to average different numbers of times
#            together at different scales, or to take into account
#            things like the Cone of Influence.
#            See discussion following Eqn (23) in Torrence & Compo.
#
#----------------------------------------------------------------------------
#   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
#   University of Colorado, Program in Atmospheric and Oceanic Sciences.
#   This software may be used, copied, or redistributed as long as it is not
#   sold and this copyright notice is reproduced on each copy made.  This
#   routine is provided as is without any express or implied warranties
#   whatsoever.
#----------------------------------------------------------------------------
function wave_signif(y,dt,scale,sigtest=-1,lag1=-1,siglvl=-1,dof=-1,mother=-1,param=-1)
    n1 = length(y)
    j1 = length(scale) - 1
    #scale[1:j1+1] = scale1   # I don't think we need this
    s0 = minimum(scale)
    dj = log(scale[2]/scale[1])/log(2.)

    variance = n1 == 1 ? y : var(y)

    if sigtest == -1; sigtest = 0; end
    if lag1 == -1;    lag1 = 0.0; end
    if siglvl == -1;  siglvl = 0.95; end
    if mother == -1;  mother = "MORLET"; end

    mother = uppercase(mother)

    # get the appropriate parameters [see Table(2)]
    if mother == "MORLET"  #----------------------------------  Morlet
        if param == -1; param = 6.; end
        k0 = param
        fourier_factor = (4*pi)/(k0 + sqrt(2 + k0^2)) # Scale-->Fourier [Sec.3h]
        empir = [2.;-1;-1;-1]
        if k0 == 6; empir[2:4]=[0.776;2.32;0.60]; end
    elseif mother == "PAUL"  #--------------------------------  Paul
        if param == -1; param = 4.; end
        m = param
        fourier_factor = 4*pi/(2*m+1)
        empir = [2.,-1,-1,-1]
        if m == 4; empir[2:4]=[1.132;1.17;1.5]; end
    elseif mother == "DOG"  #---------------------------------  DOG
        if param == -1; param = 2.; end
        m = param
        fourier_factor = 2*pi*sqrt(2 ./ (2*m+1))
        empir = [1.;-1;-1;-1]
        if m == 2; empir[2:4] = [3.541;1.43;1.4]; end
        if m == 6; empir[2:4] = [1.966;1.37;0.97]; end
    else
        error("Mother must be one of MORLET,PAUL,DOG")
    end

    period = scale.*fourier_factor
    dofmin = empir[1]     # Degrees of freedom with no smoothing
    Cdelta = empir[2]     # reconstruction factor
    gamma_fac = empir[3]  # time-decorrelation factor
    dj0 = empir[4]        # scale-decorrelation factor

    freq = dt ./ period   # normalized frequency
    fft_theor = (1 .- lag1^2) ./ (1 .- 2*lag1*cos.(freq*2*pi) .+ lag1^2)  # [Eqn(16)]
    fft_theor = variance*fft_theor  # include time-series variance
    signif = fft_theor
    if dof == -1; dof = dofmin; end

    if sigtest == 0    # no smoothing, DOF=dofmin [Sec.4]
        dof = dofmin
        #chisquare = chisquare_inv(siglvl,dof)/dof # that's Torrence & Compo ..
        chisquare = chisqinvcdf(dof,siglvl)/dof    # .. in julia we use StatsFuns
        signif = fft_theor*chisquare   # [Eqn(18)]
    elseif sigtest == 1  # time-averaged significance
        if length(dof) == 1; dof=zeros(j1+1) .+ dof; end
        #truncate = find(dof .< 1)
        #dof[truncate] = ones(length(truncate))
        map!(z->max(1.0,z), dof, dof)
        dof = dofmin*sqrt.(1 .+ (dof*dt/gamma_fac ./ scale).^2 )   # [Eqn(23)]
        #truncate = find(dof .< dofmin)
        #dof[truncate] = dofmin*ones(length(truncate))   # minimum DOF is dofmin
        map!(z->max(dofmin,z), dof, dof)
        for a1 = 1:j1+1
            #chisquare = chisquare_inv(siglvl,dof[a1])/dof[a1]
            chisquare = chisqinvcdf(dof[a1],siglvl)/dof[a1]  # use StatsFuns
            signif[a1] = fft_theor[a1]*chisquare
        end
    elseif sigtest == 2  # time-averaged significance
        if length(dof) != 2
            error("DOF must be set to [S1,S2], the range of scale-averages")
        end
        if Cdelta == -1
            error("Cdelta & dj0 not defined for " * mother *
                " with param = " * string(param))
        end
        s1 = dof[1]
        s2 = dof[2]
        avg = find(z->((z>=s1) && (z<= s2)), scale)  # scales between S1 & S2
        navg = length(avg)
        if navg == 0
            error("No valid scales between " * string(s1) * " and " * string(s2))
        end
        Savg = 1 ./ sum(1 ./ scale[avg])     # [Eqn(25)]
        Smid = exp((log(s1)+log(s2))/2.)     # power-of-two midpoint
        dof = (dofmin*navg*Savg/Smid)*sqrt(1 + (navg*dj/dj0)^2)  # [Eqn(28)]
        fft_theor = Savg*sum(fft_theor[avg] ./ scale[avg])  # [Eqn(27)]
        #chisquare = chisquare_inv(siglvl,dof) ./ dof
        chisquare = chisqinvcdf.(dof,siglvl) ./ dof     # use StatsFuns
        signif = (dj*dt/Cdelta/Savg) .* fft_theor*chisquare    # [Eqn(26)]
    else
        error("sigtest must be either 0, 1, or 2")
    end

    signif, fft_theor
end

# **********************************************************************************
# **********************************************************************************

heavi(x) = x > 0. ? 1. : 0.    # Heaviside function

function wave_bases(mother,k,scale,param)
#   Computes the wavelet function as a function of Fourier frequency,
#   used for the wavelet transform in Fourier space.
#   (This program is called automatically by WAVELET)
#
# INPUTS:
#    mother = a string, equal to 'MORLET' or 'PAUL' or 'DOG'
#    k = a vector, the Fourier frequencies at which to calculate the wavelet
#    scale = a number, the wavelet scale
#    param = the nondimensional parameter for the wavelet function
#
# OUTPUTS:
#    daughter = a vector, the wavelet function
#    fourier_factor = the ratio of Fourier period to scale
#    coi = a number, the cone-of-influence size at the scale
#    dofmin = a number, degrees of freedom for each point in the wavelet power
#             (either 2 for Morlet and Paul, or 1 for the DOG)
#             (either 2 for Morlet and Paul, or 1 for the DOG)
#
#    Copyright statement of the original Matlab version
#----------------------------------------------------------------------------
#   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
#   University of Colorado, Program in Atmospheric and Oceanic Sciences.
#   This software may be used, copied, or redistributed as long as it is not
#   sold and this copyright notice is reproduced on each copy made.  This
#   routine is provided as is without any express or implied warranties
#   whatsoever.
#----------------------------------------------------------------------------
    mother = uppercase(mother)
    n = length(k)
    if mother == "MORLET"  #-----------------------------------  Morlet
       if param == -1; param = 6.; end
       k0 = param
       expnt = -(scale*k .- k0).^2 / 2 .* heavi.(k)
       norm = sqrt(scale*k[2])*(pi^(-0.25))*sqrt(n)      # total energy=N   [Eqn(7)]
       daughter = norm*exp.(expnt)
       daughter = daughter .* heavi.(k)                  # Heaviside step function
       fourier_factor = (4.0*pi)/(k0 + sqrt(2.0 + k0^2)) # Scale-->Fourier [Sec.3h]
       coi = fourier_factor/sqrt(2.0)                    # Cone-of-influence [Sec.3g]
       dofmin = 2                                        # Degrees of freedom
     elseif mother == "PAUL"  #-----------------------------------  Paul
       if param == -1; param = 4. ; end
       m = param
       expnt = -(scale.*k).* heavi.(k)
       norm = sqrt(scale*k[2])*(2^m/sqrt(m*prod(2:(2*m-1))))*sqrt(n)
       daughter = norm*((scale.*k).^m).*exp.(expnt)
       daughter = daughter.* heavi.(k)        # Heaviside step function
       fourier_factor = 4*pi/(2.0*m+1)
       coi = fourier_factor*sqrt(2.0)
       dofmin = 2
     elseif mother == "DOG"  #-----------------------------------  DOG
       if param == -1; param = 2.; end
       m = param
       expnt = -(scale.*k).^2 ./ 2.0
       norm = sqrt(scale*k[2]/gamma(m+0.5))*sqrt(n)
       daughter = -norm*(im^m)*((scale.*k).^m).*exp.(expnt)
       fourier_factor = 2*pi*sqrt(2 ./ (2*m+1))
       coi = fourier_factor/sqrt(2)
       dofmin = 1
     else
       error("Mother must be one of MORLET,PAUL,DOG")
     end
     daughter, fourier_factor, coi, dofmin
end

# **********************************************************************************
# **********************************************************************************

#   1D contunous Wavelet transform
#
#   Computes the wavelet transform of the vector Y (length N),
#   with sampling rate DT.
#
#   By default, the Morlet wavelet (k0=6) is used.
#   The wavelet basis is normalized to have total energy=1 at all scales.
#
#
# INPUTS:
#
#    y = the time series of length N.
#    dt = amount of time between each Y value, i.e. the sampling time.
#
# OUTPUTS:
#
#    wave is the WAVELET transform of y. This is a complex array
#    of dimensions (n,j1+1). real(wave) gives the WAVELET amplitude,
#    atan(imag(wave),real(wave))  OR angle(wave) gives the WAVELET phase.
#    The WAVELET power spectrum is abs2(wave).
#    Its units are sigma^2 (the time series variance).
#
#
# OPTIONAL INPUTS:
#
# *** Note *** setting any of the following to -1 will cause the default
#               value to be used.
#
#    pad = if set to 1 (default is 0), pad time series with enough zeroes to get
#         N up to the next higher power of 2. This prevents wraparound
#         from the end of the time series to the beginning, and also
#         speeds up the FFT's used to do the wavelet transform.
#         This will not eliminate all edge effects (see COI below).
#
#    dj = the spacing between discrete scales. Default is 0.25.
#         A smaller # will give better scale resolution, but be slower to plot.
#
#    s0 = the smallest scale of the wavelet.  Default is 2*DT.
#
#    j1 = the # of scales minus one. Scales range from s0 up to s0*2^(j1*dj),
#        to give a total of (j1+1) scales. Default is j1 = (log2(n dt/s0))/dj.
#
#    mother = the mother wavelet function.
#             The choices are 'MORLET', 'PAUL', or 'DOG'
#
#    param = the mother wavelet parameter.
#            For 'MORLET' this is k0 (wavenumber), default is 6.
#            For 'PAUL' this is m (order), default is 4.
#            For 'DOG' this is m (m-th derivative), default is 2.
#
#
# OUTPUTS:
#
#    period = the vector of "Fourier" periods (in time units) that corresponds
#           to the scaels.
#
#    scale = the vector of scale indices, given by s0*2^(j*dj), j=0...j1
#            where j1+1 is the total # of scales.
#
#    coi = if specified, then return the Cone-of-Influence, which is a vector
#        of N points that contains the maximum period of useful information
#        at that particular time.
#        Periods greater than this are subject to edge effects.
#        This can be used to plot coi lines on a contour plot by doing:
#
#             IN MATLAB:
#              contour(time,log(period),log(power))
#              plot(time,log(coi),'k')
#
#----------------------------------------------------------------------------
#   Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
#
#   This software may be used, copied, or redistributed as long as it is not
#   sold and this copyright notice is reproduced on each copy made. This
#   routine is provided as is without any express or implied warranties
#   whatsoever.
#
# Notice: Please acknowledge the use of the above software in any publications:
#    ``Wavelet software was provided by C. Torrence and G. Compo,
#      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
#
# Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
#            Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78.
#
# Please send a copy of such publications to either C. Torrence or G. Compo:
#  Dr. Christopher Torrence               Dr. Gilbert P. Compo
#  Research Systems, Inc.                 Climate Diagnostics Center
#  4990 Pearl East Circle                 325 Broadway R/CDC1
#  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
#  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
#----------------------------------------------------------------------------

function wavelet(y,dt,pad=0,dj=-1,s0=-1,j1=-1,mother=-1,param=-1)

    n1 = length(y)
    if s0 == -1; s0=2*dt; end
    if dj == -1; dj = 1. / 4.; end
    if j1 == -1; j1=trunc(Int,(log(n1*dt/s0)/log(2))/dj); end
    if mother == -1; mother = "MORLET"; end

    #....construct time series to analyze, pad if necessary
    x = y .- mean(y)
    if pad == 1
        base2 = trunc(Int,log(n1)/log(2) + 0.4999)   # power of 2 nearest to N
        x = [x; zeros(2^(base2+1)-n1)]
    end
    n = length(x)

    #....construct wavenumber array used in transform [Eqn(5)]
    k = 1:trunc(Int,n/2)
    k = k .* ((2. * pi)/(n*dt))
    k = [0.; k; -k[trunc(Int,(n-1)/2):-1:1]]
    #
    #....compute FFT of the (padded) time series
    f = fft(x)    # [Eqn(3)]

    #....construct SCALE array & empty PERIOD & WAVE arrays
    scale = s0*2 .^ ((0:j1)*dj)
    period = scale
    wave = zeros(Complex,j1+1,n)   # define the wavelet array

    # loop through all scales and compute transform
    fourier_factor = 0.0  # julia seems to need a default
    for a1 = 1:j1+1
       daughter,fourier_factor,coi,dofmin = wave_bases(mother,k,scale[a1],param)
       wave[a1,:] = ifft(f.*daughter)  # wavelet transform[Eqn(4)]
    end

    period = fourier_factor*scale
    #coi = coi*dt*[1e-5; collect(1:((n1+1)/2-1)); collect((n1/2-1):-1:1); 1e-5]  # COI [Sec.3g]
    coi = coi*dt*[1e-5; (1:((n1+1)/2-1)); ((n1/2-1):-1:1); 1e-5]  # COI [Sec.3g]
    wave = wave[:,1:n1]  # get rid of padding before returning

    wave, period, scale, coi
end

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
# Usage:
#     julia benchmark/run_benchmarks.jl
using BenchmarkTools
using LocalBinaryPatterns
using ImageTransformations
using TestImages
using ImageCore
using Interpolations

include("mapwindow_impl.jl")

on_CI = haskey(ENV, "GITHUB_ACTIONS")

img = testimage("cameraman")

tst_sizes = [256, ]
tst_types = [Gray{Float32}]

const SUITE = BenchmarkGroup()

alg_list = [
    ("Original", local_binary_pattern),
    ("Original 8 bits with bilinear", img->local_binary_pattern(img, 8, 1.0, Linear())),
    ("Original 8 bits with nearest neighbours", img->local_binary_pattern(img, 8, 1.0, Constant())),
    ("Original 24 bits with bilinear", img->local_binary_pattern(img, 24, 3.0, Linear())),
    ("Modified", img->local_binary_pattern(LBP(average_mode, (1, 1)), img)),
    ("Multi-block (1x1)", img->local_binary_pattern(LBP((1, 1)), img)),
    ("Multi-block (3x3)", img->local_binary_pattern(LBP((3, 3)), img)),
    ("Modified Multi-block (3x3)", img->local_binary_pattern(LBP(average_mode, (3, 3)), img))
]

function add_algorithm_benchmark!(suite, img, alg_name, alg;
                                  tst_sizes, tst_types)
    haskey(suite, alg_name) || (suite[alg_name] = BenchmarkGroup())

    for T in tst_types
        haskey(suite[alg_name], T) || (suite[alg_name][T] = BenchmarkGroup())
        for sz in tst_sizes
            tst_img = imresize(T.(img), (sz, sz))

            suite[alg_name][T]["$sz×$sz"] = @benchmarkable $alg($tst_img)
        end
    end
end

for (alg_name, alg) in alg_list
    add_algorithm_benchmark!(SUITE, img, alg_name, alg;
                             tst_sizes=tst_sizes,
                             tst_types=tst_types)
end

# compare with the simple mapwindow version
suite = SUITE["Original (mapwindow)"] = BenchmarkGroup()
for T in tst_types
    haskey(suite, T) || (suite[T] = BenchmarkGroup())
    for sz in tst_sizes
        tst_img = imresize(T.(img), (sz, sz))
        suite[T]["$sz×$sz"] = @benchmarkable lbp_origin_mapwindow($tst_img)
    end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    SynchronousExecution

Functor object for using synchronous execution in a progressive-hedging algorithm (assuming multiple Julia cores are available). Create by supplying a [`Synchronous`](@ref) object through `execution` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SynchronousExecution{T <: AbstractFloat,
                            A <: AbstractVector,
                            PT <: AbstractPenaltyTerm} <: AbstractProgressiveHedgingExecution
    subworkers::Vector{SubWorker{T,A,PT}}

    function SynchronousExecution(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat,
                                                                           A <: AbstractVector,
                                                                           PT <: AbstractPenaltyTerm}
        return new{T,A,PT}(Vector{SubWorker{T,A,PT}}(undef, nworkers()))
    end
end

function initialize_subproblems!(ph::AbstractProgressiveHedging,
                                 execution::SynchronousExecution,
                                 scenarioproblems::DistributedScenarioProblems,
                                 penaltyterm::AbstractPenaltyTerm)
    # Create subproblems on worker processes
    initialize_subproblems!(ph,
                            execution.subworkers,
                            scenarioproblems,
                            penaltyterm)
    # Initial reductions
    update_iterate!(ph)
    update_dual_gap!(ph)
    return nothing
end

function restore_subproblems!(::AbstractProgressiveHedging, execution::SynchronousExecution)
    restore_subproblems!(execution.subworkers)
    return nothing
end

function resolve_subproblems!(ph::AbstractProgressiveHedging, execution::SynchronousExecution{T}) where T <: AbstractFloat
    partial_solutions = Vector{SubproblemSolution{T}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_solutions[i] = remotecall_fetch(resolve_subproblems!,
                                                           w,
                                                           execution.subworkers[w-1],
                                                           ph.ξ,
                                                           penalty(ph))
        end
    end
    return sum(partial_solutions)
end

function update_iterate!(ph::AbstractProgressiveHedging, execution::SynchronousExecution{T,A}) where {T <: AbstractFloat, A <: AbstractVector}
    partial_primals = Vector{A}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_primals[i] = remotecall_fetch(collect_primals,
                                                         w,
                                                         execution.subworkers[w-1],
                                                         length(ph.ξ))
        end
    end
    ξ_prev = copy(ph.ξ)
    ph.ξ .= sum(partial_primals)
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ - ξ_prev, 2) ^ 2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedging, execution::SynchronousExecution)
    # Update dual prices
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                execution.subworkers[w-1],
                ph.ξ,
                penalty(ph)) do sw, ξ, r
                    subproblems = fetch(sw)
                    if length(subproblems) > 0
                        update_subproblems!(subproblems, ξ, r)
                    end
                end
        end
    end
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedging, execution::SynchronousExecution)
    return update_dual_gap!(ph, execution.subworkers)
end

function calculate_objective_value(ph::AbstractProgressiveHedging, execution::SynchronousExecution)
    return calculate_objective_value(execution.subworkers)
end

function scalar_subproblem_reduction(value::Function, execution::SynchronousExecution{T}) where T <: AbstractFloat
    partial_results = Vector{T}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_results[i] = remotecall_fetch(w, execution.subworkers[w-1], value) do sw, value
                subproblems = fetch(sw)
                return mapreduce(+, subproblems, init = zero(T)) do subproblem
                    π = subproblem.probability
                    return π * value(subproblem)
                end
            end
        end
    end
    return sum(partial_results)
end

function vector_subproblem_reduction(value::Function, execution::SynchronousExecution{T,A}, n::Integer) where {T <: AbstractFloat, A <: AbstractVector}
    partial_results = Vector{A}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_results[i] = remotecall_fetch(w, execution.subworkers[w-1], value, n) do sw, value, n
                subproblems = fetch(sw)
                return mapreduce(+, subproblems, init = zero(T, n)) do subproblem
                    π = subproblem.probability
                    return π * value(subproblem)
                end
            end
        end
    end
    return sum(partial_results)
end

# API
# ------------------------------------------------------------
function (execution::Synchronous)(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat,
                                                                           A <: AbstractVector,
                                                                           PT <: AbstractPenaltyTerm}
    return SynchronousExecution(T,A,PT)
end

function str(::Synchronous)
    return "Synchronous "
end


text:
using PumasTutorials, Pkg
Pkg.activate("../tutorials")
PumasTutorials.weave_file(".", "test.jmd")


text:
using CrossDecomposition

# This implements the RDC as described in the paper:
#  http://arxiv.org/pdf/1304.7717v2.pdf
function rdc(x::AbstractMatrix, y::AbstractMatrix, k = 10, s = 1/6, f1 = cos, f2 = sin)
  tx = hcat(mapslices((u) -> sortperm(u)/length(u), x, 1), ones(size(x,1)))
  ty = hcat(mapslices((u) -> sortperm(u)/length(u), y, 1), ones(size(y,1)))
  wx = s/size(tx,2) * randn(size(tx,2), k)
  wy = s/size(ty,2) * randn(size(ty,2), k)
  cor(canoncor(hcat(f1(tx*wx), f2(tx*wx)), hcat(f1(tx*wx), f2(tx*wx))))[1]
end

x = randn(100,1)
y = 10x
rdc(x,y,10,1/6)

# R code:
#rdc <- function(x,y,k,s) {
#  tx <- cbind( apply(as.matrix(x),2,function(u) rank(u)/length(u)), 1 )
#  ty <- cbind( apply(as.matrix(y),2,function(u) rank(u)/length(u)), 1 )
#  wx <- matrix(rnorm(ncol(tx)*k,0,s),ncol(tx),k)
#  wy <- matrix(rnorm(ncol(ty)*k,0,s),ncol(ty),k)
#  cancor(cbind(cos(tx%*%wx),sin(tx%*%wx)), cbind(cos(ty%*%wy),sin(ty%*%wy)))$cor[1]
#}


text:
#rank-1 update of general matrices
#ger
#gerc
#geru

#rank-1 update of hermitian matrices
#her
#hpr

#rank-2 update of hermitian matrices
#her2
#hpr2

#updates of symmetric packed matrices
#spr
#spr2

#updates of symmetric matrices
#syr
#syr2

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
escapename(n::String) =
    replace(n,r"(^#)" => "_#")

fmtcof(c::Float64) =
    if     c < 0.0
        if     c < -1.0 || c > -1.0
             "- $(-c)"
        else
            "-"
        end
    elseif c > 0.0
        if c < 1.0 || c > 1.0
            "+ $c"
        else
            "+"
        end
    else
        "0"
    end

#Base.showall(f::IO, t::Task) = showlimited(f,t,0)
Base.show(f::IO, t::Task) = showlimited(f,t,20)

Base.show(io::IOContext, t::Task) =
    if get(io, :short, false)
        showlimited(io,t,20)
    else
        showlimited(io,t,0)
    end

function showlimited(f::IO, t::Mosek.Task, limit :: Int)
    numvar    = getnumvar(t)
    numbarvar = getnumbarvar(t)
    numcon    = getnumcon(t)
    numcone   = getnumcone(t)
    numanz    = getnumanz(t)
    numqobjnz = getnumqobjnz(t)
    numqconnz = if numcon > 0 sum(i -> getnumqconknz(t,i),1:numcon) else 0 end

    limitnumcon = if limit > 0 min(limit,numcon) else numcon end
    limitnumvar = if limit > 0 min(limit,numvar) else numvar end
    limitnumcone = if limit > 0 min(limit,numcone) else numcone end
    limitnumbarvar = if limit > 0 min(limit,numbarvar) else numbarvar end

    termslimit = numvar+numbarvar+numqobjnz+numqconnz
    if limit > 0 termslimit = min(limit,termslimit) end

    if numvar == 0 && numcon == 0 && numbarvar == 0
        println(f,"Empty Task")
    else
        varnames = String[ escapename(getvarname(t,i)) for i in 1:numvar]
        barvarnames = String[ escapename(getbarvarname(t,i)) for i in 1:numbarvar]
        connames = String[ escapename(getconname(t,i)) for i in 1:numcon]
        conenames = String[ escapename(getconename(t,i)) for i in 1:numcone]

        for i in 1:numvar
            if length(varnames[i]) == 0
                varnames[i] = "#x$i"
            end
        end

        for i in 1:numcon
            if length(connames[i]) == 0
                connames[i] = "#c$i"
            end
        end

        for i in 1:numcone
            if length(conenames[i]) == 0
                conenames[i] = "#k$i"
            end
        end

        for i in 1:numbarvar
            if length(barvarnames[i]) == 0
                barvarnames[i] = "#X̄$i"
            end
        end

        maxvarnamelen = if numvar > 0 max(map(n -> length(n),varnames)...) else 0 end
        maxconnamelen = if numcon > 0 max(map(n -> length(n),connames)...) else 0 end

        taskname = gettaskname(t)

        println(f,"Task '$(escape_string(taskname))'")

        _,barcidx = getbarcsparsity(t)

        _,baraidx = getbarasparsity(t)
        baraptrb = zeros(Int,numcon+1)
        for idx in baraidx
            (barai,baraj) = getbaraidxij(t,idx)
            baraptrb[barai+1] += 1
        end
        for i in 1:numcon
            baraptrb[i+1] += baraptrb[i]
        end
        baraptrb .+= 1

        baraptre = copy(baraptrb)
        baraidx_idx = Vector{Int}(undef,length(baraidx))
        for idx in baraidx
            (barai,baraj) = getbaraidxij(t,idx)
            baraidx_idx[baraptre[barai]] = idx
            baraptre[barai] += 1
        end

        # Objective
        objsense = getobjsense(t)
        objname  = escapename(getobjname(t))
        if objname == ""
            objname = "#obj"
        end

        if objsense == MSK_OBJECTIVE_SENSE_MAXIMIZE
            println(f,"    Maximize")
        else
            println(f,"    Minimize")
        end
        c = getc(t)
        print(f,"        $objname: ")
        terms = 0
        omittedterms = 0
        for j in 1:numvar
            if c[j] < 0.0 || c[j] > 0.0
               if terms < termslimit
                   print(f,"$(fmtcof(c[j])) $(varnames[j]) ")
                   terms += 1
               else
                   omittedterms += 1
               end
            end
        end

        for idx in barcidx
            if terms < termslimit
                (barcj,num,sub,w) = getbarcidx(t,idx)
                if num == 1
                    print(f,"$(fmtcof(w[1])) #MX$(sub[1]) ⋅ $(barvarnames[barcj]) ")
                elseif num > 1
                    print(f,"(")
                    printf(f,"$(fmtcof(w[1])) #MX$(sub[1])")
                    print(f,") ⋅ $(barvarnames[barcj]) ")
                end
                terms += 1
            else
                omittedterms += 1
            end
        end

        numqobjnz = getnumqobjnz(t)
        if numqobjnz > 0
            (nnz,qcsubi,qcsubj,qcval) = getqobj(t)
            for k in 1:nnz
                if terms < limitterms
                    if qcsubi[k] == qcsubj[k]
                        print(f,"$(Mosek.fmtcof(qcval[k])) $(varnames[qcsubi[k]]) ^ 2")
                    else
                        ni = varnames[qcsubi[k]]
                        nj = varnames[qcsubj[k]]
                        print(f,"$(Mosek.fmtcof(qcval[k])) $ni ⋅ $nj")
                    end
                    terms += 1
                else
                    omittedterms += 1
                end
            end
        end

        if omittedterms > 0
            print(f," ... ($omittedterms terms omitted)")
        end

        cfix = getcfix(t)
        if cfix < 0 || cfix > 0
            print(f,fmtcof(cfix))
        end
        println(f)

        # Constraints
        if numcon+numcone > 0
            println(f,"    Subject to")

            for i in 1:limitnumcon
                (nzi,subi,vali) = getarow(t,Int32(i))
                (bk,bl,bu) = getconbound(t,Int32(i))

                print(f,"        $(connames[i]): ")
                if bk == MSK_BK_LO || bk == MSK_BK_RA
                    print(f,"$bl < ")
                end
                terms = 0
                omittedterms = 0
                for k in 1:nzi
                    if terms < termslimit
                        print(f,"$(fmtcof(vali[k])) $(varnames[subi[k]]) ")
                        terms += 1
                    else
                        omittedterms += 1
                    end
                end

                for k in baraptrb[i]:baraptrb[i+1]-1
                    if terms < termslimit
                        idx = baraidx_idx[k]
                        (barai,baraj,num,sub,w) = getbaraidx(t,idx)
                        if num == 1
                            print(f,"$(fmtcof(w[1])) #MX$(sub[1]) ⋅ $(barvarnames[baraj]) ")
                        elseif num > 1
                            print(f,"(")
                            printf(f,"$(fmtcof(w[1])) #MX$(sub[1])")
                            print(f,") ⋅ $(barvarnames[baraj]) ")
                        end
                        terms += 1
                    else
                        omittedterms += 1
                    end
                end

                numqconnz = getnumqconknz(t,i)
                if numqconnz > 0
                    (nnz,qcsubi,qcsubj,qcval) = getqconk(t,i)
                    for k in 1:nnz
                        if terms < limitterms
                            if qcsubi[k] == qcsubj[k]
                                print(f,"$(Mosek.fmtcof(qcval[k])) $(varnames[qcsubi[k]]) ^ 2")
                            else
                                ni = varnames[qcsubi[k]]
                                nj = varnames[qcsubj[k]]
                                print(f,"$(Mosek.fmtcof(qcval[k])) $ni ⋅ $nj")
                            end
                            terms += 1
                        else
                            omittedterms += 1
                        end
                    end
                end

                if omittedterms > 0
                    print(f," ... ($omittedterms terms omitted)")
                end

                if bk == MSK_BK_UP || bk == MSK_BK_RA
                    println(f,"< $bu")
                elseif bk == MSK_BK_FX
                    println(f,"= $bu")
                else
                    println(f,"")
                end
            end

            numacc = getnumacc(t)
            for i in 1:numacc
                domidx = getaccdomain(t,i)
                sz = getdomainn(t,i)
                name = getaccname(t,i)
                afeidxs = getaccafeidxlist(t,i)
                domtp = getdomaintype(t,domidx)

                domname = ( if     domtp == MSK_DOMAIN_R "Unbounded"
                            elseif domtp == MSK_DOMAIN_RZERO "Zero"
                            elseif domtp == MSK_DOMAIN_RPLUS "Nonnegative"
                            elseif domtp == MSK_DOMAIN_RMINUS "Nonpositive"
                            elseif domtp == MSK_DOMAIN_INF_NORM_CONE "K_inf"
                            elseif domtp == MSK_DOMAIN_ONE_NORM_CONE "K_1"
                            elseif domtp == MSK_DOMAIN_QUADRATIC_CONE  "𝒞_q"
                            elseif domtp == MSK_DOMAIN_RQUADRATIC_CONE  "𝒞_qr"
                            elseif domtp == MSK_DOMAIN_PRIMAL_GEO_MEAN_CONE "GeoMean"
                            elseif domtp == MSK_DOMAIN_PRIMAL_POWER_CONE "𝒞_pow"
                            elseif domtp == MSK_DOMAIN_DUAL_POWER_CONE "𝒞*_pow"
                            elseif domtp == MSK_DOMAIN_PRIMAL_EXP_CONE "𝒞_exp"
                            elseif domtp == MSK_DOMAIN_DUAL_EXP_CONE "𝒞*_exp"
                            elseif domtp == MSK_DOMAIN_PSD_CONE "PSD"
                            else   "<?>"
                            end )
                println(f,"        \"$name\" $domname")
                for j in 1:sz
                    g = getafeg(t,afeidxs[j])
                    (nzi,subj,cof) = getafefrow(t,afeidxs[j])
                    (barj,termptr,midx,w) = getafebarfrow(t,afeidxs[j])
                    print(f,"            ")

                    terms = 0
                    for k in 1:nzi
                        if terms < termslimit
                            print(f,"$(fmtcof(cof[k])) $(varnames[subj[k]]) ")
                            terms += 1
                        else
                            omittedterms += 1
                        end
                    end

                    # for idx in barj
                    #     if terms < termslimit
                    #         (barcj,num,sub,w) = getbarcidx(t,idx)
                    #         if num == 1
                    #             print(f,"$(fmtcof(w[1])) #MX$(sub[1]) ⋅ $(barvarnames[barcj]) ")
                    #         elseif num > 1
                    #             print(f,"(")
                    #             printf(f,"$(fmtcof(w[1])) #MX$(sub[1])")
                    #             print(f,") ⋅ $(barvarnames[barcj]) ")
                    #         end
                    #         terms += 1
                    #     else
                    #         omittedterms += 1
                    #     end
                    # end
                    print(f,"\n")
                end
            end

            if limitnumcon < numcon
                println(f,"        ... ($(numcon-limitnumcon) constraints omitted)")
            end

            if numcone > 0
                for k in 1:limitnumcone
                    (ct,conepar,nummem,submem) = getcone(t,k)
                    dom =
                        if     ct == MSK_CT_ZERO "𝒞_0"
                        elseif ct == MSK_CT_QUAD "𝒞_q"
                        elseif ct == MSK_CT_RQUAD "𝒞_qr"
                        elseif ct == MSK_CT_PPOW "𝒞_pow{$conepar}"
                        elseif ct == MSK_CT_PEXP "𝒞_exp"
                        else "?{$ct}"
                        end

                    print(f,"        $(conenames[k]): (")

                    if nummem > 0
                        print(f,"$(varnames[submem[1]])")
                        for i in 2:nummem
                            print(f,",$(varnames[submem[i]])")
                        end
                    end
                    println(f,") ∈ $dom($nummem)")
                end

                if limitnumcone < numcone
                    println(f,"        ... ($(numcone-limitnumcone) cones omitted)")
                end
            end
        end

        # Variable bounds
        if numvar+numbarvar > 0
            println(f,"    Variables")
            if numvar > 0
                for j in 1:limitnumvar
                    bk,bl,bu = getvarbound(t,j)

                    if bk == MSK_BK_FX
                        println(f,"        $(varnames[j]) = $bl")
                    else
                        blstr =
                            if     bk == MSK_BK_LO || bk == MSK_BK_RA "[$bl"
                            else   "]-Inf"
                            end
                        bustr =
                            if     bk == MSK_BK_UP || bk == MSK_BK_RA "$bu]"
                            else   "+Inf["
                            end

                        println(f,"        $(varnames[j]) ∈ $blstr;$bustr")
                    end
                end
            end
            if limitnumvar < numvar
                println(f,"        ... ($(numvar-limitnumvar) variable bounds omitted)")
            end

            for j in 1:limitnumbarvar
                dim = getdimbarvarj(t,j)
                println(f,"        $(barvarnames[j]) ∈ S($dim)")
            end
            if limitnumbarvar < numbarvar
                println(f,"        ... ($(numbarvar-limitnumbarvar) semidefinite variables omitted)")
            end
        end

        # Matsto
        numsymmat = getnumsymmat(t)
        limitnumsymmat = if limit > 0 min(limit,numsymmat) else numsymmat end
        if numsymmat > 0
            println(f,"    Symmetric matrixes")
            for i in 1:limitnumsymmat
                dim,nz,tp = getsymmatinfo(t,i)
                (msubi,msubj,mvalij) = getsparsesymmat(t,i)
                print(f,"        #MX$i ($dim): [ " )
                if nz > 0
                    limitnz = if limit > 0 min(limit,nz) else nz end
                    print(f,"($(msubi[1]),$(msubj[1]),$(mvalij[1]))")
                    for i in 2:limitnz
                        print(f,",($(msubi[i]),$(msubj[i]),$(mvalij[i]))")
                    end
                    if limitnz < nz
                        print(f,"... ($(nz-limitnz) nonzeros omitted)")
                    end
                end
                println(f," ]")
            end
            if limitnumsymmat < numsymmat
                println(f,"        ... ($(numsymmat-limitnumsymmat) matrixes omitted)")
            end

        end
    end
end


text:
type LabeledArray{T,N,A<:AbstractArray,L<:AbstractArray} <: AbstractArray{T,N}
    data::A
    label::L
    colors::Vector{RGB}

    # Enforce that the label has Integer type
    LabeledArray{Ti<:Integer}(data::AbstractArray{T,N}, label::AbstractArray{Ti}, colors::Vector{RGB}) =
        new(data, label, colors)
end
LabeledArray{T,N}(data::AbstractArray{T,N}, label::AbstractArray, colors::Vector{RGB}) =
    LabeledArray{T,N,typeof(data),typeof(label)}(data, label, colors)

size(A::LabeledArray) = size(A.data)
size(A::LabeledArray, i::Integer) = size(A.data, i)
eltype{T}(A::LabeledArray{T}) = T
ndims{T,N}(A::LabeledArray{T,N}) = N

for N = 1:4
    @eval begin
        # All but the last of these are inside the @eval loop simply to avoid ambiguity warnings.
        # These first two are additional ones needed to avoid ambiguity warnings.
        _uint32color_gray!{T}(buf::Array{Uint32}, A::LabeledArray{T,$N}, scalei::ScaleSigned) = error("Cannot use ScaleSigned with a labeled array")
        _uint32color_gray!{T,L<:LabeledArray}(buf::Array{Uint32}, A::SubArray{T,$N,L}, scalei::ScaleSigned) = error("Cannot use ScaleSigned with a labeled array")

        function _uint32color_gray!{T}(buf::Array{Uint32}, A::LabeledArray{T,$N}, scalei::ScaleInfo = scaleinfo(Uint8, A))
            if size(buf) != size(A)
                error("Size mismatch")
            end
            dat = A.data
            label = A.label
            col = A.colors
            for i = 1:length(dat)
                gr = scale(scalei, dat[i])
                lbl = label[i]
                if lbl == 0
                    buf[i] = rgb24(gr,gr,gr)
                else
                    buf[i] = convert(RGB24, gr*col[lbl])
                end
            end
            buf
        end

        # For SubArrays, we can't efficiently use linear indexing, and in any event
        # we want to broadcast label where necessary
        function _uint32color_gray!{T,A<:LabeledArray}(buf::Array{Uint32}, S::SubArray{T,$N,A}, scalei::ScaleInfo = scaleinfo(Uint8, A))
            if size(buf) != size(S)
                error("Size mismatch")
            end
            indexes = S.indexes
            dat = slice(S.parent.data, indexes)
            plabel = S.parent.label
            newindexes = RangeIndex[size(plabel,i)==1 ? (isa(indexes[i], Int) ? 1 : (1:1)) : indexes[i]  for i = 1:ndims(plabel)]
            label = slice(plabel, newindexes...)
            col = S.parent.colors
            _uint32color_labeled(buf, dat, label, col, scalei) # type of label can't be inferred, use function boundary
        end

        function _uint32color_labeled{T}(buf, dat::AbstractArray{T,$N}, label, col, scalei)
            k = 0
            @inbounds @nloops $N i buf begin
                val = @nref $N dat i
                gr = scale(scalei, val)
                lbl = @nref $N label i
                if lbl == 0
                    buf[k+=1] = rgb24(gr,gr,gr)
                else
                    buf[k+=1] = convert(RGB24, gr*col[lbl])
                end
            end
            buf
        end
    end
end


text:
using Test
using SymDiff

@test differentiate(:(x^2), :x) == :(2 * x ^ (2 - 1) * 1)


text:
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using POMDPToolbox
using BasicPOMCP
using POMCPOW
using ParticleFilters
using CPUTime
using VDPTag
using DataFrames
using ProfileView

pomdp = VDPTagPOMDP()

solvers = Dict{String, Union{Solver,Policy}}(

    "pomcpow" => begin
        ro = ToNextML(mdp(pomdp))
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(40.0),
                               final_criterion=MaxTries(),
                               max_depth=20,
                               max_time=0.1,
                               k_action=8.0,
                               alpha_action=1/20,
                               k_observation=4.0,
                               alpha_observation=1/20,
                               estimate_value=0.0,
                               check_repeat_act=true,
                               check_repeat_obs=true,
                               rng=MersenneTwister(13)
                              )
    end,

    "pomcp" => begin
        ro = ToNextML(mdp(pomdp))
        POMCPSolver(max_depth=20,
                    max_time=0.1,
                    c=40.0,
                    tree_queries=typemax(Int),
                    # default_action=TagAction(false, 0.0),
                    estimate_value=0.0,
                    rng=MersenneTwister(17)
                   )
    end
)

n_angles = 10
dpomdp = AODiscreteVDPTagPOMDP(n_angles=n_angles, n_obs_angles=n_angles)

solver = solvers["pomcp"]

planner = solve(solver, dpomdp)
cplanner = translate_policy(planner, dpomdp, pomdp, dpomdp)

filter = ObsAdaptiveParticleFilter(deepcopy(pomdp),
                                   LowVarianceResampler(10_000),
                                   0.05, MersenneTwister(90_000))

@show simulate(RolloutSimulator(max_steps=100), pomdp, cplanner, filter)

Profile.clear()
@profile simulate(RolloutSimulator(max_steps=100), pomdp, cplanner, filter)
ProfileView.view()