text:
#===================================================================================================
  Pairwise Kernels - Consume two vectors
===================================================================================================#

abstract PairwiseKernel{T<:AbstractFloat} <: StandardKernel{T}
abstract AdditiveKernel{T<:AbstractFloat} <: PairwiseKernel{T}

function description_string(κ::PairwiseKernel, showtype::Bool = true)
    obj = typeof(κ)
    fields = fieldnames(obj)
    obj_str = string(obj.name.name) * (showtype ? string("{", eltype(κ), "}") : "")
    if length(fields) == 0
        return obj_str
    else
        fields_str = join(["$field=$(getfield(κ,field).value)" for field in fields], ",")
        *(obj_str, "(", fields_str, ")")
    end
end

function =={T<:PairwiseKernel}(κ1::T, κ2::T)
    all([getfield(κ1,field) == getfield(κ2,field) for field in fieldnames(T)])
end

#== Non-Negative Negative Definite Additive Kernel ==#

abstract NonNegNegDefAdditiveKernel{T<:AbstractFloat} <: AdditiveKernel{T}
@inline isnegdef(::NonNegNegDefAdditiveKernel)        = true
@inline attainsnegative(::NonNegNegDefAdditiveKernel) = false

doc"SquaredDistanceKernel() = (x-y)ᵀ(x-y)"
immutable SquaredDistanceKernel{T<:AbstractFloat} <: NonNegNegDefAdditiveKernel{T} end
SquaredDistanceKernel() = SquaredDistanceKernel{Float64}()
convert{T}(::Type{SquaredDistanceKernel{T}}, ::SquaredDistanceKernel) = SquaredDistanceKernel{T}()
@inline phi{T<:AbstractFloat}(κ::SquaredDistanceKernel{T}, x::T, y::T) = (x-y)^2

doc"SineSquaredKernel(p) = Σⱼsin²(p(xⱼ-yⱼ))"
immutable SineSquaredKernel{T<:AbstractFloat} <: NonNegNegDefAdditiveKernel{T}
    p::HyperParameter{T}
    SineSquaredKernel(p::Variable{T}) = new(
        HyperParameter(p, leftbounded(zero(T), :open))
    )
end
@outer_constructor(SineSquaredKernel, (π,))
@inline phi{T<:AbstractFloat}(κ::SineSquaredKernel{T}, x::T, y::T) = sin(κ.p*(x-y))^2

doc"ChiSquaredKernel() = Σⱼ(xⱼ-yⱼ)²/(xⱼ+yⱼ)"
immutable ChiSquaredKernel{T<:AbstractFloat} <: NonNegNegDefAdditiveKernel{T} end
ChiSquaredKernel() = ChiSquaredKernel{Float64}()
convert{T}(::Type{ChiSquaredKernel{T}}, ::ChiSquaredKernel) = ChiSquaredKernel{T}()
@inline function phi{T<:AbstractFloat}(κ::ChiSquaredKernel{T}, x::T, y::T)
    (x == y == zero(T)) ? zero(T) : (x-y)^2/(x+y)
end

#== Non-Negative Negative Definite Additive Kernel ==#

doc"ScalarProductKernel() = xᵀy"
immutable ScalarProductKernel{T<:AbstractFloat} <: AdditiveKernel{T} end
ScalarProductKernel() = ScalarProductKernel{Float64}()
convert{T}(::Type{ScalarProductKernel{T}}, ::ScalarProductKernel) = ScalarProductKernel{T}()
@inline phi{T<:AbstractFloat}(κ::ScalarProductKernel{T}, x::T, y::T) = x*y
@inline ismercer(::ScalarProductKernel) = true


text:
#####
##### simple comparison of various methods
#####

using MultistartOptimization
using NLopt
using Parameters

include("../test/test_functions.jl")

####
#### utility functions
####

function test(multistart_method, nlopt_method, f, N; options = ())
    n = 0
    function f̃(x)
        n += 1
        f(x)
    end
    F = MinimizationProblem(f̃, lower_bounds(f, N), upper_bounds(f, N))
    L = NLoptLocalMethod(nlopt_method; options...)
    @unpack location, value = multistart_minimization(multistart_method, L, F)
    (location = location, value = value, n = n)
end

function make_table(local_methods, functions, results; f_atol = 1e-3)
    codefence(str) = "`" * str * "`"
    table_body = map(x -> x.value ≤ f_atol ? string(x.n) : "**FAIL**", results)
    header_rows = [map(codefence ∘ string, local_methods)...]
    header_cols = [map(string ∘ Base.typename ∘ typeof, functions)...]
    table = ["" permutedims(header_cols);
             header_rows table_body]
end

function print_table(io, table)
    print_row(cells) = println(io, "| ", join(cells, " | "), " |")
    for row in axes(table, 1)
        print_row(table[row, :])
        if row == firstindex(table, 1)
            print_row(fill("----", size(table, 2)))
        end
    end
end

###
### tests
###

METHODS = (NLopt.LN_BOBYQA, NLopt.LN_NELDERMEAD, NLopt.LN_NEWUOA_BOUND, NLopt.LN_SBPLX,
           NLopt.LN_COBYLA, NLopt.LN_PRAXIS)

mm = TikTak(100)

T10 = [test(mm, method, f, 10; options = (xtol_abs = 1e-8, ))
       for method in METHODS, f in TEST_FUNCTIONS]

table = make_table(METHODS, TEST_FUNCTIONS, T10)
print_table(stdout, table)

ext: jl


text:
include("writeMat.jl")

nxs = 8
SP, X = PolynomialRing(Nemo.ZZ,["X$i" for i in 1:nxs])

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[2]+X[3];
       X[8] X[2]+X[4];
       X[1]+X[5] X[7];
       X[1]+X[6] X[8];
       X[7] X[1]+X[5];
       X[8] X[1]+X[6]]

toMaple(net, X, "M1.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[2]+X[3];
       X[8] X[2]+X[4];
       X[1]+X[5] X[7];
       X[1]+X[6] X[8]]

toMaple(net, X, "M1nonrev11.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[5] X[2]+X[3];
       X[8] X[2]+X[4];
       0 0;
       X[1]+X[6] X[8]]

toMaple(net, X, "M1nonrev01.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[7];
       0 0]

