text:
export sgd_parameter_learning

using LogicCircuits: num_nodes
using StatsFuns: logaddexp
using CUDA

function sgd_parameter_learning(pc::ProbCircuit, data; lr::Float64 = 0.01,
                                use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data),
                                reuse_values = nothing, reuse_flows = nothing,
                                reuse = (nothing, nothing))
    # Construct the low-level circuit representation
    pbc = ParamBitCircuit(pc, data)
    if use_gpu
        pbc = to_gpu(pbc)
    end

    # Main training loop
    @assert isbatched(data)
    for idx = 1 : length(data)
        sgd_parameter_learning(pbc, data[idx]; lr, use_sample_weights, use_gpu,
                               reuse_values, reuse_flows, reuse)
    end

    # Store the updated parameters back to `pc`
    estimate_parameters_cached!(pc, pbc)

    pbc.params
end
function sgd_parameter_learning(pbc::ParamBitCircuit, data; lr::Float64 = 0.01,
                                use_sample_weights::Bool = true, use_gpu::Bool = isgpu(data),
                                reuse_values = nothing, reuse_flows = nothing,
                                reuse = (nothing, nothing))
    # Extract weights from dataset
    if isweighted(data)
        # `data' is weighted according to its `weight' column
        data, weights = split_sample_weights(data)
    else
        use_sample_weights = false
    end

    # Move data to GPU if needed
    if isgpu(data)
        use_gpu = true
    elseif use_gpu && !isgpu(data)
        data = to_gpu(data)
    end

    if use_gpu && !isgpu(pbc)
        pbc.bitcircuit = to_gpu(pbc.bitcircuit)
        pbc.params = to_gpu(pbc.params)
    end

    if use_gpu
        if !isgpu(pbc)
            pbc.bitcircuit = to_gpu(pbc.bitcircuit)
            pbc.params = to_gpu(pbc.params)
        end

        if use_sample_weights
            if !isgpu(weights)
                weights = to_gpu(weights)
            end
            sgd_parameter_learning_gpu(pbc, data; lr, weights, reuse_values, reuse_flows, reuse)
        else
            sgd_parameter_learning_gpu(pbc, data; lr, reuse_values, reuse_flows, reuse)
        end
    else
        if use_sample_weights
            sgd_parameter_learning_cpu(pbc, data; lr, weights, reuse_values, reuse_flows, reuse)
        else
            sgd_parameter_learning_cpu(pbc, data; lr, reuse_values, reuse_flows, reuse)
        end
    end

    pbc.params
end

function sgd_parameter_learning_cpu(pbc::ParamBitCircuit, data; lr::Float64 = 0.01,
                                    weights = nothing, reuse_values = nothing, reuse_flows = nothing,
                                    reuse = (nothing, nothing))
    bc = pbc.bitcircuit
    params = pbc.params

    # Forward/upward pass
    values = marginal_all(pbc, data, reuse_values)

    # Compute the gradient (this part can be customized)
    log_grads = similar!(reuse[1], typeof(params), num_examples(data))
    @inbounds @views log_grads[:] .= zero(eltype(values)) .- values[:, num_nodes(pbc.bitcircuit)]

    param_grads = similar!(reuse[2], typeof(params), size(params)...)

    # Backward/downward pass
    flows = backprop_flows_down_cpu(pbc, log_grads, values, param_grads, reuse_flows; weights)

    # Apply gradients
    apply_gradients_cpu(pbc, param_grads; lr)

    pbc.params, values, flows, (log_grads, param_grads)
end

function sgd_parameter_learning_gpu(pbc::ParamBitCircuit, data; lr::Float64 = 0.01,
                                    weights = nothing, reuse_values = nothing, reuse_flows = nothing,
                                    reuse = (nothing, nothing))
    bc = pbc.bitcircuit
    params = pbc.params

    # Forward/upward pass
    values = marginal_all(pbc, data, reuse_values)

    # Compute the gradient (this part can be customized)
    log_grads = similar!(reuse[1], typeof(params), num_examples(data))
    @inbounds @views log_grads[:] .= zero(eltype(values)) .- values[:, num_nodes(pbc.bitcircuit)]

    param_grads = similar!(reuse[2], typeof(params), size(params)...)

    # Backward/downward pass
    flows = backprop_flows_down_gpu(pbc, log_grads, values, param_grads, reuse_flows; weights)

    # Apply gradients
    apply_gradients_gpu(pbc, param_grads; lr)

    pbc.params, values, flows, (log_grads, param_grads)
end

######################################
# Bit circuit backprop downward pass #
######################################

function backprop_flows_down_cpu(pbc::ParamBitCircuit, log_grads, values, param_grads, reuse = nothing; weights = nothing)
    flows = similar!(reuse, typeof(values), size(values)...)

    if weights !== nothing
        @inbounds @views log_grads .+= log.(weights)
    end

    backprop_flows_down_layers_cpu(pbc, log_grads, flows, values, param_grads)

    flows # Return for reuse
end

function backprop_flows_down_layers_cpu(pbc::ParamBitCircuit, log_grads::Vector, flows::Matrix, values::Matrix,
                                    param_grads::Vector)
    @assert flows !== values
    bc::BitCircuit = pbc.bitcircuit
    els = bc.elements
    for layer in Iterators.reverse(bc.layers)
        ## Compute gradients
        Threads.@threads for dec_id in layer
            par_start = @inbounds bc.nodes[3, dec_id]
            if iszero(par_start)
                if dec_id == num_nodes(bc)
                    # Assign (log) gradient to the root node
                    @inbounds @views flows[:, dec_id] .= log_grads
                end
                # no parents, ignore (can happen for false/true node and root)
            else
                par_end = @inbounds bc.nodes[4, dec_id]
                for j = par_start : par_end
                    par = @inbounds bc.parents[j] # Parent #j of `dec_id`
                    grandpa = @inbounds els[1, par] # Parent of `par`
                    sib_id = sibling(els, par, dec_id) # Sibling of `dec_id` w.r.t. `par`

                    single_child = has_single_child(bc.nodes, grandpa)
                    if single_child
                        if j == par_start
                            @inbounds @views @. flows[:, dec_id] = flows[:, grandpa] + values[:, sib_id]
                        else
                            @inbounds @simd for i = 1 : size(flows, 1)
                                flows[i, dec_id] = logaddexp(flows[i, dec_id], flows[i, grandpa] + values[i, sib_id])
                            end
                        end
                    else
                        θ = eltype(flows)(pbc.params[par])
                        if j == par_start
                            @inbounds @views @. flows[:, dec_id] = flows[:, grandpa] + values[:, sib_id] + θ
                        else
                            @inbounds @simd for i = 1 : size(flows, 1)
                                flows[i, dec_id] = logaddexp(flows[i, dec_id], flows[i, grandpa] + values[i, sib_id] + θ)
                            end
                        end
                    end
                    # Compute gradient only once
                    if sib_id > dec_id
                        param_grads[par] = sum(exp.(flows[:, grandpa] + values[:, dec_id] + values[:, sib_id]))
                    end
                end
            end
        end
    end
end

function apply_gradients_cpu(pbc::ParamBitCircuit, param_grads::Vector; lr::Float64 = 0.01)
    bc::BitCircuit = pbc.bitcircuit
    nodes = bc.nodes
    params = pbc.params
    for layer in Iterators.reverse(bc.layers)
        ## Apply gradients
        Threads.@threads for dec_id in layer
            if !has_single_child(nodes, dec_id)
                ele_start_id = nodes[1, dec_id]
                ele_end_id = nodes[2, dec_id]

                sum_grads = 0.0
                @inbounds for ele_id = ele_start_id : ele_end_id
                    sum_grads += param_grads[ele_id]
                end
                sum_params = -Inf
                @inbounds for ele_id = ele_start_id : ele_end_id
                    params[ele_id] += lr * param_grads[ele_id] / (sum_grads + 1e-8)
                    sum_params = logaddexp(sum_params, params[ele_id])
                end
                @inbounds for ele_id = ele_start_id : ele_end_id
                    params[ele_id] -= sum_params
                end
            end
        end
    end
end

function backprop_flows_down_gpu(pbc::ParamBitCircuit, log_grads::CuVector, values::CuMatrix,
                                 param_grads::CuVector, reuse =  nothing; dec_per_thread = 8,
                                 log2_threads_per_block = 7, weights = nothing)
    flows = similar!(reuse, typeof(values), size(values)...)
    if weights !== nothing
        @inbounds @views log_grads .+= CUDA.log.(weights)
    end

    bc::BitCircuit = pbc.bitcircuit
    @inbounds @views param_grads .= zero(eltype(param_grads))
    param_grads_device = CUDA.cudaconvert(param_grads)

    CUDA.@sync for layer in Iterators.reverse(bc.layers)
        num_examples = size(values, 1)
        num_decision_sets = length(layer)/dec_per_thread
        num_threads =  balance_threads(num_examples, num_decision_sets, log2_threads_per_block)
        num_blocks = (ceil(Int, num_examples/num_threads[1]),
                      ceil(Int, num_decision_sets/num_threads[2]))
        @cuda threads=num_threads blocks=num_blocks backprop_flows_down_layers_cuda(layer, bc.nodes, bc.elements, bc.parents, pbc.params, log_grads, flows, values, param_grads_device)
    end

    flows # Return for reuse
end

function backprop_flows_down_layers_cuda(layer, nodes, elements, parents, params, log_grads, flows, values, param_grads)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y

    for k = index_x : stride_x : size(values, 1)
        for i = index_y : stride_y : length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes, 2)
                # Assign (log) gradient to the root node
                flow = log_grads[k]
            else
                par_start = @inbounds nodes[3, dec_id]
                flow = typemin(eltype(flows)) # log(0)
                if !iszero(par_start)
                    par_end = @inbounds nodes[4, dec_id]
                    for j = par_start : par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1, par]
                        sib_id = sibling(elements, par, dec_id)

                        g_flow = @inbounds flows[k, grandpa]
                        d_value = @inbounds values[k, dec_id]
                        s_value = @inbounds values[k, sib_id]
                        if has_single_child(nodes, grandpa)
                            edge_flow = g_flow + s_value
                        else
                            θ = eltype(flows)(params[par])
                            edge_flow = g_flow + s_value + θ
                        end
                        flow = logsumexp_cuda(flow, edge_flow)

                        # Compute gradient only once
                        if sib_id > dec_id
                            grad::Float64 = CUDA.exp(g_flow + d_value + s_value)
                            CUDA.@atomic param_grads[par] += grad
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
        end
    end
    return nothing
end

function apply_gradients_gpu(pbc::ParamBitCircuit, param_grads::CuVector; lr::Float64 = 0.01)
    bc::BitCircuit = pbc.bitcircuit

    CUDA.@sync for layer in Iterators.reverse(bc.layers)
        num_threads = 2^min(ceil(Int, 2.0 * log2(length(layer))), 8)
        num_blocks = 2^ceil(Int, log2(length(layer)^2 / num_threads))
        @cuda threads=num_threads blocks=num_blocks apply_gradients_cuda(layer, bc.nodes, param_grads, pbc.params,
                                                                         lr::Float64)
    end
end

function apply_gradients_cuda(layer, nodes, param_grads, params, lr::Float64)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride_x = blockDim().x * gridDim().x

    for i = index_x : stride_x : length(layer)
        dec_id = @inbounds layer[i]
        single_child = has_single_child(nodes, dec_id)
        if !single_child
            ele_start_id = nodes[1, dec_id]
            ele_end_id = nodes[2, dec_id]

            sum_grads = zero(eltype(param_grads))
            @inbounds for ele_id = ele_start_id : ele_end_id
                sum_grads += param_grads[ele_id]
            end
            sum_params = -Inf
            @inbounds for ele_id = ele_start_id : ele_end_id
                params[ele_id] += lr * param_grads[ele_id] / (sum_grads + 1e-8)
                sum_params = logsumexp_cuda(sum_params, params[ele_id])
            end
            @inbounds for ele_id = ele_start_id : ele_end_id
                params[ele_id] -= sum_params
            end
        end
    end
end


text:
using Test
using Distributions
using Random

using EnsembleKalmanProcesses.DataStorage

@testset "DataStorage" begin
    seed=2021
    Random.seed!(seed)

    parameter_samples = rand(MvNormal(2,0.1),10) #10 samples of 4D params
    data_samples = rand(MvNormal(12,2),10) #10 samples of 12D data
    data_samples_short = data_samples[:,1:end-1]

    new_parameter_samples = rand(MvNormal(2,0.1),10) #10 samples of 4D params
    new_data_samples = rand(MvNormal(12,2),10) #10 samples of 12D data
    new_data_samples_short = new_data_samples[:,1:end-1]

    #test DataContainer
    parameter_samples_T = permutedims(parameter_samples, (2,1))
    idata = DataContainer(parameter_samples)
    idata_T = DataContainer(parameter_samples_T, data_are_columns=false)

    @test get_data(idata) == parameter_samples
    @test get_data(idata_T) == parameter_samples
    @test size(idata) == size(parameter_samples)

    #test Paired Container
    idata = DataContainer(parameter_samples)
    odata = DataContainer(data_samples)
    odata_short = DataContainer(data_samples_short)

    @test_throws DimensionMismatch PairedDataContainer(parameter_samples, data_samples_short)
    iopairs = PairedDataContainer(parameter_samples,data_samples)

    @test_throws DimensionMismatch PairedDataContainer(idata,odata_short)
    iopairs2 = PairedDataContainer(idata,odata)

    @test get_data(iopairs) == get_data(iopairs2)
    @test get_inputs(iopairs) == parameter_samples
    @test get_outputs(iopairs) == data_samples

