text:
# AUTO GENERATED FILE - DO NOT EDIT

export daterangepicker

"""
    daterangepicker(;kwargs...)

A DateRangePicker component.
Capture date input from user. For more information, see: https://mantine.dev/dates/date-range-picker/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `allowLevelChange` (Bool; optional): Allow to change level (date – month – year)
- `allowSingleDateInRange` (Bool; optional): Allow one date to be selected as range
- `amountOfMonths` (Real; optional): Amount of displayed months
- `className` (String; optional): Often used with CSS to style elements with common properties
- `clearable` (Bool; optional): Allow to clear value
- `closeCalendarOnChange` (Bool; optional): Set to false to force dropdown to stay open after date was selected
- `closeDropdownOnScroll` (Bool; optional): Set to true to disable dropdown closing on scroll
- `dates` (Array of Strings; optional): Selected date
- `description` (String; optional): Input description, displayed after label
- `disableOutsideEvents` (Bool; optional): When true dates that are outside of given month cannot be clicked or focused
- `disabled` (Bool; optional): A Datepicker can show it is currently unable to be interacted with
- `dropdownType` (a value equal to: "modal", "popover"; optional): Where to show calendar in modal or popover
- `error` (String; optional): Displays error message after input
- `firstDayOfWeek` (a value equal to: "sunday", "monday"; optional): Set first day of the week
- `fixOnBlur` (Bool; optional): call onChange with last valid value onBlur
- `focusable` (Bool; optional): Should focusable days have tabIndex={0}?
- `format` (String; optional): DatePicker display format
- `fullWidth` (Bool; optional): Set to true to make calendar take 100% of container width
- `hideOutsideDates` (Bool; optional): Remove outside dates
- `hideWeekdays` (Bool; optional): Set to false to remove weekdays row
- `initialLevel` (a value equal to: "date", "month", "year"; optional): Initial date selection level
- `initialMonth` (String; optional): Initial selected month
- `initiallyOpened` (Bool; optional): Control initial dropdown opened state
- `label` (String; optional): Input label, displayed before input
- `labelSeparator` (String; optional): Separator between dates
- `locale` (String; optional): Locale used for all labels formatting
- `maxDate` (String; optional): Maximum possible date
- `minDate` (String; optional): Minimum possible date
- `multiline` (Bool; optional): Will input have multiple lines?
- `placeholder` (String; optional): Placeholder, displayed when date is not selected
- `preventFocus` (Bool; optional): Prevent focusing upon clicking
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Input border-radius from theme or number to set border-radius in px
- `required` (Bool; optional): Adds red asterisk on the right side of label
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Input size
- `style` (Dict; optional): Inline style override
- `variant` (a value equal to: "default", "filled", "unstyled", "headless"; optional): Defines input appearance, defaults to default in light color scheme and filled in dark
- `withinPortal` (Bool; optional): Whether to render the dropdown in a Portal
- `zIndex` (Real; optional): Popper zIndex
"""
function daterangepicker(; kwargs...)
        available_props = Symbol[:id, :allowLevelChange, :allowSingleDateInRange, :amountOfMonths, :className, :clearable, :closeCalendarOnChange, :closeDropdownOnScroll, :dates, :description, :disableOutsideEvents, :disabled, :dropdownType, :error, :firstDayOfWeek, :fixOnBlur, :focusable, :format, :fullWidth, :hideOutsideDates, :hideWeekdays, :initialLevel, :initialMonth, :initiallyOpened, :label, :labelSeparator, :locale, :maxDate, :minDate, :multiline, :placeholder, :preventFocus, :radius, :required, :size, :style, :variant, :withinPortal, :zIndex]
        wild_props = Symbol[]
        return Component("daterangepicker", "DateRangePicker", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
module AsapSim

# Use this to provide default values to a bunch of fields of structs.
using Parameters
using MacroTools
using LightDES
import AutoAligns

import Base.trunc_int

export  @asap4asm,
        # AsapOpcode,
        # SrcDest,
        AsapCore,
        AsapInstruction,
        AsapInstructionKeyword,
        InstructionLabelTarget,
        Loc,
        assemble,
        # Accessors for Core types.
        dmem,
        update!,
        summarize,
        showprogram,
        # Fifo Ops
        iswriteready,
        writeupdate!,
        isreadready,
        readupdate!

include("assembler/assembler.jl")

include("cores/fifos/fifo.jl")
include("cores/fifos/dualclockfifo.jl")
include("cores/fifos/testfifo.jl")

include("cores/core/core.jl")
include("cores/core/pipeline.jl")
include("cores/core/show.jl")

include("cores/io/io.jl")

include("sim/sim.jl")

end # module

ext: jl


text:
"""
    statsrank(x::Matrix{T}, ranktol::Real=1e-8) where {T<:AbstractFloat}

Return the numerical column rank and a pivot vector.

The rank is determined from the absolute values of the diagonal of R from
a pivoted QR decomposition, relative to the first (and, hence, largest)
element of this vector.

In the full-rank case the pivot vector is `collect(axes(x, 2))`.
"""
function statsrank(x::AbstractMatrix{T}; ranktol=1e-8) where {T<:AbstractFloat}
    m, n = size(x)
    piv = collect(axes(x, 2))

    iszero(n) && return (rank=n, piv=piv)

    qrpiv = pivoted_qr(x)
    dvec = abs.(diag(qrpiv.R))
    fdv = first(dvec)
    cmp = fdv * ranktol
    (last(dvec) > cmp) && return (rank=n, piv=piv)

    rank = searchsortedlast(dvec, cmp; rev=true)
    @assert rank < n
    piv = qrpiv.p
    v1 = first(eachcol(x))
    if all(isone, v1) && first(piv) ≠ 1
        # make sure the first column isn't moved by inflating v1
        v1 .*= (fdv + one(fdv)) / sqrt(m)
        qrpiv = pivoted_qr(x)
        piv = qrpiv.p
        fill!(v1, one(T))    # restore the contents of the first column
    end

    # maintain original column order for the linearly independent columns
    sort!(view(piv, 1:rank))
    return (rank=rank, piv=piv)
end


text:
# This is an internal type, do not export!
struct SpglibDataset
    spacegroup_number::Cint
    hall_number::Cint
    international_symbol::NTuple{11,Cchar}
    hall_symbol::NTuple{17,Cchar}
    choice::NTuple{6,Cchar}
    transformation_matrix::NTuple{9,Cdouble}
    origin_shift::NTuple{3,Cdouble}
    n_operations::Cint
    rotations::Ptr{NTuple{9,Cint}}
    translations::Ptr{NTuple{3,Cdouble}}
    n_atoms::Cint
    wyckoffs::Ptr{Cint}
    site_symmetry_symbols::Ptr{NTuple{7,Cchar}}
    equivalent_atoms::Ptr{Cint}
    crystallographic_orbits::Ptr{Cint}  # Added in v1.15.0
    primitive_lattice::NTuple{9,Cdouble}  # Added in v1.15.0
    mapping_to_primitive::Ptr{Cint}
    n_std_atoms::Cint
    std_lattice::NTuple{9,Cdouble}
    std_types::Ptr{Cint}
    std_positions::Ptr{NTuple{3,Cdouble}}
    std_rotation_matrix::NTuple{9,Cdouble}
    std_mapping_to_primitive::Ptr{Cint}
    pointgroup_symbol::NTuple{6,Cchar}
end

"""
    Dataset(spacegroup_number, hall_number, international_symbol, hall_symbol, choice, transformation_matrix, origin_shift, n_operations, rotations, translations, n_atoms, wyckoffs, site_symmetry_symbols, equivalent_atoms, crystallographic_orbits, primitive_lattice, mapping_to_primitive, n_std_atoms, std_lattice, std_types, std_positions, std_rotation_matrix, std_mapping_to_primitive, pointgroup_symbol)

Represent `SpglibDataset`, see its [official documentation](https://spglib.github.io/spglib/dataset.html#spglib-dataset).

!!! note
    Fields `crystallographic_orbits` and `primitive_lattice` are added after `spglib` `v1.15.0`.
"""
struct Dataset
    spacegroup_number::Int
    hall_number::Int
    international_symbol::String
    hall_symbol::String
    choice::String
    transformation_matrix::Matrix{Float64}
    origin_shift::Vector{Float64}
    n_operations::Int
    rotations::Array{Float64,3}
    translations::Matrix{Float64}
    n_atoms::Int
    wyckoffs::Vector{Char}
    site_symmetry_symbols::Vector{String}
    equivalent_atoms::Vector{Int}
    crystallographic_orbits::Vector{Int}
    primitive_lattice::Matrix{Float64}
    mapping_to_primitive::Vector{Int}
    n_std_atoms::Int
    std_lattice::Matrix{Float64}
    std_types::Vector{Int}
    std_positions::Matrix{Float64}
    std_rotation_matrix::Matrix{Float64}
    std_mapping_to_primitive::Vector{Int}
    pointgroup_symbol::String
end

function Base.convert(::Type{Dataset}, dataset::SpglibDataset)
    r = unsafe_wrap(Vector{NTuple{9,Cint}}, dataset.rotations, dataset.n_operations)
    t = unsafe_wrap(Vector{NTuple{3,Float64}}, dataset.translations, dataset.n_operations)
    wyckoffs = unsafe_wrap(Vector{Cint}, dataset.wyckoffs, dataset.n_atoms)
    pos = unsafe_wrap(Vector{NTuple{3,Float64}}, dataset.std_positions, dataset.n_std_atoms)
    return Dataset(
        dataset.spacegroup_number,
        dataset.hall_number,
        cchars2string(dataset.international_symbol),
        cchars2string(dataset.hall_symbol),
        cchars2string(dataset.choice),
        tuple2matrix(dataset.transformation_matrix),
        collect(dataset.origin_shift),
        dataset.n_operations,
        rotsFromTuple(r, dataset.n_operations),
        transFromTuple(t, dataset.n_operations),
        dataset.n_atoms,
        [LETTERS[x+1] for x in wyckoffs],  # Need to add 1 because of C-index starts from 0
        map(
            cchars2string,
            unsafe_wrap(
                Vector{NTuple{7,Cchar}},
                dataset.site_symmetry_symbols,
                dataset.n_atoms,
            ),
        ),
        unsafe_wrap(Vector{Cint}, dataset.equivalent_atoms, dataset.n_atoms),
        unsafe_wrap(Vector{Cint}, dataset.crystallographic_orbits, dataset.n_atoms),
        transpose(tuple2matrix(dataset.primitive_lattice)),
        unsafe_wrap(Vector{Cint}, dataset.mapping_to_primitive, dataset.n_atoms),
        dataset.n_std_atoms,
        transpose(tuple2matrix(dataset.std_lattice)),
        unsafe_wrap(Vector{Cint}, dataset.std_types, dataset.n_std_atoms),
        transFromTuple(pos, dataset.n_std_atoms),
        tuple2matrix(dataset.std_rotation_matrix),
        unsafe_wrap(Vector{Cint}, dataset.std_mapping_to_primitive, dataset.n_std_atoms),
        cchars2string(dataset.pointgroup_symbol),
    )
end


text:
module HMMs

using Models
using Distributions
using Random
using Combinatorics
using LinearAlgebra
using Printf
using SpecialFunctions
using Random
using Dependencies

include("GenTDist.jl")
include("NormalCubature.jl")

mutable struct DirtyArray{A <: AbstractArray}
    data::A
    dirty::Bool
end

DirtyArray( a::A ) where {A} = DirtyArray{A}( a, true )
DirtyArray{A}() where {T,N,A <: AbstractArray{T,N}} =
    DirtyArray( convert( A, zeros( T, fill(0,N)... ) ), true )

Base.copy( da::DirtyArray ) = DirtyArray( copy( da.data ), da.dirty )

const DirtyVector{T} = DirtyArray{Vector{T}} where {T}
const DirtyMatrix{T} = DirtyArray{Matrix{T}} where {T}

mutable struct HMM{N, Dist <: Distribution, Calc <: Real, Out <: Real, T} <: Models.AbstractModel{T, Out}
    initialprobabilities::Vector{Out}
    transitionprobabilities::Matrix{Out}
    stateparameters::Matrix{Out}

    currentprobabilities::Vector{Out}

    b::DirtyMatrix{Calc}
    dlogb::DirtyArray{Array{Calc,3}}
    d2logb::DirtyArray{Array{Calc,4}}
    d2b::DirtyArray{Array{Calc,4}}

    alpha::DirtyMatrix{Calc}
    dalpha::DirtyArray{Array{Calc,3}}
    d2alpha::DirtyArray{Array{Calc,4}}

    beta::DirtyMatrix{Calc}

    xi::DirtyArray{Array{Calc,3}}

    gamma::DirtyMatrix{Calc}

    likelihood::DirtyVector{Calc}
    dlikelihood::DirtyMatrix{Calc}
    d2likelihood::DirtyArray{Array{Calc,3}}
    d2loglikelihood::DirtyArray{Array{Calc,3}}

    y::Vector{Out}

    constraintmatrix::Matrix{Out}
    constraintvector::Vector{Out}
    basis::Matrix{Out}

    inequalityconstraintmatrix::Matrix{Out}
    inequalityconstraintvector::Vector{Out}

    statecdfs::Matrix{Out}
    statedists::Vector{Distribution}

    scratch::Dict{Symbol,Any}
end

const GaussianHMM{N, Calc, Out, T} = HMM{N, Normal, Calc, Out, T}
const LaplaceHMM{N, Calc, Out, T} = HMM{N, Laplace, Calc, Out, T}

inequalityconstraintmatrix( ::Type{GenTDist} ) = [0.0 1.0 0.0; 0.0 0.0 1.0]
inequalityconstraintvector( ::Type{GenTDist} ) = [1e-4, 2.0]

function equalityconstraintmatrix( m::Int, p::Int )
    A = [0.0 + (m==1 || div(j-1,m)==i-1) for i in 1:m, j in 1:m^2]
    return [A zeros(m,m*p)]
end
equalityconstraintvector( m::Int ) = ones( m )

function HMM{N,Dist,Calc,Out,T}(
    pi::Vector{Out},
    a::Matrix{Out},
    stateparameters::Matrix{Out};
    scratch::Dict{Symbol,Any} = Dict{Symbol,Any}(),
) where {N,Dist,Calc,Out,T}
    (p,m) = size(stateparameters)

    A = equalityconstraintmatrix( m, p )

    C = [Matrix{Out}(Diagonal(ones(m^2))) zeros(m^2,m*p)]
    d = zeros(m^2)
    M = inequalityconstraintmatrix( Dist )
    v = inequalityconstraintvector( Dist )
    for i = 1:m
        for j = 1:length(v)
            constraint = zeros(m^2)
            for k = 1:p
                constraint = [constraint; zeros(i-1); M[j,k]; zeros(m-i)]
            end
            C = vcat( C, constraint' )
            push!( d, v[j] )
        end
    end

    result = HMM{N,Dist,Calc,Out,T}(
        pi, a, stateparameters,

        copy(pi),

        DirtyMatrix{Calc}(),
        DirtyArray{Array{Calc,3}}(),
        DirtyArray{Array{Calc,4}}(),
        DirtyArray{Array{Calc,4}}(),

        DirtyMatrix{Calc}(),
        DirtyArray{Array{Calc,3}}(),
        DirtyArray{Array{Calc,4}}(),

        DirtyMatrix{Calc}(),

        DirtyArray{Array{Calc,3}}(),

        DirtyMatrix{Calc}(),

        DirtyVector{Calc}(),
        DirtyMatrix{Calc}(),
        DirtyArray{Array{Calc,3}}(),
        DirtyArray{Array{Calc,3}}(),

        Vector{Out}(),

        A,
        equalityconstraintvector( m ),
        basis(A),

        C,
        d,

        zeros(0,0),
        Distribution[],

        scratch,
    )
end

Base.deepcopy( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T} =
    HMM{N,Dist,Calc,Out,T}(
        copy( hmm.initialprobabilities ),
        copy( hmm.transitionprobabilities ),
        copy( hmm.stateparameters ),

        copy( hmm.currentprobabilities ),

        copy( hmm.b ),
        copy( hmm.dlogb ),
        copy( hmm.d2logb ),
        copy( hmm.d2b ),

        copy( hmm.alpha ),
        copy( hmm.dalpha ),
        copy( hmm.d2alpha ),

        copy( hmm.beta ),
        copy( hmm.xi ),
        copy( hmm.gamma ),

        copy( hmm.likelihood ),
        copy( hmm.dlikelihood ),
        copy( hmm.d2likelihood ),
        copy( hmm.d2loglikelihood ),

        copy( hmm.y ),

        copy( hmm.constraintmatrix ),
        copy( hmm.constraintvector ),
        copy( hmm.basis ),

        copy( hmm.inequalityconstraintmatrix ),
        copy( hmm.inequalityconstraintvector ),

        copy( hmm.statecdfs ),
        copy( hmm.statedists ),

        copy( hmm.scratch ),
    )

randomparameters( ::Union{Type{Normal},Type{Laplace}}, numstates::Int ) =
    [randn( 1, numstates ); randn( 1, numstates ).^2]

function Models.rand(
    ::Type{HMM{N,D,Calc,Out,T}};
    seed::Int = -1,
) where {N, D <: Distribution, Calc <: Real, Out<:Real, T}
    if seed >= 0
        Random.seed!( seed )
    end
    initialprobabilities = Vector{Out}(rand( N ))
    initialprobabilities ./= sum( initialprobabilities )

    transitionprobabilities = zeros( Out, N, N )
    for i = 1:N
        for j = 1:N
            transitionprobabilities[i,j] = rand()
        end
    end
    transitionprobabilities ./= sum( transitionprobabilities, dims=2 )

    parameters = Matrix{Out}( randomparameters( D, N ) )
    scratch = Dict{Symbol,Any}()

    return HMM{N,D,Calc,Out,T}( initialprobabilities, transitionprobabilities, parameters, scratch=scratch )
end

function Distributions.rand!(
    hmm::HMM{N,D,Calc,Out,T};
    seed::Int = -1,
) where {N, D <: Distribution, Calc <: Real, Out<:Real, T}
    if seed >= 0
        Random.seed!( seed )
    end
    hmm.initialprobabilities[:] = Vector{Out}(rand( N ))
    hmm.initialprobabilities[:] ./= sum( initialprobabilities )

    hmm.transitionprobabilities[:] = zeros( Out, N, N )
    for i = 1:N
        for j = 1:N
            hmm.transitionprobabilities[i,j] = rand()
        end
    end
    transitionprobabilities[:,:] ./= sum( transitionprobabilities, dims=2 )

    hmm.stateparameters[:,:] = randomparameters( D, N )
    return hmm
end

function Distributions.rand!(
    hmm::HMM{N,Dist,Calc,Out,T},
    observations::AbstractVector{Out},
    n::Int = length(observations);
    seed::Int = 1,
) where {N,Dist,Calc,Out,T}
    isempty(observations) && return nothing

    state = searchsorted( cumsum( hmm.currentprobabilities ), rand() ).start
    if isempty( hmm.statedists )
        hmm.statedists = [Dist(hmm.stateparameters[:,state]...) for state in 1:length(hmm.initialprobabilities)]
    end
    observations[1] = rand( hmm.statedists[state] )

    if isempty(hmm.statecdfs)
        hmm.statecdfs = hcat( [ cumsum( hmm.transitionprobabilities[i,:] ) for i in 1:length(hmm.initialprobabilities) ]... )'
    end

    for t = 2:n
        state = searchsorted( hmm.statecdfs[state,:], rand() ).start
        observations[t] = rand( hmm.statedists[state] )
    end
end

Distributions.rand!( hmm::HMM{N,Dist,Calc,Out,T}, t::AbstractVector{T}, u::AbstractVector{Out};
                     kwargs... ) where {N,Dist,Calc,Out,T} =
                         rand!( hmm, u; kwargs... )

function clear( hmm::HMM )
    hmm.b.dirty = true
    hmm.dlogb.dirty = true
    hmm.d2logb.dirty = true
    hmm.d2b.dirty = true

    hmm.alpha.dirty = true
    hmm.dalpha.dirty = true
    hmm.d2alpha.dirty = true

    hmm.beta.dirty = true
    hmm.xi.dirty = true
    hmm.gamma.dirty = true

    hmm.likelihood.dirty = true
    hmm.dlikelihood.dirty = true
    hmm.d2likelihood.dirty = true
    hmm.d2loglikelihood.dirty = true
end

function free( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    for field in [:b, :dlogb, :d2logb, :d2b, :alpha, :dalpha, :d2alpha, :beta, :xi, :gamma,
                  :likelihood, :dlikelihood, :d2likelihood, :d2loglikelihood]
        value = getfield(hmm, field)
        if !isempty(value.data)
            setfield!( hmm, field, typeof(getfield(hmm,field))() )
        end
    end
end

function writearray( io::IO, v::Array{T,N} ) where {T,N}
    for i in size(v)
        Base.write( io, i )
    end
    for x in v
        Base.write( io, x )
    end
end

function readarray( io::IO, ::Type{Array{T,N}} ) where {T,N}
    size = Int[]
    for i in 1:N
        push!( size, Base.read( io, Int ) )
    end
    v = zeros( T, size... )
    for i = 1:prod(size)
        v[i] = Base.read( io, T )
    end
    return v
end

function Base.write( io::IO, hmm::HMM )
    writearray( io, hmm.initialprobabilities )
    writearray( io, hmm.transitionprobabilities )
    writearray( io, hmm.stateparameters )
    writearray( io, hmm.y )
end

function Base.read( io::IO, ::Type{HMM{N,Dist,Calc,Out}} ) where {N,Dist,Calc,Out}
    initialprobabilities = readarray( io, Vector{Out} )
    transitionprobabilities = readarray( io, Matrix{Out} )
    parameters = readarray( io, Matrix{Out} )
    y = readarray( io, Vector{Out} )
    hmm = HMM{N,Dist,Calc,Out}( initialprobabilities, transitionprobabilities, parameters )
    setobservations( hmm, y )
    return hmm
end

function getparameters!( hmm::HMM{N,Dist,Calc,Out,T}, parameters::Vector{Out} ) where {N,Dist,Calc,Out,T}
    parameters[:] = [hmm.transitionprobabilities'[:]; hmm.stateparameters'[:]]
end

getparameters( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T} =
    [hmm.transitionprobabilities'[:]; ; hmm.stateparameters'[:]]

function setparameters!( hmm::HMM{N,Dist,Calc,Out,T}, parameters::AbstractVector{Out} ) where {N,Dist,Calc,Out,T}
    oldparameters = getparameters( hmm )
    delta = parameters - oldparameters

    C = hmm.inequalityconstraintmatrix
    d = hmm.inequalityconstraintvector
    active = C * parameters .< d
    alpha = sum(active) == 0 ? 1.0 : minimum((d[active] - C[active,:]*oldparameters) ./ (C[active,:] * delta))
    newparameters = oldparameters + alpha * delta

    (p, m) = size(hmm.stateparameters)

    hmm.transitionprobabilities'[:] = newparameters[1:m*m]
    sums = sum(hmm.transitionprobabilities, dims=2)
    hmm.transitionprobabilities ./= sums
    hmm.stateparameters'[:] = newparameters[m*m+1:m*(m+p)]
    clear( hmm )
end

function setobservations( hmm::HMM{N,Dist,Calc,T}, y::Vector{U} ) where {N,Dist,Calc,T, U <: Real}
    m = length(hmm.initialprobabilities)

    free( hmm )
    hmm.y = y
end

observations( hmm::HMM ) = hmm.y

function probabilities( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    # makes the assumption that the first parameter is location and second is scale
    if hmm.b.dirty
        y = observations( hmm )
        t = length(y)
        m = length(hmm.initialprobabilities)
        if isempty(hmm.b.data)
            hmm.b.data = zeros( Calc, (t,m) )
        end

        for i in 1:length(hmm.initialprobabilities)
            if isnan(hmm.stateparameters[1,i]) || isnan(hmm.stateparameters[2,i])
                for j = 1:size(hmm.b.data[:,1],1)
                    hmm.b.data[j,i] = NaN
                end
            elseif hmm.stateparameters[2,i] == 0.0
                for j = 1:size(hmm.b.data[:,1],1)
                    hmm.b.data[j,i] = y == hmm.stateparameters[1,i] ? Inf : 0
                end
            else
                dist = Dist( hmm.stateparameters[:,i]... )
                hmm.b.data[:,i] = pdf.( dist, y )
            end
        end

        hmm.b.dirty = false
    end
    return hmm.b.data
end

function dlogprobabilities( hmm::GaussianHMM{N,Calc,Out,T} ) where {N,Calc,Out,T}
    if hmm.dlogb.dirty
        b = probabilities( hmm )
        (t,m) = size(b)
        if isempty(hmm.dlogb.data)
            hmm.dlogb.data = zeros( Calc, m*(m+2), m, t )
        end
        y = observations( hmm )
        for i = 1:m
            for t = 1:t
                z = (y[t] - hmm.stateparameters[1,i])/hmm.stateparameters[2,i]
                hmm.dlogb.data[m^2 + i,i,t] = z / hmm.stateparameters[2,i]
                hmm.dlogb.data[m*(m+1) + i,i,t] = (z^2 - 1)/hmm.stateparameters[2,i]
            end
        end
        hmm.dlogb.dirty = false
    end
    return hmm.dlogb.data
end

function dlogprobabilities( hmm::HMM{N,GenTDist,Calc,Out,T} ) where {N,Calc,Out,T}
    if hmm.dlogb.dirty
        b = probabilities( hmm )
        (t,m) = size(b)
        if isempty(hmm.dlogb.data)
            hmm.dlogb.data = zeros( Calc, m*(m+3), m, t )
        end

        y = observations( hmm )

        for i = 1:m
            (mu,sigma,nu) = hmm.stateparameters[:,i]
            for t = 1:t
                centeredy = y[t] - mu
                normalysq = (centeredy/sigma)^2
                hmm.dlogb.data[m^2 + i, i, t] = (nu+1)*centeredy/(nu*sigma^2 + centeredy^2)
                hmm.dlogb.data[m*(m+1) + i, i, t] = -1/sigma + (nu+1)* normalysq / (sigma*(nu + normalysq))

                hmm.dlogb.data[m*(m+2) + i, i, t] = digamma((nu+1)/2)/2 - 1/(2*nu) - digamma(nu/2)/2
                hmm.dlogb.data[m*(m+2) + i, i, t] += -log(1 + normalysq/nu)/2 + (nu+1)/(2*nu) * normalysq/(nu + normalysq)
            end
        end
        hmm.dlogb.dirty = false
    end
    return hmm.dlogb.data
end

function d2logprobabilities( hmm::GaussianHMM{N,Calc,Out,T} ) where {N,Calc,Out,T}
    if hmm.d2logb.dirty
        dlogb = dlogprobabilities( hmm )
        (p,m,t) = size(dlogb)
        if isempty(hmm.d2logb.data)
            hmm.d2logb.data = zeros( Calc, p, p, m, t )
        end
        y = observations( hmm )
        for i = 1:m
            for t = 1:t
                sigma = hmm.stateparameters[2,i]
                z = (y[t] - hmm.stateparameters[1,i])/sigma

                mui = m^2+i
                hmm.d2logb.data[mui,mui,i,t] = -1/sigma^2

                sigmai = m*(m+1)+i
                dmusigma = -2*z/sigma^2
                hmm.d2logb.data[mui,sigmai,i,t] = dmusigma
                hmm.d2logb.data[sigmai,mui,i,t] = dmusigma

                hmm.d2logb.data[sigmai,sigmai,i,t] = (1 - 3z^2)/sigma^2
            end
        end
        hmm.d2logb.dirty = false
    end
    return hmm.d2logb.data
end

function d2logprobabilities( hmm::HMM{N,GenTDist,Calc,T} ) where {N,Calc,T}
    if hmm.d2logb.dirty
        dlogb = dlogprobabilities( hmm )
        (p,m,t) = size(dlogb)
        if isempty(hmm.d2logb.data)
            hmm.d2logb.data = zeros( Calc, p, p, m, t )
        end
        y = observations( hmm )
        for i = 1:m
            (mu,sigma,nu) = hmm.stateparameters[:,i]
            centeredy = y .- mu
            centeredysq = centeredy .^ 2
            normalysq = (centeredy./sigma).^2
            denom = (nu * sigma^2 .+ centeredysq).^2

            mui = m^2 + i
            sigmai = m*(m+1) + i
            nui = m*(m+2) + i
            hmm.d2logb.data[mui,mui,i,:] = (nu+1) .* (centeredysq .- nu*sigma^2) ./ denom
            hmm.d2logb.data[mui,sigmai,i,:] = hmm.d2logb.data[sigmai,mui,i,:] = -2*(nu+1)*nu*sigma .* centeredy ./ denom
            hmm.d2logb.data[mui,nui,i,:] = hmm.d2logb.data[nui,mui,i,:] = (centeredy.^3 .- sigma^2 .* centeredy) ./ denom

            hmm.d2logb.data[sigmai,sigmai,i,:] =
                1/sigma^2 .- (nu+1) .* centeredysq .* (3*nu*sigma^2 .+ centeredysq) ./ (sigma^2 .* denom)
            hmm.d2logb.data[sigmai,nui,i,:] = hmm.d2logb.data[nui,sigmai,i,:] =
                sigma .* centeredysq .* (centeredysq .- sigma^2) ./ (sigma^2 .* denom)

            hmm.d2logb.data[nui,nui,i,:] =
                polygamma(1,(nu+1)/2)/4 + 1/(2*nu^2) - polygamma(1,nu/2)/4 .+
                ((nu-1) .* normalysq .^ 2 ./ 2 .- normalysq .* nu) ./ (nu .* (nu .+ normalysq)).^2
        end
    end
    return hmm.d2logb.data
end

function d2probabilities( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    if hmm.d2b.dirty
        b = probabilities( hmm )
        dlogb = dlogprobabilities( hmm )
        d2logb = d2logprobabilities( hmm )
        (p,m,t) = size(dlogb)
        if isempty(hmm.d2b.data)
            hmm.d2b.data = zeros( Calc, p, p, m, t )
        end
        for i = 1:m
            for t = 1:t
                hmm.d2b.data[:,:,i,t] = b[t,i] .* (d2logb[:,:,i,t] + dlogb[:,i,t] * dlogb[:,i,t]')
            end
        end

        hmm.d2b.dirty = false
    end
    return hmm.d2b.data
end

function stationary( hmm::HMM )
    P = Matrix{Float64}( hmm.transitionprobabilities )
    N = length(hmm.initialprobabilities)
    I = one(P)
    P -= I
    P[:,1] = ones(N)
    return I[1,:]'*pinv(P)
end

function variance( hmm )
    pi = stationary( hmm )
    mus = hmm.stateparameters[1,:]
    mu = dot( pi, mus )
    result = dot(pi, (hmm.stateparameters[2,:].^2 + (mus .- mu).^2))
    return result
end

function autocovariance( hmm, lag::Int )
    if lag == 0
        return variance( hmm )
    else
        pi = stationary( hmm )
        P = hmm.transitionprobabilities
        mus = hmm.stateparameters[1,:]
        mu = dot(pi, mus)
        jointprobabilities = pi' .* P^lag
        sum( jointprobabilities .* ((mus .- mu) * (mus .- mu)') )
    end
end

function forwardprobabilities( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    if hmm.alpha.dirty
        b = probabilities( hmm )
        (t,m) = size(b)
        if isempty( hmm.alpha.data )
            hmm.alpha.data = zeros( Calc, t, m )
        end

        hmm.alpha.data[1,:] = hmm.initialprobabilities .* b[1,:]
        hmm.alpha.data[2:t,:] = zeros( Calc, (t-1,m) )
        for i = 2:t
            for from = 1:N
                for to = 1:N
                    hmm.alpha.data[i,to] += hmm.transitionprobabilities[from, to] * hmm.alpha.data[i-1,from] * b[i,to]
                end
            end
        end
        hmm.alpha.dirty = false
    end
    return hmm.alpha.data
end

function dforwardprobabilities( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    if hmm.dalpha.dirty
        b = probabilities( hmm )
        dlogb = dlogprobabilities( hmm )
        alpha = forwardprobabilities( hmm )

        (p,m,t) = size(dlogb)

        if isempty( hmm.dalpha.data )
            hmm.dalpha.data = zeros( Calc, p, m, t )
        end

        hmm.dalpha.data[:,:,1] = hmm.initialprobabilities' .* (dlogb[:,:,1] .* b[1,:]')
        hmm.dalpha.data[:,:,2:t] = zeros( p, m, t-1 )
        for i = 2:t
            for from = 1:N
                for to = 1:N
                    paramindex = (from-1)*m + to
                    hmm.dalpha.data[paramindex,to,i] += alpha[i-1,from] * b[i,to]

                    hmm.dalpha.data[:,to,i] += hmm.transitionprobabilities[from, to] * hmm.dalpha.data[:,from,i-1] * b[i,to]

                    hmm.dalpha.data[:,to,i] +=
                        hmm.transitionprobabilities[from, to] * alpha[i-1,from] * (dlogb[:,to,i] .* b[i,to])
                end
            end
        end
        hmm.dalpha.dirty = false
    end
    return hmm.dalpha.data
end

function d2forwardprobabilities( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    if hmm.d2alpha.dirty
        b = probabilities( hmm )
        dlogb = dlogprobabilities( hmm )
        d2b = d2probabilities( hmm )

        alpha = forwardprobabilities( hmm )
        dalpha = dforwardprobabilities( hmm )

        (p,m,t) = size(dlogb)
        if isempty( hmm.d2alpha.data )
            hmm.d2alpha.data = zeros( Calc, p, p, m, t )
        end

        for i = 1:m
            hmm.d2alpha.data[:,:,i,1] = hmm.initialprobabilities[i] .* d2b[:,:,i,1]
        end
        hmm.d2alpha.data[:,:,:,2:t] = zeros( p, p, m, t-1 )
        for i = 2:t
            for from = 1:N
                for to = 1:N
                    paramindex = (from-1)*m + to

                    hmm.d2alpha.data[:,:,to,i] +=
                        hmm.transitionprobabilities[from,to] * hmm.d2alpha.data[:,:,from,i-1] * b[i,to]

                    hmm.d2alpha.data[paramindex,:,to,i] += dalpha[:,from,i-1] * b[i,to]
                    hmm.d2alpha.data[:,paramindex,to,i] += dalpha[:,from,i-1] * b[i,to]

                    hmm.d2alpha.data[paramindex,:,to,i] += alpha[i-1,from] * (dlogb[:,to,i] .* b[i,to])
                    hmm.d2alpha.data[:,paramindex,to,i] += alpha[i-1,from] * (dlogb[:,to,i] .* b[i,to])

                    # adding these together all at once preserves symmetry
                    hmm.d2alpha.data[:,:,to,i] +=
                        dalpha[:,from,i-1] * (dlogb[:,to,i]' .* b[i,to]) * hmm.transitionprobabilities[from,to] +
                        (dlogb[:,to,i] .* b[i,to]) * dalpha[:,from,i-1]' * hmm.transitionprobabilities[from,to]

                    hmm.d2alpha.data[:,:,to,i] +=
                        alpha[i-1,from] * hmm.transitionprobabilities[from,to] * d2b[:,:,to,i]
                end
            end
        end
        hmm.d2alpha.dirty = false
    end
    return hmm.d2alpha.data
end

function backwardprobabilities( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    if hmm.beta.dirty
        b = probabilities( hmm )
        (t,m) = size(b)
        if isempty( hmm.beta.data )
            hmm.beta.data = zeros( Calc, t, m )
        end

        hmm.beta.data[end,:] = ones(Calc,length(hmm.initialprobabilities))
        hmm.beta.data[1:t-1,:] = zeros( Calc, (t-1,m) )
        for i = t:-1:2
            for from = 1:N
                for to = 1:N
                    hmm.beta.data[i-1,from] += hmm.transitionprobabilities[from,to] * hmm.beta.data[i,to] * b[i,to]
                end
            end
        end
        hmm.beta.dirty = false
    end
    return hmm.beta.data
end

function likelihood( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    if hmm.likelihood.dirty
        alpha = forwardprobabilities( hmm )
        (t,m) = size(alpha)
        if isempty( hmm.likelihood.data )
            hmm.likelihood.data = zeros( Calc, t )
        end

        hmm.likelihood.data[:] = sum(alpha, dims=2)

        hmm.likelihood.dirty = false
    end
    return hmm.likelihood.data
end

function dlikelihood( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    if hmm.dlikelihood.dirty
        dalpha = dforwardprobabilities( hmm )
        (p,m,t) = size(dalpha)
        if isempty( hmm.dlikelihood.data )
            hmm.dlikelihood.data = zeros( Calc, p, t )
        end

        hmm.dlikelihood.data[:,:] = reshape(sum(dalpha[:,:,:],dims=2), (p,t))

        hmm.dlikelihood.dirty = false
    end
    return hmm.dlikelihood.data
end

function d2likelihood( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    if hmm.d2likelihood.dirty
        d2alpha = d2forwardprobabilities( hmm )
        (p,p,m,t) = size(d2alpha)
        if isempty( hmm.d2likelihood.data )
            hmm.d2likelihood.data = zeros( Calc, p, p, t )
        end

        hmm.d2likelihood.data[:,:,:] = reshape(sum(d2alpha[:,:,:,:],dims=3), (p,p,t))

        hmm.d2likelihood.dirty = false
    end
    return hmm.d2likelihood.data
end

function d2loglikelihood( hmm::HMM{N,Dist,Calc,T} ) where {N,Dist,Calc,T}
    if hmm.d2loglikelihood.dirty
        l = likelihood( hmm )
        dl = dlikelihood( hmm )
        d2l = d2likelihood( hmm )
        (p,t) = size(dl)

        if isempty( hmm.d2loglikelihood.data )
            hmm.d2loglikelihood.data = zeros(Calc,p,p,t)
        end

        for i = 1:t
            hmm.d2loglikelihood.data[:,:,i] = d2l[:,:,i] ./ l[i] - dl[:,i] * dl[:,i]' ./ l[i]^2
        end

        hmm.d2loglikelihood.dirty = false
    end
    return hmm.d2loglikelihood.data
end

function basis( A::Matrix{Out}; epsilon::Out = 1e-12 ) where {Out <: Number}
    (m, n) = size(A)
    (U,S,V) = svd( [A; zeros(n - m, n)] )
    indices = abs.(S) .< epsilon
    result = V[:,indices]
    return result
end

function addconstraints( hmm::HMM{N,Dist,Calc,Out,T}, A::AbstractMatrix{Out}, b::AbstractVector{Out} ) where {N,Dist,Calc,Out,T}
    hmm.constraintmatrix = [hmm.constraintmatrix; A]
    hmm.constraintvector = [hmm.constraintvector; b]
    hmm.basis = basis( hmm.constraintmatrix )
end

dcollapse( hmm::HMM ) = hmm.basis'

dexpand( hmm::HMM ) = hmm.basis

function Models.sandwich( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    # let's set transition probabilities
    (p,m) = size( hmm.stateparameters )
    hmm.constraintmatrix = equalityconstraintmatrix( m, p )
    hmm.constraintvector = equalityconstraintvector( m )
    parameters = getparameters( hmm )
    indices = findall(abs.(parameters[1:m^2]) .< 1e-8)
    if !isempty( indices )
        A = Matrix{Out}(hcat([[zeros(index-1); 1.0; zeros(length(parameters)-index)] for index in indices]...)')
        b = zeros(length(indices))
        addconstraints( hmm, A, b )
    end

    dc = dcollapse( hmm )

    l = likelihood( hmm )
    dl = dc * dlikelihood( hmm )

    (p,t) = size(dl)
    dlogl = [zeros(1,p); dl'./l]
    dlogln = diff( dlogl, dims=1 )
    V = dlogln' * dlogln

    J = dc * d2loglikelihood( hmm )[:,:,end] * dc'

    result = inv(J) * V * inv(J)
    return Matrix{Out}((result + result')/2)
end

function testle( hmm::HMM{N,Dist,Calc,Out,T}, A::AbstractMatrix{Out}, b::AbstractVector{Out} ) where {N,Dist,Calc,Out,T}
    dc = dcollapse( hmm )
    C = convert( Matrix{Out}, sandwich( hmm ) )
    Adc = A * dc'
    cov = Adc * C * Adc'
    p = getparameters( hmm )
    if length(b) == 1
        return cdf( Normal(0, sqrt(cov[1])), (A * p - b)[1] )
    else
        @assert( "Cubature not yet implemented" )
    end

end

function conditionaljointstateprobabilities( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    if hmm.xi.dirty
        alpha = forwardprobabilities( hmm )
        beta = backwardprobabilities( hmm )
        proby = likelihood( hmm )[end]
        b = probabilities( hmm )
        (t,m) = size(alpha)

        if isempty( hmm.xi.data )
            hmm.xi.data = zeros( Calc, t-1, m, m )
        end

        hmm.xi.data[1:t-1,:,:] = zeros( Calc, (t-1,m,m) )
        for i = 1:t-1
            for from = 1:N
                for to = 1:N
                    hmm.xi.data[i,from,to] +=
                        hmm.transitionprobabilities[from,to] * alpha[i,from] * beta[i+1,to] * b[i+1,to]
                end
            end
        end
        hmm.xi.data[:,:,:] /= proby

        hmm.xi.dirty = false
    end
    return hmm.xi.data
end

function conditionalstateprobabilities( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    if hmm.gamma.dirty
        xi = conditionaljointstateprobabilities( hmm )
        (Tm1, m, trash) = size(xi)
        if isempty( hmm.gamma.data )
            hmm.gamma.data = zeros( Calc, Tm1 + 1, m )
        end

        hmm.gamma.data[1:end-1,:] = sum(xi, dims=3)
        hmm.gamma.data[end,:] = sum(xi[end,:,:], dims=1)

        hmm.gamma.dirty = false
    end
    return hmm.gamma.data
end

function fit_mle!(
    ::Type{Normal},
    parameters::AbstractVector{Out},
    x::Vector{Out},
    w::Vector{Calc},
    scratch::Dict{Symbol,Any};
    kwargs...
) where {Calc,Out}
    mu = dot( w, x )
    sigma = sqrt( dot( w, (x .- mu).^2 ) )
    parameters[1] = mu
    parameters[2] = sigma
end

function fit_mle!(
    ::Type{Laplace},
    parameters::AbstractVector{Out},
    x::Vector{Out},
    w::Vector{Calc},
    scratch::Dict{Symbol,Any};
    kwargs...
) where {Calc,Out}
    if !haskey( scratch, :sortperm )
        scratch[:sortperm] = sortperm(x)
    end
    perm = scratch[:sortperm]

    cumw = cumsum( w[perm] )
    mu = x[perm[findfirst(cumw.>=0.5*cumw[end])]]
    sigma = dot( w, abs.(x .- mu) )
    parameters[1] = mu
    parameters[2] = sigma
end

function emstep(
    hmm::HMM{N,Dist,Calc,Out,T},
    nexthmm::HMM{N,Dist,Calc,Out,T};
    max_iter::Int = 0,
    print_level::Int = 0,
) where {N,Dist,Calc,Out,T}
    y = observations( hmm )
    t = length(y)

    gamma = conditionalstateprobabilities( hmm )
    occupation = sum(gamma[1:end-1,:],dims=1)

    xi = conditionaljointstateprobabilities( hmm )

    m = length(hmm.initialprobabilities)
    nexthmm.transitionprobabilities = reshape(sum(xi, dims=1), (m,m))./occupation'
    nexthmm.initialprobabilities = gamma[1,:]./sum(gamma[1,:])

    for i = 1:m
        nexthmm.stateparameters[:,i] = hmm.stateparameters[:,i]
        if max_iter > 0
            fit_mle!( Dist, view( nexthmm.stateparameters, :, i ), y, gamma[:,i]/occupation[i], hmm.scratch,
                      max_iter = max_iter, print_level=print_level )
        else
            fit_mle!( Dist, view( nexthmm.stateparameters, :, i ), y, gamma[:,i]/occupation[i], hmm.scratch,
                      print_level=print_level )
        end
    end

    clear( nexthmm )
end

function em(
    hmm::HMM{N,Dist,Calc,Out,T};
    epsilon::Float64 = 0.0,
    debug::Int = 0,
    maxiterations::Iter = Inf,
    keepintermediates = false,
    acceleration = Inf,
    accelerationlinestart = 500,
    accelerationmaxhalves = 5,
    max_iter::Int = 0,
    print_level::Int = 0,
    timefractiontowardszero = 0.5,
    observations = 10,
    finishbig = 0,
) where {N,Dist, Calc, Out, T, Iter <: Number}
    if acceleration < Inf
        @assert( keepintermediates )
    end

    t0 = Base.time()
    nexthmm = deepcopy( hmm )
    hmms = [hmm, nexthmm]
    oldlikelihood = zero(Calc)
    newlikelihood = likelihood( hmm )[end]
    done = false
    i = 1
    nextacceleration = acceleration
    m = length(hmm.initialprobabilities)

    if keepintermediates
        intermediates = [getparameters( hmm )]
    end

    iterations = 1
    while !done
        if debug >= 2
            println( "Iteration $iterations, likelihood = $newlikelihood" )
        end
        emstep( hmms[i], hmms[3-i], max_iter=max_iter, print_level=print_level )
        oldlikelihood = newlikelihood
        done = any(isnan.(hmms[3-i].initialprobabilities)) ||
            any(isnan.(hmms[3-i].transitionprobabilities)) ||
            any(isnan.(hmms[3-i].stateparameters)) || any(hmms[3-i].stateparameters[2,:].<=0) ||
            iterations >= maxiterations
        if !done
            newlikelihood = likelihood( hmms[3-i] )[end]
            done = newlikelihood / oldlikelihood - 1 <= epsilon
            if !done
                i = 3-i
            end
        else
            i = 3-i
        end
        iterations += 1
        if keepintermediates
            push!( intermediates, getparameters( hmms[i] ) )
        end
        if iterations >= nextacceleration
            y = hcat( intermediates[end-observations+1:end]... )'
            t = -collect(reverse(0:observations-1))
            X = [ones(observations) t t.^2]
            beta = inv(X'*X)*X'*y

            t = accelerationlinestart

            # let's bound t to where the probabilities go <0 or >1
            m = length(hmms[3-i].initialprobabilities)
            for j = 1:m^2
                (c,b,a) = beta[:,j]

                for target=0:1
                    discriminant = b^2 - 4 * a * (c - target) * timefractiontowardszero
                    if discriminant >= 0
                        srd = sqrt(discriminant)
                        roots = (-b .+ [-srd,srd])/(2*a)
                        roots = roots[roots .> 0]
                        if !isempty(roots)
                            t = min( t, minimum(roots) )
                        end
                    end
                end
            end

            newerlikelihood = 0.0
            optimalparameters = nothing
            optimalt = 0
            for j = 1:accelerationmaxhalves
                x = [1.0, t, t^2]
                p = x' * beta
                p = reshape( p, length(p) )
                setparameters!( hmms[3-i], p )

                hmms[3-i].transitionprobabilities ./= sum( hmms[3-i].transitionprobabilities, dims=2 )

                if any(hmms[3-i].stateparameters[2,:] .< 0 )
                    continue
                end

                evennewerlikelihood = likelihood( hmms[3-i] )[end]
                if evennewerlikelihood > newerlikelihood
                    newerlikelihood = evennewerlikelihood
                    optimalparameters = getparameters( hmms[3-i] )
                    optimalt = t
                end
                t = t/2
            end
            if newerlikelihood > newlikelihood
                if debug >= 2
                    println( "Accepting acceleration with $(newerlikelihood) and t=$optimalt" )
                end
                newlikelihood = newerlikelihood
                i = 3-i
                setparameters!( hmms[i], optimalparameters )
                push!( intermediates, optimalparameters )
            else
                if debug >= 2
                    println( "Acceleration not accepted with $(newerlikelihood)" )
                end
            end
            nextacceleration += acceleration
        end
    end

    if finishbig > 0
        bighmm = HMM{N,Dist,BigFloat,Out}(
            hmm.initialprobabilities,
            hmm.transitionprobabilities,
            hmm.stateparameters,
        )
        setobservations( bighmm, hmm.y )

        # no acceleration
        em( bighmm,
            epsilon=epsilon,
            debug=debug,
            maxiterations=finishbig,
            max_iter=max_iter,
            print_level=print_level,
            )

        setparameters!( hmm, getparameters( bighmm ) )
        clear( hmm )
    end

    hmm.initialprobabilities = hmms[i].initialprobabilities
    hmm.transitionprobabilities = hmms[i].transitionprobabilities
    hmm.stateparameters = hmms[i].stateparameters
    hmm.scratch = hmms[i].scratch

    hmm.scratch[:iterations] = iterations
    hmm.scratch[:time] = Base.time() - t0
    if keepintermediates
        hmm.scratch[:intermediates] = intermediates
    end

    if debug >= 1
        println( "Final likelihood = $(HMMs.likelihood(hmm)[end]); iterations = $iterations, time = $(hmm.scratch[:time])" )
        flush(stdout)
    end

    # the Model infrastructure requires this
    alpha = hmm.alpha.data[end,:]
    hmm.currentprobabilities = convert( Vector{Out}, alpha./sum(alpha) )
    reorder!( hmm )
    clear( hmm )

    return hmm
end

time( hmm::HMM ) = hmm.scratch[:time]

iterations( hmm::HMM ) = hmm.scratch[:iterations]

function permutederror( hmm1::GaussianHMM, hmm2::GaussianHMM )
    # use transition matrix to find best permutation; apply to means and stds
    n = length(hmm1.initialprobabilities)
    @assert( n == length(hmm2.initialprobabilities) )
    minerror = Inf
    minperm = nothing
    tp1 = convert( Matrix{Float64}, hmm1.transitionprobabilities )
    tp2 = convert( Matrix{Float64}, hmm2.transitionprobabilities )
    for perm in permutations(1:n)
        error = norm( tp1[perm, perm] - tp2, Inf )
        if error < minerror
            minerror = error
            minperm = perm
        end
    end
    meanerror = norm( hmm1.stateparameters[1,minperm] - hmm2.stateparameters[1,:], Inf )
    stderror = norm( hmm1.stateparameters[2,minperm] - hmm2.stateparameters[2,:], Inf )
    return (transitionprobabilities=minerror, means=meanerror, stds=stderror)
end

function applyperm!( hmm, perm )
    for i = 1:size(hmm.stateparameters,1)
        hmm.stateparameters[i,:] = hmm.stateparameters[i,perm]
    end
    hmm.initialprobabilities = hmm.initialprobabilities[perm]
    hmm.transitionprobabilities = hmm.transitionprobabilities[perm,perm]
    clear( hmm )
    return hmm
end

reorder!( hmm ) = applyperm!( hmm, sortperm(hmm.stateparameters[1,:]) )

ppv( io, s, v ) = println( io, rpad( s, 32 ), join([@sprintf("%8.2f", 100*convert(Float64,x)) for x in v]) )

function Base.show( io::IO, hmm::HMM{N,T} ) where {N,T <: Union{Normal,Laplace}}
    println( io )
    ppv( io, "initial probabilities:", hmm.initialprobabilities )
    println( io )
    ppv( io, "transition probabilities:", hmm.transitionprobabilities[1,:] )
    for i = 2:size(hmm.transitionprobabilities,1)
        ppv( io, "", hmm.transitionprobabilities[i,:] )
    end
    println( io )
    ppv( io, "locations:", hmm.stateparameters[1,:] )
    ppv( io, "scales:", hmm.stateparameters[2,:] )
end

function Base.show( io::IO, hmm::HMM{N,GenTDist} ) where {N}
    println( io )
    ppv( io, "initial probabilities:", hmm.initialprobabilities )
    println( io )
    ppv( io, "transition probabilities:", hmm.transitionprobabilities[1,:] )
    for i = 2:size(hmm.transitionprobabilities,1)
        ppv( io, "", hmm.transitionprobabilities[i,:] )
    end
    println( io )
    ppv( io, "locations:", hmm.stateparameters[1,:] )
    ppv( io, "scales:", hmm.stateparameters[2,:] )
    ppv( io, "nus:", hmm.stateparameters[3,:] )
end

function draw( outputfile::String, hmm::HMM )
    inputfile = tempname()

    reorder!( hmm )

    open( inputfile, "w" ) do io
        println( io, "digraph G {" )
        for index in CartesianIndices(hmm.transitionprobabilities)
            (i,j) = Tuple(index)
            weight = sqrt(1 - hmm.transitionprobabilities[index])
            value = @sprintf( "%02x", Int(round(convert(Float64,255*weight))) )
            if value != "00"
                println( io, "$i -> $j [color=\"#$value$value$value\"];" )
            end
        end
        println( io, "}" )
    end
    run( `dot -Tpdf $inputfile -o $outputfile` )
    rm( inputfile )
end

Models.initializationcount( ::Type{HMM{N,Dist,Calc,Out,T}} ) where {N,Dist,Calc,Out,T} = 0

function Models.initialize( hmm::HMM{N,Dist,Calc,Out,T}, args... ) where {N,Dist,Calc,Out,T}
    hmm.currentprobabilities[:] = hmm.initialprobabilities
    free( hmm )
end

function roll( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T}
    hmm.currentprobabilities[:] = hmm.currentprobabilities' * hmm.transitionprobabilities
    clear( hmm )
end

function Models.update( hmm::HMM{N,Dist,Calc,Out,T}, t::T, u::Out ) where {N,Dist,Calc,Out,T}
    push!( hmm.y, u )
    free( hmm )
    roll( hmm )
    probabilities = [pdf( Dist( hmm.stateparameters[:,i]... ), u ) for i in 1:length(hmm.initialprobabilities)]
    alpha = hmm.currentprobabilities .* probabilities
    hmm.currentprobabilities[:] = alpha/sum(alpha)
end

Dependencies.compress( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T} = free( hmm )

Models.state( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T} = hmm.currentprobabilities

Models.rootmodel( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T} = hmm

Models.getcompressedparameters( hmm::HMM{N,Dist,Calc,Out,T} ) where {N,Dist,Calc,Out,T} =
    hmm.basis' * getparameters( hmm )

function Models.setcompressedparameters!( hmm::HMM{N,Dist,Calc,Out,T}, p::AbstractVector{Float64} ) where {N,Dist,Calc,Out,T}
    A = hmm.constraintmatrix
    b = hmm.constraintvector
    B = hmm.basis
    setparameters!( hmm, A'*inv(A*A')*b + B*p )
end

end # module


text:
#TODO
function hclust(data, dist)
end

function random_init(features,k)
  features = shuffle!(features)
  centroids = Array(FeatureVector,k)
  while k > 0
    centroids[k] = features[k]
    k-=1
  end
  return centroids
end
random_init(features,k,dist_func) = random_init(features,k)

function max_min_init(features,k,dist_func)
  rand_num = (abs(rand(Int64)%Base.length(features)) + 1)
  orig_cent = features[rand_num]
  centroids = vcat(orig_cent)
  cents_to_be_found = k-1
  while cents_to_be_found > 0
    max_min_dist = 0
    next_cent = FeatureVector()
    for fv in features
      min_dist = Inf
      for centroid in centroids
        cent_dist = dist_func(centroid,fv)
        if cent_dist < min_dist
          min_dist = cent_dist
        end
      end
      if min_dist > max_min_dist
        max_min_dist = min_dist
        next_cent = fv
      end
    end
    centroids = vcat(centroids,next_cent)
    cents_to_be_found -= 1
  end
  return centroids
end

function kmeans(clust::Dict, cents::Array=[], k=iceil(sqrt(length(clust)/2)), init_cent_func=max_min_init, dist_func=dist_cos, max_iter=10000)
  # find initial k centroids
  features = collect(Base.values(clust))
  clust_keys = collect(Base.keys(clust))

  j = 1
  if k == 1
    single_clust = Cluster()
    for fv in features
      single_clust[clust_keys[j]] = fv
      j += 1
    end
    return [single_clust]
  end

  # check if user sent in own array of centroids
  if Base.length(cents) == 0
    centroids = init_cent_func(features,k,dist_func)
  else
    centroids = cents
    length_array = Base.length(centroids)
    if k > length_array
      Base.warn("The k($k) you entered is bigger than the amount of centroids in the array, reverting k to $length_array")
    end
  end

  # make Array of k clusters
  new_clusters = []
  for centroid in centroids
    new_clusters = vcat(new_clusters, Cluster())
  end

  # find distance between fv and centroid
  iteration = 1
  changed = true
  while changed && iteration < max_iter
    println("Start iteration: $iteration")
    i = 1
    for fv in features
      dist = Inf
      j = 1
      min_dist_cluster = Cluster()
      for cluster in new_clusters
        current_dist = dist_func(centroids[j],fv)
        if current_dist < dist
          dist = current_dist
          min_dist_cluster = cluster
        end
        j += 1
      end
      min_dist_cluster[clust_keys[i]] = fv
      i+=1
    end

    # recompute new centroids
    old_centroids = centroids
    new_centroids = Array(Any,length(centroids))
    x = 1
    for cluster in new_clusters
      new_cent = centroid(cluster)
      new_centroids[x] = new_cent
      x+=1
    end

    # checking if centroids moved.
    changed = false
    i = 1
    for centroid in old_centroids
      dist = dist_func(centroid,new_centroids[i])
      if dist > .000001
        changed = true
        centroids = new_centroids
        break
      end
      i += 1
    end

    # reset clusters if there are no changed.
    if changed
      new_clusters = []
      for centroid in centroids
        new_clusters = vcat(new_clusters, Cluster())
      end
    end
    iteration += 1
  end

  return new_clusters
end
kmeans(clust::Dict, k, init_cent_func, dist_func) = kmeans(clust,[],k,init_cent_func,dist_func)
kmeans(clust::Dict, k, init_cent_func) = kmeans(clust,[],k,init_cent_func)
kmeans(clust::Dict, k) = kmeans(clust,[],k)

function elbow_method(clust::Dict, dist_func::Function, low_bound, high_bound)
  if high_bound > length(clust)
    k = high_bound
    high_bound = length(clust)
    Base.warn("The high_bound($k) you entered is bigger than the amount of centroids in the array, reverting k to $high_bound")
  end
  if low_bound < 1
    k = low_bound
    low_bound = 1
    Base.warn("The low_bound($k) you entered is bigger than the amount of centroids in the array, reverting k to $low_bound")
  end

  temp_low = copy(low_bound)
  distances = []
  elbow_array = []
  while temp_low <= high_bound
    println("Start Elbow cluster $temp_low")
    clusters = kmeans(clust,temp_low,max_min_init,dist_func)
    avg_dist = 0
    for cluster in clusters
      features = values(cluster.vectors)
      cluster_avg_dist = 0
      cent = centroid(cluster)
      for fv in features
        cluster_avg_dist += dist_func(cent,fv)
      end
       avg_dist += (cluster_avg_dist/length(features))
    end
    distances = vcat(distances,(avg_dist/length(clusters)))
    elbow_array = vcat(elbow_array,(temp_low,(avg_dist/length(clusters))))
    temp_low += 1
  end
  println(scatterplot(collect(low_bound:high_bound),distances,sym='*'))

  return elbow_array
end

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
abstract type PropertyPrimitive end

# Meaningless isless function used to sort in optimize_batching
function Base.isless(a::T, b::T) where T <: PropertyPrimitive
    for field in fieldnames(T)
        x = getfield(a, field)
        y = getfield(b, field)
        if isa(x, Colorant)
            if color_isless(x, y)
                return true
            elseif color_isless(y, x)
                return false
            end
        else
            if x < y
                return true
            elseif x > y
                return false
            end
        end
    end
    return false
end

struct Property{P <: PropertyPrimitive} <: ComposeNode
    primitives::Vector{P}
end

isempty(p::Property) = isempty(p.primitives)

isscalar(p::Property) = length(p.primitives) == 1

# Some properties can be applied multiple times, most cannot.
isrepeatable(p::Property) = false

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::Property{T}) where T =
        Property{T}([resolve(box, units, t, primitive) for primitive in p.primitives])

# Property primitive catchall: most properties don't need measure transforms
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, primitive::PropertyPrimitive) = primitive

# Stroke
# ------

struct StrokePrimitive <: PropertyPrimitive
	color::RGBA{Float64}
end

const Stroke = Property{StrokePrimitive}

stroke(c::Nothing) = Stroke([StrokePrimitive(RGBA{Float64}(0, 0, 0, 0))])
stroke(c::Union{Colorant, AbstractString}) = Stroke([StrokePrimitive(parse_colorant(c))])
stroke(cs::AbstractArray) = Stroke([StrokePrimitive(c == nothing ?
        RGBA{Float64}(0, 0, 0, 0) : parse_colorant(c)) for c in cs])

prop_string(::Stroke) = "s"

# Fill
# ----

struct FillPrimitive <: PropertyPrimitive
	color::RGBA{Float64}
end

const Fill = Property{FillPrimitive}

fill(c::Nothing) = Fill([FillPrimitive(RGBA{Float64}(0.0, 0.0, 0.0, 0.0))])

"""
    fill(c)

Define a fill color, where `c` can be a `Colorant` or `String`.
"""
fill(c::Union{Colorant, AbstractString}) = Fill([FillPrimitive(parse_colorant(c))])

"""
    fill(cs::AbstractArray)

Arguments can be passed in arrays in order to perform multiple drawing operations at once.
"""
fill(cs::AbstractArray) = Fill([FillPrimitive(c == nothing ?
        RGBA{Float64}(0.0, 0.0, 0.0, 0.0) : parse_colorant(c)) for c in cs])

prop_string(::Fill) = "f"

# StrokeDash
# ----------

struct StrokeDashPrimitive <: PropertyPrimitive
    value::Vector{Measure}
end

const StrokeDash = Property{StrokeDashPrimitive}

strokedash(values::AbstractArray) = StrokeDash([StrokeDashPrimitive(collect(Measure, values))])
strokedash(values::AbstractArray{<:AbstractArray}) =
        StrokeDash([StrokeDashPrimitive(collect(Measure, value)) for value in values])

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, primitive::StrokeDashPrimitive) =
        StrokeDashPrimitive([resolve(box, units, t, v) for v in primitive.value])

prop_string(::StrokeDash) = "sd"

# StrokeLineCap
# -------------

abstract type LineCap end
struct LineCapButt <: LineCap end
struct LineCapSquare <: LineCap end
struct LineCapRound <: LineCap end

struct StrokeLineCapPrimitive <: PropertyPrimitive
    value::LineCap
end

StrokeLineCapPrimitive(value::Type{LineCap}) = StrokeLineCapPrimitive(value())

const StrokeLineCap = Property{StrokeLineCapPrimitive}

strokelinecap(value::Union{LineCap, Type{LineCap}}) =
        StrokeLineCap([StrokeLineCapPrimitive(value)])
strokelinecap(values::AbstractArray) =
        StrokeLineCap([StrokeLineCapPrimitive(value) for value in values])

prop_string(::StrokeLineCap) = "slc"

# StrokeLineJoin
# --------------

abstract type LineJoin end
struct LineJoinMiter <: LineJoin end
struct LineJoinRound <: LineJoin end
struct LineJoinBevel <: LineJoin end

struct StrokeLineJoinPrimitive <: PropertyPrimitive
    value::LineJoin
end

StrokeLineCapPrimitive(value::Type{LineJoin}) = new(value())

const StrokeLineJoin = Property{StrokeLineJoinPrimitive}

strokelinejoin(value::Union{LineJoin, Type{LineJoin}}) =
        StrokeLineJoin([StrokeLineJoinPrimitive(value)])
strokelinejoin(values::AbstractArray) =
        StrokeLineJoin([StrokeLineJoinPrimitive(value) for value in values])

prop_string(::StrokeLineJoin) = "slj"

# LineWidth
# ---------

struct LineWidthPrimitive <: PropertyPrimitive
    value::Measure

    function LineWidthPrimitive(value)
        return new(size_measure(value))
    end
end

const LineWidth = Property{LineWidthPrimitive}

linewidth(value::Union{Measure, Number}) = LineWidth([LineWidthPrimitive(value)])
linewidth(values::AbstractArray) =
        LineWidth([LineWidthPrimitive(value) for value in values])

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, primitive::LineWidthPrimitive) =
        LineWidthPrimitive(resolve(box, units, t, primitive.value))

prop_string(::LineWidth) = "lw"

# Visible
# -------

struct VisiblePrimitive <: PropertyPrimitive
    value::Bool
end

const Visible = Property{VisiblePrimitive}

visible(value::Bool) = Visible([VisiblePrimitive(value)])
visible(values::AbstractArray) = Visible([VisiblePrimitive(value) for value in values])

prop_string(::Visible) = "v"

# FillOpacity
# -----------

struct FillOpacityPrimitive <: PropertyPrimitive
    value::Float64

    function FillOpacityPrimitive(value_::Number)
        value = Float64(value_)
        (value < 0.0 || value > 1.0) && error("Opacity must be between 0 and 1.")
        return new(value)
    end
end

const FillOpacity = Property{FillOpacityPrimitive}

"""
    fillopacity(value)

Define a fill opacity, where 0≤value≤1.  For svg, nested contexts  will inherit from parent contexts e.g. `(context(), fillopacity(a), (context(), fill(c::String), circle()))`.
"""
fillopacity(value::Float64) = FillOpacity([FillOpacityPrimitive(value)])

"""
    fillopacity(values::AbstractArray)

Arguments can be passed in arrays in order to perform multiple drawing operations at once.
"""
fillopacity(values::AbstractArray) =
        FillOpacity([FillOpacityPrimitive(value) for value in values])

prop_string(::FillOpacity) = "fo"

# StrokeOpacity
# -------------

struct StrokeOpacityPrimitive <: PropertyPrimitive
    value::Float64

    function StrokeOpacityPrimitive(value_::Number)
        value = Float64(value_)
        (value < 0.0 || value > 1.0) && error("Opacity must be between 0 and 1.")
        return new(value)
    end
end

const StrokeOpacity = Property{StrokeOpacityPrimitive}

strokeopacity(value::Float64) = StrokeOpacity([StrokeOpacityPrimitive(value)])
strokeopacity(values::AbstractArray) =
        StrokeOpacity([StrokeOpacityPrimitive(value) for value in values])

prop_string(::StrokeOpacity) = "so"

# Clip
# ----

struct ClipPrimitive{P <: Vec} <: PropertyPrimitive
    points::Vector{P}
end

const Clip = Property{ClipPrimitive}

clip() = Clip([ClipPrimitive(Array{Vec}(undef, 0))])

"""
    clip(points::AbstractArray)

`clip()` is a property.  Only forms inside the clip shape will be visible.
"""
function clip(points::AbstractArray{T}) where T <: XYTupleOrVec
    XM, YM = narrow_polygon_point_types(Vector[points])
    if XM == Any
        XM = Length{:cx, Float64}
    end
    if YM == Any
        YM = Length{:cy, Float64}
    end
    VecType = Tuple{XM, YM}

    prim = ClipPrimitive(VecType[(x_measure(point[1]), y_measure(point[2])) for point in points])
    return Clip(typeof(prim)[prim])
end

"""
    clip(point_arrays::AbstractArray...)

Arguments can be passed in arrays in order to perform multiple clipping operations at once.
"""
function clip(point_arrays::AbstractArray...)
    XM, YM = narrow_polygon_point_types(point_arrays)
    VecType = XM == YM == Any ? Vec : Vec{XM, YM}
    PrimType = XM == YM == Any ? ClipPrimitive : ClipPrimitive{VecType}

    clipprims = Array{PrimType}(undef, length(point_arrays))
    for (i, point_array) in enumerate(point_arrays)
        clipprims[i] = ClipPrimitive(VecType[(x_measure(point[1]), y_measure(point[2]))
                                             for point in point_array])
    end
    return Property{PrimType}(clipprims)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, primitive::ClipPrimitive) =
        ClipPrimitive{AbsoluteVec2}( AbsoluteVec2[
            resolve(box, units, t, point) for point in primitive.points])

prop_string(::Clip) = "clp"

# Font
# ----

struct FontPrimitive <: PropertyPrimitive
    family::AbstractString
end

const Font = Property{FontPrimitive}

font(family::AbstractString) = Font([FontPrimitive(family)])
font(families::AbstractArray) =
        Font([FontPrimitive(family) for family in families])

prop_string(::Font) = "fnt"

Base.hash(primitive::FontPrimitive, h::UInt) = hash(primitive.family, h)

==(a::FontPrimitive, b::FontPrimitive) = a.family == b.family

# FontSize
# --------

struct FontSizePrimitive <: PropertyPrimitive
    value::Measure

    function FontSizePrimitive(value)
        return new(size_measure(value))
    end
end

const FontSize = Property{FontSizePrimitive}

fontsize(value::Union{Number, Measure}) = FontSize([FontSizePrimitive(value)])

fontsize(values::AbstractArray) = FontSize([FontSizePrimitive(value) for value in values])

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, primitive::FontSizePrimitive) =
        FontSizePrimitive(resolve(box, units, t, primitive.value))

prop_string(::FontSize) = "fsz"

# SVGID
# -----

struct SVGIDPrimitive <: PropertyPrimitive
    value::AbstractString
end

const SVGID = Property{SVGIDPrimitive}

svgid(value::AbstractString) = SVGID([SVGIDPrimitive(value)])
svgid(values::AbstractArray) = SVGID([SVGIDPrimitive(value) for value in values])

prop_string(::SVGID) = "svgid"

Base.hash(primitive::SVGIDPrimitive, h::UInt) = hash(primitive.value, h)

==(a::SVGIDPrimitive, b::SVGIDPrimitive) = a.value == b.value

# SVGClass
# --------

struct SVGClassPrimitive <: PropertyPrimitive
    value::String
end

const SVGClass = Property{SVGClassPrimitive}

svgclass(value::AbstractString) = SVGClass([SVGClassPrimitive(value)])
svgclass(values::AbstractArray) =
        SVGClass([SVGClassPrimitive(value) for value in values])

function prop_string(svgc::SVGClass)
    if isscalar(svgc)
        return string("svgc(", svgc.primitives[1].value, ")")
    else
        return string("svgc(", svgc.primitives[1].value, "...)")
    end
end

Base.hash(primitive::SVGClassPrimitive, h::UInt) = hash(primitive.value, h)

==(a::SVGClassPrimitive, b::SVGClassPrimitive) = a.value == b.value

# SVGAttribute
# ------------

struct SVGAttributePrimitive <: PropertyPrimitive
    attribute::String
    value::String
end

const SVGAttribute = Property{SVGAttributePrimitive}

svgattribute(attribute::AbstractString, value) =
        SVGAttribute([SVGAttributePrimitive(attribute, string(value))])
svgattribute(attribute::AbstractString, values::AbstractArray) =
        SVGAttribute([SVGAttributePrimitive(attribute, string(value)) for value in values])

svgattribute(attributes::AbstractArray, values::AbstractArray) =
    SVGAttribute( @makeprimitives SVGAttributePrimitive,
            (attribute in attributes, value in values),
            SVGAttributePrimitive(attribute, string(value)))

prop_string(::SVGAttribute) = "svga"

function Base.hash(primitive::SVGAttributePrimitive, h::UInt)
    h = hash(primitive.attribute, h)
    h = hash(primitive.value, h)
    return h
end

==(a::SVGAttributePrimitive, b::SVGAttributePrimitive) =
        a.attribute == b.attribute && a.value == b.value

# JSInclude
# ---------

struct JSIncludePrimitive <: PropertyPrimitive
    value::AbstractString
    jsmodule::Union{Nothing, Tuple{AbstractString, AbstractString}}
end

const JSInclude = Property{JSIncludePrimitive}

jsinclude(value::AbstractString, module_name=nothing) =
        JSInclude([JSIncludePrimitive(value, module_name)])

# Don't bother with a vectorized version of this. It wouldn't really make #
# sense.

prop_string(::JSInclude) = "jsip"

# JSCall
# ------

struct JSCallPrimitive <: PropertyPrimitive
    code::AbstractString
    args::Vector{Measure}
end

const JSCall = Property{JSCallPrimitive}

jscall(code::AbstractString, arg::Vector{Measure}=Measure[]) = JSCall([JSCallPrimitive(code, arg)])

jscall(codes::AbstractArray, args::AbstractArray{Vector{Measure}}=Vector{Measure}[Measure[]]) =
        JSCall( @makeprimitives JSCallPrimitive,
            (code in codes, arg in args),
            JSCallPrimitive(code, arg))

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform, primitive::JSCallPrimitive)
    # we are going to build a new string by scanning across "code" and
    # replacing %x with translated x values, %y with translated y values
    # and %s with translated size values.
    newcode = IOBuffer()

    i = 1
    validx = 1
    while true
        j = findnext(primitive.code, "%", i)

        if j === nothing
            write(newcode, primitive.code[i:end])
            break
        end

        write(newcode, primitive.code[i:j-1])
        if j == length(primitive.code)
            write(newcode, '%')
            break
        elseif primitive.code[j+1] == '%'
            write(newcode, '%')
        elseif primitive.code[j+1] == 'x'
            val = resolve(box, units, t, (primitive.args[validx], 0mm))
            write(newcode, svg_fmt_float(val[1].value))
            validx += 1
        elseif primitive.code[j+1] == 'y'
            val = resolve(box, units, t, (0mm, primitive.args[validx]))
            write(newcode, svg_fmt_float(val[2].value))
            validx += 1
        elseif primitive.code[j+1] == 's'
            val = resolve(box, units, t, primitive.args[validx])
            write(newcode, svg_fmt_float(val.value))
            validx += 1
        else
            write(newcode, '%', primitive.code[j+1])
        end

        i = j + 2
    end

    return JSCallPrimitive(String(take!(newcode)), Measure[])
end

isrepeatable(p::JSCall) = true

Base.isless(a::FillPrimitive, b::FillPrimitive) = color_isless(a.color, b.color)
Base.isless(a::StrokePrimitive, b::StrokePrimitive) = color_isless(a.color, b.color)

prop_string(::JSCall) = "jsc"

# Arrow Property
# -------

struct ArrowPrimitive <: PropertyPrimitive
    value::Bool
end

const Arrow = Property{ArrowPrimitive}

"""
    arrow()

    `arrow() = arrow(true)`
"""
arrow() = Arrow([ArrowPrimitive(true)])

"""
    arrow(value::Bool)

    `arrow()` is a property of arcs, lines and curves. The color of the arrowhead is the same as `stroke()`, but for svg the results will be browser-dependent.
"""
arrow(value::Bool) = Arrow([ArrowPrimitive(value)])

"""
    arrow(values::AbstractArray)

Arguments can be passed in arrays in order to perform multiple drawing operations at once.
"""
arrow(values::AbstractArray) = Arrow([ArrowPrimitive(value) for value in values])

prop_string(::Arrow) = "arrow"


text:
# call with
# ARGS[1] = T
# ARGS[2] = frac err
# ARGS[3] = nameOutput
# ARGS[4] = port1
# ARGS[5] = port2
@show ARGS

using Caesar, IncrementalInference, RoME, HDF5, JLD
d = jldopen("data/VicPrkBasic.jld", "r") do file
  read(file, "dBasic")
end
f = jldopen("data/VicPrkBasic.jld", "r") do file
  read(file, "fBasic")
end
# MM = jldopen("test/data/VicPrkBasic.jld", "r") do file
#   read(file, "MMBasic")
# end
include("data/Vic120MM.jl")
MMr = MMwithRedirects(MM);
frac = parse(Float64,ARGS[2])
MMe = MMRandErrs(MMr,frac=frac)
T=parse(Int,ARGS[1])

isamdict, fgu = doISAMSolve(d,f,toT=T, savejld=false, MM=MMe, retfg=true, port1=parse(Int,ARGS[4]),port2=parse(Int,ARGS[5]))

# isamdict = jldopen("data/fgWithISAMrefERR01.jld", "r") do file
#   read(file, "isamdict")
# end
# fgu = jldopen("data/fgWithISAMrefERR01.jld", "r") do file
#   read(file, "fgu")
# end

fg = emptyFactorGraph();
idx = appendFactorGraph!(fg, d, f, toT=T, lcmode=:mmodal, MM=MMe);
tree = prepBatchTree!(fg, ordering=:qr,drawpdf=true);
inferOverTree!(fg,tree, N=100);
fg1 = deepcopy(fg)
@save "results/fgT$(T)_$(ARGS[3])_perErr$(ARGS[2]).jld" fg1 isamdict fgu MM MMr MMe T frac

tic()
inferOverTree!(fg,tree, N=100);
comptime = toc()
fg2 = deepcopy(fg)
@save "results/fgT$(T)_$(ARGS[3])_perErr$(ARGS[2]).jld" fg1 fg2 isamdict fgu MM MMr MMe T frac comptime

tic()
inferOverTree!(fg,tree, N=100);
comptime = toc()
fg3 = fg
@save "results/fgT$(T)_$(ARGS[3])_perErr$(ARGS[2]).jld" fg1 fg2 fg3 isamdict fgu MM MMr MMe T frac comptime


text:
export HDF5Array, extractsparse, extractdense, sparse
import HDF5
import SparseArrays

"""
The `HDF5Array` is an abstract type that describes the concept of a HDF5-backed array abstraction.
In other words, data is stored inside a HDF5 file and is retrieved on demand rather than being loaded into memory.
`T` is the type of the data while `N` is the dimensionality.

Concrete subtypes include the [`DenseHDF5Array`](@ref) and the [`SparseHDF5Matrix`](@ref).
Subtypes are expected to implement the [`extractdense`](@ref) and [`extractsparse`](@ref) methods.
They may also override the `SparseArrays.issparse` method to indicate whether they contain sparse data.

We provide conversion functions to quickly create in-memory `Array`s or `SparseMatrixCSC` objects from any `HDF5Array`s.
In addition, subsetting operations will automatically views on the original array rather than immediately loading data from file.
This enables lazy evaluation for memory-efficient operations on a subset of the dataset.

Currently, all `HDF5Array`s are read-only; calling `setindex!` will fail.
"""
abstract type HDF5Array{T, N} <: AbstractArray{T, N} end

"""
    getindex(x::HDF5Array{T,N}, I...)

Create a view into a `HDF5Array{T,N}` at the indices `I`.
This does not read any data from file, only acting as a delayed subsetting operation.

This view-returning behavior continues when `x` is already a view on a `HDF5Array`.
In that case, a new view is returned containing recomputed indices based on `I` and the `parentindices(x)`.

Calling [`extractdense`](@ref) or [`extractsparse`](@ref) on this view will use the indices to extract the desired subset.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> exampledense(tmp, "stuff", (20, 10))

julia> x = DenseHDF5Array(tmp, "stuff");

julia> y = x[1:5, :];

julia> isa(y, SubArray)
true

julia> z = y[1:2,1:5];

julia> isa(z, SubArray)
true
```
"""
function Base.getindex(x::HDF5Array{T,N}, I...) where {T,N}
    return view(x, I...)
end

function Base.getindex(x::SubArray{T,N,P,I,L}, indices...) where {T,N,P<:HDF5Array{T,N},I,L}
    sub = sub_indices(x, indices...)
    return view(parent(x), sub...)
end

"""
    setindex!(x::HDF5Array{T,N}, v, I::Vararg{Integer,N})

All `HDF5Array`s are strictly read-only, so any attempt to call this function will throw an error.
"""
function Base.setindex!(x::HDF5Array{T,N}, v, I::Vararg{Integer,N}) where {T,N}
    throw(ErrorException(typeof(x) * " instances are strictly read-only"))
end

"""
    Array(x::HDF5Array{T,N}})

Convert a `HDF5Array` into an in-memory `Array` of the same data type and dimension.
This is equivalent to using [`extractdense`](@ref) while requesting the full extent of each dimension.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> exampledense(tmp, "stuff", (20, 10))

julia> x = DenseHDF5Array(tmp, "stuff");

julia> y = Array(x);

julia> size(y)
(20, 10)
```
"""
function Array{T,N}(x::HDF5Array{T,N}) where {T,N}
    colons = Vector{Any}(undef, N)
    fill!(colons, :)
    return extractdense(x, colons...)
end

"""
    sparse(x::HDF5Array{T,N})

Convert a 2-dimensional `HDF5Array` into an in-memory sparse matrix of the same data type and dimension,
assuming that `T` is some numeric or boolean type.
This is equivalent to using [`extractsparse`](@ref) while requesting the full extent of each dimension.
Note that this only makes sense if `x` contains a high proportion of zeros.

# Examples
```jldoctest
julia> using HDF5Arrays, SparseArrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.2)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> y = sparse(x);

julia> typeof(y)
SparseMatrixCSC{Float64, Int64}
```
"""
function SparseArrays.sparse(x::HDF5Array{T,2}) where {T<:Union{Number,Bool}}
    return extractsparse(x, :, :)
end

function sub_indices(x::SubArray{T,N,P,I,L}, indices...) where {T,N,P<:HDF5Array{T,N},I,L}
    sub = Vector{Any}(undef, N)
    idx = to_indices(x, indices)
    par = parentindices(x)
    for i in 1:N
        sub[i] = par[i][idx[i]]
    end
    return sub
end

"""
    extractdense(x::SubArray{T,N,P<:HDF5Array{T,N},I,L}, indices...; blockdim = nothing)

Extract an in-memory dense `Array` from a `SubArray` `x` of a parent `HDF5Array`.
The returned matrix contains the same values as `x[indices...]`, in addition to any subsetting used to create the view on the parent.

`blockdim` is forwarded to the method for the parent `HDF5Array`.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> exampledense(tmp, "stuff", (20, 10))

julia> x = DenseHDF5Array(tmp, "stuff");

julia> y = x[2:10, [2,3,5,9,10]];

julia> z = extractdense(y, :, :);

julia> size(z)
(9, 5)
```
"""
function extractdense(x::SubArray{T,N,P,I,L}, indices...; blockdim = nothing) where {T,N,P<:HDF5Array{T,N},I,L}
    sub = sub_indices(x, indices...)
    return extractdense(parent(x), sub...; blockdim = blockdim)
end

"""
    extractsparse(x::SubArray{T,N,P<:HDF5Array{T,N},I,L}, i, j; blockdim = nothing)

Extract a sparse matrix from a `SubArray` `x` of a parent `HDF5Array`.
The returned matrix contains the same values as `x[i, j]`, in addition to any subsetting used to create the view on the parent.

`blockdim` is forwarded to the method for the parent `HDF5Array`.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.1)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> y = x[2:10, [2,3,5,9,10]];

julia> z = extractsparse(y, :, :);

julia> size(z)
(9, 5)
```
"""
function extractsparse(x::SubArray{T,2,P,I,L}, i, j; blockdim = nothing) where {T,P<:HDF5Array{T,2},I,L}
    sub = sub_indices(x, i, j)
    return extractsparse(parent(x), sub...; blockdim = blockdim)
end

"""
    Array(x::SubArray{T,N,P<:HDF5Array{T,N},I,L})

Convert a view of a `HDF5Array` into an in-memory `Array` of the same type and dimension.
This is equivalent to using [`extractdense`](@ref) while requesting the full extent of each dimension.

# Examples
```jldoctest
julia> using HDF5Arrays

julia> tmp = tempname();

julia> exampledense(tmp, "stuff", (20, 10))

julia> x = DenseHDF5Array(tmp, "stuff");

julia> y = Array(x[1:10,2:5]);

julia> size(y)
(10, 4)
```
"""
function Array{T,N}(x::SubArray{T,N,P,I,L}) where {T,N,P<:HDF5Array{T,N},I,L}
    colons = Vector{Any}(undef, N)
    fill!(colons, :)
    return extractdense(x, colons...)
end

"""
    sparse(x::SubArray{T,2,P<:HDF5Array{T,2},I,L})

Convert a view of a 2-dimensional `HDF5Array` into an in-memory sparse matrix of the same type and dimension.
This is equivalent to using [`extractsparse`](@ref) while requesting the full extent of each dimension.

# Examples
```jldoctest
julia> using HDF5Arrays, SparseArrays

julia> tmp = tempname();

julia> examplesparse(tmp, "stuff", (20, 10), 0.2)

julia> x = SparseHDF5Matrix(tmp, "stuff");

julia> y = sparse(x[1:10,2:5]);

julia> typeof(y)
SparseMatrixCSC{Float64, Int64}
```
"""
function SparseArrays.sparse(x::SubArray{T,2,P,I,L}) where {T,P<:HDF5Array{T,2},I,L}
    return extractsparse(x, :, :)
end


text:
t=0:1/100:1
y = 2*sin.(5*t)+3*sin.(10*t)+2*sin.(30*t)
pks, vals = findmaxima(y)
_, proms = peakprom(pks, y; strict=false)
promlinesx = vec(vcat(reshape(repeat(t[pks]; inner=2),2,length(pks)), fill(NaN, (1,length(pks)))))
promlinesy = vec([vals - proms vals fill(NaN, (length(pks),1))]');

_, widths, lower, upper = peakwidth(pks, y, proms; strict=true, relheight=prevfloat(1.))
fullwidthlinesx = vec([(lower .- 1)./fs (upper .- 1)./fs fill(NaN, (length(pks),1))]')
fullwidthlinesy = vec(vcat(reshape(repeat(vals - proms.*1; inner=2),2,length(pks)), fill(NaN, (1,length(pks)))))

_, widths, lower, upper = peakwidth(pks, y, proms; strict=true)
halfwidthlinesx = vec([(lower .- 1)./fs (upper .- 1)./fs fill(NaN, (length(pks),1))]')
halfwidthlinesy = vec(vcat(reshape(repeat(vals - proms.*0.5; inner=2),2,length(pks)), fill(NaN, (1,length(pks)))))

p = plot([
    scatter(;x=t, y, line_color="black", name="Signal", line_shape="spline", legendrank=1),
    scatter(;x=promlinesx, y=promlinesy, mode="lines", line_color="#4063d8", name="Prominence", legendrank=3),
    scatter(;x=fullwidthlinesx, y=fullwidthlinesy, mode="lines", line_color="#4063d8", showlegend=false),
    scatter(;x=t[pks], y=vals, marker=attr(;color="#cb3c33", size=8), mode="markers", name="Maxima", legendrank=2),
    scatter(;x=halfwidthlinesx, y=halfwidthlinesy, line_dash="dot", line_color="grey", name="Width", legendrank=4),
    ], Layout(;
    width=600,
    height=400,
    margin=attr(;
        t=5,
        r=5,
        b=25,
        l=30,
        pad=5,
        autoexpand=false,
    ),
    font_size="16",

))
savefig(p, "../docs/src/assets/images/peaks_prom_width.png"; width=600, height=400, scale=4)

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
using DFTK
using Test
using LinearAlgebra: norm

include("testcases.jl")

function get_scf_energies(testcase, supersampling, functionals)
    Ecut=3
    grid_size=15
    scf_tol=1e-12  # Tolerance in total enengy
    n_bands = 10
    kcoords = [[.2, .3, .4]]

    # force symmetries to false because the symmetrization is weird at low ecuts
    model = model_DFT(testcase.lattice, testcase.atoms, testcase.positions, functionals; symmetries=false)
    fft_size = compute_fft_size(model, Ecut, kcoords;
                                supersampling, ensure_smallprimes=false, algorithm=:precise)

    kweights = [1]
    basis = PlaneWaveBasis(model, Ecut, kcoords, kweights; fft_size)
    scfres = self_consistent_field(basis; tol=scf_tol)
    values(scfres.energies)
end

@testset "Energy is exact for supersampling>2 without XC" begin
    energies = [get_scf_energies(silicon, supersampling, []) for supersampling in (1, 2, 3)]

    @test abs(sum(energies[1]) - sum(energies[2])) > 1e-10

    # supersampling == 2 is exact and going beyond has no effect on the total energy
    @test abs(sum(energies[2]) - sum(energies[3])) < 1e-10

    # Individual terms are not variational and only converged to the amount
    # the density and the orbitals are converged, which is about sqrt(energy)
    @test norm(energies[2] .- energies[3]) < 1e-5
end

@testset "Energy is not exact for supersampling>2 with XC" begin
    energies = [get_scf_energies(silicon, supersampling, [:lda_x, :lda_c_vwn])
                for supersampling in (1, 2, 3)]

    @test abs(sum(energies[1]) - sum(energies[2])) > 1e-10

    # supersampling == 2 is not exact in total energy for XC
    @test abs(sum(energies[2]) - sum(energies[3])) > 1e-10
end


text:
export evaluate

# NOTE: Point are sparse vectors, hence non referenced variables are assumed to be null.
function evaluate(p::Polynomial, pt::Point; partial=false)
    res=0
    expldeg = conjdeg = 0
    for (expo, λ) in p
        res += λ*evaluate(expo, pt; partial=partial)
    end
    typeof(res)<:Polynomial && update_degree!(res)

    # Fully evaluated polynomial
    if typeof(res) == Polynomial && res.degree == Degree(0,0)
        return first(res)[2]
    end
    return res
end

function evaluate(expo::Exponent, pt::Point; partial=false)
    res=1
    for (var, deg) in expo.expo
        res *= (evaluate(var, pt; partial=partial)^deg.explvar) * (conj(evaluate(var, pt; partial=partial))^deg.conjvar)
    end
    return res
end

function evaluate(x::Variable, pt::Point; partial=false)
    if !haskey(pt, x)
        return partial ? x : 0.0
    elseif x.kind<:Complex
        return pt[x]
    else
        return real(pt[x])
    end
end


text:
#### General Assumptions of these TNEP Models ####
#
#

export run_tnep

""
function run_tnep(file, model_constructor, solver; kwargs...)
    return run_generic_model(file, model_constructor, solver, post_tnep; solution_builder = get_tnep_solution, kwargs...)
end

"the general form of the tnep optimization model"
function post_tnep(pm::GenericPowerModel)
    variable_branch_ne(pm)
    variable_voltage(pm)
    variable_voltage_ne(pm)
    variable_generation(pm)
    variable_branch_flow(pm)
    variable_dcline_flow(pm)
    variable_branch_flow_ne(pm)

    objective_tnep_cost(pm)

    constraint_voltage(pm)
    constraint_voltage_ne(pm)

    for i in ids(pm, :ref_buses)
        constraint_theta_ref(pm, i)
    end

    for i in ids(pm, :bus)
        constraint_kcl_shunt_ne(pm, i)
    end

    for i in ids(pm, :branch)
        constraint_ohms_yt_from(pm, i)
        constraint_ohms_yt_to(pm, i)

        constraint_voltage_angle_difference(pm, i)

        constraint_thermal_limit_from(pm, i)
        constraint_thermal_limit_to(pm, i)
    end

    for i in ids(pm, :ne_branch)
        constraint_ohms_yt_from_ne(pm, i)
        constraint_ohms_yt_to_ne(pm, i)

        constraint_voltage_angle_difference_ne(pm, i)

        constraint_thermal_limit_from_ne(pm, i)
        constraint_thermal_limit_to_ne(pm, i)
    end

    for i in ids(pm, :dcline)
        constraint_dcline(pm, i)
    end
end

""
function get_tnep_solution(pm::GenericPowerModel, sol::Dict{String,Any})
    add_bus_voltage_setpoint(sol, pm)
    add_generator_power_setpoint(sol, pm)
    add_branch_flow_setpoint(sol, pm)
    add_branch_flow_setpoint_ne(sol, pm)
    add_branch_ne_setpoint(sol, pm)
end


text: group_by_evtno(table::TypedTables.Table) = TypedTables.Table(consgroupedview(table[sortperm(table.evtno)].evtno, Tables.columns(table)))


text:
module Widgets

using Observables, OrderedCollections
using Dates
using Colors: Colorant

import OrderedCollections: OrderedDict
import Observables: off, on, Observable, AbstractObservable, observe, ObservablePair, @map, @map!, @on

export Widget, widget, @layout!, node, scope, scope!

include("widget.jl")
include("utils.jl")
include("layout.jl")
include("backend.jl")
include("defaults.jl")
include("manipulate.jl")
include("modifiers.jl")

end # module


text:
abstract Layer{A<:AbstractArray, B<:AbstractArray}

# Pointwise layers always use the same dimensionality for input and output
abstract Pointwise{A<:AbstractArray} <: Layer{A, A}

abstract Container{A<:AbstractArray, B<:AbstractArray} <: Layer{A, B}

# By default, accGradParameters! does nothing
function accGradParameters!{A, B}(self :: Layer{A, B}, input :: A,
                                  gradOutput :: B, scale :: eltype(A))
end

# Call the scale version if scale unspecified
function accGradParameters!{A, B}(self :: Layer{A, B}, input :: A,
                                  gradOutput :: B)
    return accGradParameters!(self, input, one(eltype(A)))
end

function forward!{A, B}(self :: Layer{A, B}, input :: A)
    return updateOutput!(self, input)
end
export forward!

function backward!{A, B}(self :: Layer{A, B}, input :: A, gradOutput :: B,
                         scale :: eltype(A))
    updateGradInput!(self, input, gradOutput)
    accGradParameters!(self, input, gradOutput, scale)
    return self.gradInput
end

# Call the scale version if scale unspecified
function backward!{A, B}(self :: Layer{A, B}, input :: A, gradOutput :: B)
    return backward!(self, input, gradOutput, one(eltype(A)))
end
export backward!

# Add a layer to the container
function add_layer!{A, B}(self :: Container{A, B}, layer :: Layer)
    push!(self.modules, layer)
end
export add_layer!

function validate{A, B}(self :: Layer{A, B})
end

function validate{A, B}(self :: Container{A, B})
    for layer in self.modules
        validate(layer)
    end
end
export validate

# Iteration; layers return themselves; containers return themselves and
# children

import Base.start
import Base.done
import Base.next

function start(self :: Layer)
    return false
end

function done(self :: Layer, state)
    return state
end

function next(self :: Layer, state)
    return self, true
end

# state is (0, false) for initial setup (must return self next)
#          (idx, state_for_modules[idx]) afterwards
function start(self :: Container)
    return 0, false
end

function done(self :: Container, state)
    local idx, o_state
    idx, o_state = state
    if idx == 0 then
        return false
    end
    if idx > length(self.modules)
        return true
    end

    # Note that this check is simpler than in the general case, as we know
    # each layer will return at least one element (self). So, as long
    # as idx < length(self.modules), we know we're not done; this would no
    # longer be true if we allowed empty layers.
    return done(self.modules[idx], o_state) && idx == length(self.modules)
end

function next(self :: Container, state)
    local idx, o_state
    idx, o_state = state
    if idx == 0
        return self, (1, isempty(self.modules) ? true : start(self.modules[1]))
    end

    if done(self.modules[idx], o_state)
        idx += 1
        o_state = start(self.modules[idx])
    end

    item, o_state = next(self.modules[idx], o_state)
    return item, (idx, o_state)
end

function zeroGradParameters!(self :: Layer)
    for layer in self
        for field in (:gradWeight, :gradBias)
            if field in names(layer)
                local a = getfield(layer, field)
                fill!(a, zero(eltype(a)))
            end
        end
    end
end
export zeroGradParameters!

using Distributions
function initWeights!(self :: Layer)
    for layer in self
        if :weight in names(layer)
            rand!(Distributions.Uniform(-0.1, 0.1), layer.weight)
        end
        if :bias in names(layer)
            rand!(Distributions.Uniform(0, 0.1), layer.bias)
        end
    end
end
export initWeights!

function getParameters(self :: Layer)
    local params = Array((AbstractArray, AbstractArray), 0)
    for layer in self
        if :weight in names(layer) && :gradWeight in names(layer)
            push!(params, (layer.weight, layer.gradWeight))
        end
        if :bias in names(layer) && :gradBias in names(layer)
            push!(params, (layer.bias, layer.gradBias))
        end
    end
    return params
end
export getParameters

abstract Criterion{A<:AbstractArray, B<:AbstractArray}

function forward!{A,B}(c :: Criterion{A,B}, x :: A, t :: A)
    return updateOutput!(c, x, t)
end

function backward!{A,B}(c :: Criterion{A,B}, x :: A, t :: A)
    return updateGradInput!(c, x, t)
end

#######

type Parallel <: Container
    # TODO: Define my fields here...
end

type SoftMax <: Layer
    # TODO: Define my fields here...
end

type LogSoftMax <: Layer
    # TODO: Define my fields here...
end

type Affine <: Layer
    # TODO: Define my fields here...
end

type SpatialMaxPool <: Layer
    # TODO: Define my fields here...
end

type SpatialConv{T, I <: Integer} <: Layer
    nInputPlane::I
    nOutputPlane::I
    kH::I
    kW::I
    dH::I
    dW::I
    weight::Array{T}
    bias::Array{T}
    grad_weight::Array{T}
    grad_bias::Array{T}
    output::Array{T}
    gradInput::Array{T}
end


text:
"""
    instantiate(gl::GitLink; verbose = false)

Create the file structure and clone the `GitLink` locally.
Returns `true` if the action was successful.
"""
function instantiate(gl::GitLink; verbose = true)

    # Sync GitLink
    return sync_link(gl; verbose, force = true)
end


text:
module Dash
import HTTP, JSON2, CodecZlib, MD5
using Sockets
using MacroTools
const ROOT_PATH = realpath(joinpath(@__DIR__, ".."))
include("Components.jl")
include("Front.jl")
import .Front
using .Components

export dash, Component, Front, <|, @callid_str, CallbackId, callback!,
enable_dev_tools!, ClientsideFunction,
run_server, PreventUpdate, no_update, @var_str

#ComponentPackages.@reg_components()
include("env.jl")
include("utils.jl")
include("app.jl")
include("resources/registry.jl")
include("resources/application.jl")
include("handlers.jl")

@doc """
    module Dash

Julia backend for [Plotly Dash](https://github.com/plotly/dash)

# Examples
```julia

using Dash
app = dash(external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"]) do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        dcc_graph(id="graph",
            figure = (
                data = [(x = [1,2,3], y = [3,2,8], type="bar")],
                layout = Dict(:title => "Graph")
            )
        )

    end
end
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children") do type, value
    "You've entered: '\$(value)' into a '\$(type)' input control"
end
callback!(app, callid"graphTitle.value => graph.figure") do value
    (
        data = [
            (x = [1,2,3], y = abs.(randn(3)), type="bar"),
            (x = [1,2,3], y = abs.(randn(3)), type="scatter", mode = "lines+markers", line = (width = 4,))
        ],
        layout = (title = value,)
    )
end
handle = make_handler(app, debug = true)
run_server(handle, HTTP.Sockets.localhost, 8050)
```

""" Dashboards

"""
    run_server(app::DashApp, host = HTTP.Sockets.localhost, port = 8050; debug::Bool = false)

Run Dash server

#Arguments
- `app` - Dash application
- `host` - host
- `port` - port
- `debug::Bool = false` - Enable/disable all the dev tools

#Examples
```jldoctest
julia> app = dash() do
    html_div() do
        html_h1("Test Dashboard")
    end
end
julia>
julia> run_server(handler,  HTTP.Sockets.localhost, 8050)
```

"""
function run_server(app::DashApp, host = HTTP.Sockets.localhost, port = 8050;
            debug = nothing,
            dev_tools_ui = nothing,
            dev_tools_props_check = nothing,
            dev_tools_serve_dev_bundles = nothing,
            dev_tools_hot_reload = nothing,
            dev_tools_hot_reload_interval = nothing,
            dev_tools_hot_reload_watch_interval = nothing,
            dev_tools_hot_reload_max_retry = nothing,
            dev_tools_silence_routes_logging = nothing,
            dev_tools_prune_errors = nothing
            )
    @env_default!(debug, Bool, false)
    enable_dev_tools!(app,
        debug = debug,
        dev_tools_ui = dev_tools_ui,
        dev_tools_props_check = dev_tools_props_check,
        dev_tools_serve_dev_bundles = dev_tools_serve_dev_bundles,
        dev_tools_hot_reload = dev_tools_hot_reload,
        dev_tools_hot_reload_interval = dev_tools_hot_reload_interval,
        dev_tools_hot_reload_watch_interval = dev_tools_hot_reload_watch_interval,
        dev_tools_hot_reload_max_retry = dev_tools_hot_reload_max_retry,
        dev_tools_silence_routes_logging = dev_tools_silence_routes_logging,
        dev_tools_prune_errors = dev_tools_prune_errors
    )
    main_func = () -> begin
        ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)
        handler = make_handler(app);
        try
            task = @async HTTP.serve(handler, host, port)
            @info string("Running on http://", host, ":", port)
            wait(task)
        catch e
            if e isa InterruptException
                @info "exited"
            else
                rethrow(e)
            end

        end
    end
    start_server = () -> begin
        handler = make_handler(app);
        server = Sockets.listen(get_inetaddr(host, port))
        task = @async HTTP.serve(handler, host, port; server = server)
        @info string("Running on http://", host, ":", port)
        return (server, task)
    end

    if get_devsetting(app, :hot_reload) && !is_hot_restart_available()
        @warn "Hot reloading is disabled for interactive sessions. Please run your app using julia from the command line to take advantage of this feature."
    end

    if get_devsetting(app, :hot_reload) && is_hot_restart_available()
        hot_restart(start_server, check_interval = get_devsetting(app, :hot_reload_watch_interval))
    else
        (server, task) = start_server()
        try
            wait(task)
            println(task)
        catch e
            close(server)
            if e isa InterruptException
                println("finished")
                return
            else
                rethrow(e)
            end
        end
    end
end
get_inetaddr(host::String, port::Integer) = Sockets.InetAddr(parse(IPAddr, host), port)
get_inetaddr(host::IPAddr, port::Integer) = Sockets.InetAddr(host, port)

end # module


text:
"""
    prune_tree(root::T, node_names::Vector{String})::T where T<:AbstractNode

This function removes specific nodes, including their descendants, from a tree.

* `root` : root Node of tree to prune.

* `node_names` : vector of strings, used to specify nodes to remove.
"""
function prune_tree(root::T, node_names::Vector{String})::T where T<:AbstractNode
    # copy the tree and call the inplace version of the function on the copy
    if root.name in node_names
        throw(ArgumentError("trying to prune root, please set root to nothing instead"))
    end
    copyroot = deepcopy(root)
    prune_tree!(copyroot, node_names)
    return copyroot
end

"""
    prune_tree!(root::T, node_names::Vector{String})::Nothing where T<:AbstractNode

In-place version of prune_tree.

* `root` : root Node of tree to prune.

* `node_names` : vector of strings, used to specify nodes to remove.
"""
function prune_tree!(root::T, node_names::Vector{String})::Nothing where T<:AbstractNode
    nodes = post_order(root)
    names = [node.name for node in nodes]
    nodes_to_prune = Vector{Node}()
    for name in node_names
        indeces = findall(x->x == name, names)
        if length(indeces) > 1
            throw(ArgumentError("Multiple nodes are named \"$name\""))
        elseif length(indeces) == 0
            print("Warning: No node named \"$name\" in tree")
        else
            push!(nodes_to_prune, nodes[indeces[1]])
        end # if/else
    end # for
    if length(nodes_to_prune) == 0
        throw(ArgumentError("None of the node names correspond to a node in the tree"))
    end # if
    prune_tree!(root, nodes_to_prune)
end

"""
    prune_tree!(root::T, node_names::Vector{T})::Nothing where T<:AbstractNode

In-place version of prune_tree.

* `root` : root node of tree to prune.

* `node_names`: vector of Node objects to be removed from tree.
"""
function prune_tree!(root::T, nodes::Vector{T})::Nothing where T<:AbstractNode
    if root in nodes
        throw(ArgumentError("trying to prune root, please set root to nothing instead"))
    else
        for node in nodes
            remove_child!(node.mother, node)
        end
    end
end