toMaple(net, X, "M1nonrev10.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0]

toMaple(net, X, "M1nonrev00.mpl")

max_stars_repo_name: LauraBMo/CRNTscripts

max_issues_repo_name: LauraBMo/CRNTscripts

max_forks_repo_name: LauraBMo/CRNTscripts


text:
include("../temperature_render.jl")
include("../iron_line_profile.jl")
include("../colourbar.jl")
using Measures
gr()

function combined_plot(;
                        mass = 1,
                        spin = 0.998,
                        obs_angle = 85.0,
                        tolerance = 1e-8,
                        dtmax = 1000.0,
                        size_multiplier::Int64 = 1,
                        resolution = 400,
                        fov = 3.0,
                        nbins=200
                        )
    plt, hmap, title, cache = iron_line_profile(
                                    mass = mass,
                                    spin = spin,
                                    obs_angle = obs_angle,
                                    tolerance = tolerance,
                                    dtmax = dtmax,
                                    size_multiplier = size_multiplier,
                                    resolution = resolution,
                                    fov = fov,
                                    nbins=nbins
                                    )

    verticle_line_x = fill(6.4, 100)
    verticle_line_y = LinRange(0, 1.25, 100)
    plot!(
        plt,
        verticle_line_x,
        verticle_line_y,
        linestyle=:dash,
        linecolor=:black,
        label=false,
        # xlims=(4,7.5),
        # ylims=(0,maximum(last(line_profiles*1.1))),
        lw=2
    )

    hmap = plot(hmap, grid=false,
    ticks=false,
    cbar=false,
    foreground_color_subplot=colorant"white"
    )
    plt = plot(plt, grid=false, legend=false, bottom_margin=5mm, framestyle=:box, ylims=(0,1.25))
    cbar = colourbar()
    # combined plot
    # subplot just for title
    titleplot = plot(title=title, grid = false, showaxis = false, bottom_margin = -50Plots.px, ticks=false)

    l = @layout[A{0.1h}; [B{};C{0.05h} ] D{0.4w} ]
    combined_plot = plot(titleplot, hmap, cbar, plt, layout=l, size=(1000,500))
    return combined_plot
end

# comb = combined_plot(
#                 mass=10,
#                 spin=0.998,
#                 obs_angle=15.0,
#                 tolerance=1e-12,
#                 size_multiplier=1,
#                 fov=10,
#                 dtmax=20
# )
# display(comb)


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: Example

import ..Example2
import FlatBuffers

FlatBuffers.@UNION(Any_, (
    Nothing,
    Monster,
    TestSimpleTableWithEnum,
    Example2.Monster,
))


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

@doc raw"""StatusCause provides more information about an api.Status failure, including cases when multiple errors are encountered.

    IoK8sApimachineryPkgApisMetaV1StatusCause(;
        field=nothing,
        message=nothing,
        reason=nothing,
    )

    - field::String : The field of the resource that has caused this error, as named by its JSON serialization. May include dot and postfix notation for nested attributes. Arrays are zero-indexed.  Fields may appear more than once in an array of causes due to fields having multiple errors. Optional.  Examples:   \"name\" - the field \"name\" on the current resource   \"items[0].name\" - the field \"name\" on the first array entry in \"items\"
    - message::String : A human-readable description of the cause of the error.  This field may be presented as-is to a reader.
    - reason::String : A machine-readable description of the cause of the error. If this value is empty there is no information available.
"""
mutable struct IoK8sApimachineryPkgApisMetaV1StatusCause <: SwaggerModel
    field::Any # spec type: Union{ Nothing, String } # spec name: field
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason

    function IoK8sApimachineryPkgApisMetaV1StatusCause(;field=nothing, message=nothing, reason=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusCause, Symbol("field"), field)
        setfield!(o, Symbol("field"), field)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusCause, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusCause, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1StatusCause

const _property_map_IoK8sApimachineryPkgApisMetaV1StatusCause = Dict{Symbol,Symbol}(Symbol("field")=>Symbol("field"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"))
const _property_types_IoK8sApimachineryPkgApisMetaV1StatusCause = Dict{Symbol,String}(Symbol("field")=>"String", Symbol("message")=>"String", Symbol("reason")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1StatusCause))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1StatusCause[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1StatusCause }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1StatusCause[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1StatusCause)
    true
end

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


text:
include("./imports.jl")

@testset ExtendedTestSet "ProgressBarLogger" begin
    learner = testlearner(coeff = 1, callbacks = [ProgressBarLogger()])
    @test_nowarn fit!(learner, 1)
end

#= Need to suppress stdout
@testset ExtendedTestSet "MetricsLogger" begin
    learner = testlearner(coeff = 1, callbacks = [MetricsLogger()])
    @test_nowarn fit!(learner, 1)
end
=#

@testset ExtendedTestSet "StopOnNanLoss" begin
    # is used by default
    learner = testlearner(coeff = NaN, callbacks = [StopOnNaNLoss()])
    # Epoch will be cancelled
    @test_nowarn fit!(learner, 1)
    @test learner.state.history.epochs == 0
end

@testset ExtendedTestSet "CustomCallback" begin
    cb = CustomCallback{EpochEnd,TrainingPhase}((learner) -> CancelFittingException("test"))
    learner = testlearner(coeff = 3, callbacks = [cb])
    fit!(learner, TrainingPhase())
    @test learner.state.history.epochs == 1
end


text:
module opesciSLIM

using PyCall

# function to prepend Python paths
function prependMyPyPath(d::String)
    mypath=Pkg.dir("opesciSLIM")
    myd=joinpath(mypath,d)
    unshift!(PyVector(pyimport("sys")["path"]), myd)
end

# prepend Python paths
prependMyPyPath("Python/operators/acoustic")
prependMyPyPath("Python/containers")

# submodule TimeModeling
include("TimeModeling/TimeModeling.jl")

# submodule LeastSquaresMigration
include("LeastSquaresMigration/LeastSquaresMigration.jl")

# submodule Optimization
include("Optimization/SLIM_optim.jl")

end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
module PhysicalConstants

using Measurements, Unitful

import Measurements: value, uncertainty
import Unitful: AbstractQuantity

# The type

"""
    PhysicalConstant{name,T,D,U} <: AbstractQuantity{T,D,U}

A type representing a physical constant, with units and -- optionally -- with
error as standard deviation.

Each type is a singleton and is parametrised by

* `name`: a `Symbol` representing its name
* `T`: the numerical type of the constant, e.g. `Float64`
* `D`: the physical dimension, from [`Unitful.jl`](https://painterqubits.github.io/Unitful.jl/stable/)
* `U`: the physical unit, from [`Unitful.jl`](https://painterqubits.github.io/Unitful.jl/stable/)

See [`@constant`](@ref) and [`@derived_constant`](@ref) for how to define a new
physical constant.
"""
struct PhysicalConstant{name,T,D,U} <: AbstractQuantity{T,D,U} end

_name(::PhysicalConstant{name,T,D,U}) where {name,T,D,U} = name

# Functions composing the building blocks of the macros

function _constant_preamble(name, sym, unit, def)
    ename = esc(name)
    qname = esc(Expr(:quote, name))
    esym = esc(sym)
    qsym = esc(Expr(:quote, sym))
    eunit = esc(unit)
    _bigconvert = isa(def,Symbol) ? quote
        function _big(::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U}
            c = BigFloat()
            ccall(($(string("mpfr_const_", def)), :libmpfr),
                  Cint, (Ref{BigFloat}, Int32), c, MPFR.ROUNDING_MODE[])
            return c
        end
    end : quote
        _big(::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} = $(esc(def))
    end
    return ename, qname, esym, qsym, eunit, _bigconvert
end

function _constant_begin(qname, ename, esym, eunit, val, _bigconvert)
    quote
        const $ename = PhysicalConstant{gensym($qname),Float64,dimension($eunit),typeof($eunit)}()
        export $ename
        const $esym = $ename

        Base.float(::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} = $val * $eunit
        Base.float(FT::DataType, ::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} =
            FT($val) * $eunit
        $_bigconvert
        Base.big(x::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} = _big(x) * $eunit
        Base.float(::Type{BigFloat}, x::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} = big(x)
    end
end

function _constant_end(qname, ename, qsym, descr, val, reference, eunit)
    quote
        Measurements.measurement(::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} =
            measurement(Float64, $ename)

        Unitful.unit(::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U}      = $eunit
        Unitful.dimension(::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} = D

        function Base.show(io::IO, x::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U}
            unc = uncertainty(ustrip(measurement($ename)))
            println(io, $descr, " (", $qsym, ")")
            println(io, "Value                         = ", float($ename))
            println(io, "Standard uncertainty          = ",
                    iszero(unc) ? "(exact)" : unc * $eunit)
            println(io, "Relative standard uncertainty = ",
                    iszero(unc) ? "(exact)" : round(unc / $val, sigdigits=2))
            print(io,   "Reference                     = ", $reference)
        end

        @assert isa(ustrip(float($ename)), Float64)
        @assert isa(ustrip(big($ename)), BigFloat)
        @assert isa(ustrip(measurement($ename)), Measurement{Float64})
        @assert isa(ustrip(measurement(Float32, $ename)), Measurement{Float32})
        @assert ustrip(float(Float64, $ename)) == Float64(ustrip(big($ename)))
        @assert ustrip(float(Float32, $ename)) == Float32(ustrip(big($ename)))
        @assert isapprox(Float64(value(ustrip(measurement(BigFloat, $ename)))),
                         value(ustrip(measurement($ename))), rtol = 1e-15)
        @assert Float64(uncertainty(ustrip(measurement(BigFloat, $ename)))) ==
            uncertainty(ustrip(measurement($ename)))
        @assert ustrip(big($ename)) == value(ustrip(measurement(BigFloat, $ename)))
    end
end

# Currently AbstractQuantity defines some operations in terms of the inner field `val`.

function Base.getproperty(c::PhysicalConstant{s,T,D,U}, sym::Symbol) where {s,T,D,U}
    if sym === :val
        return ustrip(float(T, c))
    else # fallback to getfield
        return getfield(c, sym)
    end
end

# @constant and @derived_constant macros

"""
    @constant(name, sym, descr, val, def, unit, unc, bigunc, reference) -> PhysicalConstant

Macro to define a new [`PhysicalConstant`](@ref).

The arguments are:

* `name`: a symbol with the long name of the constant.  This symbol is automatically eported
* `sym`: a non-exported alias for the constant
* `descr`: a description of the constant
* `val`: the numerical `Float64` value of the constant, in the reference units
* `def`:  the expression to use to compute the constant in any precision, in the reference units
* `unit`: the reference units
* `unc`: the numerical `Float64` value of the uncertainty, in the reference units
* `bigunc`: the expression to use to compute the uncertainty in any precision, in the reference units
* `reference`: the name of the reference

```jldoctest
julia> using PhysicalConstants, Unitful

julia> PhysicalConstants.@constant(MyConstant, mc, "A custom constant",
           12.34, BigFloat(1234) / BigFloat(100), u"m/s",
           0.56, BigFloat(56) / BigFloat(100), "My lab notebook")

julia> MyConstant
A custom constant (mc)
Value                         = 12.34 m s^-1
Standard uncertainty          = 0.56 m s^-1
Relative standard uncertainty = 0.045
Reference                     = My lab notebook
```
"""
macro constant(name, sym, descr, val, def, unit, unc, bigunc, reference)
    ename, qname, esym, qsym, eunit, _bigconvert = _constant_preamble(name, sym, unit, def)
    tag = Measurements.tag_counters[Base.Threads.threadid()] += 1
    quote
       $(_constant_begin(qname, ename, esym, eunit, val, _bigconvert))

        function Measurements.measurement(FT::DataType,
                                          ::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U}
            vl = FT($val)
            newder = Measurements.empty_der2(vl)
            if iszero($unc)
                return Measurement{FT}(vl, FT($unc), UInt64(0), newder) * $eunit
            else
                return Measurement{FT}(vl, FT($unc), $tag,
                                       Measurements.Derivatives(newder,
                                                                (vl, $unc, $tag)=>one(FT))) * $eunit
            end
        end
        function Measurements.measurement(::Type{BigFloat},
                                          x::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U}
            vl = _big(x)
            unc = BigFloat($bigunc)
            newder = Measurements.empty_der2(vl)
            if iszero($unc)
                return Measurement{BigFloat}(vl, unc, UInt64(0), newder) * $eunit
            else
                return Measurement{BigFloat}(vl, unc, $tag,
                                             Measurements.Derivatives(newder,
                                                                      (vl, unc, $tag)=>one(BigFloat))) * $eunit
            end
        end

        $(_constant_end(qname, ename, qsym, esc(descr), val, reference, eunit))
    end
end

"""
    @derived_constant(name, sym, descr, val, def, unit, unc, bigunc, reference) -> PhysicalConstant

Macro to define a new [`PhysicalConstant`](@ref) derived from another existing `PhysicalConstant`.

The arguments are:

* `name`: a symbol with the long name of the constant.  This symbol is automatically eported
* `sym`: a non-exported alias for the constant
* `descr`: a description of the constant
* `val`: the numerical `Float64` value of the constant, in the reference units
* `def`:  the expression to use to compute the constant in any precision, in the reference units
* `unit`: the reference units
* `measure64`: the numerical `Measurement{Float64}` value of the constant, in the reference units
* `measurebig`: the expression to use to compute the comstamt as a `Measurement{BigFloat}`, in the reference units
* `reference`: the name of the reference

```jldoctest
julia> using PhysicalConstants, Unitful, Measurements

julia> PhysicalConstants.@constant(MyConstant, mc, "A custom constant",
           12.34, BigFloat(1234) / BigFloat(100), u"m/s",
           0.56, BigFloat(56) / BigFloat(100), "My lab notebook")

julia> MyConstant
A custom constant (mc)
Value                         = 12.34 m s^-1
Standard uncertainty          = 0.56 m s^-1
Relative standard uncertainty = 0.045
Reference                     = My lab notebook

julia> PhysicalConstants.@derived_constant(MyDerivedConstant, mdc, "A custom derived constant",
           96.252, ustrip(big(mc)) * BigFloat(78) / BigFloat(10), u"m/s",
           measurement(mc) * 7.8, measurement(BigFloat, mc)  * BigFloat(78) / BigFloat(10),
           "My lab notebook")

julia> MyDerivedConstant
A custom derived constant (mdc)
Value                         = 96.252 m s^-1
Standard uncertainty          = 4.368 m s^-1
Relative standard uncertainty = 0.045
Reference                     = My lab notebook
```
"""
macro derived_constant(name, sym, descr, val, def, unit, measure64, measurebig, reference)
    ename, qname, esym, qsym, eunit, _bigconvert = _constant_preamble(name, sym, unit, def)
    quote
        $(_constant_begin(qname, ename, esym, eunit, val, _bigconvert))

        Measurements.measurement(::Type{Float64},
                                 ::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} =
                                     $(esc(measure64))
        Measurements.measurement(::Type{BigFloat},
                                 ::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} =
                                     $(esc(measurebig))
        Measurements.measurement(FT::DataType,
                                 x::PhysicalConstant{_name($ename),T,D,U}) where {T,D,U} =
                                     convert(Measurement{FT}, ustrip(measurement(x))) * $eunit

        $(_constant_end(qname, ename, qsym, esc(descr), val, reference, eunit))
    end
end

"""
    float(::PhysicalConstant{name,T,D,U}) where {T,D,U}
    float(FloatType, ::PhysicalConstant{name,T,D,U}) where {T,D,U}

Return the physical constant as a `Quantity` with the floating type optionally specified by
`FloatType`, `Float64` by default.

```jldoctest
julia> using PhysicalConstants.CODATA2018: G

julia> G
Newtonian constant of gravitation (G)
Value                         = 6.6743e-11 m^3 kg^-1 s^-2
Standard uncertainty          = 1.5e-15 m^3 kg^-1 s^-2
Relative standard uncertainty = 2.2e-5
Reference                     = CODATA 2018

julia> float(G)
6.6743e-11 m^3 kg^-1 s^-2

julia> float(Float32, G)
6.6743f-11 m^3 kg^-1 s^-2
```
"""
float(::PhysicalConstant)

"""
    measurement(::PhysicalConstant{name,T,D,U}) where {T,D,U}
    measurement(FloatType, ::PhysicalConstant{name,T,D,U}) where {T,D,U}

Return the physical constant as a `Quantity` with standard uncertainty.  The floating-point
precision can be optionally specified with the `FloatType`, `Float64` by default.

```jldoctest
julia> using PhysicalConstants.CODATA2018, Measurements

julia> import PhysicalConstants.CODATA2018: μ_0

julia> μ_0
Vacuum magnetic permeability (μ_0)
Value                         = 1.25663706212e-6 N A^-2
Standard uncertainty          = 1.9e-16 N A^-2
Relative standard uncertainty = 1.5e-10
Reference                     = CODATA 2018

julia> measurement(μ_0)
1.25663706212e-6 ± 1.9e-16 N A^-2

julia> measurement(Float32, μ_0)
1.256637e-6 ± 1.9e-16 N A^-2
```
"""
measurement(::PhysicalConstant)

include("promotion.jl")
include("codata2014.jl")
include("codata2018.jl")

end # module

max_issues_repo_name: JuliaPhysics/Constants

max_forks_repo_name: JuliaPhysics/Constants


text:
using Printf
using Oceananigans.Grids: size_summary, scalar_summary
using Oceananigans.Utils: prettysummary
using Oceananigans.BoundaryConditions: bc_str

location_str(::Type{Face})    = "Face"
location_str(::Type{Center})  = "Center"
location_str(::Type{Nothing}) = "⋅"
show_location(LX, LY, LZ) = "($(location_str(LX)), $(location_str(LY)), $(location_str(LZ)))"
show_location(field::AbstractField) = show_location(location(field)...)

function Base.summary(field::Field)
    LX, LY, LZ = location(field)
    prefix = string(size_summary(size(field)), " Field{$LX, $LY, $LZ}")

    grid_name = typeof(field.grid).name.wrapper
    reduced_dims = reduced_dimensions(field)

    suffix = reduced_dims === () ?
        string(" on ", grid_name, " on ", summary(architecture(field))) :
        string(" reduced over dims = ", reduced_dims,
               " on ", grid_name, " on ", summary(architecture(field)))

    return string(prefix, suffix)
end

data_summary(field) = string("max=", prettysummary(maximum(field)), ", ",
                             "min=", prettysummary(minimum(field)), ", ",
                             "mean=", prettysummary(mean(field)))

function Base.show(io::IO, field::Field)

    bcs = field.boundary_conditions

    prefix =
        string("$(summary(field))\n",
               "├── grid: ", summary(field.grid), '\n',
               "├── boundary conditions: ", summary(bcs), '\n',
               "│   └── west: ", bc_str(bcs.west), ", east: ", bc_str(bcs.east),
                     ", south: ", bc_str(bcs.south), ", north: ", bc_str(bcs.north),
                     ", bottom: ", bc_str(bcs.bottom), ", top: ", bc_str(bcs.top),
                     ", immersed: ", bc_str(bcs.immersed), '\n')

    middle = isnothing(field.operand) ? "" :
        string("├── operand: ", summary(field.operand), '\n',
               "├── status: ", summary(field.status), '\n')

    suffix = string("└── data: ", summary(field.data), '\n',
                    "    └── ", data_summary(field))

    print(io, prefix, middle, suffix)
end

Base.summary(status::FieldStatus) = "time=$(status.time)"

Base.summary(::ZeroField{N}) where N = "ZeroField{$N}"
Base.show(io::IO, z::ZeroField) = print(io, summary(z))

Base.show(io::IO, ::MIME"text/plain", f::AbstractField) = show(io, f)

const FieldTuple = Tuple{Field, Vararg{Field}}
const NamedFieldTuple = NamedTuple{S, <:FieldTuple} where S

function Base.show(io::IO, ft::NamedFieldTuple)
    names = keys(ft)
    N = length(ft)

    grid = first(ft).grid
    all_same_grid = true
    for field in ft
        if field.grid !== grid
            all_same_grid = false
        end
    end

    print(io, "NamedTuple with ", N, " Fields ")

    if all_same_grid
        print(io, "on ", summary(grid), ":\n")
    else
        print(io, "on different grids:", '\n')
    end

    for name in names[1:end-1]
        field = ft[name]
        print(io, "├── $name: ", summary(field), '\n')

        if !all_same_grid
            print(io, "│   └── grid: ", summary(field.grid), '\n')
        end
    end

    name = names[end]
    field = ft[name]
    print(io, "└── $name: ", summary(field))

    if !all_same_grid
        print(io, '\n')
        print(io, "    └── grid: ", summary(field.grid))
    end
end


text:
using MzPlots
using MzXML
using AbstractPlotting
using Test

@testset "MzPlots AbstractPlotting" begin
    mzxmldir = dirname(dirname(pathof(MzXML)))
    testfile = joinpath(mzxmldir, "test", "test64.mzXML")
    scans, info = MzXML.load(testfile)
    scene, ax = mzplot(scans)
    @test isa(scene, Scene)
end

import Plots
using MzPlots.MzCore
using MzCore.AxisArrays: Axis

@testset "MzPlots Plots" begin
    mzxmldir = dirname(dirname(pathof(MzXML)))
    testfile = joinpath(mzxmldir, "test", "test64.mzXML")
    scans, info = MzXML.load(testfile)
    lims = limits(scans)
    sz = (10, 2)
    asrange(axi::Axis{name}, n) where name = Axis{name}(range(minimum(axi.val), stop=maximum(axi.val), length=n))
    aa = AxisArray(zeros(sz), map(asrange, lims, sz)...)
    copyto!(aa, scans)
    p = Plots.plot(aa)
    @test isa(p, Plots.Plot)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

using .Optim

function _kliep_coeffs(K_nu, K_de, dre::KLIEP, optlib::Type{OptimLib})
  # retrieve parameters
  σ, b = dre.σ, size(K_de, 2)

  # number of numerator and denominator samples
  n_nu, n_de = size(K_nu, 1), size(K_de, 1)

  # constants for objective function
  K = K_nu

  # constants for equality constraints
  A = sum(K_de, dims=1)
  lc = uc = [n_de]

  # constants for inequality constraints
  T  = eltype(K_de)
  lx = fill(zero(T), b)
  ux = fill(Inf, b)

  # objective
  f(α) = -sum(log, K*α)
  function ∇f!(g, α)
    p = K*α
    for l in 1:b
      g[l] = -sum(K[j,l] / p[j] for j in 1:n_nu)
    end
  end
  function ∇²f!(h, α)
    p = K*α
    for k in 1:b, l in 1:b
      h[k,l] = sum(view(K,:,k) .* view(K,:,l) ./ p)
    end
  end

  # equality constraint
  c!(c, α)    = c  .= A*α
  J!(J, α)    = J  .= A
  H!(H, α, λ) = H .+= 0

  # initial guess
  αₒ = fill(n_de/sum(A), b)

  # optimization problem
  objective   = TwiceDifferentiable(f, ∇f!, ∇²f!, αₒ)
  constraints = TwiceDifferentiableConstraints(c!, J!, H!, lx, ux, lc, uc)
  initguess   = αₒ

  # solve problem with interior-point primal-dual Newton
  solution = optimize(objective, constraints, initguess, IPNewton())

  # optimal coefficients
  solution.minimizer
end


text:
using GPFlux
using Zygote
using Calculus
using Test, Random

@testset "power function (scalar)" begin
	Random.seed!(4)
	xs = 10*rand(4)
	αs = [-0.5, 3.5]

	for x in xs, α in αs
		dfdx, dfdα = gradient(GPFlux.power, x, α)
		Finit_dfdx, Finit_dfdα = Calculus.gradient(z->GPFlux.power(z[1], z[2]), [x, α])
		@test dfdx ≈ Finit_dfdx
		@test dfdα ≈ Finit_dfdα
	end
end

@testset "power function (Array)" begin
	Random.seed!(4)
	x = rand(3, 3)
	w = rand(3)
	f(x, α) = w'*GPFlux.power(x, α)*w

	dfdx, dfdα = gradient(f, x, 2.5)
	Finit_dfdα = Calculus.gradient(z->f(x, first(z)), [2.5])
	@test dfdα ≈ Finit_dfdα[1]
end

@testset "positive constraint" begin
	Random.seed!(4)
	x = randn(3, 5)
	w = rand(3)
	v = rand(5)
	f(x) = w'*GPFlux.positive(x)*v
	g(y) = w'*GPFlux.positive(reshape(y, 3, 5))*v

	dfdx, = gradient(f, x)
	Finit_dfdx = reshape(Calculus.gradient(g, vec(x)), 3, 5)
	@test dfdx ≈ Finit_dfdx
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
macro subset(d, con, verbose = false)
  # @show d, typeof(d)
  dname = string(d)
  con = string(con)
  # @show dname, con

  expr = con_dt_transform(con; dname=dname)
  verbose && @show expr
  esc(Meta.parse(expr))
end

export @subset


text:
using RAI
using Test
import HTTP, Arrow
using Mocking

Mocking.activate()

# -----------------------------------
# v2 transactions

make_patch(response) = @patch RAI.request(ctx::Context, args...; kw...) = response

const v2_async_response = HTTP.Response(200, [
        "Content-Type" => "application/json",
    ],
    body = """{"id":"1fc9001b-1b88-8685-452e-c01bc6812429","state":"CREATED"}""")

const v2_get_results_response() = join([
        "--8a89e52be8efe57f0b68ea75388314a3",
        "Content-Disposition: form-data; name=\"/:output/Int64\"; filename=\"/:output/Int64\"",
        "Content-Type: application/vnd.apache.arrow.stream",
        "",
        "\xff\xff\xff\xffx\0\0\0\x10\0\0\0\0\0\n\0\f\0\n\0\b\0\x04\0\n\0\0\0\x10\0\0\0\x01\0\x04\0\b\0\b\0\0\0\x04\0\b\0\0\0\x04\0\0\0\x01\0\0\0\x14\0\0\0\x10\0\x14\0\x10\0\0\0\x0e\0\b\0\0\0\x04\0\x10\0\0\0\x10\0\0\0\x18\0\0\0\0\0\x02\0\x1c\0\0\0\0\0\0\0\b\0\f\0\b\0\a\0\b\0\0\0\0\0\0\x01@\0\0\0\x02\0\0\0v1\0\0\xff\xff\xff\xff\x88\0\0\0\x14\0\0\0\0\0\0\0\f\0\x16\0\x14\0\x12\0\f\0\x04\0\f\0\0\0\b\0\0\0\0\0\0\0\x14\0\0\0\0\0\x03\0\x04\0\n\0\x18\0\f\0\b\0\x04\0\n\0\0\0\x14\0\0\08\0\0\0\x01\0\0\0\0\0\0\0\0\0\0\0\x02\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\b\0\0\0\0\0\0\0\0\0\0\0\x01\0\0\0\x01\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\x04\0\0\0\0\0\0\0\xff\xff\xff\xff\0\0\0\0",
    ], "\r\n")

const v2_get_relation_count_response() = join([
        "--8a89e52be8efe57f0b68ea75388314a3",
        "Content-Disposition: form-data; name=\"relation-count\"; filename=\"\"",
        "Content-Type: application/json",
        "",
        "{\"relation_count\":1}",
    ], "\r\n")

const v2_get_transaction_results_response = HTTP.Response(200, [
        "Content-Type" => "Content-Type: multipart/form-data; boundary=8a89e52be8efe57f0b68ea75388314a3",
        "Transfer-Encoding" => "chunked",
    ],
    body = join([
    "",
    v2_get_results_response(),
    v2_get_relation_count_response(),
    "--8a89e52be8efe57f0b68ea75388314a3--",
    "",
], "\r\n"))

const v2_fastpath_response = HTTP.Response(200, [
        "Content-Type" => "Content-Type: multipart/form-data; boundary=8a89e52be8efe57f0b68ea75388314a3",
        "Transfer-Encoding" => "chunked",
    ],
    body = join([
    "",
    "--8a89e52be8efe57f0b68ea75388314a3",
    "Content-Disposition: form-data; name=\"transaction\"; filename=\"\"",
    "Content-Type: application/json",
    "",
    """{"id":"a3e3bc91-0a98-50ba-733c-0987e160eb7d","results_format_version":"2.0.1","state":"COMPLETED"}""",
    "--8a89e52be8efe57f0b68ea75388314a3",
    "Content-Disposition: form-data; name=\"metadata\"; filename=\"\"",
    "Content-Type: application/json",
    "",
    """[{"relationId":"/:output/Int64","types":[":output","Int64"]}]""",
    "--8a89e52be8efe57f0b68ea75388314a3",
    "Content-Disposition: form-data; name=\"problems\"; filename=\"\"",
    "Content-Type: application/json",
    "",
    """[]""",
    v2_get_results_response(),
    "--8a89e52be8efe57f0b68ea75388314a3--",
    "",
], "\r\n"))

function make_arrow_table(vals)
    io = IOBuffer()
    Arrow.write(io, (v1=vals,))
    seekstart(io)
    return Arrow.Table(io)
end

@testset "exec_async" begin
    ctx = Context("region", "scheme", "host", "2342", nothing)

    @testset "async response" begin
        patch = make_patch(v2_async_response)

        apply(patch) do
            rsp = RAI.exec_async(ctx, "engine", "database", "2+2")
            @test rsp == Dict(
                "transaction" => JSON3.read("""{"id":"1fc9001b-1b88-8685-452e-c01bc6812429","state":"CREATED"}""")
            )
        end
    end

    @testset "sync response" begin
        patch = make_patch(v2_fastpath_response)

        apply(patch) do
            rsp = RAI.exec_async(ctx, "engine", "database", "2+2")
            @test rsp["transaction"] == JSON3.read("""{
                    "id": "a3e3bc91-0a98-50ba-733c-0987e160eb7d",
                    "results_format_version": "2.0.1",
                    "state": "COMPLETED"
                }""")
            @test rsp["metadata"] == [JSON3.read("""{
                "relationId": "/:output/Int64",
                    "types": [
                                ":output",
                                "Int64"
                            ]
            }""")]
            @test rsp["problems"] == Union{}[]

            # Test for the expected arrow data:
            expected_data = make_arrow_table([4])
            # Arrow.Tables can't be compared via == (https://github.com/apache/arrow-julia/issues/310)
            @test length(rsp["results"]) == 1
            @test rsp["results"][1][1] == "/:output/Int64"
            @test collect(rsp["results"][1][2]) == collect(expected_data)
        end
    end

    @testset "get_transaction_results" begin
        patch = make_patch(v2_get_transaction_results_response)

        apply(patch) do
            rsp = RAI.get_transaction_results(ctx, "fake-txn-id")
            @test !isempty(rsp)
            @test rsp[1][2] isa Arrow.Table
        end
    end
end

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
using PMDPs
using Test
using POMDPSimulators

using StaticArrays
using POMDPs
using MCTS, DiscreteValueIteration
using StatsBase
using Plots
using POMDPSimulators
using D3Trees
using POMDPPolicies
using POMDPLinter
using Random
using DataFrames
using POMDPSimulators

using DrWatson
using Distributions
using BeliefUpdaters

println("Processing")
mdp_params = Dict(pairs( (n_edges = 3, c_init = 2, demand = Float64[1,1,1], selling_horizon_end = [40,45,50], actions = 15:5:90, objective=:revenue)))
#mdp_params = Dict(pairs( (n_edges = 2, c_init = 1, demand = Float64[1,1], selling_horizon_end = [20,25], actions = 15:5:30)))
# mcts_params = Dict(solver= MCTSSolver, n_iterations=1000, depth=30, exploration_constant=40.0, reuse_tree=true)
mcts_params = Dict(pairs( (solver= DPWSolver, n_iterations=50, depth=30, exploration_constant=40.0, enable_state_pw = true, keep_tree=true, show_progress=false)))
exp_params = Dict(pairs( (n_runs = 20, vi=true, save=:stats) ))
params = Dict(:mdp=>mdp_params, :mcts=>mcts_params, :exp=>exp_params)

result, filepath =  makesim(params);

r = result[:r]
u = result[:u]
result[:t]
println()
display("rewards: $r")
display("utilization: $u")

# function get_trace(h)
#     [rec for rec in collect(eachstep(h, "s, a, r, info")) if (sum(rec.s.p)>0 || rec.s.t==length(h)-1)]
# end

# h_i = nothing
# if h_i ≠ nothing
#     h_mc = hs_mc[h_i]
#     run_vi ? h_vi = hs_vi[h_i] : nothing

#     display((rs[h_i], us[h_i]))
#     display(get_trace(h_mc))
#     run_vi ? display(get_trace(h_vi)) : nothing
# end

# # Load results
# using Distributions
# using BeliefUpdaters
# res_foldername = readdir(datadir("sims", "exp_small_u"))[1]
# exp_file = readdir(datadir("sims", "tst"))[1]
# res = wload(datadir("sims", "tst", "tst.bson"))

# # # analyze results
# using DataFrames

# df = collect_results(
#     datadir("sims", "exp_small_u"),
#     white_list = [:r, :u, :t, :params], subfolders=true
# )

# rng = MersenneTwister(123)
# RolloutEstimator(RandomSolver(rng))


text:
# using Revise
# using Neo4j

# Debug logging
# using DistributedFactorGraphs
using Test
using Logging
using Neo4j
using DistributedFactorGraphs
using IncrementalInference
logger = SimpleLogger(stdout, Logging.Debug)
global_logger(logger)

dfg = CloudGraphsDFG{SolverParams}("localhost", 7474, "neo4j", "test",
                            "testUser", "testRobot", "testSession",
                            nothing,
                            nothing,
                            IncrementalInference.decodePackedType,
                            IncrementalInference.rebuildFactorMetadata!,
                            solverParams=SolverParams())
user = User(Symbol(dfg.userId), "Bob Zack", "Description", Dict{Symbol, String}())
@test_throws Exception createUser(dfg, user) == user

isFullyConnected(dfg)
hasOrphans(dfg)
dfgSubgraph = getSubgraphAroundNode(dfg, getVariable(dfg, :x1), 2)


text:
using Pkg

# this script installs all the packages necessary to reproduce the paper
# except for GMT.jl, which requires a local copy of GMT
# add following packages
pkgs = ["AWSCore", "AWSS3", "DataFrames", "CSV", "Glob", "JLD2",
        "Plots", "CuArrays", "CUDAnative", "CUDAdrv", "BenchmarkTools",
        "PyCall", "Conda","StatsPlots", "HTTP"]
for p in pkgs
      Pkg.add(p)
end
Pkg.add(PackageSpec(name="SeisIO",rev="master"))
Pkg.add(PackageSpec(name="SeisNoise",rev="master"))

using Conda
Conda.add("numpy")
Conda.add("scipy")
Conda.add("obspy",channel="conda-forge")


text:
import LinearAlgebra

function _pseudoinverse(m::AbstractMatrix)
    try
        return LinearAlgebra.inv(m)
    catch ex
        @debug("", exception = (ex, catch_backtrace()))
    end
    return LinearAlgebra.pinv(m)
end


text:
abstract type AbstractCIFAR <: Name end
struct CIFAR10 <: AbstractCIFAR end
struct CIFAR20 <: AbstractCIFAR end
struct CIFAR100 <: AbstractCIFAR end

# download options
task(::Type{<:AbstractCIFAR}) = MultiClass
format(::Type{<:AbstractCIFAR}) = ColorImages
source(::Type{<:AbstractCIFAR}) = "https://www.cs.toronto.edu/~kriz/cifar.html"
make_datadep(::Type{<:AbstractCIFAR}) = nothing
datadepname(::Type{CIFAR10}) = "CIFAR10"
datadepname(::Type{<:AbstractCIFAR}) = "CIFAR100"

loadraw(::Type{CIFAR10}, type) = loadraw_mldatasets(MLDatasets.CIFAR10, type)

function loadraw(::Type{CIFAR20}, type, T = Float32)
    if type == :train
        x, y, ~ = MLDatasets.CIFAR100.traindata(T)
    elseif type == :test
        x, y, ~ = MLDatasets.CIFAR100.testdata(T)
    else
        error("$(type) dataset not provided for $(nameof(N))")
    end
    return x, y
end

function loadraw(::Type{CIFAR100}, type, T = Float32)
    if type == :train
        x, ~, y = MLDatasets.CIFAR100.traindata(T)
    elseif type == :test
        x, ~, y = MLDatasets.CIFAR100.testdata(T)
    else
        error("$(type) dataset not provided for $(nameof(N))")
    end
    return x, y
end

nclasses(::Type{CIFAR10}) = 10
nclasses(::Type{CIFAR20}) = 20
nclasses(::Type{CIFAR100}) = 100
classes(::Type{CIFAR10}) = collect(0:9)
classes(::Type{CIFAR20}) = collect(0:19)
classes(::Type{CIFAR100}) = collect(0:69)
classes_orig(::Type{CIFAR10}) = collect(0:9)
classes_orig(::Type{CIFAR20}) = collect(0:19)
classes_orig(::Type{CIFAR100}) = collect(0:99)
nattributes(::Type{<:AbstractCIFAR}) = (32, 32, 3)
ninstances(::Type{<:AbstractCIFAR}) = (60000, 0, 10000)

# dataset description
name(N::Type{<:AbstractCIFAR}) = "$(nameof(N))"
author(::Type{<:AbstractCIFAR}) = ["Alex Krizhevsky", "Vinod Nair", "Geoffrey Hinton"]

function citation(::Type{<:AbstractCIFAR})
    return """
    @article{krizhevsky2009learning,
        title = {Learning multiple layers of features from tiny images},
        author = {Krizhevsky, Alex and Hinton, Geoffrey and others},
        year = {2009},
        publisher = {Citeseer}
    }
    """
end

max_stars_repo_name: VaclavMacha/DatasetProvider

max_issues_repo_name: VaclavMacha/DatasetProvider

max_forks_repo_name: VaclavMacha/DatasetProvider


text:
# Hofbauer invariant measure h with C=1
heaviside(x::AbstractFloat) = ifelse(x < 0, zero(x), ifelse(x > 0, one(x), oftype(x,0.5)))
hofbauer(x,β,N=100) = heaviside(1-x) + [β^-j for j in 1:N]'*heaviside.(bernoulli_map(1e0,β,N)[2:end] .- x)

max_stars_repo_name: milankl/BernoulliMap

max_issues_repo_name: milankl/BernoulliMap

max_forks_repo_name: milankl/BernoulliMap


text:
"""
$(TYPEDEF)

Message prior on all dimensions of a variable node in the factor graph.

Notes
- Only temporary existance during CSM operations.
"""
struct MsgPrior{T <: SamplableBelief} <: AbstractPrior
  Z::T
  inferdim::Float64
end

# MsgPrior{T}() where {T} = new{T}()
# MsgPrior{T}(z::T, infd::R) where {T <: SamplableBelief, R <: Real} = new{T}(z, infd)
# function MsgPrior(z::T, infd::R) where {T <: SamplableBelief, R <: Real}
#     MsgPrior{T}(z, infd)
# end

getSample(cf::CalcFactor{<:MsgPrior}, N::Int=1) = (reshape(rand(cf.factor.Z,N),:,N), )

# this is a developmental type, will be standardized after conclusion of #1010
# TODO resolve type instability
const MsgRelativeType = Vector{NamedTuple{(:variables, :likelihood), Tuple{Vector{Symbol},DFG.AbstractRelative}}}

const MsgPriorType = Dict{Symbol, MsgPrior{BallTreeDensity}}

function (cfo::CalcFactor{<:MsgPrior})(z, x1)
  return  z .- x1
end

struct PackedMsgPrior <: PackedInferenceType
  Z::String
  inferdim::Float64
end

function convert(::Type{PackedMsgPrior}, d::MsgPrior)
  PackedMsgPrior(convert(PackedSamplableBelief, d.Z), d.inferdim)
end
function convert(::Type{<:MsgPrior}, d::PackedMsgPrior)
  MsgPrior(convert(SamplableBelief, d.Z), d.inferdim)
end

"""
    $TYPEDEF

Internal development types used during consolidation.  Stores relative and prior information making up a joint likelihood
message passed upward on the Bayes tree.
"""
mutable struct _MsgJointLikelihood
  relatives::IIF.MsgRelativeType
  priors::IIF.MsgPriorType
end

_MsgJointLikelihood(;relatives::MsgRelativeType=MsgRelativeType(),priors::MsgPriorType=MsgPriorType() ) = _MsgJointLikelihood(relatives, priors)

function Base.show(io::IO, x::_MsgJointLikelihood)
  println(io, )
  printstyled(io, " _MsgJointLikelihood:\n", color=:blue)
  print(io, "  .relatives: ")
  for tp in x.relatives
    print(io, tp.variables, "::", typeof(tp.likelihood).name)
    print(io, "; ")
  end
  println(io)
  print(io, "  .priors: ")
  for k in keys(x.priors)
    print(io, k, ", ")
  end
  println(io)
end

Base.show(io::IO, ::MIME"text/plain", x::_MsgJointLikelihood) = show(io, x)

"""
  $(TYPEDEF)
Belief message for message passing on the tree.  This should be considered an incomplete joint probility.

Notes:
- belief -> Dictionary of [`TreeBelief`](@ref)
- variableOrder -> Ordered variable id list of the seperators in cliqueLikelihood
- cliqueLikelihood -> marginal distribution (<: `SamplableBelief`) over clique seperators.
- Older names include: productFactor, Fnew, MsgPrior, LikelihoodMessage

DevNotes:
- Used by both nonparametric and parametric.
- Objective for parametric case: `MvNormal(μ=[:x0;:x2;:l5], Σ=[+ * *; * + *; * * +])`.
- Part of the consolidation effort, see #459.
- Better conditioning for joint structure in the works using deconvolution, see #579, #635.
  - TODO confirm why <: Singleton.

$(TYPEDFIELDS)
"""
mutable struct LikelihoodMessage{T <: MessageType} <: AbstractPrior
  status::CliqStatus
  belief::Dict{Symbol, TreeBelief} # will eventually be deprecated
  variableOrder::Vector{Symbol}
  cliqueLikelihood::Union{Nothing,SamplableBelief}
  msgType::T
  hasPriors::Bool
  # this is different from belief[].inferdim, as the total available infer dims remaining during down msgs -- see #910
  childSolvDims::Dict{Int, Float64}
  # calc differential factors for joint in the child clique
  jointmsg::_MsgJointLikelihood
  # diffJoints::Vector{NamedTuple{(:variables, :likelihood), Tuple{Vector{Symbol},DFG.AbstractRelative}}}
end

LikelihoodMessage(; status::CliqStatus=NULL,
                    beliefDict::Dict=Dict{Symbol, TreeBelief}(),
                    variableOrder::Vector{Symbol}=Symbol[],
                    cliqueLikelihood=nothing,
                    msgType::T=NonparametricMessage(),
                    hasPriors::Bool=true,
                    childSolvDims::Dict{Int, Float64}=Dict{Int, Float64}(),
                    jointmsg::_MsgJointLikelihood=_MsgJointLikelihood(),
                  ) where {T <: MessageType} =
        LikelihoodMessage{T}(status, beliefDict, variableOrder, cliqueLikelihood, msgType, hasPriors, childSolvDims, jointmsg)
#

function Base.show(io::IO, ::MIME"text/plain", msg::LikelihoodMessage)
  t = typeof(msg)
  fields = fieldnames(t)
  nf = nfields(msg)
  println(io, "LikelihoodMessage:")
  for f in fields
    printstyled(io, f,": ", color=:blue)
    show(io, getproperty(msg, f))
    println(io)
  end

end

function compare( l1::LikelihoodMessage,
                  l2::LikelihoodMessage;
                  skip::Vector{Symbol}=[] )
  #
  TP = true
  TP = TP && l1.status == l2.status
  TP = TP && l1.variableOrder == l2.variableOrder
  TP = TP && l1.msgType == l2.msgType
  TP = TP && l1.cliqueLikelihood |> typeof == l2.cliqueLikelihood |> typeof
  for (k,v) in l1.belief
    TP = TP && haskey(l2.belief, k)
    TP = TP && compare(v, l2.belief[k])
  end
end

==(l1::LikelihoodMessage,l2::LikelihoodMessage) = compare(l1,l2)

## =========================================================================================
## DEPRECATE BELOW AS ABLE
## =========================================================================================

# figure out how to deprecate (not critical at the moment)
# used in RoMEPlotting
const TempUpMsgPlotting = Dict{Symbol,Vector{Tuple{Symbol, Int, BallTreeDensity, Float64}}}


text:
include("Synapses.jl")
include("DBFs.jl")
include("PointArrays.jl")
include("SWCs.jl")
include("NodeNets.jl")
#include("Manifests.jl")

include("Segments.jl")
include("Neurons.jl")
include("NBLASTs.jl")

include("Utils/include.jl")


text:
type Mul <: Op; alpha; beta; Mul(;alpha=1,beta=1,o...)=new(alpha,beta); end

ninputs(::Mul)=2
canoverwrite(::Mul)=false
back_reads_x(::Mul)=true
back_reads_y(::Mul)=false

function forw(l::Mul, x1, x2, y; o...)
    x2 == nothing || size(x2) == size(y) || throw(DimensionMismatch("x2 and y should have the same size in mul"))
    x1 == nothing || size(x1) == size(y) || (x1 = reshape_to_match(x1,y))
    if x1==nothing # nothing represents zero
        l.alpha == 1 || error("mul does not support nothing^p for p!=1")
        nothing
    elseif x2==nothing
        l.beta  == 1 || error("mul does not support nothing^p for p!=1")
        nothing
    else
        bmul!(l.alpha,x1,l.beta,x2,y)
    end
end

# z = mul(x,y; alpha,beta) = x^alpha * y^beta
# dJ/dx = dJ/dz dz/dx = dJ/dz * y^beta * x^(alpha-1) * alpha
# dJ/dy = dJ/dz dz/dy = dJ/dz * x^alpha * y^(beta-1) * beta

# When broadcasting:
# zi = xj^alpha * yi^beta   where multiple i map to the same j
# dJ/dyi = dJ/dzi dzi/dyi = dJ/dzi * xj^alpha * yi(beta-1) * beta   same as before, bmul! is sufficient.
# dJ/dxj = sum_i dJ/dzi dzi/dxj = sum_i dJ/dzi * yi^beta * xj^(alpha-1) * alpha
# For dJ/dzi * yi^beta mul2 is sufficient
# For sum_i we need row sums, baddback?
# For * xj^(alpha-1) * alpha we need regular mul2 and scaling

function back(l::Mul, dy, dx1, dx2; x=nothing, o...)
    (x1,x2) = (length(x)==2 ? x : error("back(mul) needs inputs x1 and x2"))
    dx2 == nothing || size(dx2) == size(dy) || throw(DimensionMismatch("x2 and y should have the same size in mul"))
    x2  == nothing || size(x2)  == size(dy) || throw(DimensionMismatch("x2 and y should have the same size in mul"))
    dx1 == nothing || size(dx1) == size(dy) || (dx1 = reshape_to_match(dx1,dy))
    x1  == nothing || size(x1)  == size(dy) || (x1  = reshape_to_match(x1,dy))
    x1  != nothing || l.alpha == 1 || error("mul does not support nothing^p for p!=1")
    x2  != nothing || l.beta  == 1 || error("mul does not support nothing^p for p!=1")

    if dx2 != nothing  # dJ/dy = dJ/dz dz/dy = dJ/dz * x^alpha * y^(beta-1) * beta
        if x1 == nothing      # representing zero
            fillsync!(dx2, 0)
        else
            bmul!(l.alpha,x1,1,dy,dx2)
            l.beta == 1 || scale!(l.beta, bmul!(l.beta-1,x2,1,dx2,dx2))
        end
    end

    if dx1 != nothing  # dJ/dxj = sum_i dJ/dzi dzi/dxj = sum_i dJ/dzi * yi^beta * xj^(alpha-1) * alpha
        if x2 == nothing
            fillsync!(dx1, 0)
        elseif size(dx1) == size(dy)
            bmul!(l.beta,x2,1,dy,dx1)
            l.alpha == 1 || scale!(l.alpha, bmul!(1,dx1,l.alpha-1,x1,dx1))
        else
            tmp = similar(dy)
            bmul!(l.beta,x2,1,dy,tmp)
            baddback!(tmp,dx1)
            free(tmp)
            l.alpha == 1 || scale!(l.alpha, bmul!(l.alpha-1,x1,1,dx1,dx1))
        end
    end
    gpusync()
end

### bmul! broadcasting multiplication: c=a^alpha * b^beta

@gpu function bmul!{T}(alpha::Number,a::CudaArray{T},beta::Number,b::CudaArray{T},c::CudaArray{T})
    size(b) == size(c) || throw(DimensionMismatch("b and c should have the same size in bmul!"))
    if size(a) == size(b)
        T <: Float32 ? ccall((:mul32,libknet),Void,(Cint,Cfloat,Ptr{Cfloat},Cfloat,Ptr{Cfloat},Ptr{Cfloat}),length(a),T(alpha),a,T(beta),b,c) :
        T <: Float64 ? ccall((:mul64,libknet),Void,(Cint,Cdouble,Ptr{Cdouble},Cdouble,Ptr{Cdouble},Ptr{Cdouble}),length(a),T(alpha),a,T(beta),b,c) :
        error("$T not supported")
    else
        ndims(b) <= 8 || error("mul kernel supports dimensions up to 8")
        a = reshape_to_match(a,b)
        T <: Float32 ? ccall((:bmul32,libknet),Void,(Cint,Cfloat,Ptr{Cint},Ptr{Cfloat},Cfloat,Ptr{Cint},Ptr{Cfloat},Ptr{Cfloat}),ndims(a),T(alpha),cudadims(a),a,T(beta),cudadims(b),b,c) :
        T <: Float64 ? ccall((:bmul64,libknet),Void,(Cint,Cdouble,Ptr{Cint},Ptr{Cdouble},Cdouble,Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble}),ndims(a),T(alpha),cudadims(a),a,T(beta),cudadims(b),b,c) :
        error("$T not supported")
    end
    gpusync()
    return c
end

function cudadims(a)
    n = ndims(a)
    s = zeros(Cint,n)
    for i=1:n; s[i]=size(a,i); end
    CudaArray(s)
end

# infersize for mul is exactly the same as add except for 1-D x1 with CUDNN_ADD_SAME_C trick:
function infersize(a::Mul, x1, x2, y)
    if x1==x2==y==nothing
        nothing
    elseif x2==y==nothing
        (x1, nothing, nothing)
    elseif x2==nothing
        infersize(a, x1, y, y)
    elseif y==nothing
        infersize(a, x1, x2, x2)
    else
        ydims = map(x2, y) do xdim, ydim
            xdim == ydim ? xdim :
            xdim == 0 ? ydim :
            ydim == 0 ? xdim :
            throw(DimensionMismatch())
        end
        x1==nothing && (return (ydims, ydims, ydims))
        length(x1) <= length(ydims) || throw(DimensionMismatch())
        xdims = [x1...]; ydims = [ydims...]
        for i=1:length(xdims)
            xdims[i] == ydims[i] ? continue :
            xdims[i] == 1 ? continue :
            xdims[i] == 0 ? continue :
            ydims[i] == 0 ? (ydims[i] = xdims[i]) :
            throw(DimensionMismatch())
        end
        xdims = tuple(xdims...); ydims = tuple(ydims...)
        return (xdims, ydims, ydims)
    end
end

# function infersize(m::Mul, x1, x2, y)
#     if x1==x2==y==nothing
#         nothing
#     elseif x1==nothing || x2==nothing || y==nothing
#         n = length(x1!=nothing ? x1 : x2!=nothing ? x2 : y!=nothing ? y : error())
#         x1 == nothing && (x1 = ntuple(i->0,n))
#         x2 == nothing && (x2 = ntuple(i->0,n))
#         y == nothing && (y = ntuple(i->0,n))
#         infersize(m,x1,x2,y)
#     else
#         length(x1)==length(x2)==length(y) || throw(DimensionMismatch())
#         dims = map(x1,x2,y) do a,b,c
#             n = 0
#             a==0 || a==n ? nothing : n==0 ? n=a : throw(DimensionMismatch())
#             b==0 || b==n ? nothing : n==0 ? n=b : throw(DimensionMismatch())
#             c==0 || c==n ? nothing : n==0 ? n=c : throw(DimensionMismatch())
#             n
#         end
#         (dims, dims, dims)
#     end
# end

# ### mul2 element-wise multiplication:

# mul2!(alpha::Number,a::Array,beta::Number,b::Array,c::Array)=(for i=1:length(c); c[i] = a[i]^alpha*b[i]^beta; end; c)
# mul2!(alpha::Number,a::CudaArray{Float32},beta::Number,b::CudaArray{Float32},c::CudaArray{Float32})=(ccall((:mul2_32,libknet),Void,(Cint,Cfloat,Ptr{Cfloat},Cfloat,Ptr{Cfloat},Ptr{Cfloat}),length(a),Cfloat(alpha),a,Cfloat(beta),b,c); gpusync(); c)
# mul2!(alpha::Number,a::CudaArray{Float64},beta::Number,b::CudaArray{Float64},c::CudaArray{Float64})=(ccall((:mul2_64,libknet),Void,(Cint,Cdouble,Ptr{Cdouble},Cdouble,Ptr{Cdouble},Ptr{Cdouble}),length(a),Cdouble(alpha),a,Cdouble(beta),b,c); gpusync(); c)

# mul(x1,x2,y; alpha=1, beta=1)=(Mul(alpha,beta),x1,x2,y)


text:
@doc raw"""
    conv_v2h(w, v)

Internal function used to compute inputs from a visible configurations `v` to the hidden
layer, where `w` are the convolutional RBM weights.

```math
I_\mu^{k_1,\dots,k_n,b} = \sum_{c,j_1,\dots,j_n} w_{c,j_1,\dots,j_n,\mu} v_{c,j_1+k_1-1,\dots,j_n+k_n-1,b}
```

Assumes that:

* `w` is of size `(C,J₁,...,Jₙ,M)`
* `v` is of size `(C,N₁,...,Nₙ,B)`

Here `C` is the flattened channel dimension,
`M` is the number of hidden units and `B` is the batch size.
Therefore the hidden and batch dimensions must be flattened before calling `conv_v2h`.
The output `I` is of size `(M, N₁ - J₁ + 1, ..., Nₙ - Jₙ + 1, B)`.

!!! warning
    This is an internal function and is not part of the public API.

!!! warning
    Only works for `n = 1, 2, 3` due to a technical limitations.
"""
function conv_v2h end

function conv_v2h(w::AbstractTensor{3}, v::AbstractTensor{3})
    _conv_v2h_check_size(w, v)
    @tullio I[μ,k,b] := w[c,j,μ] * v[c,j+k-1,b]
end

function conv_v2h(w::AbstractTensor{4}, v::AbstractTensor{4})
    _conv_v2h_check_size(w, v)
    @tullio I[μ,k₁,k₂,b] := w[c,j₁,j₂,μ] * v[c,j₁+k₁-1,j₂+k₂-1,b]
end

function conv_v2h(w::AbstractTensor{5}, v::AbstractTensor{5})
    _conv_v2h_check_size(w, v)
    @tullio I[μ,k₁,k₂,k₃,b] := w[c,j₁,j₂,j₃,μ] * v[c,j₁+k₁-1,j₂+k₂-1,j₃+k₃-1,b]
end

function conv_v2h(w::AbstractTensor{6}, v::AbstractTensor{6})
    _conv_v2h_check_size(w, v)
    @tullio I[μ,k₁,k₂,k₃,k₄,b] := w[c,j₁,j₂,j₃,j₄,μ] * v[c,j₁+k₁-1,j₂+k₂-1,j₃+k₃-1,j₄+k₄-1,b]
end

function _conv_v2h_check_size(w::AbstractTensor{N}, v::AbstractTensor{N}) where {N}
    @assert size(w, 1) == size(v, 1) # channel size
    @assert all((size(w) .≤ size(v))[2:(end - 1)])
end

@doc raw"""
    conv_h2v(w, v)

Internal function used to compute inputs from a hidden configurations `h` to the visible
layer, where `w` are the convolutional RBM weights.

```math
I_{i_1,\dots,i_n}^{k_1,\dots,k_n,b} = \sum_{c,j_1,\dots,j_n} w_{c,j_1,\dots,j_n,\mu} v_{c,j_1+k_1-1,\dots,j_n+k_n-1,b}
```

Assumes that:

* `w` is of size `(C,J₁,...,Jₙ,M)`
* `h` is of size `(M,K₁,...,Kₙ,B)`

Here `C` is the channel dimension, `M` is the number of hidden units and `B` is the batch size.
These three dimensions are flat, so that the hidden and batch dimensions must be
flattened before calling `conv_h2v`.
The output `I` is of size `(C, K₁ + J₁ - 1, ..., Kₙ + Jₙ - 1, B)`.

!!! warning
    This is an internal function and is not part of the public API.

!!! warning
    Only works for `n = 1, 2, 3, 4`.
    Due to a technical limitation of Tullio.jl
    ([#129](https://github.com/mcabbott/Tullio.jl/issues/129)),
    `conv_h2v` is defined by hand for different values of `n`.
    Therefore only a small number can be supported.
    If you need some higher value of `n`, consider doing a PR to define the corresponding
    `conv_h2v` method.
"""
function conv_h2v(w::AbstractTensor{N}, h::AbstractTensor{N}) where {N}
    Isz = (size(w,1), (size(w) .+ size(h) .- 1)[2:(end - 1)]..., size(h)[end])
    I = zeros(promote_type(eltype(w), eltype(h)), Isz)
    return conv_h2v!(I, w, h)
end

function conv_h2v!(I::AbstractTensor{N}, w::AbstractTensor{N}, h::AbstractTensor{N}) where {N}
    kernel = size(w)[(begin + 1):(end - 1)]
    output = size(h)[(begin + 1):(end - 1)]
    inputs = kernel .+ output .- 1
    @assert size(I) == (size(w, 1), inputs..., size(h, N))
    @assert size(w) == (size(w, 1), kernel..., size(h, 1))
    @assert size(h) == (size(h, 1), output..., size(h, N))
    for b in axes(h, N), k in CartesianIndices(output), j in CartesianIndices(kernel), μ in axes(h, 1), c in axes(w, 1)
        i = j + k - oneunit(j)
        I[c,i,b] += w[c,j,μ] * h[μ,k,b]
    end
    return I
end

function ChainRulesCore.rrule(::typeof(conv_h2v), w::AbstractTensor{N}, h::AbstractTensor{N}) where {N}
    function conv_h2v_pullback(δI)
        @tullio δw[c,j,μ] = δI[c,k + j - 1,b] * h[μ,k,b]
        return (NoTangent(), δw, NoTangent()) # note we drop gradients of `h`!
    end
    return conv_h2v(w, h), conv_h2v_pullback
end

function _conv_h2v_check_size(I::AbstractTensor{N}, w::AbstractTensor{N}, h::AbstractTensor{N}) where {N}
    input_size = (size(w) .+ size(h) .- 1)[(begin + 1):(end - 1)]
    @assert size(I) == (size(w,1), input_size..., size(h)[end])
    @assert size(w)[end] == size(h, 1) # number of hidden units
end

function conv_h2v!(I::AbstractTensor{3}, w::AbstractTensor{3}, h::AbstractTensor{3}) where {N}
    _conv_h2v_check_size(I, w, h)
    for k ∈ axes(h, 2), μ ∈ axes(h, 1), j ∈ axes(w, 2), c in axes(w, 1)
        i = j + k - 1
        I[c,i,b] = w[c,j,μ] * h[μ,k,b]
    end
end

function conv_h2v!(I::AbstractTensor{4}, w::AbstractTensor{4}, h::AbstractTensor{4}) where {N}
    _conv_h2v_check_size(I, w, h)
    for b ∈ , k ∈ axes(h, 2), j ∈ axes(w, 2)
        i = j + k - 1
        for μ ∈ axes(h, 1), c in axes(w, 1)
            I[c,i,b] += w[c,j,μ] * h[μ,k,b]
        end
    end
end

function _conv_h2v!(I::AbstractTensor{3}, w::AbstractTensor{3}, h::AbstractTensor{3})
    _conv_h2v_check_size(I, w, h)
    @tullio I[c,j+k-1,b] = w[c,j,μ] * h[μ,k,b]
end

function _conv_h2v!(I::AbstractTensor{4}, w::AbstractTensor{4}, h::AbstractTensor{4})
    _conv_h2v_check_size(I, w, h)
    @tullio I[c,j₁+k₁-1,j₂+k₂-1,b] = w[c,j₁,j₂,μ] * h[μ,k₁,k₂,b]
end

function _conv_h2v!(I::AbstractTensor{5}, w::AbstractTensor{5}, h::AbstractTensor{5})
    _conv_h2v_check_size(I, w, h)
    @tullio I[c,j₁+k₁-1,j₂+k₂-1,j₃+k₃-1,b] = w[c,j₁,j₂,j₃,μ] * h[μ,k₁,k₂,k₃,b]
end

function _conv_h2v!(I::AbstractTensor{6}, w::AbstractTensor{6}, h::AbstractTensor{6})
    _conv_h2v_check_size(I, w, h)
    @tullio I[c,j₁+k₁-1,j₂+k₂-1,j₃+k₃-1,j₄+k₄-1,b] = w[c,j₁,j₂,j₃,j₄,μ] * h[μ,k₁,k₂,k₃,k₄,b]
end


text:
module BoundaryConditions

export
    BCType, Flux, Gradient, Value, Open,
    BoundaryCondition, getbc, setbc!,
    PeriodicBoundaryCondition, OpenBoundaryCondition, NoFluxBoundaryCondition,
    FluxBoundaryCondition, ValueBoundaryCondition, GradientBoundaryCondition,
    FieldBoundaryConditions,
    apply_x_bcs!, apply_y_bcs!, apply_z_bcs!,
    fill_halo_regions!

using CUDA
using KernelAbstractions

using Oceananigans.Architectures: device
using Oceananigans.Utils: work_layout, launch!
using Oceananigans.Operators: Ax, Ay, Az, volume
using Oceananigans.Grids

include("boundary_condition_classifications.jl")
include("boundary_condition.jl")
include("discrete_boundary_function.jl")
include("continuous_boundary_function.jl")
include("field_boundary_conditions.jl")
include("show_boundary_conditions.jl")

include("fill_halo_regions.jl")
include("fill_halo_regions_value_gradient.jl")
include("fill_halo_regions_open.jl")
include("fill_halo_regions_periodic.jl")
include("fill_halo_regions_flux.jl")
include("fill_halo_regions_nothing.jl")

include("apply_flux_bcs.jl")

end


text:
## distributions.jl (c) 2015 David A. van Leeuwen

## conversion to GMM
function GMM(m::MixtureModel{Multivariate,Continuous,FullNormal})
    Σ = eltype(FullCov{eltype(m)})[cholinv(c.Σ.mat) for c in components(m)]
    μ = hcat([c.μ for c in components(m)]...)'
    w = probs(m)
    n, d = size(μ)
    h = [History(@sprintf("Initialization from MixtureModel n=%d, d=%d, kind=full", n, d))]
    GMM(w, μ, Σ, h, 0)
end

function GMM(m::MixtureModel{Multivariate,Continuous,DiagNormal})
    Σ = hcat([c.Σ.diag for c in components(m)]...)'
    μ = hcat([c.μ for c in components(m)]...)'
    w = probs(m)
    n, d = size(μ)
    h = [History(@sprintf("Initialization from MixtureModel n=%d, d=%d, kind=diag", n, d))]
    GMM(w, μ, Σ, h, 0)
end

## conversion to MixtureModel
function Distributions.MixtureModel(gmm::GMM{T}) where {T<:AbstractFloat}
    if gmm.d == 1
        mixtures = [Normal(gmm.μ[i,1], sqrt(gmm.Σ[i,1])) for i=1:gmm.n]
    elseif kind(gmm) == :full
        mixtures = [MvNormal(vec(gmm.μ[i,:]), covar(gmm.Σ[i])) for i=1:gmm.n]
    else
        mixtures = [MvNormal(vec(gmm.μ[i,:]), sqrt.(vec(gmm.Σ[i,:]))) for i=1:gmm.n]
    end
    MixtureModel(mixtures, gmm.w)
end


text:
# Wrapper for level-1 BLIS matrix-diagonal routines.
#

# Level1d API of common forms are put together.
#
macro blis_group_level1d_form1(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          diagoffa, BliDoff,
                          diaga,    BliDiag,
                          transa,   BliTrans,
                          m,        BliDim,
                          n,        BliDim,
                          a,        Ptr{xType}, rsa, BliInc, csa, BliInc,
                          b,        Ptr{xType}, rsb, BliInc, csb, BliInc)
    end
end

@blis_group_level1d_form1 addd
@blis_group_level1d_form1 copyd
@blis_group_level1d_form1 subd

macro blis_group_level1d_form2(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          diagoffa, BliDoff,
                          diaga,    BliDiag,
                          transa,   BliTrans,
                          m,        BliDim,
                          n,        BliDim,
                          α,        Ptr{xType},
                          a,        Ptr{xType}, rsa, BliInc, csa, BliInc,
                          b,        Ptr{xType}, rsb, BliInc, csb, BliInc)
    end
end

@blis_group_level1d_form2 axpyd
@blis_group_level1d_form2 scal2d

macro blis_group_level1d_form3(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          conjα,    BliConj,
                          diagoffa, BliDoff,
                          m,        BliDim,
                          n,        BliDim,
                          α,        Ptr{xType},
                          a,        Ptr{xType}, rsa, BliInc, csa, BliInc)
    end
end

@blis_group_level1d_form3 scald
@blis_group_level1d_form3 setd

# Forms with only 1 member are defined directly.
#
@blis_ccall_group(invertd,
                  Cvoid,
                  diagoffa, BliDoff,
                  m,        BliDim,
                  n,        BliDim,
                  a,        Ptr{xType}, rsa, BliInc, csa, BliInc)
@blis_ccall_group(setid,
                  Cvoid,
                  diagoffa, BliDoff,
                  m,        BliDim,
                  n,        BliDim,
                  α,        Ptr{rType},
                  a,        Ptr{xType}, rsa, BliInc, csa, BliInc)
@blis_ccall_group(shiftd, Cvoid,
                  diagoffa, BliDoff,
                  m,        BliDim,
                  n,        BliDim,
                  α,        Ptr{rType},
                  a,        Ptr{xType}, rsa, BliInc, csa, BliInc)
@blis_ccall_group(xpbyd, Cvoid,
                  diagoffa, BliDoff,
                  diaga,    BliDiag,
                  transa,   BliTrans,
                  m,        BliDim,
                  n,        BliDim,
                  a,        Ptr{xType}, rsa, BliInc, csa, BliInc,
                  β,        Ptr{xType},
                  b,        Ptr{xType}, rsb, BliInc, csb, BliInc)


text:
uniform(u::Real) = 1//2 * Float64(abs(u) <= 1.0)
uniform_unnormalized(u::Real) = 1//2 * Float64(abs(u) <= 1.0)

triangular(u::Real) = (1 - abs(u)) * Float64(abs(u) <= 1.0)
triangular_unnormalized(u::Real) = (1 - abs(u)) * Float64(abs(u) <= 1.0)

epanechnikov(u::Real) = 3//4 * (1 - u^2) * Float64(abs(u) <= 1.0)
epanechnikov_unnormalized(u::Real) = (1 - u^2) * Float64(abs(u) <= 1.0)

biweight(u::Real) = 15//16 * (1 - u^2)^2 * Float64(abs(u) <= 1.0)
biweight_unnormalized(u::Real) = (1 - u^2)^2 * Float64(abs(u) <= 1.0)

triweight(u::Real) = 35//32 * (1 - u^2)^3 * Float64(abs(u) <= 1.0)
triweight_unnormalized(u::Real) = (1 - u^2)^3 * Float64(abs(u) <= 1.0)

tricube(u::Real) = 70//81 * (1 - abs(u)^3)^3 * Float64(abs(u) <= 1.0)
tricube_unnormalized(u::Real) = (1 - abs(u)^3)^3 * Float64(abs(u) <= 1.0)

gaussian(u::Real) = (1 / sqrt(2 * pi)) *  exp(-1//2 * u^2)
gaussian_unnormalized(u::Real) = exp(-1//2 * u^2)

cosine(u::Real) = (pi / 4) * cos((pi / 2) * u) * Float64(abs(u) <= 1.0)
cosine_unnormalized(u::Real) = cos((pi / 2) * u) * Float64(abs(u) <= 1.0)

logistic(u::Real) = 1 / (exp(u) + 2 + exp(-u))
logistic_unnormalized(u::Real) = 1 / (exp(u) + 2 + exp(-u))

kernels = Dict(
            :uniform => uniform,
            :triangular => triangular,
            :epanechnikov => epanechnikov,
            :biweight => biweight,
            :triweight => triweight,
            :tricube => tricube,
            :gaussian => gaussian,
            :cosine => cosine,
            :logistic => logistic
          )

unnormalized_kernels = Dict(
                        :uniform => uniform_unnormalized,
                        :triangular => triangular_unnormalized,
                        :epanechnikov => epanechnikov_unnormalized,
                        :biweight => biweight_unnormalized,
                        :triweight => triweight_unnormalized,
                        :tricube => tricube_unnormalized,
                        :gaussian => gaussian_unnormalized,
                        :cosine => cosine_unnormalized,
                        :logistic => logistic_unnormalized
                      )