end

ext: jl


text:
# Conversions between Julia and Python types for the PyCall module.

#########################################################################
# Conversions of simple types (numbers and nothing)

# conversions from Julia types to PyObject:

PyObject(i::Unsigned) = PyObject(@pycheckn ccall(pyint_from_size_t::Ptr{Void},
                                                 PyPtr, (Uint,), i))
PyObject(i::Integer) = PyObject(@pycheckn ccall(pyint_from_ssize_t::Ptr{Void},
                                                PyPtr, (Int,), i))

PyObject(b::Bool) = PyObject(@pycheckn ccall((@pysym :PyBool_FromLong),
                                             PyPtr, (Clong,), b))

PyObject(r::Real) = PyObject(@pycheckn ccall((@pysym :PyFloat_FromDouble),
                                             PyPtr, (Cdouble,), r))

PyObject(c::Complex) = PyObject(@pycheckn ccall((@pysym :PyComplex_FromDoubles),
                                                PyPtr, (Cdouble,Cdouble),
                                                real(c), imag(c)))

PyObject(n::Nothing) = begin @pyinitialize; pyerr_check("PyObject(nothing)", pyincref(pynothing.o)); end

# conversions to Julia types from PyObject

convert{T<:Integer}(::Type{T}, po::PyObject) =
  convert(T, @pycheck ccall(pyint_as_ssize_t::Ptr{Void}, Int, (PyPtr,), po))

convert(::Type{Bool}, po::PyObject) =
  convert(Bool, @pycheck ccall(pyint_as_ssize_t::Ptr{Void}, Int, (PyPtr,), po))

convert{T<:Real}(::Type{T}, po::PyObject) =
  convert(T, @pycheck ccall((@pysym :PyFloat_AsDouble), Cdouble, (PyPtr,), po))

convert{T<:Complex}(::Type{T}, po::PyObject) =
  convert(T,
    begin
        re = @pycheck ccall((@pysym :PyComplex_RealAsDouble),
                            Cdouble, (PyPtr,), po)
        complex128(re, ccall((@pysym :PyComplex_ImagAsDouble),
                             Cdouble, (PyPtr,), po))
    end)

convert(::Type{Nothing}, po::PyObject) = nothing

#########################################################################
# String conversions (both bytes arrays and unicode strings)

PyObject(s::UTF8String) =
  PyObject(@pycheckn ccall(PyUnicode_DecodeUTF8::Ptr{Void},
                           PyPtr, (Ptr{Uint8}, Int, Ptr{Uint8}),
                           bytestring(s), length(s.data), C_NULL))

function PyObject(s::String)
    @pyinitialize
    if pyunicode_literals::Bool
        sb = bytestring(s)
        PyObject(@pycheckn ccall(PyUnicode_DecodeUTF8::Ptr{Void},
                                 PyPtr, (Ptr{Uint8}, Int, Ptr{Uint8}),
                                 sb, length(sb.data), C_NULL))
    else
        PyObject(@pycheckn ccall(pystring_fromstring::Ptr{Void},
                                 PyPtr, (Ptr{Uint8},), bytestring(s)))
    end
end

function convert{T<:String}(::Type{T}, po::PyObject)
    @pyinitialize
    if pyisinstance(po, @pysym :PyUnicode_Type)
        convert(T, PyObject(@pycheckni ccall(PyUnicode_AsUTF8String::Ptr{Void},
                                             PyPtr, (PyPtr,), po)))
    else
        convert(T, bytestring(@pycheckni ccall(pystring_asstring::Ptr{Void},
                                               Ptr{Uint8}, (PyPtr,), po)))
    end
end

# TODO: should symbols be converted to a subclass of Python strings/bytes,
#       so that PyAny conversion can convert it back to a Julia symbol?
PyObject(s::Symbol) = PyObject(string(s))
convert(::Type{Symbol}, po::PyObject) = symbol(convert(String, po))

#########################################################################
# Pointer conversions, using ctypes, PyCObject, or PyCapsule

PyObject(p::Ptr) = begin @pyinitialize; (py_void_p::Function)(p); end

function convert(::Type{Ptr{Void}}, po::PyObject)
    @pyinitialize
    if pyisinstance(po, c_void_p_Type::PyObject)
        v = po["value"]
        # ctypes stores the NULL pointer specially, grrr
        pynothing_query(v) == Nothing ? C_NULL :
          convert(Ptr{Void}, convert(Uint, po["value"]))
    elseif pyisinstance(po, PyCObject_Type::Ptr{Void})
        @pychecki ccall((@pysym :PyCObject_AsVoidPtr), Ptr{Void}, (PyPtr,), po)
    elseif pyisinstance(po, PyCapsule_Type::Ptr{Void})
        @pychecki ccall((@pysym :PyCapsule_GetPointer),
                        Ptr{Void}, (PyPtr,Ptr{Uint8}),
                        po, ccall((@pysym :PyCapsule_GetName),
                                  Ptr{Uint8}, (PyPtr,), po))
    else
        convert(Ptr{Void}, convert(Uint, po))
    end
end

pyptr_query(po::PyObject) = pyisinstance(po, c_void_p_Type::PyObject) || pyisinstance(po, PyCObject_Type::Ptr{Void}) || pyisinstance(po, PyCapsule_Type::Ptr{Void}) ? Ptr{Void} : None

#########################################################################
# for automatic conversions, I pass Vector{PyAny}, NTuple{PyAny}, etc.,
# but since PyAny is an abstract type I need to convert this to Any
# before actually creating the Julia object

# I want to use a union, but this seems to confuse Julia's method
# dispatch for the convert function in some circumstances
# typealias PyAny Union(PyObject, Int, Bool, Float64, Complex128, String, Function, Dict, Tuple, Array)
abstract PyAny

# I originally implemented this via multiple dispatch, with
#   pyany_toany(::Type{PyAny}), pyany_toany(x::Tuple), and pyany_toany(x),
# but Julia seemed to get easily confused about which one to call.
pyany_toany(x) = isa(x, Type{PyAny}) ? Any : (isa(x, Tuple) ?
                                              map(pyany_toany, x) : x)

# no-op conversions
for T in (:PyObject, :Int, :Bool, :Float64, :Complex128, :String,
          :Function, :Dict, :Tuple, :Array)
    @eval convert(::Type{PyAny}, x::$T) = x
end

#########################################################################
# Function conversion (see callback.jl for conversion the other way)

convert(::Type{Function}, po::PyObject) =
    (args...) -> pycall(po, PyAny, args...)

#########################################################################
# Tuple conversion

function PyObject(t::Tuple)
    o = PyObject(@pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,),
                                 length(t)))
    for i = 1:length(t)
        oi = PyObject(t[i])
        @pycheckzi ccall((@pysym :PyTuple_SetItem), Cint, (PyPtr,Int,PyPtr),
                         o, i-1, oi)
        pyincref(oi) # PyTuple_SetItem steals the reference
    end
    return o
end

function convert(tt::NTuple{Type}, o::PyObject)
    len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
    if len != length(tt)
        throw(BoundsError())
    end
    ntuple(len, i ->
           convert(tt[i], PyObject(ccall((@pysym :PySequence_GetItem), PyPtr,
                                         (PyPtr, Int), o, i-1))))
end

#########################################################################
# PyVector: no-copy wrapping of a Julia object around a Python sequence

type PyVector{T} <: AbstractVector{T}
    o::PyObject
    function PyVector(o::PyObject)
        if o.o == C_NULL
            throw(ArgumentError("cannot make PyVector from NULL PyObject"))
        elseif pysequence_query(o) == None
            throw(ArgumentError("only List and Sequence objects can be converted to PyVector"))
        end
        new(o)
    end
end

PyVector(o::PyObject) = PyVector{PyAny}(o)
convert(::Type{PyVector}, o::PyObject) = PyVector(o)
convert{T}(::Type{PyVector{T}}, o::PyObject) = PyVector{T}(o)
convert(::Type{PyPtr}, a::PyVector) = a.o.o

# when a PyVector is copied it is converted into an ordinary Julia Vector
similar(a::PyVector, T, dims::Dims) = Array(T, dims)
similar{T}(a::PyVector{T}) = similar(a, pyany_toany(T), size(a))
similar{T}(a::PyVector{T}, dims::Dims) = similar(a, pyany_toany(T), dims)
similar{T}(a::PyVector{T}, dims::Int...) = similar(a, pyany_toany(T), dims)
eltype{T}(::PyVector{T}) = pyany_toany(T)
eltype{T}(::Type{PyVector{T}}) = pyany_toany(T)

size(a::PyVector) = ((@pycheckzi ccall((@pysym :PySequence_Size), Int, (PyPtr,), a)),)

ref(a::PyVector) = ref(a, 1)
ref{T}(a::PyVector{T}, i::Integer) = convert(T, PyObject(@pycheckni ccall((@pysym :PySequence_GetItem), PyPtr, (PyPtr, Int), a, i-1)))

assign(a::PyVector, v) = assign(a, v, 1)
assign(a::PyVector, v, i::Integer) = @pycheckzi ccall((@pysym :PySequence_SetItem), Cint, (PyPtr, Int, PyPtr), a, i-1, PyObject(v))

function delete!(a::PyVector, i::Integer)
    v = a[i]
    @pycheckzi ccall((@pysym :PySequence_DelItem), Cint, (PyPtr, Int), a, i-1)
    v
end

pop!(a::PyVector) = delete!(a, length(a))

summary{T}(a::PyVector{T}) = string(Base.dims2string(size(a)), " ",
                                   string(pyany_toany(T)), " PyVector")

#########################################################################
# Lists and 1d arrays.

# recursive conversion of A to a list of list of lists... starting
# with dimension dim and index i in A.
function array2py{T, N}(A::AbstractArray{T, N}, dim::Integer, i::Integer)
    if dim > N
        return PyObject(A[i])
    elseif dim == N # special case last dim to coarsen recursion leaves
        len = size(A, dim)
        s = N == 1 ? 1 : stride(A, dim)
        o = PyObject(@pycheckn ccall((@pysym :PyList_New), PyPtr, (Int,), len))
        for j = 0:len-1
            oi = PyObject(A[i+j*s])
            @pycheckzi ccall((@pysym :PyList_SetItem), Cint, (PyPtr,Int,PyPtr),
                             o, j, oi)
            pyincref(oi) # PyList_SetItem steals the reference
        end
        return o
    else # dim < N: store multidimensional array as list of lists
        len = size(A, dim)
        s = stride(A, dim)
        o = PyObject(@pycheckn ccall((@pysym :PyList_New), PyPtr, (Int,), len))
        for j = 0:len-1
            oi = array2py(A, dim+1, i+j*s)
            @pycheckzi ccall((@pysym :PyList_SetItem), Cint, (PyPtr,Int,PyPtr),
                             o, j, oi)
            pyincref(oi) # PyList_SetItem steals the reference
        end
        return o
    end
end

array2py(A::AbstractArray) = array2py(A, 1, 1)

PyObject(A::AbstractArray) =
   ndims(A) <= 1 || method_exists(stride,(typeof(A),Int)) ? array2py(A) :
   pyjlwrap_new(A)

function py2array{TA,N}(T, A::Array{TA,N}, o::PyObject,
                        dim::Integer, i::Integer)
    if dim > N
        A[i] = convert(T, o)
        return A
    elseif dim == N
        len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
        if len != size(A, dim)
            error("dimension mismatch in py2array")
        end
        s = stride(A, dim)
        for j = 0:len-1
            A[i+j*s] = convert(T, PyObject(ccall((@pysym :PySequence_GetItem),
                                                 PyPtr, (PyPtr, Int), o, j)))
        end
        return A
    else # dim < N: recursively extract list of lists into A
        len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
        if len != size(A, dim)
            error("dimension mismatch in py2array")
        end
        s = stride(A, dim)
        for j = 0:len-1
            py2array(T, A, PyObject(ccall((@pysym :PySequence_GetItem),
                                       PyPtr, (PyPtr, Int), o, j)),
                     dim+1, i+j*s)
        end
        return A
    end
end

# figure out if we can treat o as a multidimensional array, and return
# the dimensions
function pyarray_dims(o::PyObject)
    len = ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
    if len < 0 || # not a sequence
       len+1 < 0 || # object pretending to be a sequence of infinite length
       pystring_query(o) != None ||
       pyisinstance(o, @pysym :PyTuple_Type)
        pyerr_clear()
        return ()
    elseif len == 0
        return (0,)
    end
    dims0 = pyarray_dims(PyObject(ccall((@pysym :PySequence_GetItem),
                                        PyPtr, (PyPtr, Int), o, 0)))
    if length(dims0) == 0 # not a nested sequence
        return (len,)
    end
    for j = 1:len-1
        dims = pyarray_dims(PyObject(ccall((@pysym :PySequence_GetItem),
                                           PyPtr, (PyPtr, Int), o, j)))
        if dims != dims0
            # elements don't have equal lengths, cannot
            # treat as multidimensional array
            return (len,)
        end
    end
    return tuple(len, dims0...)
end

function py2array(T, o::PyObject)
    dims = pyarray_dims(o)
    A = Array(pyany_toany(T), dims)
    py2array(T, A, o, 1, 1)
end

function convert{T}(::Type{Vector{T}}, o::PyObject)
    len = ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
    if len < 0 || # not a sequence
       len+1 < 0  # object pretending to be a sequence of infinite length
        throw(ArgumentError("expected Python sequence"))
    end
    py2array(T, Array(pyany_toany(T), len), o, 1, 1)
end

convert(::Type{Array}, o::PyObject) = py2array(PyAny, o)

# NumPy conversions (multidimensional arrays)
include("numpy.jl")

#########################################################################
# PyDict: no-copy wrapping of a Julia object around a Python dictionary

type PyDict{K,V} <: Associative{K,V}
    o::PyObject
    isdict::Bool # whether this is a Python Dict (vs. generic Mapping object)

    function PyDict(o::PyObject)
        if o.o == C_NULL
            throw(ArgumentError("cannot make PyDict from NULL PyObject"))
        elseif pydict_query(o) == None
            throw(ArgumentError("only Dict and Mapping objects can be converted to PyDict"))
        end
        new(o, pyisinstance(o, @pysym :PyDict_Type))
    end
    function PyDict()
        @pyinitialize
        new(PyObject(@pycheckni ccall((@pysym :PyDict_New), PyPtr, ())), true)
    end
end

PyDict(o::PyObject) = PyDict{PyAny,PyAny}(o)
PyDict() = PyDict{PyAny,PyAny}()
PyDict{K,V}(d::Associative{K,V}) = PyDict{K,V}(PyObject(d))
convert(::Type{PyDict}, o::PyObject) = PyDict(o)
convert{K,V}(::Type{PyDict{K,V}}, o::PyObject) = PyDict{K,V}(o)
convert(::Type{PyPtr}, d::PyDict) = d.o.o

has(d::PyDict, key) = 1 == ccall(d.isdict ? (@pysym :PyDict_Contains)
                                          : (@pysym :PyMapping_HasKey),
                                  Cint, (PyPtr, PyPtr), d, PyObject(key))

pyobject_call(d::PyDict, vec::String) = PyObject(@pycheckni ccall((@pysym :PyObject_CallMethod), PyPtr, (PyPtr,Ptr{Uint8},Ptr{Uint8}), d, bytestring(vec), C_NULL))

keys{T}(::Type{T}, d::PyDict) = convert(Vector{T}, d.isdict ? PyObject(@pycheckni ccall((@pysym :PyDict_Keys), PyPtr, (PyPtr,), d)) : pyobject_call(d, "keys"))

values{T}(::Type{T}, d::PyDict) = convert(Vector{T}, d.isdict ? PyObject(@pycheckni ccall((@pysym :PyDict_Values), PyPtr, (PyPtr,), d)) : pyobject_call(d, "values"))

similar{K,V}(d::PyDict{K,V}) = Dict{pyany_toany(K),pyany_toany(V)}()
eltype{K,V}(a::PyDict{K,V}) = (pyany_toany(K),pyany_toany(V))

function assign(d::PyDict, v, k)
    @pycheckzi ccall((@pysym :PyObject_SetItem), Cint, (PyPtr, PyPtr, PyPtr),
                     d, PyObject(k), PyObject(v))
    d
end

function get{K,V}(d::PyDict{K,V}, k, default)
    vo = ccall((@pysym :PyObject_GetItem), PyPtr, (PyPtr,PyPtr), d, PyObject(k))
    if vo == C_NULL
        pyerr_clear()
        return default
    else
        return convert(V, PyObject(vo))
    end
end

function delete!(d::PyDict, k)
    v = d[k]
    @pycheckzi ccall(d.isdict ? (@pysym :PyDict_DelItem)
                              : (@pysym :PyObject_DelItem),
                     Cint, (PyPtr, PyPtr), d, PyObject(k))
    return v
end

function delete!(d::PyDict, k, default)
    try
        return delete!(d, k)
    catch
        return default
    end
end

function empty!(d::PyDict)
    if d.isdict
        @pychecki ccall((@pysym :PyDict_Clear), Void, (PyPtr,), d)
    else
        # for generic Mapping items we must delete keys one by one
        for k in keys(d)
            delete!(d, k)
        end
    end
    return d
end

length(d::PyDict) = @pycheckz ccall(d.isdict ? (@pysym :PyDict_Size)
                                             : (@pysym :PyObject_Size),
                                    Int, (PyPtr,), d)
isempty(d::PyDict) = length(d) == 0

type PyDict_Iterator
    # arrays to pass key, value, and pos pointers to PyDict_Next
    ka::Array{PyPtr}
    va::Array{PyPtr}
    pa::Vector{Int}

    items::PyObject # items list, for generic Mapping objects

    i::Int # current position in items list (0-based)
    len::Int # length of items list
end

function start(d::PyDict)
    if d.isdict
        PyDict_Iterator(Array(PyPtr,1), Array(PyPtr,1), zeros(Int,1),
                        PyObject(), 0, length(d))
    else
        items = convert(Vector{PyObject}, pyobject_call(d, PyObject, "items"))
        PyDict_Iterator(Array(PyPtr,0), Array(PyPtr,0), zeros(Int,0),
                        items, 0,
                        @pycheckz ccall((@pysym :PySequence_Size),
                                        Int, (PyPtr,), items))
    end
end

done(d::PyDict, itr::PyDict_Iterator) = itr.i >= itr.len

function next{K,V}(d::PyDict{K,V}, itr::PyDict_Iterator)
    if itr.items.o == C_NULL
        # Dict object, use PyDict_Next
        if 0 == ccall((@pysym :PyDict_Next), Cint,
                      (PyPtr, Ptr{Int}, Ptr{PyPtr}, Ptr{PyPtr}),
                      d, itr.pa, itr.ka, itr.va)
            error("unexpected end of PyDict_Next")
        end
        ko = pyincref(itr.ka[1]) # PyDict_Next returns
        vo = pyincref(itr.va[1]) #   borrowed ref, so incref
        ((convert(K,ko), convert(V,vo)),
         PyDict_Iterator(itr.ka, itr.va, itr.pa, itr.items, itr.i+1, itr.len))
    else
        # generic Mapping object, use items list
        (convert((K,V), PyObject(@pycheckni ccall((@pysym :PySequence_GetItem),
                                                  PyPtr, (PyPtr,Int),
                                                  itr.items, itr.i))),
         PyDict_Iterator(itr.ka, itr.va, itr.pa, itr.items, itr.i+1, itr.len))
    end
end

function filter!(f::Function, d::PyDict)
    # We must use items(d) here rather than (k,v) in d,
    # because PyDict_Next does not permit changing the set of keys
    # during iteration.
    for (k,v) in items(d)
        if !f(k,v)
            delete!(d,k)
        end
    end
    return d
end

#########################################################################
# Dictionary conversions (copies)

function PyObject(d::Associative)
    o = PyObject(@pycheckn ccall((@pysym :PyDict_New), PyPtr, ()))
    for k in keys(d)
        @pycheckzi ccall((@pysym :PyDict_SetItem), Cint, (PyPtr,PyPtr,PyPtr),
                         o, PyObject(k), PyObject(d[k]))
    end
    return o
end

function convert{K,V}(::Type{Dict{K,V}}, o::PyObject)
    @pyinitialize
    copy(PyDict{K,V}(o))
end

#########################################################################
# Inferring Julia types at runtime from Python objects:
#
# [Note that we sometimes use the PyFoo_Check API and sometimes we use
#  PyObject_IsInstance(o, PyFoo_Type), since sometimes the former API
#  is a macro (hence inaccessible in Julia).]

# A type-query function f(o::PyObject) returns the Julia type
# for use with the convert function, or None if there isn't one.

pyint_query(o::PyObject) = pyisinstance(o, pyint_type::Ptr{Void}) ?
  (pyisinstance(o, @pysym :PyBool_Type) ? Bool : Int) : None

pyfloat_query(o::PyObject) = pyisinstance(o, @pysym :PyFloat_Type) ? Float64 : None

pycomplex_query(o::PyObject) =
  pyisinstance(o, @pysym :PyComplex_Type) ? Complex128 : None

pystring_query(o::PyObject) = pyisinstance(o, pystring_type::Ptr{Void}) ? String : pyisinstance(o, @pysym :PyUnicode_Type) ? UTF8String : None

pyfunction_query(o::PyObject) = pyisinstance(o, @pysym :PyFunction_Type) || pyisinstance(o, BuiltinFunctionType::PyObject) || pyisinstance(o, ufuncType::PyObject) || pyisinstance(o, TypeType::PyObject) || pyisinstance(o, MethodType::PyObject) || pyisinstance(o, MethodWrapperType::PyObject) ? Function : None

pynothing_query(o::PyObject) = o.o == pynothing.o ? Nothing : None

# we check for "items" attr since PyMapping_Check doesn't do this (it only
# checks for __getitem__) and PyMapping_Check returns true for some
# scipy scalar array members, grrr.
pydict_query(o::PyObject) = pyisinstance(o, @pysym :PyDict_Type) || (pyquery((@pysym :PyMapping_Check), o) && ccall((@pysym :PyObject_HasAttrString), Cint, (PyPtr,Array{Uint8}), o, "items") == 1) ? Dict{PyAny,PyAny} : None

function pysequence_query(o::PyObject)
    # pyquery(:PySequence_Check, o) always succeeds according to the docs,
    # but it seems we need to be careful; I've noticed that things like
    # scipy define "fake" sequence types with intmax lengths and other
    # problems
    if pyisinstance(o, @pysym :PyTuple_Type)
        len = @pycheckzi ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
        return ntuple(len, i ->
                      pytype_query(PyObject(ccall((@pysym :PySequence_GetItem),
                                                  PyPtr, (PyPtr,Int), o,i-1)),
                                   PyAny))
    else
        try
            otypestr = PyObject(@pycheckni ccall((@pysym :PyObject_GetItem), PyPtr, (PyPtr,PyPtr,), o["__array_interface__"], PyObject("typestr")))
            typestr = convert(String, otypestr)
            T = npy_typestrs[typestr[2:end]]
            return Array{T}
        catch
            # only handle PyList for now
            return pyisinstance(o, @pysym :PyList_Type) ? Array : None
        end
    end
end

macro return_not_None(ex)
    quote
        T = $ex
        if T != None
            return T
        end
    end
end

function pytype_query(o::PyObject, default::Type)
    @pyinitialize
    # Would be faster to have some kind of hash table here, but
    # that seems a bit tricky when we take subclasses into account
    @return_not_None pyint_query(o)
    @return_not_None pyfloat_query(o)
    @return_not_None pycomplex_query(o)
    @return_not_None pystring_query(o)
    @return_not_None pyfunction_query(o)
    @return_not_None pydict_query(o)
    @return_not_None pysequence_query(o)
    @return_not_None pyptr_query(o)
    @return_not_None pynothing_query(o)
    return default
end

pytype_query(o::PyObject) = pytype_query(o, PyObject)

function convert(::Type{PyAny}, o::PyObject)
    if (o.o == C_NULL)
        return o
    end
    try
        T = pytype_query(o)
        if T == PyObject && is_pyjlwrap(o)
            return unsafe_pyjlwrap_to_objref(o.o)
        end
        convert(pytype_query(o), o)
    catch
        pyerr_clear() # just in case
        o
    end
end

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


text:
# TODO: when building constraints, check whether they have already been built to avoid having them twice in the formulation.

function objective_edge_expression(::RoutingModel, edge_obj::EdgeWiseObjectiveFunction, ::Edge, ::Dict{Edge{Int}, Float64})
    error("Edge-wise objective function not implemented: $edge_obj")
end

function objective_edge_expression(rm::RoutingModel, edge_obj::EdgeWiseObjectiveFunction, e::Edge)
    if length(rm.data.traffic_matrix) == 0
        error("Creating an edge-wise objective term with an empty default traffic matrix. " *
              "Set one in the corresponding RoutingData object or use the `dm` argument.")
    end

    return objective_edge_expression(rm, edge_obj, e, rm.data.traffic_matrix)
end

function objective_edge_expression(rm::RoutingModel, ::Load, e::Edge, dm::Dict{Edge{Int}, Float64})
    return total_flow_in_edge(rm, e, dm) / capacity(rm, e)
end

function objective_edge_expression(rm::RoutingModel, k::KleinrockLoad, e::Edge,
                                   dm::Dict{Edge{Int}, Float64})
    # Model the function load / (1 - load).
    kleinrock = @variable(rm.model, base_name="kleinrock[$e]", lower_bound=0)

    load_e = objective_edge_expression(rm, Load(), e, dm)

    # The constraint is a rotated SOCP:
    #     load / (1 - load) ≤ kl    ⟺    (√2)² ≤ 2 (1 - load) (1 + kl)
    c = @constraint(rm.model,
                    [1 - load_e, 1 + kleinrock, sqrt(2)]
                    in RotatedSecondOrderCone())
    set_name(c, "kleinrock[$(e)]")

    if k.use_nonconvex_bilinear_formulation_for_equality
        # To get the equality, also model the reversed inequality (which is
        # nonconvex):
        #     kl ≤ load / (1 - load)    ⟺    kl - kl ⋅ load ≤ load
        c = @constraint(rm.model, kleinrock - load_e * kleinrock <= load_e)
        set_name(c, "kleinrock_bilin[$(e)]")
    end

    return kleinrock
end

function objective_edge_expression(rm::RoutingModel, ::FortzThorupLoad, e::Edge,
                                   dm::Dict{Edge{Int}, Float64})
    fortzthorup = @variable(rm.model, base_name="fortzthorup[$e]", lower_bound=0)
    load_e = objective_edge_expression(rm, Load(), e, dm)

    c1 = @constraint(rm.model, fortzthorup >= load_e)
    set_name(c1, "fortzthorup_1[$(e)]")

    c2 = @constraint(rm.model, fortzthorup >= 3 * load_e - 2 / 3)
    set_name(c2, "fortzthorup_2[$(e)]")

    c3 = @constraint(rm.model, fortzthorup >= 10 * load_e - 16 / 3)
    set_name(c3, "fortzthorup_3[$(e)]")

    c4 = @constraint(rm.model, fortzthorup >= 70 * load_e - 178 / 3)
    set_name(c4, "fortzthorup_4[$(e)]")

    c5 = @constraint(rm.model, fortzthorup >= 500 * load_e - 1486 / 3)
    set_name(c5, "fortzthorup_5[$(e)]")

    c6 = @constraint(rm.model, fortzthorup >= 5000 * load_e - 16318 / 3)
    set_name(c6, "fortzthorup_6[$(e)]")

    return fortzthorup
end

function objective_edge_expression(rm::RoutingModel, edge_obj::AlphaFairness, e::Edge, dm::Dict{Edge{Int}, Float64})
    # Model the function load^α / (1 - α). Depending on the value of α,
    # this constraint has very different models. They are all reformulation
    # of the inequality load^α / (1 - α) ≥ af.

    α = edge_obj.α
    @assert edge_obj.α >= 0.0

    alphafairness = @variable(rm.model, base_name="alphafairness[$e]")

    # Implement basic bounds on the value of α-fairness. For α = 1, both signs
    # are possible. Add bounds on the absolute value to improve numerics.
    LARGE_NUMBER = 20_000 # TODO: parameter?
    if α < 1.0
        @constraint(rm.model, alphafairness >= 0)
        @constraint(rm.model, alphafairness <= LARGE_NUMBER)
    elseif α > 1.0
        @constraint(rm.model, alphafairness >= - LARGE_NUMBER)
        @constraint(rm.model, alphafairness <= 0)
    else
        @constraint(rm.model, alphafairness >= - LARGE_NUMBER)
        @constraint(rm.model, alphafairness <= LARGE_NUMBER)
    end

    # Actual implementation of the objective functions.
    load_e = objective_edge_expression(rm, Load(), e, dm)

    c = if α == 0.0
        # Very simple case, actually linear (it can even be minimised):
        #     af == load
        @constraint(rm.model, load_e == alphafairness)
    elseif α == 0.5 && ! edge_obj.force_power_cone
        # This is a very specific case where the power cone is not really
        # needed, one second-order cone is sufficient. This cone is more
        # widely tolerated.
        #     af ≤ 2 √load
        #   ⟺ af² ≤ 4 load = 2 ⋅ 2 ⋅ load
        @constraint(rm.model, [2, load_e, alphafairness] in RotatedSecondOrderCone())
    elseif α == 1.0
        # The function becomes a logarithm:
        #     af ≤ log(load)
        #     1 × exp(af / 1) ≤ load
        #     (af, 1, l) ∈ EXP
        @constraint(rm.model, [alphafairness, 1, load_e] in MOI.ExponentialCone())
    elseif α == 1.5 && ! edge_obj.force_power_cone
        # Another very specific case:
        #     af ≤ -2 / √load
        #     af ≥ 2 / √load
        # Use the DCP to propose a reformulation of this:
        #     ϕ ≤ √load,    af ≥ 2 / ϕ
        #                   (af, ϕ, 2) ∈ RSOC
        phi = @variable(rm.model, lower_bound=0.0)
        set_name(phi, "alphafairness_phi[$(e)]")

        c = @constraint(rm.model, [1/2, load_e, phi] in RotatedSecondOrderCone())
        set_name(c, "alphafairness_phi[$(e)]")

        @constraint(rm.model, [- alphafairness, phi, 2] in RotatedSecondOrderCone())
    elseif α == 2.0 && ! edge_obj.force_power_cone
        # Another very specific case:
        #     af ≤ -1 / load
        # Use the DCP to propose a reformulation of this:
        #     ϕ ≥ 1 / load,    af <= - ϕ
        phi = @variable(rm.model, lower_bound=0.0)
        set_name(phi, "alphafairness_phi[$(e)]")

        c = @constraint(rm.model, [phi / sqrt(2), load_e / sqrt(2), 1] in RotatedSecondOrderCone())
        set_name(c, "alphafairness_phi[$(e)]")

        @constraint(rm.model, alphafairness <= - phi)
    elseif 0.0 < α < 1.0
        # Generic case: 0.0 < α < 1.0.
        #     af ≤ load^(1-α) / (1 - α),     with α < 1, i.e. 1 - α > 0
        #   ⟺ af ≤ load^(1-α) × 1^α / (1 - α)
        #   ⟺ (1 - α) × af ≤ load^(1-α) × 1^α
        # Caution: with α = 0, this is not a power cone!
        @constraint(rm.model, [load_e, 1, alphafairness * (1 - α)] in MOI.PowerCone(1.0 - α))
    else
        @assert α > 1.0
        # Last case: α > 1.0.
        #     af ≤ load^(1-α) / (1 - α),     with α > 1, i.e. 1-α < 0
        #   ⟺ 1 ≤ load^(1-α) × t^-1 / (1 - α)
        #   ⟺ 1 - α ≤ load^(1-α) × t^-1
        #   ⟺ (1 - α)^(-α)  ≤ load^(1-α)/(-α) × t^-1/(-α)
        #   ⟺ 1 / (1 - α)^α ≤ load^(α-1)/α × 1^(1/α)
        neg_af = @variable(rm.model, lower_bound=0.0)
        set_name(neg_af, "alphafairness_neg[$(e)]")

        c = @constraint(rm.model, alphafairness == -neg_af)
        set_name(c, "alphafairness_neg[$(e)]")

        β = α - 1 # > 0
        @constraint(rm.model, [neg_af, load_e, (1.0 / β) ^ (1.0 / (1.0 + β))] in MOI.PowerCone(1.0 / (1.0 + β)))
        # @constraint(rm.model, [1.0 / (α - 1)^α, load_e, neg_af] in MOI.PowerCone((α - 1.0) / α))
    end

    set_name(c, "alphafairness[$(e)]")
    return alphafairness
end


text:
type RNNLayer # a single layer of a multi-layer RNN
    wxh::NNMatrix # input (X or prev layer) to hidden weights & gradients
    whh::NNMatrix # hidden to hidden weights & gradients
    bhh::NNMatrix # bias of hidden to hidden
    function RNNLayer(prevsize::Int, hiddensize::Int, std::Float64)
        wxh = randNNMat(hiddensize, prevsize, std)
        whh = randNNMat(hiddensize, hiddensize, std)
        wbh = NNMatrix(hiddensize, 1, zeros(hiddensize,1), zeros(hiddensize,1))
        new(wxh, whh, wbh)
    end
end

type RNN <: Model
    hdlayers::Array{RNNLayer,1} # holds variable number of hidden layers
    whd::NNMatrix # hidden to decoder weights & gradients
    bd::NNMatrix  # bias of hidden to decoder layer
    matrices::Array{NNMatrix,1} # used by solver - holds references to each of matrices in model
    hiddensizes::Array{Int,1}
    function RNN(inputsize::Int, hiddensizes::Array{Int,1}, outputsize::Int, std::Float64=0.08)
        hdlayers = Array(RNNLayer, length(hiddensizes))
        matrices = Array(NNMatrix, 0)
        for d in 1:length(hiddensizes)
            prevsize = d == 1? inputsize : hiddensizes[d-1]
            layer = RNNLayer(prevsize, hiddensizes[d], std)
            hdlayers[d] = layer
            # store a reference to all the matrices in each
            push!(matrices, layer.wxh) # input (X or prev layer) to hidden weights & gradients
            push!(matrices, layer.whh) # hidden to hidden weights & gradients
            push!(matrices, layer.bhh) # bias of hidden to hidden
        end
        whd = randNNMat(outputsize, hiddensizes[end], std)
        bd = NNMatrix(outputsize, 1, zeros(outputsize,1), zeros(outputsize,1))

        # store a reference to these matrices
        push!(matrices, whd)
        push!(matrices, bd)
#         println("bh.w size($(size(bd.w)))")
        new(hdlayers, whd, bd, matrices, hiddensizes)
    end
end

function forwardprop(g::Graph, model::RNN, x, prev)

    # forward prop for a single tick of RNN
    # G is graph to append ops to
    # model contains RNN parameters
    # x is 1D column vector with observation
    # prev is a struct containing hidden activations and output from last step

    prevhd, _ = prev # extract previous hidden from the tuple
    hiddenprevs = Array(NNMatrix,0)
    if length(prevhd) == 0
        for hdsize in model.hiddensizes
            push!(hiddenprevs, NNMatrix(hdsize,1))
        end
    else
        hiddenprevs = prevhd
    end

    hidden = Array(NNMatrix,0)
    for d in 1:length(model.hiddensizes) # for each hidden layer

        input = d == 1 ? x : hidden[d-1]
        hdprev = hiddenprevs[d]
        wxh = model.hdlayers[d].wxh
        whh = model.hdlayers[d].whh
        bhh = model.hdlayers[d].bhh

#         println((d,size(wxh.w),(size(input.w))))
        h0 = mul(g, wxh, input)
        h1 = mul(g, whh, hdprev)
#         println((typeof(h0),typeof(h1),typeof(bhh)))
        hidden_d = relu(g, add(g, add(g, h0, h1), bhh))

        push!(hidden,hidden_d)
    end

    # one decoder to outputs at end
#     println((size(model.whd.w), size(hidden[end].w), size(model.bd.w)))
    output = add(g, mul(g, model.whd, hidden[end]),model.bd)

    # return hidden representation and output
    return hidden, output
end

max_issues_repo_name: Andy-P/RecurrentNN

max_forks_repo_name: Andy-P/RecurrentNN


text:
module TestCollectors

include("preamble.jl")

@testset "stable eltype" begin
    @testset for unsafe in [false, true]
        n = 10
        @testset "$label" for (label, src) in [
            "vector" => 1:n,
            "generator" => (x + 1 for x in 1:n),
            "filter" => (x for x in 1:n if isodd(x)),
        ]
            data = Vector{Int}(undef, n)
            c = collector(data, Val(unsafe))
            @test finish!(append!!(c, src)) === data
            @test data == collect(src)
        end
    end
end

@testset "unstable eltype" begin
    @testset for unsafe in [false, true]
        n = 10
        @testset "$label" for (label, src) in [
            "generator" => (x > 5 ? nothing : x for x in 1:n),
            "filter" => (x > 5 ? nothing : x for x in 1:n if isodd(x)),
        ]
            data = Vector{Int}(undef, n)
            c = collector(data, Val(unsafe))
            result = finish!(append!!(c, src))
            @test result !== data
            @test result == collect(src)
            @test data[1] == 1
        end
    end
end

@testset "combine" begin
    c1 = append!!(collector(), [1, 2])
    c2 = append!!(collector(), [3, 4, 5])
    @test finish!(append!!(c1, c2)) == 1:5

    @test finish!(append!!(collector(), collector())) ==ₜ Union{}[]
end

@testset "Empty" begin
    @test finish!(append!!(collector(), Union{}[])) ==ₜ Union{}[]
    @test finish!(append!!(collector(), ())) ==ₜ Union{}[]
end

end  # module


text:
"""
    VertexWeightedType

Whether a cycle/path constraint has a weight for vertices:

* the constraint has a fixed weight for each vertex: `FIXED_WEIGHT_VERTEX`
* the constraint has a variable weight for each vertex:
  `VARIABLE_WEIGHT_VERTEX` (typically, such a constraint has a dimension for
  each vertex to give its weight)
* the constraint does not have a notion of per-vertex weight:
  `UNWEIGHTED_VERTEX`

Typically, constraints with a type other than `UNWEIGHTED_VERTEX` have a
dimension for the total weight due to vertices.
"""
@enum VertexWeightedType begin
    FIXED_WEIGHT_VERTEX
    VARIABLE_WEIGHT_VERTEX
    UNWEIGHTED_VERTEX
end

"""
    EdgeWeightedType

Whether a cycle/path constraint has a weight for edges:

* the constraint has a fixed weight for each edge: `FIXED_WEIGHT_EDGE`
* the constraint has a variable weight for each edge:
  `VARIABLE_WEIGHT_EDGE` (typically, such a constraint has a dimension for
  each edge to give its weight)
* the constraint does not have a notion of per-edge weight:
  `UNWEIGHTED_EDGE`

Typically, constraints with a type other than `UNWEIGHTED_EDGE` have a
dimension for the total weight due to edges.
"""
@enum EdgeWeightedType begin
    FIXED_WEIGHT_EDGE
    VARIABLE_WEIGHT_EDGE
    UNWEIGHTED_EDGE
end

"""
    WalkType

The type of walk implemented by the constraint. In graph theory, a walk is a
sequence of edges that join vertices (also called nodes): a walk starts at some
vertex, then uses an edge to get to a second vertex, etc.

The possible values are:

* a trail (walk whose edges are distinct, but not necessarily the vertices):
  `TRAIL_WALK`
* a cycle (walk whose edges are distinct, starting and ending at the same
  vertex): `CYCLE_WALK`
* a path (walk whose edges and vertices are distinct; sometimes called a
  "simple path"): `PATH_WALK`

Although each walk is supposed to have a first vertex `s` (often called the
source vertex) and a last vertex `t` (often called the destination vertex),
this distinction does not make sense for cycles.
"""
@enum WalkType begin
    TRAIL_WALK
    CYCLE_WALK
    PATH_WALK
end

"""
    WalkSubType

The specific type of walk (e.g., a path or a cycle) implemented by the
constraint:

* a simple walk: `NO_SPECIFIC_WALK`
* a Eulerian walk (every edge is visited exactly once): `EULERIAN_WALK`
* a Hamiltonian walk (every vertex is visited exactly once): `HAMILTONIAN_WALK`

Subtypes can be thought of as adjectives for the walk types of `WalkType`.
"""
@enum WalkSubType begin
    EULERIAN_WALK
    HAMILTONIAN_WALK
    NO_SPECIFIC_WALK
end

"""
    WalkSourceType

The way the source vertex is given:

* no source node (typically, only for cycles): `NO_SOURCE_VERTEX`
* a fixed source vertex: `FIXED_SOURCE_VERTEX`
* a variable source vertex: `VARIABLE_SOURCE_VERTEX`
"""
@enum WalkSourceType begin
    NO_SOURCE_VERTEX
    FIXED_SOURCE_VERTEX
    VARIABLE_SOURCE_VERTEX
end

"""
    WalkDestinationType

The way the destination vertex is given:

* no destination node (typically, only for cycles): `NO_DESTINATION_VERTEX`
* a fixed destination vertex: `FIXED_DESTINATION_VERTEX`
* a variable destination vertex: `VARIABLE_DESTINATION_VERTEX`
"""
@enum WalkDestinationType begin
    NO_DESTINATION_VERTEX
    FIXED_DESTINATION_VERTEX
    VARIABLE_DESTINATION_VERTEX
end

"""
    Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int)

A walk in an undirected graph.

If the vector `x` describes the walk within a `Walk` constraint, each `x[i]`
denotes the next vertex in the `n`-vertex graph, for `i ∈ [1, n]`.

The considered graph is an undirected complete graph with `n` vertices.
To model a walk in a noncomplete graph, you can add constraints on the
variables: if the vertex `i` only has edges towards `j` and `k`, then `x[i]`
should only have the possible values `j` and `k`.

The dimensions of this set are as follows:

* first, the description of the walk, typically denoted by `x`
* second, the source vertex, depending on `WsT`
* third, the destination vertex, depending on `WtT`
* fourth, the individual weights, depending on `VWT` and `EWT` -- vertices
  (`VWT`) come before edges (`EWT`)
* fifth, the total weights, depending on `VWT` and `EWT` -- vertices (`VWT`)
  come before edges (`EWT`)
* sixth, the total weight, depending on `VWT` and `EWT` (sum of the weight over
  the vertices [`VWT`] and the edges [`EWT`])

For cycles, all the variables describing the cycle are implied to have an
integer value between `1` and `n`. For other walks, the walk-description
variables

Some variants are called `circuit` or `atour`.

GCC: https://sofdem.github.io/gccat/gccat/Ccircuit.html

## Unweighted walk

`Walk{UNWEIGHTED_VERTEX, UNWEIGHTED_EDGE, WalkType, WalkSubType,
WalkSourceType, WalkDestinationType, T}` considers an unweighted walk, for all
`WalkType`, `WalkSubType`, `WalkSourceType`, `WalkDestinationType`, and
`T <: Real`.

`x`-in-`Walk{UNWEIGHTED_VERTEX, UNWEIGHTED_EDGE, WT, WST, WsT, WtT, T}(n)`:
a walk in the complete graph of `n` vertices. `x[i]` is the index of the next
vertex after `i` in the walk.

## Fixed-edge-weight cycle

`Walk{UNWEIGHTED_VERTEX, FIXED_WEIGHT_EDGE, CYCLE_WALK, NO_SPECIFIC_WALK,
NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}` considers a cycle whose edge
weights are fixed: having an edge in the cycle increases the total weight.
Vertices are unweighted, because all of them must be included in a cycle.

`[x, tw]`-in-`Walk{UNWEIGHTED_VERTEX, FIXED_WEIGHT_EDGE, CYCLE_WALK,
NO_SPECIFIC_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}(n,
edge_weights)` where the elements of `edge_weights` have type `T`:

* `x` is a cycle in the complete graph of `n` vertices,
  `x[i]` is the index of the next vertex in the cycle
* `tw` is the total weight of the cycle, with `edge_weights` being indexed by
  the vertices:

  ``\\mathtt{tw} = \\sum_{i=1}^{n} \\mathtt{edge_weights[i, x[i]]}``

## Variable-vertex-weight path

`Walk{VARIABLE_WEIGHT_VERTEX, UNWEIGHTED_EDGE, PATH_WALK, NO_SPECIFIC_WALK,
FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}` considers a path whose
vertex weights are variable. Having a vertex in the cycle increases the total
weight, but edges do not contribute in this case (although there is no reason
why they could not).

`[x, w, tw]`-in-`Walk{VARIABLE_WEIGHT_VERTEX, UNWEIGHTED_EDGE, PATH_WALK,
NO_SPECIFIC_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}(n)`:

* `x` is a path in the complete graph of `n` vertices,
  `x[i]` is the index of the next vertex in the path
* `w` is the weight of each vertex (a vector indexed by the vertex indices)
* `tw` is the total vertex weight of the path:

  ``\\mathtt{tw} = \\sum_{i=1}^{n} \\mathtt{w[x[i]]}``
"""
struct Walk{VWT, EWT, WT, WST, WsT, WtT, T <: Real} <: MOI.AbstractVectorSet
    n_nodes::Int
    s::Int
    t::Int
    vertex_weights::AbstractVector{T}
    edge_weights::AbstractMatrix{T}

    # TODO: add a constructor to check if the parameters make sense.
    # - VARIABLE_WEIGHT_* or UNWEIGHTED_*: no vertex_weights, edge_weights
    # - cycle: no s or t
    # - VARIABLE_*_VERTEX: no s or t
end

Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, 0, 0, zeros(T, 0), zeros(T, 0, 0))
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, vertex_weights::AbstractVector{T}) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, 0, 0, vertex_weights, zeros(T, 0, 0))
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, edge_weights::AbstractMatrix{T}) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, 0, 0, zeros(T, 0), edge_weights)
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, vertex_weights::AbstractVector{T}, edge_weights::AbstractMatrix{T}) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, 0, 0, vertex_weights, edge_weights)
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, s::Int, t::Int) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, s, t, zeros(T, 0), zeros(T, 0, 0))
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, vertex_weights::AbstractVector{T}, s::Int, t::Int) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, s, t, vertex_weights, zeros(T, 0, 0))
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, edge_weights::AbstractMatrix{T}, s::Int, t::Int) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, s, t, zeros(T, 0), edge_weights)
Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes::Int, vertex_weights::AbstractVector{T}, edge_weights::AbstractMatrix{T}, s::Int, t::Int) where {VWT, EWT, WT, WST, WsT, WtT, T} = Walk{VWT, EWT, WT, WST, WsT, WtT, T}(n_nodes, s, t, vertex_weights, edge_weights)

function MOI.dimension(set::Walk{VWT, EWT, WT, WST, WsT, WtT, T}) where {VWT, EWT, WT, WST, WsT, WtT, T}
    dim = set.n_nodes

    if VWT == FIXED_WEIGHT_VERTEX
        dim += 1
    elseif VWT == VARIABLE_WEIGHT_VERTEX
        dim += set.n_nodes + 1
    end

    if EWT == FIXED_WEIGHT_EDGE
        dim += 1
    elseif EWT == VARIABLE_WEIGHT_EDGE
        dim += set.n_nodes ^ 2 + 1
    end

    if WsT == VARIABLE_SOURCE_VERTEX
        dim += 1
    end
    if WsT == VARIABLE_DESTINATION_VERTEX
        dim += 1
    end

    return dim
end

function copy(set::Walk{VWT, EWT, WT, WST, WsT, WtT, T}) where {VWT, EWT, WT, WST, WsT, WtT, T}
    return Walk{VWT, EWT, WT, WST, WsT, WtT, T}(set.n_nodes, set.s, set.t, copy(set.vertex_weights), copy(set.edge_weights))
end

function Base.:(==)(x::Walk{VWT, EWT, WT, WST, WsT, WtT, T}, y::Walk{VWT, EWT, WT, WST, WsT, WtT, T}) where {VWT, EWT, WT, WST, WsT, WtT, T}
    return x.n_nodes == y.n_nodes && x.s == y.s && x.t == y.t && x.vertex_weights == y.vertex_weights && x.edge_weights == y.edge_weights
end

# Some shortcuts.
const HamiltonianCycle{T} = Walk{UNWEIGHTED_VERTEX, UNWEIGHTED_EDGE, CYCLE_WALK, HAMILTONIAN_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}
const HamiltonianPath{T} = Walk{UNWEIGHTED_VERTEX, UNWEIGHTED_EDGE, PATH_WALK, HAMILTONIAN_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}

const FixedWeightHamiltonianCycle{T} = Walk{UNWEIGHTED_VERTEX, FIXED_WEIGHT_EDGE, CYCLE_WALK, HAMILTONIAN_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}
const FixedWeightHamiltonianPath{T} = Walk{UNWEIGHTED_VERTEX, FIXED_WEIGHT_EDGE, PATH_WALK, HAMILTONIAN_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}

const VariableWeightHamiltonianCycle{T} = Walk{UNWEIGHTED_VERTEX, VARIABLE_WEIGHT_EDGE, CYCLE_WALK, HAMILTONIAN_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}
const VariableWeightHamiltonianPath{T} = Walk{UNWEIGHTED_VERTEX, VARIABLE_WEIGHT_EDGE, PATH_WALK, HAMILTONIAN_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}

const EulerianCycle{T} = Walk{UNWEIGHTED_VERTEX, UNWEIGHTED_EDGE, CYCLE_WALK, EULERIAN_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}
const EulerianPath{T} = Walk{UNWEIGHTED_VERTEX, UNWEIGHTED_EDGE, PATH_WALK, EULERIAN_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}

const FixedWeightEulerianCycle{T} = Walk{FIXED_WEIGHT_VERTEX, UNWEIGHTED_EDGE, CYCLE_WALK, EULERIAN_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}
const FixedWeightEulerianPath{T} = Walk{FIXED_WEIGHT_VERTEX, UNWEIGHTED_EDGE, PATH_WALK, EULERIAN_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}

const VariableWeightEulerianCycle{T} = Walk{VARIABLE_WEIGHT_VERTEX, UNWEIGHTED_EDGE, CYCLE_WALK, EULERIAN_WALK, NO_SOURCE_VERTEX, NO_DESTINATION_VERTEX, T}
const VariableWeightEulerianPath{T} = Walk{VARIABLE_WEIGHT_VERTEX, UNWEIGHTED_EDGE, PATH_WALK, EULERIAN_WALK, FIXED_SOURCE_VERTEX, FIXED_DESTINATION_VERTEX, T}


text:
## (c) 2013--2014 David A. van Leeuwen, (c) 2005--2012, Dan Ellis

## Freely adapted from Dan Ellis's rastamat package.
## We've kept routine names the same, but the interface has changed a bit.

## we haven't implemented rasta filtering, yet, in fact.  These routines are a minimum for
## encoding HTK-style mfccs

# powspec tested against octave with simple vectors

using DSP
using FFTW
using LinearAlgebra

function powspec(x::Vector{T}, sr::Real=8000.0; wintime=0.025, steptime=0.01, dither=true) where {T<:AbstractFloat}
    nwin = round(Integer, wintime * sr)
    nstep = round(Integer, steptime * sr)

    nfft = 2 .^ Integer((ceil(log2(nwin))))
    window = hamming(nwin)      # overrule default in specgram which is hanning in Octave
    noverlap = nwin - nstep

    y = spectrogram(x .* (1<<15), nwin, noverlap, nfft=nfft, fs=sr, window=window, onesided=true).power
    ## for compability with previous specgram method, remove the last frequency and scale
    y = y[1:end-1, :] ##  * sumabs2(window) * sr / 2
    y .+= dither * nwin / (sum(abs2, window) * sr / 2) ## OK with julia 0.5, 0.6 interpretation as broadcast!

    return y
end

# audspec tested against octave with simple vectors for all fbtypes
function audspec(x::Matrix{T}, sr::Real=16000.0; nfilts=ceil(Int, hz2bark(sr/2)), fbtype::FilterBank=Bark(),
                 minfreq=0., maxfreq=sr/2, sumpower=true, bwidth=1.0) where {T<:AbstractFloat}
    nfreqs, nframes = size(x)
    nfft = 2(nfreqs-1)
    wts = fft2band(fbtype, nfft, nfilts, sr=sr, width=bwidth, minfreq=minfreq, maxfreq=maxfreq)

    wts = wts[:, 1:nfreqs]
    if sumpower
        return wts * x
    else
        return (wts * sqrt.(x)).^2
    end
end

function fft2band(fbtype::Bark, nfft::Int, nfilts::Int; sr=8000.0, width=1.0, minfreq=0., maxfreq=sr/2)
    hnfft = nfft >> 1
    minbark = hz2band(minfreq, fbtype)
    nyqbark = hz2band(maxfreq, fbtype) - minbark
    wts=zeros(nfilts, nfft)
    stepbark = nyqbark/(nfilts-1)
    binbarks = hz2band.((0:hnfft) * sr / nfft, fbtype)
    for i in 1:nfilts
        midbark = minbark + (i-1) * stepbark
        lof = (binbarks .- midbark) / width .- 0.5
        hif = (binbarks .- midbark) / width .+ 0.5
        logwts = min.(0, hif, -2.5lof)
        wts[i, 1:1+hnfft] = 10.0.^logwts
    end
    return wts
end

scale_channel!(::FilterBank, wts::AbstractArray, nfilts::Int, binfreqs) = nothing
## unclear what this does...
## Slaney-style mel is scaled to be approx constant E per channel
scale_channel!(::HTKMel, wts::AbstractArray, nfilts::Int, binfreqs) = wts = broadcast(*, 2 ./ ((binfreqs[3:nfilts+2]) - binfreqs[1:nfilts]), wts)

function fft2band(fbtype::MelFilterBank, nfft::Int, nfilts::Int; sr=8000.0, width=1.0, minfreq=0.0, maxfreq=sr/2)
    wts = zeros(nfilts, nfft)
    # Center freqs of each DFT bin
    fftfreqs = collect(0:nfft-1) / nfft * sr;
    # 'Center freqs' of mel bands - uniformly spaced between limits
    minmel = hz2band(minfreq, fbtype)
    maxmel = hz2band(maxfreq, fbtype)
    binfreqs = band2hz(minmel .+ collect(0:(nfilts+1)) / (nfilts + 1) * (maxmel - minmel), fbtype);
    ## binbin = iround(binfrqs/sr*(nfft-1));

    for i in 1:nfilts
        fs = binfreqs[i .+ (0:2)]
        # scale by width
        fs = fs[2] .+ (fs .- fs[2])width
        # lower and upper slopes for all bins
        loslope = (fftfreqs .- fs[1]) / (fs[2] - fs[1])
        hislope = (fs[3] .- fftfreqs) / (fs[3] - fs[2])
        # then intersect them with each other and zero
        wts[i,:] = max.(0, min.(loslope, hislope))
    end

    scale_channel!(fbtype, wts, nfilts, binfreqs)

    # Make sure 2nd half of DFT is zero
    wts[:, (nfft>>1)+1:nfft] .= 0.
    return wts
end

hz2band(::T, ::F) where {F <: FilterBank, T <: AbstractFloat }= error("hz2band not implemented for filter bank $F")
hz2band(f::T, ::Bark) where {T<:AbstractFloat} = 6asinh(f / 600)
hz2band(f::T, ::HTKMelFilterBank) where {T<:AbstractFloat} = 2595 * log10(1 + f / 700)
function hz2band(f::T, ::Mel) where {T<:AbstractFloat}
    f0 = 0.0
    fsp = 200/3
    brkfrq = 1000.0
    brkpt = (brkfrq - f0) / fsp
    logstep = exp(log(6.4) / 27)
    f < brkfrq ? f / brkfrq / log(logstep) : brkpt + log(f / brkfrq) / log(logstep)
end

function hz2band(v::Vector{T}, fbtype::F) where {F <: FilterBank, T<:AbstractFloat}
    hz2band_scalar(x::T) = hz2band(x, fbtype)
    hz2band_scalar.(v)
end

band2hz(::T, ::F) where {F <: FilterBank, T <: AbstractFloat  }= error("band2hz not implemented for filter bank $F")
band2hz(b::T, ::Bark) where {T<:AbstractFloat} = 600 * sinh(b / 6)
band2hz(b::T, ::HTKMelFilterBank) where {T<:AbstractFloat} = 700 * (10 ^ (b / 2595) - 1)
function band2hz(b::T, ::Mel) where {T<:AbstractFloat}
    f0 = 0.0
    fsp = 200/3
    brkfrq = 1000.0
    brkpt = (brkfrq - f0) / fsp
    logstep = exp(log(6.4) / 27)
    b < brkpt ? f0 + fsp * b :  brkfrq * exp(log(logstep) * (b - brkpt))
end

function band2hz(v::Vector{T}, fbtype::F) where {F <: FilterBank, T<:AbstractFloat}
    band2hz_scalar(x::T) = band2hz(x, fbtype)
    band2hz_scalar.(v)
end

function postaud(x::Matrix{T}, fmax::Real, fbtype::FilterBank=Bark(), broaden=false) where {T<:AbstractFloat}
    (nbands,nframes) = size(x)
    nfpts = nbands + 2broaden

    bandcfhz = band2hz(range(0, hz2band(fmax, fbtype), length=nfpts), fbtype)

    # Remove extremal bands (the ones that will be duplicated)
    bandcfhz = bandcfhz[1+broaden:nfpts-broaden];
    # Hynek's magic equal-loudness-curve formula
    fsq = bandcfhz.^2
    ftmp = fsq + 1.6e5
    eql = ((fsq ./ ftmp).^2) .* ((fsq + 1.44e6) ./ (fsq + 9.61e6))
    # weight the critical bands
    z = broadcast(*, eql, x)
    # cube root compress
    z .^= 0.33
    # replicate first and last band (because they are unreliable as calculated)
    if broaden
        z = z[[1, 1:nbands, nbands], :];
    else
        z = z[[2, 2:(nbands-1), nbands-1], :]
    end
    return z
end

function dolpc(x::Array{T}, modelorder::Int=8) where {T<:AbstractFloat}
    nbands, nframes = size(x)
    r = real(ifft(vcat(x, x[collect(nbands-1:-1:2), :]), 1)[1:nbands, :])
    # Find LPC coeffs by durbin
    y, e = levinson(r, modelorder)
    # Normalize each poly by gain
    y ./= e
end

function lpc2cep(a::Array{T}, ncep::Int=0) where {T<:AbstractFloat}
    nlpc, nc = size(a)
    order = nlpc - 1
    if ncep == 0
        ncep = nlpc
    end
    c = zeros(ncep, nc)
    # Code copied from HSigP.c: LPC2Cepstrum
    # First cep is log(Error) from Durbin
    c[1,:] = -log.(a[1, :])
    # Renormalize lpc A coeffs
    a ./= a[1, :]
    for n in 2:ncep
        sum = zero(T)
        for m in 2:n
            sum += (n - m) * a[m, :] .* c[n - m + 1, :]
        end
        c[n, :] = -a[n, :] - sum / (n - 1)
    end
    return c
end

function spec2cep(spec::Array{T}, ncep::Int=13, dcttype::Int=2) where {T<:AbstractFloat}
    # no discrete cosine transform option
    dcttype == -1 && return log(spec)

    (nr, nc) = size(spec)
    dctm = zeros(typeof(spec[1]), ncep, nr)
    if 1 < dcttype < 4          # type 2,3
        for i in 1:ncep
            dctm[i, :] = cos.((i - 1) * collect(1:2:2nr-1)π / (2nr)) * √(2/nr)
        end
        if dcttype == 2
            dctm[1, :] /= √2
        end
    elseif dcttype == 4           # type 4
        for i in 1:ncep
            dctm[i, :] = 2cos.((i-1) * collect(1:nr)π/(nr+1))
            dctm[i, 1] += 1
            dctm[i, nr] += (-1)^(i-1)
        end
        dctm /= 2(nr + 1)
    else                        # type 1
        for i in 1:ncep
            dctm[i, :] = cos.((i-1) * collect(0:nr-1)π / (nr - 1)) / (nr - 1)
        end
        dctm[:, [1, nr]] /= 2
    end
    return dctm * log.(spec)
end

function lifter(x::Array{T}, lift::Real=0.6, invs=false) where {T<:AbstractFloat}
    (ncep, nf) = size(x)
    if lift == 0
        return x
    end
    if lift > 0
        if lift > 10
            error("Too high lift number")
        end
        liftw = [1; collect(1:ncep-1).^lift]
    else
        # Hack to support HTK liftering
        if !isa(lift, Integer)
            error("Negative lift must be interger")
        end
        lift = -lift            # strictly speaking unnecessary...
        liftw = vcat(1, (1 .+ lift/2 * sin.(collect(1:ncep-1)π / lift)))
    end
    if invs
        liftw = 1 ./ liftw
    end
    y = broadcast(*, x, liftw)
    return y
end

## Freely after octave's implementation, by Paul Kienzle <pkienzle@users.sf.net>
## Permission granted to usee this in a MIT license on 20 dec 2013 by the author Paul Kienzle:
## "You are welcome to move my octave code from GPL to MIT like core Julia."
## untested
## only returns a, v
function levinson(acf::Vector{T}, p::Int) where {T<:Real}
    if length(acf) < 1
        error("empty autocorrelation function")
    end
    if p < 0
        error("negative model order")
    end
    if p < 100
        ## direct solution [O(p^3), but no loops so slightly faster for small p]
        ## Kay & Marple Eqn (2.39)
        R = toeplitz(acf[1:p], conj(acf[1:p]))
        a = R \ -acf[2:p+1]
        pushfirst!(a, 1)
        v = real(a.*conj(acf[1:p+1]))
    else
        ## durbin-levinson [O(p^2), so significantly faster for large p]
        ## Kay & Marple Eqns (2.42-2.46)
        ref = zeros(p)
        g = -acf[2] / acf[1]
        a = [g]
        v = real((1 - abs2(g)) * acf[1])
        ref[1] = g
        for t=2:p
            g = -(acf[t+1] + a ⋅ acf[t:-1:2]) / v
            a = [a + g * conj(a[t-1:-1:1]), g]
            v *= 1 - abs2(g)
            ref[t] = g
        end
        pushfirst!(a, 1)
    end
    return (a, v)
end

function levinson(acf::Array{T}, p::Int) where {T<:Real}
    (nr,nc) = size(acf)
    a = zeros(p + 1, nc)
    v = zeros(p + 1, nc)
    for i in 1:nc
        a[:,i], v[:,i] = levinson(acf[:,i], p)
    end
    return (a, v)
end

## Freely after octave's implementation, ver 3.2.4, by jwe && jh
## skipped sparse implementation
function toeplitz(c::Vector{T}, r::Vector{T}=c) where {T<:Real}
    nc = length(r)
    nr = length(c)
    res = zeros(typeof(c[1]), nr, nc)
    if nc == 0 || nr == 0
        return res
    end
    if r[1] != c[1]
        ## warn
    end
    if typeof(c) <: Complex
        conj!(c)
        c[1] = conj(c[1])       # bug in julia?
    end
    ## if issparse(c) && ispsparse(r)
    data = [r[end:-1:2], c]
    for (i, start) in zip(1:nc, nc:-1:1)
        res[:,i] = data[start:start+nr-1]
    end
    return res
end


text:
using Flux, CuArrays, Metalhead

struct Bottleneck
  layer
end

Flux.treelike(Bottleneck)

Bottleneck(in_planes, growth_rate) = Bottleneck(
                                           Chain(BatchNorm(in_planes, relu),
                                           Conv((1, 1), in_planes=>4growth_rate),
                                           BatchNorm(4growth_rate, relu),
                                           Conv((3, 3), 4growth_rate=>growth_rate, pad = (1, 1))))

(b::Bottleneck)(x) = cat(3, b.layer(x), x)

Transition(chs::Pair{<:Int, <:Int}) = Chain(BatchNorm(chs[1], relu),
                                            Conv((1, 1), chs),
                                            x -> meanpool(x, (2, 2)))

function _make_dense_layers(block, in_planes, growth_rate, nblock)
  local layers = []
  for i in 1:nblock
    push!(layers, block(in_planes, growth_rate))
    in_planes += growth_rate
  end
  Chain(layers...)
end

function _DenseNet(nblocks; in_chs = 1, block = Bottleneck, growth_rate = 12, reduction = 0.5, num_classes = 1)
  num_planes = 2growth_rate
  layers = []
  push!(layers, Conv((7, 7), in_chs=>num_planes, stride = (2, 2), pad = (3, 3)))
  push!(layers, BatchNorm(num_planes, relu))
  push!(layers, x -> maxpool(x, (3, 3), stride = (2, 2), pad = (1, 1)))

  for i in 1:3
    push!(layers, _make_dense_layers(block, num_planes, growth_rate, nblocks[i]))
    num_planes += nblocks[i] * growth_rate
    out_planes = Int(floor(num_planes * reduction))
    push!(layers, Transition(num_planes=>out_planes))
    num_planes = out_planes
  end

  push!(layers, _make_dense_layers(block, num_planes, growth_rate, nblocks[4]))
  num_planes += nblocks[4] * growth_rate
  push!(layers, BatchNorm(num_planes, relu))

  Chain(layers..., x -> meanpool(x, (7, 7)),
        x -> reshape(x, :, size(x, 4)),
        Dense(num_planes, num_classes), x -> σ.(x))
end

"""
  get_densenet_model(depth, pretrained::Bool = false)

Creates the Densenet model as per the arguments passed.

Arguments:
1. `depth`: The depth of the Densenet Model. The current allowed
            values are `121`, `169`, `201` and `264`.
2. `pretrained`: Gets the network with pretrained Imagenet weights.
                 `Metalhead.jl` currently has only the weights for
                 Densenet 121 so sending `depth` anything other than
                 `121` results in an **ERROR**.
"""

function get_densenet_model(depth, pretrained::Bool = false)
  if pretrained && depth != 121
    error("Pretrained weights are only available for Densenet 121")
  end
  if depth == 121
    if pretrained
      m = DenseNet() |> gpu
      Flux.testmode!(m, false)
      m
    else
      _DenseNet([6, 12, 24, 16]) |> gpu
    end
  elseif depth == 169
    _DenseNet([6, 12, 32, 32]) |> gpu
  elseif depth == 201
    _DenseNet([6, 12, 48, 32]) |> gpu
  elseif depth == 264
    _DenseNet([6, 12, 64, 48]) |> gpu
  else
    error("No such model is available")
  end
end


text:
using BenchmarkTools
using MPI
using Trilinos
using Test

MPI.Init()

#Kokkos.initialize(Kokkos.OpenMP(), 4)

comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))

my_rank = Teuchos.getRank(comm)
if my_rank != 0
  redirect_stdout(open("/dev/null", "w"))
  redirect_stderr(open("/dev/null", "w"))
end

function indefinite(n)
  # Generate an indefinite "hard" matrix
  srand(1)
  A = 20 * sparse(1.0I, n, n) + sprand(n, n, 20.0 / n)
  A = (A + A') / 2
  x = ones(n)
  b = A * x

  return A, b, ones(n)
end

# The following Laplace 2D implementation is based on:
# https://github.com/lruthotto/KrylovMethods.jl/blob/master/benchmarks/benchmark2DLaplacian.ipynb

"""
ibc,iin = getBoundaryIndices(n::Int)

returns indices of boundary and interior nodes of regular mesh with n^2 cells.
"""
function getBoundaryIndices(n::Int)
    ids = reshape(collect(1:(n+1)^2),n+1,n+1)
    iin = vec(ids[2:end-1,2:end-1])
    ibc = setdiff(vec(ids),vec(iin))
    return ibc,iin
end

"""
L,Lin,Lib = getLaplacian(n::Int)

returns discrete Laplacian, interior part, and boundary part, on regular mesh with n^2 cells.
"""
function getLaplacian(n::Int)
    h  = 1/n
    d2x = spdiagm(-1 => fill(-1.0/h,n), 0 => fill(2.0/h,n+1), 1 => fill(-1.0/h,n))
    d2x[1] /= 2
    d2x[end] /= 2
    L  = kron(sparse(1.0I, n+1, n+1),d2x) + kron(d2x,sparse(1.0I, n+1, n+1))

    # split into boundary and interior part
    ibc,iin = getBoundaryIndices(n)
    Lin = L[iin,iin]
    Lib = L[iin,ibc]
    return L,Lin,Lib
end
"""
g = getBoundaryConditions(n::Int,fctn=(x,y)->5.*y*sin(pi*x))

discretizes function on regular grid with n^2 cells and returns boundary and interior values as well as the source term value
"""
function getBoundaryConditions(n::Int,dirichlet = (x,y) -> (1-x^2)*(1-y^2), sourceterm = (x, y) -> 2*((1-x^2)+(1-y^2)))
    x  = linspace(0,1,n+1)
    g = map(Float64,[dirichlet(x1,x2) for x1 in x, x2 in x]')
    s = map(Float64,[sourceterm(x1,x2) for x1 in x, x2 in x]')
    ibc,iin = getBoundaryIndices(n)
    return vec(g)[ibc], vec(g)[iin], vec(s)[iin]
end

function laplace2d(n=1000)
  L,Lin,Lib = getLaplacian(n)
  ubc,refsol,f      = getBoundaryConditions(n)
  return Lin,vec(-Lib*ubc)+f,refsol
end

function tpetra_system(S)
  A,b,x = S
  n = length(b)
  rowmap = Tpetra.Map(n, 0, comm) #, Kokkos.KokkosOpenMPWrapperNode)
  At = Tpetra.CrsMatrix(A,rowmap)

  rmap = Tpetra.getRangeMap(At)
  bt = Tpetra.Vector(rmap)
  xt = Tpetra.Vector(rmap)

  num_my_elements = Tpetra.getNodeNumElements(rowmap)

  bv = Tpetra.device_view(bt)
  xv = Tpetra.device_view(xt)
  for local_row = 0:num_my_elements-1
    global_row = Tpetra.getGlobalElement(rmap, local_row)
    bv[local_row] = b[global_row+1]
    xv[local_row] = x[global_row+1]
  end

  return At, bt, xt
end

"""
Builds a Belos / Tpetra linear system, solving it once and checking the result against the reference solution
"""
function build_checked_system(f, solver_name, prec_type, n, tol, restart::Int, maxiter::Int)
  A, b, refsol = tpetra_system(f(n))

  println("Tolerance = ", tol, "; restart = ", restart, "; max #iterations = ", maxiter)

  params = Trilinos.default_parameters(solver_name)
  solver_params = params["Linear Solver Types"]["Belos"]["Solver Types"][solver_name]
  solver_params["Convergence Tolerance"] = tol
  solver_params["Verbosity"] = Belos.StatusTestDetails + Belos.FinalSummary + Belos.TimingDetails
  solver_params["Maximum Restarts"] = Int32(restart)
  solver_params["Block Size"] = Int32(maxiter)
  solver_params["Maximum Iterations"] = Int32(maxiter)
  params["Preconditioner Type"] = prec_type

  solver = TpetraSolver(A, params)

  sol = solver \ b
  v = Tpetra.device_view(sol)
  refv = Tpetra.device_view(refsol)
  maxerr = 0
  for i in LinearIndices(v)
    maxerr = max(abs(v[i] - refv[i]), maxerr)
  end
  @test maxerr < 100*tol

  return solver, b
end

sol_ndef(solver_name; n = 500000, tol = 1e-5, restart::Int = 15, maxiter::Int = 1500) = build_checked_system(indefinite, solver_name, "None", n, tol, restart, maxiter)

"""
Note: n is not the size of the matrix, but n in an n × n grid!
"""
sol_posdef(solver_name; n = 1000, tol = 1e-10, maxiter::Int = 1500) = build_checked_system(laplace2d, solver_name, "Ifpack2", n, tol, 0, maxiter)

suite = BenchmarkGroup()

suite["GMRES"] = @benchmarkable x=A\b setup=((A,b) = sol_ndef("GMRES"))
suite["BICGSTAB"] = @benchmarkable x=A\b setup=((A,b) = sol_ndef("BICGSTAB"))
suite["CG"] = @benchmarkable A\b setup=((A,b) = sol_posdef("CG",n=300))

result = run(suite)
if(my_rank == 0)
  display(result)
end

MPI.Finalize()


text:
export Random, random
"""
    random(x;args)

Randomly generate a vector of size `x`, based on provided (optional arguments) :
- `L` is the typical wavelength (default is `L=1`),
- `s` is the (real) Sobolev index regularity (default is `s=∞`, smooth data),
- `λ` is the length of spatial localization (default is `λ=∞`, no localization),
- `a` is the amplitude of the returned vector (default is `a=1`).

The vector is generated through randomly chosen Fourier coefficients,
multiplied with weigth `w=10^(-|k|L/(2π))` if `s=∞`, or `w=1/(1+9(|k|L/(2π))^(s+1/2))` otherwise.
If `λ≠∞`, the function in spatial variables is multiplied by `exp(-|x/λ|^2)`,
and in any case normalized to have maximum absolute value 1.

"""
function random( x;L=1,s=Inf,λ=Inf,a=1 )

    k = Mesh( x ).k
    if s == Inf
        w = 10 .^(-abs.(k*L/(2*π)))
    else
        w = 1 ./( 1 .+ 9*abs.(k*L/(2*π)).^(s+1/2) )
    end
    if λ == Inf
        φ = zero(x).+1
    else
        φ = exp.(-abs.((x/λ).^2))
    end

    θ = 2*π*rand(Float64,length(x))
    r = rand(Float64,length(x))
    Fourier = r.*exp.(-1im.*θ).*w
    Physic = real.(ifft(Fourier)).*φ
    return Physic./(maximum(Physic)-minimum(Physic))*a
end

"""
    Random(param;args)

Randomly generated initial data, based on provided (optional arguments) :
- `L` is the typical wavelength (default is `L=1`),
- `s` is the (real) Sobolev index regularity (default is `s=∞`),
- `λ` is the length of spatial localization (default is `λ=∞`, no localization),
- `a` is the couple of amplitudes of the surface deformation, and velocity (default is `a=(1,1)`).

Return an initial data `init::InitialData`,  to be used in initial-value problems `Problem(model, init, param)`,
with `init.η` and `init.v` two samples of `random(x;L,s,λ)`,
where `x` is the set of collocation points generated by `Mesh(param)`.

"""
struct Random <: InitialData

    η
    v
    label :: String
    info  :: String

    function Random(param;L=1,s=Inf,λ=Inf,a=(1,1))
        mesh=Mesh(param);x=mesh.x;
        η = random( x; L=L,s=s,λ=λ,a=a[1] )
        v = random( x; L=L,s=s,λ=λ,a=a[2] )
        init=Init(x,η,v)
        if s == Inf s="∞" end
        if λ == Inf λ="∞" end
        label = "randomly generated"
        info = "Randomly generated initial data with \n\
        ├───typical wavelength: L = $L,\n\
        ├─────regularity index: s = $s,\n\
        ├─spatial localization: λ = $λ,\n\
        ├────surface amplitude: $(a[1]),\n\
        └───velocity amplitude: $(a[2])."
    	new( init.η,init.v,label,info )

    end

end


text:
module Swarms

include("Solvers.jl")

using TensorOperations

export Swarm, step!

"""
    Swarm(func; [size, dim, tol, step_size, positions])

Structure storing all information needed to describe a swarm of particles.
"""
struct Swarm
    func::Function
    size::Int64
    dim::Int8
    tol::Float64
    positions::Matrix{Float64}
    step_size::Vector{Float64}
    error_history::Array{Float64}

    function Swarm(func::Function;
        size::Int64 = 1000,
        dim::Int64 = 3,
        tol::Float64 = 1e-8,
        step_size::Vector{Float64} = repeat([1e-2], size),
        positions::Matrix{Float64} = rand(Float64, dim, size) .*2 .-1)

        if dim != 3 error("Number of dimensions not supported.") end

        new(func, size, dim, tol, positions, step_size, Array{Float64}(undef, dim, size, 0))
    end
end

function Swarm(func::Function, default_function_params::NamedTuple)

    # if dim != 3 error("Number of dimensions not supported") end

    Swarm(x -> func(x; default_function_params...))
end

"""
    step!(swarm, solver)

Step the swarm one step using the provided solver.
"""
function step!(swarm::Swarm, solver::RungeKuttaMethod)

    h = swarm.step_size

    dx⃗ = zeros(Float64, size(swarm.positions)..., length(solver.β))
    y_step = Matrix{Float64}(undef, size(swarm.positions)...)

    for i in 1:length(solver.β)
        weights = solver.γ[i,:]
        @tensor y_step[i,j] = dx⃗[i,j,k]*weights[k]
        dx⃗[:,:,i] = swarm.func(swarm.positions + h'.*y_step)'
    end

    @tensor y_step[i,j] = dx⃗[i,j,k]*solver.α[k]

    swarm.positions .+= h'.*y_step
end

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Graphics (Modia3D/Graphics/_module.jl)
#

# Shadow definitions
@enum ShadowMask NoShadows=0 CastsShadows=1 ReceivesShadows=2 CastsAndReceivesShadows=3

# Material for visual elements
"""
    material = Modia3D.Material(;color=defaultColor(), wireframe=false, transparency=0.0,
                                 reflectslight=true, shininess=0.7,
                                 shadowMask=CastsAndReceivesShadows))

Return a material object that defines attributes for the visualization of an Object3D
that has visual or solid properties.

# Arguments
- `color`: This argument is passed to function [`Modia3D.rgb`](@ref)`(color)` to return the RGB color value in
           form of a vector. E.g. `color="Red"` or `color=[255,0,0]`.
- `wireframe`: = false, if solid, otherwise wireframe representation.
- `transparency`: = 0.0 (opaque) ... 1.0 (fully transparent).
- `reflectslight`: = true if it reflects light and false, if it does not reflect light.
- `shininess`: = 0.0 (matte surface) ... 1.0 (very shiny surface).
- `shadowMask`: defines whether or not an object casts or receives shadows. Possible values:
   NoShadows, CastsShadows, ReceivesShadows, CastsAndReceivesShadows.
"""
mutable struct Material
  color::MVector{3,Cint}  # RGB color
  wireframe::Bool         # if true, the wireframe of object is displayed
  transparency::Float64   # [0,1]... 0.0 is opaque, 1.0 is transparent
  reflectslight::Bool     # if true, specular Highlights are active
  shininess::Float64      # [0,1]... 0.0 is matte surface, 1.0 is very shiny
  shadowMask::ShadowMask  # defines whether or not an object casts or receives shadows 0: No Shadows, 1: Casts Shadows, 2: Receives Shadows, 3: Casts and Recieves Shadows

  function Material(;color=defaultColor(), wireframe=false, transparency=0.0, reflectslight=true, shininess=0.7, shadowMask=CastsAndReceivesShadows)
    color2 = rgb(color)
    @assert(color2[1]>=0 && color2[1]<=255)
    @assert(color2[2]>=0 && color2[2]<=255)
    @assert(color2[3]>=0 && color2[3]<=255)
    @assert(transparency >= 0.0 && transparency <= 1.0)
    @assert(shininess    >= 0.0 && shininess    <= 1.0)
    new(color2, wireframe, transparency, reflectslight, shininess, shadowMask)
  end
end


text:
module DiffEqBenchmarks

end


text:
using VectorEngine
using GPUCompiler
using LinearAlgebra
using LLVM, LLVM.Interop
using Test

using Random
Random.seed!(1)

include("setup.jl")

@testset "VectorEngine" begin
    @testset "Memory Operations" begin
        include("memory.jl")
    end
    @testset "Pointer Operations" begin
        include("pointer.jl")
    end
    @testset "Arrays" begin
        include("array.jl")
    end
    @testset "Execution" begin
        include("execution.jl")
    end
    @testset "Vectorize" begin
        include("vectorize.jl")
    end
    #@testset "Examples" begin
    #    include("examples.jl")
    #end

    #@testset "Device Tests" begin
    #    include("device/globals.jl")
    #end
end


text:
struct LeastActionResult
    path
    lam # Least action energy
    res
end

@with_kw struct curve_ind
    integral_func
    func! # ODE function (du, u,p, t)
    DiffusionMatrix
    dt
    dim
    n_points
    p
    t = nothing # for ode: placeholder
end

function (self::curve_ind)(x)
    points = reshape(x, (self.dim, self.n_points) )

    intd = self.integral_func(points, self.func!, self.DiffusionMatrix, self.dt, self.p; t=self.t )
    return intd
end

"""
Least Action Method

Arguements
----------
- `n_points` {Integer}: Number of time points
- `tmax` {Float}: maximum time
- `point_start` {Array}: initial variables
- `point_end` {Array}: final state
- `ode!` {Function}: ODE function (non-allocation type) f(du,u,p,t)
- `p` {Array}: parameters `p` of `ode!`

Return
------
- `op` {Optim.MultivariateOptimizationResults} : keywords. `ls_success` (success of optimization). Use `Optim.minimizer` function to get the arg minimum.
"""
function action(n_points, tmax, point_start, point_end, ode!, p, Jacobian, DiffusionMatrix)

    dim = length(point_start)
    dt = tmax/n_points

    # A straight line from A → B
    initpath = LeastAction.linearpath(point_start, point_end, n_points)
    init_vec = reshape(initpath, 1, :)

    # Function of curve integration
    curve_func = curve_ind(integral, ode!, DiffusionMatrix, dt, dim, n_points, p, nothing)
    g! = grad!(Jacobian, dt, dim, n_points)

    op = optimize(curve_func, g!, init_vec, ConjugateGradient())

    path = reshape(Optim.minimizer(op), size(initpath))
    lam = op.minimum

    return LeastActionResult(path, lam, op)

end

"""
Integration along the path
"""
function integral(points, func!, DiffusionMatrix, dt, p; t=nothing)
    dim = size(points)[1]
    integral = 0
    for k = 2: size(points)[end]
        du = get_dxdt(func!, points[:,k-1], p; t=nothing)

        Tmp = Δ(points, k, dt) .- du #eq18.
        integral = integral + Tmp' * DiffusionMatrix(dim)^(-1) * Tmp * dt
    end
    integral = integral/4

    return integral
end

"""
Gradient for optimzation
"""
struct grad!
    Jacobian
    dt
    dim
    N
end

function (self::grad!)(G, x)
    x_ = reshape(x, self.dim, self.N)

    m_jac = self.Jacobian(self.dt, x_)

    G[:] = reshape(m_jac, size(x))
end


text:
module BernoulliTest

using Test
using ReactiveMP
using Distributions
using Random

@testset "Bernoulli" begin

    # Bernoulli comes from Distributions.jl and most of the things should be covered there
    # Here we test some extra ReactiveMP.jl specific functionality

    @testset "vague" begin
        d = vague(Bernoulli)

        @test typeof(d) <: Bernoulli
        @test mean(d) === 0.5
        @test succprob(d) === 0.5
        @test failprob(d) === 0.5
    end

    @testset "prod" begin
        @test prod(ProdAnalytical(), Bernoulli(0.5), Bernoulli(0.5))   ≈ Bernoulli(0.5)
        @test prod(ProdAnalytical(), Bernoulli(0.1), Bernoulli(0.6))   ≈ Bernoulli(0.14285714285714285)
        @test prod(ProdAnalytical(), Bernoulli(0.78), Bernoulli(0.05)) ≈ Bernoulli(0.1572580645161291)
    end

    @testset "probvec" begin
        @test probvec(Bernoulli(0.5)) === (0.5, 0.5)
        @test probvec(Bernoulli(0.3)) === (0.3, 0.7)
        @test probvec(Bernoulli(0.6)) === (0.6, 0.4)
    end

end

end


text:
# push!(LOAD_PATH, "Y:\\")
# using Pkg
# Pkg.add("NLopt")
# Pkg.add("RecipesBase")

using LikelihoodProfiler, Test

include("./cases_func.jl")
methods_list = [:CICO_ONE_PASS, :LIN_EXTRAPOL, :QUADR_EXTRAPOL]

println("Starting tests for get_interval")
@testset "get_interval" begin include("test_get_interval.jl") end

println("Starting tests for get_endpoint")
@testset "get_endpoint" begin include("test_get_endpoint.jl") end
@testset "get_endpoint CICO_ONE_PASS" begin include("test_get_endpoint.jl") end
@testset "get_endpoint errors and warnings" begin include("test_get_endpoint_errors.jl") end

println("Starting tests for all all methods of get_right_endpoint")
@testset "get_right_endpoint methods" begin include("test_methods.jl") end
@testset "get_right_endpoint CICO_ONE_PASS" begin include("test_cico_one_pass.jl") end
@testset "get_right_endpoint methods errors and warnings" begin include("test_methods_errors.jl") end

println("Starting tests for profile")
@testset "profile" begin include("test_profile.jl") end

println("Starting tests for different fitting alg")
@testset "local_alg" begin include("test_algorithms.jl") end

println("Starting tests for Plot @recipe")
@testset "PLOT_INTERVAL" begin include("test_plots.jl") end

println("Starting tests for loss error")
@testset "LOSS_ERROR" begin include("test_loss_error.jl") end


text:
function codegen_dynamics(model)
    # dimensions
    nq = model.nq

    # variables
    @variables q[1:nq] q̇[1:nq]

    # Lagrangian
    L = lagrangian(model, q, q̇)

    ddL = Symbolics.hessian(L, [q; q̇])
    dLq = Symbolics.gradient(L, q)
    ddLq̇q = ddL[nq .+ (1:nq), 1:nq]

    M = ddL[nq .+ (1:nq), nq .+ (1:nq)]
    C = ddLq̇q * q̇ - dLq

    # mass matrix and dynamics bias
    mass_matrix = eval(Symbolics.build_function(M, q)[1])
    dynamics_bias = eval(Symbolics.build_function(C, q, q̇)[1])

    return mass_matrix, dynamics_bias
end


text:
using SurrealNumbers
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

out_dir = "Data/"

#### PRESUMES THAT DOT (from GraphVis) is installed

# recursive construction and one and zero
x0 = SurrealFinite("0", ϕ,   ϕ)
x1 = SurrealFinite("1", [x0], ϕ)
x11 = SurrealFinite("-1", ϕ,  [x0])
x4 = SurrealFinite( [x11],  [x1,x0])
x5 = SurrealFinite( [x11,x0,x4],  [x1])
x21  = SurrealFinite("2", [x1],  ϕ)
x22  = SurrealFinite( [x0, x1],  ϕ)
x23  = SurrealFinite( [x11,x1],  ϕ)
x24  = SurrealFinite( [x11,x0,x1],  ϕ)
x25 = [x11,x0,x1] ≀ ϕ
z = zero(x1)
z = zero(SurrealFinite)
z = one(x1)
s1 = convert(SurrealFinite, 1//2)
s2 = convert(SurrealFinite, 3//4)
s3 = convert(SurrealFinite, 1//2) + convert(SurrealFinite, 1//4)

x3 = convert(SurrealFinite, 3)
x41 = convert(SurrealFinite, 4)
x42 = SurrealFinite( [x3],  ϕ)
x43 = convert(SurrealFinite, 2) * convert(SurrealFinite, 2)

x6 = convert(SurrealFinite, 3) * convert(SurrealFinite, 2)
x8 = convert(SurrealFinite, 2) * convert(SurrealFinite, 4)
x9 = convert(SurrealFinite, 3) * convert(SurrealFinite, 3)

x03 = convert(SurrealFinite, 3) - convert(SurrealFinite, 3)

x00 = SurrealFinite("0", [-x1],   [x1])
x000 = x00 + x00

surreal2dot(stdout, x1)

file = "$(out_dir)test_dot_x41.dot"
FID = open(file, "w")
surreal2dot(FID, x41)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x41.dot"
FID = open(file, "w")
surreal2dag(FID, x41)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dot_s2.dot"
FID = open(file, "w")
surreal2dot(FID, s2; direction="back")
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_s2.dot"
FID = open(file, "w")
surreal2dag(FID, s2)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dot_s3.dot"
FID = open(file, "w")
surreal2dot(FID, s3)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_s3.dot"
FID = open(file, "w")
surreal2dag(FID, s3)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dot_x43.dot"
FID = open(file, "w")
surreal2dot(FID, x43)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x43.dot"
FID = open(file, "w")
surreal2dag(FID, x43)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dot_x1.dot"
FID = open(file, "w")
surreal2dot(FID,  convert(SurrealFinite, 1))
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dot_x00.dot"
FID = open(file, "w")
surreal2dot(FID,  x00; direction="back")
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

# temp1 = values(ExistingSurreals)
# temp1 = sort(collect(temp1))

file = "$(out_dir)test_dot_x5.dot"
FID = open(file, "w")
surreal2dot(FID, x5)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x5.dot"
FID = open(file, "w")
surreal2dag(FID, x5)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

# file = "$(out_dir)test_dot_x6.dot"
# FID = open(file, "w")
# surreal2dot(FID, x6)
# close(FID)
# run(`dot -Tpdf -O $file`)
# run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x6.dot"
FID = open(file, "w")
surreal2dag(FID, x6)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x8.dot"
FID = open(file, "w")
surreal2dag(FID, x8)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x9.dot"
FID = open(file, "w")
surreal2dag(FID, x9)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x1.dot"
FID = open(file, "w")
surreal2dag(FID,  convert(SurrealFinite, 1))
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x00.dot"
FID = open(file, "w")
surreal2dag(FID,  x00; direction="back")
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

file = "$(out_dir)test_dag_x000.dot"
FID = open(file, "w")
surreal2dag(FID,  x000)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

temp2 = values(ExistingSurreals)
temp2 = sort(collect(temp2))


text:
using Mocking

# Note: Explicitly setting Mocking.ENABLE should only be needed on Julia 0.4
if VERSION < v"0.5-"
    Mocking.enable()
end

using Base.Test
import Mocking: apply

include("expr.jl")
include("patch.jl")

include("concept.jl")
include("scope.jl")
include("closure.jl")
include("import.jl")
include("real-open.jl")
include("real-isfile.jl")
include("real-nested.jl")
include("readme.jl")
include("optional.jl")


text:
"""
	minrequired(::Type{XXX})

Returns the minimum required properties.  Other classes implement this to check whether a Spectrum or Dict has all the
necessary properties for the specified algorithm or data structure.
"""
minproperties(::Type{Any}) =
    @assert "minproperties(...) needs to be specialized for this algorithm."

"""
    hasminrequired(ty::Type, item::Union{Spectrum,Dict{Symbol,Any}})

Does this spectrum have the minimal set of required properties?
"""
hasminrequired(ty::Type, dict::Dict{Symbol,Any}) =
    all(haskey(dict, a) for a in minproperties(ty))

"""
    requiredbutmissing(ty::Type, item::Union{Spectrum,Dict{Symbol,Any}})

List any required but missing properties.
"""
requiredbutmissing(ty::Type, dict::Dict{Symbol,Any}) =
    filter(a -> !haskey(dict, a), minproperties(ty))