text:
using .Instantiation
using .ModelElaboration

using Unitful
using .StructuralTransform
using .Synchronous: sample, Clock, previous, hold, positive, positiveChange, positiveEdge

"""
Shortcut for `Variable`
"""
Var(; args...) = Variable(; args...)

"""
Create a floating-point `Variable`
"""
Float(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, Float64, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

"""
Create a boolean `Variable`
"""
Boolean(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, Bool, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

"""
Create an integer `Variable`
"""
Integ(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, Int, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

"""
Create a string `Variable`
"""
Str(value=nothing; info="", size=nothing, unit=NoUnits, displayUnit=NoUnits,
    min=nothing, max=nothing, start=nothing, fixed::Bool=false, nominal=nothing,
    variability=continuous,
    flow::Bool=false, state::Bool=true) = Variable(variability, String, size, value,
    unit, displayUnit, min, max, start, fixed, nominal, info, flow, state, general)

#=
Float(; args...) = Var(T=Float64; args...)
Float0(; args...) = Float(size=(); args...)

Boolean(; args...) = Var(T=Bool; args...)
Boolean0(; args...) = Boolean(size=(); args...)

Integ(; args...) = Var(T=Int; args...)
Integ0(; args...) = Integ(size=(); args...)

Str(; args...) = Var(T=String; args...)
Str0(; args...) = Str(size=(); args...)
=#

"""
Create a `Variable` with `parameter` variability, meaning it
is an input variable that is constant with time
"""
Parameter(; args...) = Variable(variability=parameter; args...)
Parameter(value; args...) = Variable(variability=parameter, value=value; args...)

"""
Shortcut for `Parameter`
"""
Par(; args...) = Variable(variability=parameter; args...)
Par(value; args...) = Variable(variability=parameter, value=value; args...)

"""
A value meant to be filled in later
"""
const undefined = nothing


text:
module SIMDIndexingMappedArrays

using Base: @propagate_inbounds

using ..SIMDIndexing: Index
using ..SIMDIndexing.MappedArrays
using ..SIMDIndexing.MappedArrays: _getindex, _setindex!,
    AbstractMultiMappedArray, MultiMappedArray

# IndexLinear implementations
@propagate_inbounds Base.getindex(A::AbstractMappedArray, i::Index) =
    A.f(A.data[i])
@propagate_inbounds Base.getindex(M::AbstractMultiMappedArray, i::Index) =
    M.f(_getindex(i, M.data...)...)

@propagate_inbounds function Base.setindex!(A::MappedArray{T},
                                            val,
                                            i::Index) where {T}
    A.data[i] = A.finv(convert(T, val)::T)
end
@propagate_inbounds function Base.setindex!(A::MultiMappedArray{T},
                                            val,
                                            i::Index) where {T}
    vals = A.finv(convert(T, val)::T)
    _setindex!(A.data, vals, i)
    return vals
end

# IndexCartesian implementations
@propagate_inbounds function Base.getindex(A::AbstractMappedArray{T,N},
                                           i::Vararg{Index,N}) where {T,N}
    A.f(A.data[i...])
end
@propagate_inbounds function Base.getindex(A::AbstractMultiMappedArray{T,N},
                                           i::Vararg{Index,N}) where {T,N}
    A.f(_getindex(CartesianIndex(i), A.data...)...)
end

@propagate_inbounds function Base.setindex!(A::MappedArray{T,N},
                                            val,
                                            i::Vararg{Index,N}) where {T,N}
    A.data[i...] = A.finv(convert(T, val)::T)
end
@propagate_inbounds function Base.setindex!(A::MultiMappedArray{T,N},
                                            val,
                                            i::Vararg{Index,N}) where {T,N}
    vals = A.finv(convert(T, val)::T)
    _setindex!(A.data, vals, i...)
    return vals
end

end # module

ext: jl


text:
using ModelingToolkit
using Graphs
using DiffEqBase
using Test
using UnPack

# Define some variables
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) xˍˍt(t) yˍˍt(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T*x,
        D(D(y)) ~ T*y - g,
        0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)
lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2)

lowered_eqs = [D(xˍt) ~ T*x,
               D(yˍt) ~ T*y - g,
               D(x) ~ xˍt,
               D(y) ~ yˍt,
               0 ~ x^2 + y^2 - L^2,]
@test ODESystem(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g], name=:pendulum) == lowered_sys
@test isequal(equations(lowered_sys), lowered_eqs)

# Simple pendulum in cartesian coordinates
eqs = [D(x) ~ w,
       D(y) ~ z,
       D(w) ~ T*x,
       D(z) ~ T*y - g,
       0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)

state = TearingState(pendulum)
@unpack graph, var_to_diff = state.structure
@test StructuralTransformations.maximal_matching(graph, eq->true, v->var_to_diff[v] === nothing) == map(x -> x == 0 ? StructuralTransformations.unassigned : x, [1, 2, 3, 4, 0, 0, 0, 0, 0])

using ModelingToolkit
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t)
D = Differential(t)
idx1_pendulum = [D(x) ~ w,
                 D(y) ~ z,
                 #0 ~ x^2 + y^2 - L^2,
                 D(w) ~ T*x,
                 D(z) ~ T*y - g,
                 # intermediate 1: 0 ~ 2x*D(x) + 2y*D(y) - 0,
                 # intermediate 2(a): 0 ~ 2x*w + 2y*z - 0, (substitute D(x) and D(y))
                 #0 ~ 2x*w + 2y*z,
                 # D(D(x)) ~ D(w) and substitute the rhs
                 D(xˍt) ~ T*x,
                 # D(D(y)) ~ D(z) and substitute the rhs
                 D(yˍt) ~ T*y - g,
                 # 2x*D(D(x)) + 2*D(x)*D(x) + 2y*D(D(y)) + 2*D(y)*D(y) and
                 # substitute the rhs
                 0 ~ 2x*(T*x) + 2*xˍt*xˍt + 2y*(T*y - g) + 2*yˍt*yˍt]
@named idx1_pendulum = ODESystem(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g])
first_order_idx1_pendulum = ode_order_lowering(idx1_pendulum)

using OrdinaryDiffEq
using LinearAlgebra
prob = ODEProblem(ODEFunction(first_order_idx1_pendulum),
        #  [x, y, w, z, xˍt, yˍt, T]
           [1, 0, 0, 0,   0,   0, 0.0],# 0, 0, 0, 0],
           (0, 10.0),
           [1, 9.8],
           mass_matrix=calculate_massmatrix(first_order_idx1_pendulum))
sol = solve(prob, Rodas5());
#plot(sol, vars=(1, 2))

new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2))

prob_auto = ODEProblem(new_sys,
                       [D(x)=>0,
                        D(y)=>0,
                        x=>1,
                        y=>0,
                        T=>0.0],
                       (0, 100.0),
                       [1, 9.8])
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(x, y))

# Define some variables
@parameters t L g
@variables x(t) y(t) T(t)
D = Differential(t)

eqs2 = [D(D(x)) ~ T*x,
        D(D(y)) ~ T*y - g,
        0 ~ x^2 + y^2 - L^2]
pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum)

# Turn into a first order differential equation system
first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2)

# Perform index reduction to get an Index 1 DAE
new_sys = dae_index_lowering(first_order_sys)

u0 = [
  D(x)    => 0.0,
  D(y)    => 0.0,
  x       => 1.0,
  y       => 0.0,
  T       => 0.0
]

p = [
    L => 1.0,
    g => 9.8
]

prob_auto = ODEProblem(new_sys,u0,(0.0,10.0),p)
sol = solve(prob_auto, Rodas5());
#plot(sol, vars=(D(x), y))

let pss_pendulum2 = partial_state_selection(pendulum2)
    # This currently selects `T` rather than `x` at top level. Needs tearing priorities to fix.
    @test_broken length(equations(pss_pendulum2)) == 3
    @test_broken length(equations(ModelingToolkit.ode_order_lowering(pss_pendulum2))) == 4
end

eqs = [D(x) ~ w,
       D(y) ~ z,
       D(w) ~ T*x,
       D(z) ~ T*y - g,
       0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)

let pss_pendulum = partial_state_selection(pendulum)
    # This currently selects `T` rather than `x` at top level. Needs tearing priorities to fix.
    @test_broken length(equations(pss_pendulum)) == 3
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
using EponymTuples
using EponymTuples: var_and_type # also test internals
using Test

@testset "internals" begin
    @test var_and_type(:Foo) == (:Foo => Any)
    @test var_and_type(:(bar::Int)) == (:bar => :Int)
    @test_throws ArgumentError var_and_type(:(1+1))
end

module TestModule               # wrap in a module to check escaping

using EponymTuples

foo(@eponymargs(a, b)) = (:generic, a, b)
foo(@eponymargs(a::Int, b::Int)) = (:ints, a, b)

struct Bar{T}
    x::T
end

bar(@eponymargs(z::Bar{T})) where {T <: Real} = z.x

pack1(a) = @eponymtuple(a, b = 2)
pack2(a, b) = @eponymtuple(a, b)

end

@testset "dispatch" begin
    TestModule.foo((a = 1, b = 2.0)) ≡ (:generic, 1, 2.0)
    TestModule.foo((a = 1, b = 2)) ≡ (:ints, 1, 2)
    @test_throws MethodError TestModule.foo((1, 2))
    @test_throws MethodError TestModule.foo((a = 1, c = 2))

    @test TestModule.bar((z = TestModule.Bar(2), )) == 2
    @test_throws MethodError TestModule.bar((z = TestModule.Bar("a fish"), ))
end

@testset "eponymtuple expansion" begin
    @test :((a = a, b = b)) == macroexpand(Main, :(@eponymtuple(a, b)))
    @test :((a = a, b = 2)) == macroexpand(Main, :(@eponymtuple(a, b = 2)))
    @test_throws LoadError macroexpand(Main, :(@eponymtuple(a, bogus::T)))
end

@testset "eponymtuple return values" begin
    @test TestModule.pack1(9.0) ≡ (a = 9.0, b = 2)
    @test TestModule.pack2(42, "a fish") ≡ (a = 42, b = "a fish")
end

@testset "non-unique varnames" begin
    # just checking that these are caught by Julia
    @test_throws ErrorException @eval let a=1; @eponymtuple(a, a) end
    @test_throws ErrorException @eval let a=1; @eponymtuple(a, a = 1) end
    # these pass through the language at the moment
    @test_throws LoadError @eval f(@eponymargs(a, a)) = 2*a
end


text:
using LinearAlgebra
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, HermOrSym

import LinearAlgebra: lmul!, mul!

export rankUpdate!

# Rank one update

## General
### BLAS
rankUpdate!(A::StridedMatrix{T}, x::StridedVector{T}, y::StridedVector{T}, α::T) where {T<:BlasReal} = BLAS.ger!(α, x, y, A)

### Generic
function rankUpdate!(A::StridedMatrix, x::StridedVector, y::StridedVector, α::Number)
    m, n = size(A, 1), size(A, 2)
    m == length(x) || throw(DimensionMismatch("x vector has wrong length"))
    n == length(y) || throw(DimensionMismatch("y vector has wrong length"))
    for j = 1:n
        yjc = y[j]'
        for i = 1:m
            A[i,j] += x[i]*α*yjc
        end
    end
end

# Deprecated 11 October 2018
Base.@deprecate rankUpdate!(α::Number, x::StridedVector, y::StridedVector, A::StridedMatrix) rankUpdate!(A, x, y, α)

## Hermitian
rankUpdate!(A::HermOrSym{T,S}, a::StridedVector{T}, α::T) where {T<:BlasReal,S<:StridedMatrix} = BLAS.syr!(A.uplo, α, a, A.data)
rankUpdate!(A::HermOrSym{T,S}, a::StridedVector{T}) where {T<:BlasReal,S<:StridedMatrix} = rankUpdate!(one(T), a, A)

### Generic
function rankUpdate!(A::Hermitian, a::StridedVector, α::Real)
    n = size(A, 1)
    n == length(a) || throw(DimensionMismatch("a vector has wrong length"))
    @inbounds for j in 1:n
        ajc = a[j]'
        for i in ((A.uplo == 'L') ? (j:n) : (1:j))
            A.data[i,j] += a[i]*α*ajc
        end
    end
    return A
end

# Deprecated 11 October 2018
Base.@deprecate rankUpdate!(α::Real, a::StridedVector, A::Hermitian) rankUpdate!(A, a, α)

# Rank k update
## Real
rankUpdate!(C::HermOrSym{T,S}, A::StridedMatrix{T}, α::T, β::T) where {T<:BlasReal,S<:StridedMatrix} = syrk!(C.uplo, 'N', α, A, β, C.data)

## Complex
rankUpdate!(C::Hermitian{T,S}, A::StridedMatrix{Complex{T}}, α::T, β::T) where {T<:BlasReal,S<:StridedMatrix} = herk!(C.uplo, 'N', α, A, β, C.data)

### Generic
function rankUpdate!(C::Hermitian, A::StridedVecOrMat, α::Real)
    n = size(C, 1)
    n == size(A, 1) || throw(DimensionMismatch("first dimension of A has wrong size"))
    @inbounds if C.uplo == 'L' # branch outside the loop to have larger loop to optimize
        for k in 1:size(A, 2)
            for j in 1:n
                ajkc = A[j,k]'
                for i in j:n
                    C.data[i,j] += A[i,k]*α*ajkc
                end
            end
        end
    else
        for k in 1:size(A, 2)
            for j in 1:n
                ajkc = A[j,k]'
                for i in 1:j
                    C.data[i,j] += A[i,k]*α*ajkc
                end
            end
        end
    end
    return C
end

# Deprecated 11 October 2018
Base.@deprecate rankUpdate!(α::Real, A::StridedVecOrMat, C::Hermitian) rankUpdate!(C, A, α)
Base.@deprecate rankUpdate!(α::Real, A::StridedVecOrMat, β::Real, C::Hermitian) rankUpdate!(C, A, α, β)

# BLAS style mul!
## gemv
mul!(y::StridedVector{T}, A::StridedMatrix{T}, x::StridedVector{T}, α::T, β::T) where {T<:BlasFloat} = gemv!('N', α, A, x, β, y)
mul!(y::StridedVector{T}, A::Adjoint{T,<:StridedMatrix{T}}, x::StridedVector{T}, α::T, β::T) where {T<:BlasFloat} = gemv!('C', α, parent(adjA), x, β, y)

## gemm
mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasFloat} = BLAS.gemm!('N', 'N', α, A, B, β, C)
mul!(C::StridedMatrix{T}, adjA::Adjoint{T,<:StridedMatrix{T}}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasFloat} = BLAS.gemm!('C', 'N', α, parent(adjA), B, β, C)
# Not optimized since it is a generic fallback. Can probably soon be removed when the signatures in base have been updated.
function mul!(C::StridedVecOrMat,
              A::StridedMatrix,
              B::StridedVecOrMat,
              α::Number,
              β::Number)

    m, n = size(C, 1), size(C, 2)
    k = size(A, 2)

    if β != 1
        if β == 0
            fill!(C, 0)
        else
            rmul!(C, β)
        end
    end
    for j = 1:n
        for i = 1:m
            for l = 1:k
                C[i,j] += α*A[i,l]*B[l,j]
            end
        end
    end
    return C
end
function mul!(C::StridedVecOrMat,
              adjA::Adjoint{<:Number,<:StridedMatrix},
              B::StridedVecOrMat,
              α::Number,
              β::Number)

    A = parent(adjA)
    m, n = size(C, 1), size(C, 2)
    k = size(A, 1)

    if β != 1
        if β == 0
            fill!(C, 0)
        else
            rmul!(C, β)
        end
    end
    for j = 1:n
        for i = 1:m
            for l = 1:k
                C[i,j] += α*A[l,i]'*B[l,j]
            end
        end
    end
    return C
end

## trmm like
### BLAS versions
mul!(A::UpperTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'U', 'N', 'N', α, A.data, B)
mul!(A::LowerTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'L', 'N', 'N', α, A.data, B)
mul!(A::UnitUpperTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'U', 'N', 'U', α, A.data, B)
mul!(A::UnitLowerTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'L', 'N', 'U', α, A.data, B)
mul!(A::Adjoint{T,UpperTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'U', 'C', 'N', α, parent(A).data, B)
mul!(A::Adjoint{T,LowerTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'L', 'C', 'N', α, parent(A).data, B)
mul!(A::Adjoint{T,UnitUpperTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'U', 'C', 'U', α, parent(A).data, B)
mul!(A::Adjoint{T,UnitLowerTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:BlasFloat,S} = trmm!('L', 'L', 'C', 'U', α, parent(A).data, B)

### Generic fallbacks
function lmul!(A::UpperTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = A.data
    m, n = size(B)
    for i = 1:m
        for j = 1:n
            B[i,j] = α*AA[i,i]*B[i,j]
            for l = i + 1:m
                B[i,j] += α*AA[i,l]*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::LowerTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = A.data
    m, n = size(B)
    for i = m:-1:1
        for j = 1:n
            B[i,j] = α*AA[i,i]*B[i,j]
            for l = 1:i - 1
                B[i,j] += α*AA[i,l]*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::UnitUpperTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = A.data
    m, n = size(B)
    for i = 1:m
        for j = 1:n
            B[i,j] = α*B[i,j]
            for l = i + 1:m
                B[i,j] = α*AA[i,l]*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::UnitLowerTriangular{T,S}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = A.data
    m, n = size(B)
    for i = m:-1:1
        for j = 1:n
            B[i,j] = α*B[i,j]
            for l = 1:i - 1
                B[i,j] += α*AA[i,l]*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::Adjoint{T,UpperTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = parent(A).data
    m, n = size(B)
    for i = m:-1:1
        for j = 1:n
            B[i,j] = α*AA[i,i]*B[i,j]
            for l = 1:i - 1
                B[i,j] += α*AA[l,i]'*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::Adjoint{T,LowerTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = parent(A).data
    m, n = size(B)
    for i = 1:m
        for j = 1:n
            B[i,j] = α*AA[i,i]*B[i,j]
            for l = i + 1:m
                B[i,j] += α*AA[l,i]'*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::Adjoint{T,UnitUpperTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = parent(A).data
    m, n = size(B)
    for i = m:-1:1
        for j = 1:n
            B[i,j] = α*B[i,j]
            for l = 1:i - 1
                B[i,j] += α*AA[l,i]'*B[l,j]
            end
        end
    end
    return B
end
function lmul!(A::Adjoint{T,UnitLowerTriangular{T,S}}, B::StridedMatrix{T}, α::T) where {T<:Number,S}
    AA = parent(A).data
    m, n = size(B)
    for i = 1:m
        for j = 1:n
            B[i,j] = α*B[i,j]
            for l = i + 1:m
                B[i,j] = α*AA[l,i]'*B[l,j]
            end
        end
    end
    return B
end


text:
mutable struct NeuralNetworkRegressor{B,O,L} <: MLJModelInterface.Deterministic
    builder::B
    optimiser::O    # mutable struct from Flux/src/optimise/optimisers.jl
    loss::L         # can be called as in `loss(yhat, y)`
    epochs::Int     # number of epochs
    batch_size::Int # size of a batch
    lambda::Float64 # regularization strength
    alpha::Float64  # regularizaton mix (0 for all l2, 1 for all l1)
    optimiser_changes_trigger_retraining::Bool
    acceleration::AbstractResource       # To use GPU
end

function NeuralNetworkRegressor(; builder::B   = Linear()
                                , optimiser::O = Flux.Optimise.ADAM()
                                , loss::L      = Flux.mse
                                , epochs       = 10
                                , batch_size   = 1
                                , lambda       = 0
                                , alpha        = 0
                                , optimiser_changes_trigger_retraining=false
                                , acceleration  = CPU1()
                                ) where {B,O,L}

    model = NeuralNetworkRegressor{B,O,L}(builder
                                          , optimiser
                                          , loss
                                          , epochs
                                          , batch_size
                                          , lambda
                                          , alpha
                                          , optimiser_changes_trigger_retraining
                                          , acceleration)

   message = clean!(model)
   isempty(message) || @warn message

    return model
end

mutable struct MultitargetNeuralNetworkRegressor{B,O,L} <: MLJModelInterface.Deterministic
    builder::B
    optimiser::O    # mutable struct from Flux/src/optimise/optimisers.jl
    loss::L         # can be called as in `loss(yhat, y)`
    epochs::Int          # number of epochs
    batch_size::Int # size of a batch
    lambda::Float64 # regularization strength
    alpha::Float64  # regularizaton mix (0 for all l2, 1 for all l1)
    optimiser_changes_trigger_retraining::Bool
    acceleration::AbstractResource
end

function MultitargetNeuralNetworkRegressor(; builder::B   = Linear()
                                           , optimiser::O = Flux.Optimise.ADAM()
                                           , loss::L      = Flux.mse
                                           , epochs       = 10
                                           , batch_size   = 1
                                           , lambda       = 0
                                           , alpha        = 0
                                           , optimiser_changes_trigger_retraining=false
                                           , acceleration = CPU1()
                                           ) where {B,O,L}

    model = MultitargetNeuralNetworkRegressor{B,O,L}(builder
                                                     , optimiser
                                                     , loss
                                                     , epochs
                                                     , batch_size
                                                     , lambda
                                                     , alpha
                                                     , optimiser_changes_trigger_retraining
                                                     , acceleration)

    message = clean!(model)
    isempty(message) || @warn message

    return model
end

const Regressor =
    Union{NeuralNetworkRegressor, MultitargetNeuralNetworkRegressor}

function MLJModelInterface.fit(model::Regressor, verbosity::Int, X, y)

    # (assumes  no categorical features)
    n_input = Tables.schema(X).names |> length
    data = MLJFlux.collate(model, X, y)

    target_is_multivariate = Tables.istable(y)
    if target_is_multivariate
        target_column_names = Tables.schema(y).names
    else
        target_column_names = [""] # We won't be using this
    end

    n_output = length(target_column_names)
    chain = build(model.builder, n_input, n_output)

    optimiser = deepcopy(model.optimiser)

    chain, history = fit!(chain,
                          optimiser,
                          model.loss,
                          model.epochs,
                          model.lambda,
                          model.alpha,
                          verbosity,
                          #data,
                          model.acceleration,
                          data[1],
                          data[2])

    # note: "state" part of `optimiser` is now mutated!

    cache = (deepcopy(model), data, history, n_input, n_output, optimiser)
    fitresult = (chain, target_is_multivariate, target_column_names)
    report = (training_losses=history,)

    return fitresult, cache, report

end

function MLJModelInterface.update(model::Regressor,
                                  verbosity::Int,
                                  old_fitresult,
                                  old_cache,
                                  X,
                                  y)

    # note: the `optimiser` in `old_cache` stores "state" (eg,
    # momentum); the "state" part of the `optimiser` field of `model`
    # and of `old_model` play no role

    old_model, data, old_history, n_input, n_output, optimiser = old_cache
    old_chain, target_is_multivariate, target_column_names = old_fitresult

    optimiser_flag = model.optimiser_changes_trigger_retraining &&
        model.optimiser != old_model.optimiser

    keep_chain = !optimiser_flag && model.epochs >= old_model.epochs &&
        MLJModelInterface.is_same_except(model, old_model, :optimiser, :epochs)

    if keep_chain
        chain = old_chain
        epochs = model.epochs - old_model.epochs
    else
        chain = build(model.builder, n_input, n_output)
        data = collate(model, X, y)
        epochs = model.epochs
    end

    # we only get to keep the optimiser "state" carried over from
    # previous training if we're doing a warm restart and the user has not
    # changed the optimiser hyper-parameter:
    if !keep_chain ||
        !MLJModelInterface._equal_to_depth_one(model.optimiser,
                                              old_model.optimiser)
        optimiser = deepcopy(model.optimiser)
    end

    chain, history = fit!(chain,
                          optimiser,
                          model.loss,
                          epochs,
                          model.lambda,
                          model.alpha,
                          verbosity,
                          #data,
                          model.acceleration,
                          data[1],
                          data[2])
    if keep_chain
        # note: history[1] = old_history[end]
        history = vcat(old_history[1:end-1], history)
    end

    fitresult = (chain, target_is_multivariate, target_column_names)
    cache = (deepcopy(model), data, history, n_input, n_output, optimiser)
    report = (training_losses=history,)

    return fitresult, cache, report

end

function MLJModelInterface.predict(model::Regressor, fitresult, Xnew_)

    chain , target_is_multivariate, target_column_names = fitresult

    Xnew_ = MLJModelInterface.matrix(Xnew_)

    if target_is_multivariate
        ypred = [chain(values.(tomat(Xnew_[i, :])))
                 for i in 1:size(Xnew_, 1)]
        return MLJModelInterface.table(reduce(hcat, y for y in ypred)',
                                       names=target_column_names)
    else
        return [chain(values.(tomat(Xnew_[i, :])))[1]
                for i in 1:size(Xnew_, 1)]
    end
end

MLJModelInterface.fitted_params(::Regressor, fitresult) = (chain=fitresult[1],)

MLJModelInterface.metadata_model(NeuralNetworkRegressor,
               input=Table(Continuous),
               target=AbstractVector{<:Continuous},
               path="MLJFlux.NeuralNetworkRegressor",
               descr="A neural network model for making "*
                     "deterministic predictions of a "*
                     "`Continuous` target, given a table of "*
                     "`Continuous` features. ")

MLJModelInterface.metadata_model(MultitargetNeuralNetworkRegressor,
               input=Table(Continuous),
               target=Table(Continuous),
               path="MLJFlux.MultitargetNeuralNetworkRegressor",
               descr = "A neural network model for making "*
                       "deterministic predictions of a "*
                       "`Continuous` multi-target, presented "*
                       "as a table, given a table of "*
                       "`Continuous` features. ")


text:
"""
Convert a graph game to a simple game based on Myerson (1977). Namely, assume v(S) > 0 iff S is connected subgraph.
"""
function graphtoMyerson(G::GraphGame)::SimpleGame
    function vᴸ(S::Players)
        subgraph, vmap = induced_subgraph(G.g, S)
        comps = [
            [vmap[i] for i in comp]
            for comp in connected_components(subgraph)]

        return ∑(G.v(T) for T in comps)
    end

    return SimpleGame(G.N, vᴸ)
end


text:
include("clang/gdal_cdefs.jl")
include("clang/gdal_cfuns.jl")

include("gdal/type.jl")
include("gdal/utilities.jl")
include("gdal/drivers.jl")
include("gdal/fileio.jl") # depends: drivers.jl
include("gdal/dataset.jl")
include("gdal/rasterband.jl") # depends: gdal_type.jl
include("gdal/gcp.jl")
include("gdal/color.jl")
include("gdal/statistics.jl")


text:
include("logistic.jl")
include("bayesiansvm.jl")

const ClassificationLikelihood = BernoulliLikelihood

function (l::BernoulliLikelihood)(y::Real, f::Real)
    return pdf(l(f), y)
end

function init_local_vars(::BernoulliLikelihood, batchsize::Int, T::DataType=Float64)
    return (; c=rand(T, batchsize), θ=zeros(T, batchsize))
end

function compute_proba(
    l::BernoulliLikelihood{L}, μ::AbstractVector{T}, σ²::AbstractVector{T}
) where {L,T<:Real}
    N = length(μ)
    pred = zeros(T, N)
    σ²_pred = zeros(T, N)
    for i in 1:N
        x = pred_nodes .* sqrt(max(σ²[i], zero(T))) .+ μ[i]
        pred[i] = dot(pred_weights, l.invlink.(x))
        σ²_pred[i] = max(dot(pred_weights, abs2.(l.invlink.(x))) - pred[i]^2, zero(T))
    end
    return pred, σ²_pred
end

# Return the labels in a vector of vectors for multiple outputs
function treat_labels!(y::AbstractVector{<:Real}, ::BernoulliLikelihood)
    labels = unique(y)
    # y isa AbstractVector{<:Union{Int,Bool}} || error("y labels should be Int")
    if sort(Int64.(labels)) == [0; 1]
        return (y .- 0.5) * 2
    elseif sort(Int64.(labels)) == [-1; 1]
        return y
    else
        throw(ArgumentError("Labels of y should be binary {-1,1} or {0,1}"))
    end
end

function treat_labels!(::AbstractVector, ::BernoulliLikelihood)
    return error(
        "For classification target(s) should be real valued (Bool, Integer or Float)"
    )
end

predict_y(::BernoulliLikelihood, μ::AbstractVector{<:Real}) = μ .> 0
predict_y(::BernoulliLikelihood, μ::Tuple{<:AbstractVector}) = first(μ) .> 0


text:
# OpenCL.Kernel

type Kernel <: CLObject
    id :: CL_kernel

    function Kernel(k::CL_kernel, retain=false)
        if retain
            @check api.clRetainKernel(k)
        end
        kernel = new(k)
        finalizer(kernel, _finalize)
        return kernel
    end
end

function _finalize(k::Kernel)
    if k.id != C_NULL
        @check api.clReleaseKernel(k.id)
        k.id = C_NULL
    end
end

Base.pointer(k::Kernel) = k.id

Base.show(io::IO, k::Kernel) = begin
    print(io, "OpenCL.Kernel(\"$(k[:name])\" nargs=$(k[:num_args]))")
end

Base.getindex(k::Kernel, kinfo::Symbol) = info(k, kinfo)

function Kernel(p::Program, kernel_name::String)
    for (dev, status) in info(p, :build_status)
        if status != CL_BUILD_SUCCESS
            msg = "OpenCL.Program has to be built before Kernel constructor invoked"
            throw(ArgumentError(msg))
        end
    end
    err_code = Array(CL_int, 1)
    kernel_id = api.clCreateKernel(p.id, kernel_name, err_code)
    if err_code[1] != CL_SUCCESS
        throw(CLError(err_code[1]))
    end
    return Kernel(kernel_id)
end

immutable LocalMem{T}
    nbytes::Csize_t
end

LocalMem{T}(::Type{T}, len::Integer) = begin
    @assert len > 0
    nbytes = sizeof(T) * len
    return LocalMem{T}(convert(Csize_t, nbytes))
end

Base.ndims(l::LocalMem) = 1
Base.eltype{T}(l::LocalMem{T}) = T
Base.sizeof{T}(l::LocalMem{T}) = l.nbytes
Base.length{T}(l::LocalMem{T}) = Int(l.nbytes ÷ sizeof(T))

function set_arg!(k::Kernel, idx::Integer, arg::Void)
    @assert idx > 0
    @check api.clSetKernelArg(k.id, cl_uint(idx-1), sizeof(CL_mem), C_NULL)
    return k
end

function set_arg!(k::Kernel, idx::Integer, arg::Ptr{Void})
    if arg != C_NULL
        throw(AttributeError("set_arg! for void pointer $arg is undefined"))
    end
    set_arg!(k, idx, nothing)
end

function set_arg!(k::Kernel, idx::Integer, arg::CLMemObject)
    @assert idx > 0
    @check api.clSetKernelArg(k.id, cl_uint(idx-1), sizeof(CL_mem), [arg.id])
    return k
end

function set_arg!(k::Kernel, idx::Integer, arg::LocalMem)
    @assert idx > 0
    @check api.clSetKernelArg(k.id, cl_uint(idx-1), arg.nbytes, C_NULL)
    return k
end

#TODO: vector types...
#TODO: type safe calling of set args for kernel (with clang)

# set scalar/vector kernel args
for cl_type in [:CL_char, :CL_uchar, :CL_short, :CL_ushort,
                :CL_int,  :CL_uint,  :CL_long,  :CL_ulong,
                :CL_half, :CL_float, :CL_double]
    @eval begin
        function set_arg!(k::Kernel, idx::Integer, arg::$cl_type)
            @assert idx > 0
            boxed_arg = $cl_type[arg,]
            @check api.clSetKernelArg(k.id, cl_uint(idx-1), sizeof($cl_type), boxed_arg)
            return k
        end
    end
end

function set_args!(k::Kernel, args...)
    for (i, a) in enumerate(args)
        set_arg!(k, i, a)
    end
end

function work_group_info(k::Kernel, winfo::CL_kernel_work_group_info, d::Device)
    if (winfo == CL_KERNEL_LOCAL_MEM_SIZE ||
        winfo == CL_KERNEL_PRIVATE_MEM_SIZE)
        result = CL_ulong[0]
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo,
                                            sizeof(CL_ulong), result, C_NULL)
        return Int(result[1])
    elseif winfo == CL_KERNEL_COMPILE_WORK_GROUP_SIZE
        size = Csize_t[0]
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo, 0, C_NULL, size)
        result = Array(Csize_t, size[1])
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo, sizeof(result), result, C_NULL)
        return map(Int, result)
    else
        result = Csize_t[0]
        @check api.clGetKernelWorkGroupInfo(k.id, d.id, winfo,
                                            sizeof(CL_ulong), result, C_NULL)
        return Int(result[1])
    end
end

function work_group_info(k::Kernel, winfo::Symbol, d::Device)
    if winfo == :size
        work_group_info(k, CL_KERNEL_WORK_GROUP_SIZE, d)
    elseif winfo == :compile_size
        work_group_info(k, CL_KERNEL_COMPILE_WORK_GROUP_SIZE, d)
    elseif winfo == :local_mem_size
        work_group_info(k, CL_KERNEL_LOCAL_MEM_SIZE, d)
    elseif winfo == :private_mem_size
        work_group_info(k, CL_KERNEL_PRIVATE_MEM_SIZE, d)
    elseif winfo == :prefered_size_multiple
        work_group_info(k, CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE, d)
    else
        throw(ArgumentError(("Unknown work_group_info flag: :$winfo")))
    end
end

# produce a cl.call thunk with kernel queue, global/local sizes
Base.getindex(k::Kernel, args...) = begin
    if length(args) < 2 || length(args) > 3
        throw(ArgumentError("kernel must be called with a queue & global size as arguments"))
    end
    if !(isa(args[1], CmdQueue))
        throw(ArgumentError("kernel first argument must a a CmdQueue"))
    end
    if !(isa(args[2], Dims)) || length(args[2]) > 3
        throw(ArgumentError("kernel global size must be of Dims type (dim <= 3)"))
    end
    if length(args) == 3 && (!(isa(args[3], Dims)) || length(args[3]) > 3)
        throw(ArgumentError("kernel local size must be of Dims type (dim <= 3)"))
    end
    queue = args[1]
    global_size = args[2]
    local_size  = length(args) == 3 ? args[3] : nothing
    # TODO: we cannot pass keywords in anon functions yet, return kernel call thunk
    return (args...) -> queue(k, global_size, local_size, args...)
end

# blocking kernel call that finishes queue
function (q::CmdQueue)(k::Kernel, global_work_size, local_work_size,
                      args...; global_work_offset=nothing,
                      wait_on::Union{Void,Vector{Event}}=nothing)
    set_args!(k, args...)
    evt = enqueue_kernel(q, k,
                         global_work_size,
                         local_work_size,
                         global_work_offset=global_work_offset,
                         wait_on=wait_on)
    finish(q)
    return evt
end

function enqueue_kernel(q::CmdQueue, k::Kernel, global_work_size)
    enqueue_kernel(q, k, global_work_size, nothing)
end

function enqueue_kernel(q::CmdQueue,
                                k::Kernel,
                                global_work_size,
                                local_work_size;
                                global_work_offset=nothing,
                                wait_on::Union{Void,Vector{Event}}=nothing)
    device = q[:device]
    max_work_dim = device[:max_work_item_dims]
    work_dim     = length(global_work_size)
    if work_dim > max_work_dim
        throw(ArgumentError("global_work_size has max dim of $max_work_dim"))
    end
    gsize = Array(Csize_t, work_dim)
    for (i, s) in enumerate(global_work_size)
        gsize[i] = s
    end

    goffset = C_NULL
    if global_work_offset !== nothing
        if length(global_work_offset) > max_work_dim
            throw(ArgumentError("global_work_offset has max dim of $max_work_dim"))
        end
        if length(global_work_offset) != work_dim
            throw(ArgumentError("global_work_size and global_work_offset have differing dims"))
        end
        goffset = Array(Csize_t, work_dim)
        for (i, o) in enumerate(global_work_offset)
            goffset[i] = o
        end
    end

    lsize = C_NULL
    if local_work_size !== nothing
        if length(local_work_size) > max_work_dim
            throw(ArgumentError("local_work_offset has max dim of $max_work_dim"))
        end
        if length(local_work_size) != work_dim
            throw(ArgumentError("global_work_size and local_work_size have differing dims"))
        end
        lsize = Array(Csize_t, work_dim)
        for (i, s) in enumerate(local_work_size)
            lsize[i] = s
        end
    end

    if wait_on !== nothing
        n_events = cl_uint(length(wait_on))
        wait_event_ids = [evt.id for evt in wait_on]
    else
        n_events = cl_uint(0)
        wait_event_ids = C_NULL
    end

    ret_event = Array(CL_event, 1)
    @check api.clEnqueueNDRangeKernel(q.id, k.id, cl_uint(work_dim), goffset, gsize, lsize,
                                      n_events, wait_event_ids, ret_event)
    return Event(ret_event[1], retain=false)
end

function enqueue_task(q::CmdQueue, k::Kernel; wait_for=nothing)
    n_evts  = 0
    evt_ids = C_NULL
    #TODO: this should be split out into its own function
    if wait_for !== nothing
        if isa(wait_for, Event)
            n_evts = 1
            evt_ids = [wait_for.id]
        else
            @assert all([isa(evt, Event) for evt in wait_for])
            n_evts = length(wait_for)
            evt_ids = [evt.id for evt in wait_for]
        end
    end
    ret_event = Array(CL_event, 1)
    @check api.clEnqueueTask(q.id, k.id, n_evts, evt_ids, ret_event)
    return ret_event[1]
end

let name(k::Kernel) = begin
        size = Ref{Csize_t}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_FUNCTION_NAME,
                                   0, C_NULL, size)
        result = Array(Cchar, size[])
        @check api.clGetKernelInfo(k.id, CL_KERNEL_FUNCTION_NAME,
                                   size[], result, size)
        return CLString(result)
    end

    num_args(k::Kernel) = begin
        ret = Ref{CL_uint}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_NUM_ARGS,
                                   sizeof(CL_uint), ret, C_NULL)
        return ret[]
    end

    reference_count(k::Kernel) = begin
        ret = Ref{CL_uint}()
        @check api.clGetKernelInfo(k.id, CL_KERNEL_REFERENCE_COUNT,
                                   sizeof(CL_uint), ret, C_NULL)
        return ret[]
    end

    program(k::Kernel) = begin
        ret = Array(CL_program, 1)
        @check api.clGetKernelInfo(k.id, CL_KERNEL_PROGRAM,
                                   sizeof(CL_program), ret, C_NULL)
        return Program(ret[1], retain=true)
    end

    attributes(k::Kernel) = begin
        size = Ref{Csize_t}()
        api.clGetKernelInfo(k.id, CL_KERNEL_ATTRIBUTES,
                            0, C_NULL, size)
        if size[] <= 1
            return ""
        end
        result = Array(CL_char, size[])
        @check api.clGetKernelInfo(k.id, CL_KERNEL_ATTRIBUTES,
                                   size[], result, size)
        return CLString(result)
    end

    const info_map = Dict{Symbol, Function}(
        :name => name,
        :num_args => num_args,
        :reference_count => reference_count,
        :program => program,
        :attributes => attributes
    )

    function info(k::Kernel, kinfo::Symbol)
        try
            func = info_map[kinfo]
            func(k)
        catch err
            if isa(err, KeyError)
                error("OpenCL.Kernel has no info for: $kinfo")
            else
                throw(err)
            end
        end
    end
end

#TODO set_arg sampler...
# OpenCL 1.2 function
#TODO: get_arg_info(k::Kernel, idx, param)
#TODO: enqueue_async_kernel()


text:
using Odoo
using Base.Test

# write your own tests here
demo = Odoo.demo_server()
@show demo
@show Odoo.version(demo)

@show Odoo.execute_kw(demo,
    "res.partner", "check_access_rights",
    ["read"], Dict("raise_exception" => false))


text:
module plugin

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Compat
using Compat.Test # Base.Test on 0.6, and Test on 0.7
using Compat.LinearAlgebra # dot

using Einsum

## Test that vars in Main aren't overwritten by einsum
let
    i = -1
    y = randn(10)
    @einsum x[i] := y[i]
    @test i == -1
end

## Test that B is overwritten by := operator
let
    B = randn(10, 10)
    A = randn(5, 10)
    @einsum B[i, j] := A[i, j] # this should run without a problem
    @test size(B) == size(A)
end

## CP decomposition test case ##
let

    # preallocated test case
    A = zeros(5, 6, 7)
    B = similar(A)
    C = similar(A)

    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)

    @einsum A[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    @einsimd B[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    @vielsum C[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]

    for i = 1:5
        for j = 1:6
            for k = 1:7
                s = 0.0
                for r = 1:2
                    s += X[i, r] * Y[j, r] * Z[k, r]
                end
                @test isapprox(A[i, j, k], s)
                @test isapprox(B[i, j, k], s)
                @test isapprox(C[i, j, k], s)
            end
        end
    end

    # without preallocation
    @einsum A2[i, j, k] := X[i, r] * Y[j, r] * Z[k, r]
    @test isapprox(A, A2)

end

# Interesting test case, can throw an error that
# local vars are declared twice.
let
    A = zeros(5, 6, 7)
    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)
    if true
        @einsum A[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    else
        @einsum A[i, j, k] = X[i, r] * Y[j, r] * Z[k, r]
    end
end

# From #21: local `T` does not interfer with internal T
let
    function test(x::Vector{T}, y::Vector{T}) where T
        @einsum z := x[i] * y[i]
        return z
    end
    @test_nowarn test(rand(3), rand(3))
end
# From #20: local `s` does not interfere with internal s
let
    x = rand(2, 3)
    @test_nowarn @einsum y[i] := x[i, s]
end

# At one point this threw an error because the lhs
# had no indices/arguments
let
    x = randn(10)
    y = randn(10)
    @einsum k := x[i] * y[i]
    @test isapprox(k, dot(x, y))
end

# Elementwise multiplication (this should create nested loops with no
# no summation.)
let
    x = randn(10)
    y = randn(10)
    @einsum k[i] := x[i] * y[i]
    @einsimd k2[i] := x[i] * y[i]
    @vielsum k3[i] := x[i] * y[i]
    @test isapprox(k, x .* y)
    @test isapprox(k2, x .* y)
    @test isapprox(k3, x .* y)
end

# Transpose a block matrix
let
    z = [rand(2, 2) for i = 1:2, j = 1:2]
    @einsum t[i, j] := transpose(z[j, i])
    @test isapprox(z[1, 1],  t[1, 1]')
    @test isapprox(z[2, 2],  t[2, 2]')
    @test isapprox(z[1, 2],  t[2, 1]')
    @test isapprox(z[2, 1],  t[1, 2]')
end

# Mapping functions
let
    A = randn(10, 10)
    @einsum B[i, j] := exp(A[i, j])
    @test isapprox(exp.(A), B)
end

# Example from numpy
let
    A = reshape(collect(1:25), 5, 5)
    @einsum B[i] := A[i, i]
    @test all(B .== [1, 7, 13, 19, 25])
end

# TODO: consider adding support for this:
# @einsum A[i, j] = A[i,j] + 50

## Test in-place operations ##
let
    A = randn(5, 6, 7)
    B = randn(5, 6, 7)
    A1 = copy(A)
    B1 = copy(B)

    X = randn(5, 2)
    Y = randn(6, 2)
    Z = randn(7, 2)

    @einsum A[i, j, k] += X[i, r] * Y[j, r] * Z[k, r]
    @einsimd B[i, j, k] += X[i, r] * Y[j, r] * Z[k, r]

    for i = 1:5
        for j = 1:6
            for k = 1:7
                s = 0.0
                for r = 1:2
                    s += X[i, r] * Y[j, r] * Z[k, r]
                end
                @test isapprox(A[i, j, k], A1[i, j, k] + s)
                @test isapprox(B[i, j, k], B1[i, j, k] + s)
            end
        end
    end

    x = randn(10)
    y = randn(10)
    k0 = randn()
    k = k0
    @einsum k += x[i] * y[i]
    @test isapprox(k, k0 + dot(x, y))

    # test multiplication

    A1[:] = A[:]
    B1[:] = B[:]

    @einsum A[i, j, k]  *= X[i, r] * Y[j, r] * Z[k, r]
    @einsimd B[i, j, k] *= X[i, r] * Y[j, r] * Z[k, r]

    for i = 1:5
        for j = 1:6
            for k = 1:7
                s = 0.0
                for r = 1:2
                    s += X[i, r] * Y[j, r] * Z[k, r]
                end
                @test isapprox(A[i, j, k], A1[i, j, k] * s)
                @test isapprox(B[i, j, k], B1[i, j, k] * s)
            end
        end
    end

    x = randn(10)
    y = randn(10)
    k0 = randn()
    k = k0
    @einsum k *= x[i] * y[i]
    @test isapprox(k, k0 * dot(x, y))
end

# Test offsets
let
    X = randn(10)

    # without preallocation
    @einsum A[i] := X[i + 5]
    @test size(A) == (5,)
    @test all(A .== X[6:end])

    # with preallocation
    B = zeros(10)
    @einsum B[i] = X[i + 5]
    @test size(B) == (10,)
    @test all(B[1:5] .== X[6:end])
end

# Test symbolic offsets
let
    offset = 5
    X = randn(10)

    # without preallocation
    # @einsum A[i] := X[i + :offset] # error on 1.0
    # @test size(A) == (5,)
    # @test all(A .== X[6:end])

    # with preallocation
    B = zeros(10)
    # @einsum B[i] = X[i + :offset] # error on 1.0
    # @test size(B) == (10,)
    # @test all(B[1:5] .== X[6:end])
end

# Test adding/subtracting constants
let
    k = 5
    X = randn(10)

    # without preallocation
    @einsum A[i] := X[i] + k
    @einsum B[i] := X[i] - k
    @test isapprox(A, X .+ k)
    @test isapprox(B, X .- k)

    # with preallocation
    C, D = zeros(10), zeros(10)
    @einsum C[i] = X[i] + k
    @einsum D[i] = X[i] - k
    @test isapprox(C, X .+ k)
    @test isapprox(D, X .- k)
end

# Test multiplying/dividing constants
let
    k = 5
    X = randn(10)

    # without preallocation
    @einsum A[i] := X[i] * k
    @einsum B[i] := X[i] / k
    @test isapprox(A, X .* k)
    @test isapprox(B, X ./ k)

    # with preallocation
    C, D = zeros(10), zeros(10)
    @einsum C[i] = X[i] * k
    @einsum D[i] = X[i] / k
    @test isapprox(C, X .* k)
    @test isapprox(D, X ./ k)
end

# Test indexing with a constant
let
    A = randn(10, 2)
    j = 2
    # @einsum B[i] := A[i, :j] # error on 1.0
    # @test all(B .== A[:, j])
    @einsum C[i] := A[i, 1]
    @test all(C .== A[:, 1])

    D = zeros(10, 3)
    # @einsum D[i, 1] = A[i, :j]
    # @test isapprox(D[:, 1], A[:, j])
    # @einsum D[i, :j] = A[i, :j]
    # @test isapprox(D[:, j], A[:, j])
end

# Better type inference on allocating arrays
let
    B1 = ones(Int, 5)
    B2 = ones(Float32, 5)
    B3 = ones(5)
    C = randn(5)
    @einsum A1[i, j] := B1[i] * C[j]
    @einsum A2[i, j] := B2[i] * C[j]
    @einsum A3[i, j] := B3[i] * C[j]

    @test eltype(A1) == Float64
    @test eltype(A2) == Float64
    @test eltype(A3) == Float64
    @test isapprox(A1, A3)
    @test isapprox(A2, A3)
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.

include("../src/sudoku.jl") # Defines Sudoku module

# Puzzle example from
# https://en.wikipedia.org/wiki/Sudoku accessed 1/1/2022
# Puzzle image license is CC0 (Public Domain Dedication) per
# https://en.wikipedia.org/wiki/Sudoku#/media/File:Sudoku_Puzzle_by_L2G-20050714_standardized_layout.svg
puzzle_values = [
    5 3 0 0 7 0 0 0 0;
    6 0 0 1 9 5 0 0 0;
    0 9 8 0 0 0 0 6 0;
    8 0 0 0 6 0 0 0 3;
    4 0 0 8 0 3 0 0 1;
    7 0 0 0 2 0 0 0 6;
    0 6 0 0 0 0 2 8 0;
    0 0 0 4 1 9 0 0 5;
    0 0 0 0 8 0 0 7 9]
puzzle = Sudoku.SolvablePuzzle(3)
Sudoku.assign_values!(puzzle,puzzle_values)
print("\nUnsolved puzzle\n")
print(Sudoku.as_text_grid(puzzle))
iterations, uncertainty = Sudoku.iterative_solve!(puzzle,1)
print("\n\nSolution after $iterations steps\n")
print(Sudoku.as_text_grid(puzzle))

# Solution image is CC BY-SA 3.0 (Attribution-ShareAlike) per
# https://en.wikipedia.org/wiki/Sudoku#/media/File:Sudoku_Puzzle_by_L2G-20050714_solution_standardized_layout.svg
#
# It is an exercise for reader to compare the calculated solution to the reference
if uncertainty != 0
    println("Failure: puzzle solution is incomplete")
else
    result_values = Sudoku.as_values(puzzle)
    if Sudoku.satisfies(puzzle_values,result_values)
        println("Verified that solution satisfies initial puzzle")
    else
        println("Solution does not match original puzzle inputs")
    end
    if Sudoku.valid_puzzle(result_values)
        println("Verified that solution satisfies Suduko rules")
    else
        println("Solution is not a valid Sudoku")
    end
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
# mri/z-list.jl

include("kspace.jl")
include("sensemap-sim.jl")


text:
"""
`Fields = Dict{String, String}`. Stores the fields `name => value` of an entry.
"""
const Fields = Dict{String,String}

"""
    const entries = [
        :article,
        :book,
        :booklet,
        :inbook,
        :incollection,
        :inproceedings,
        :manual,
        :masterthesis,
        :misc,
        :phdthesis,
        :proceedings,
        :techreport,
        :unpublished,
    ]
List of possible entries (currently based on bibtex). Keep it sorted for readability.
"""
const entries = [
    :article,
    :book,
    :booklet,
    :inbook,
    :incollection,
    :inproceedings,
    :manual,
    :masterthesis,
    :misc,
    :phdthesis,
    :proceedings,
    :techreport,
    :unpublished,
]

"""
    const fields = [
        :address,
        :annote,
        :archivePrefix,
        :author,
        :booktitle,
        :chapter,
        :crossref,
        :edition,
        :editor,
        :eprint,
        :howpublished,
        :institution,
        :journal,
        :key,
        :month,
        :note,
        :number,
        :organization,
        :pages,
        :primaryClass,
        :publisher,
        :school,
        :series,
        :title,
        :type,
        :volume,
        :year
    ]
List of possible fields (currently based on bibtex). Keep it sorted for readability
"""
const fields = [
    :address,
    :annote,
    :archivePrefix,
    :author,
    :booktitle,
    :chapter,
    :crossref,
    :edition,
    :editor,
    :eprint,
    :howpublished,
    :institution,
    :journal,
    :key,
    :month,
    :note,
    :number,
    :organization,
    :pages,
    :primaryClass,
    :publisher,
    :school,
    :series,
    :title,
    :type,
    :volume,
    :year
]

"""
    const maxfieldlength
For output formatting purpose, for instance, export to BibTeX format.
"""
const maxfieldlength = maximum(map(s -> length(string(s)), fields))


text:
using OrdinaryDiffEq

# Auto callback merging

do_nothing = DiscreteCallback(
     (u,t,integrator) -> true,
     integrator -> nothing
 )
problem = ODEProblem(
   (u,p,t)->-u,
   1.0, (0.0,1.0),
   callback = do_nothing
)
solve(
   problem, Euler(),
   dt = 0.1,
   callback = do_nothing
)


text:
function transformIndexes(bP__::AbstractVector{C}, aTb) where {C <: CartesianIndex{2}}
  bV = zeros(3)
  bV[3] = 1.0

  aV = zeros(3)
  aP__ = Vector{Vector{Float64}}(undef, length(bP__))

  for (i,kp) in enumerate(bP__)
    bV[1] = kp[1]
    bV[2] = kp[2]

    aV[:] = aTb*bV
    aP__[i] = aV[1:2]
  end

  return aP__
end

function cost(matches, aTb_)
  pts1 = (x->x[1]).(matches) .|> x-> [x[1];x[2]]
  _pts2 = (x->x[2]).(matches)
  pts2 = transformIndexes(_pts2, aTb_)

  total = 0.0
  for (i,p2) in enumerate(pts2)
    total += norm( pts1[i] .- p2 )
  end

  return total
end

function cost_xyr!(xyr, matches, aTb = [1. 0 0; 0 1 0; 0 0 1])
  aTb[1:2,1:2] .= _Rot.RotMatrix(xyr[3])
  aTb[1:2, 3] .= xyr[1:2]
  cost(matches, aTb)
end

# draw matches

function mosaicMatches(img1, img2, matches)
  grid = hcat(img1, img2)
  offset = CartesianIndex(0, size(img1, 2))
  map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)
  grid
end

function getFeaturesDescriptors(img::AbstractMatrix)
  # img = parent(img_)
  brisk_params = BRISK()
  features = Features(Keypoints(imcorner(img, method=harris)))
  desc, ret_features = create_descriptor(Gray.(img), features, brisk_params)
end

function getCornersDescriptors( img::AbstractMatrix;
                                kparg1 = 12, kparg2=0.3,
                                size = 256, window = 20, seed = 123 )
  #
  brief_params = BRIEF(;size, window, seed)
  keypoints = Keypoints(fastcorners(img, kparg1, kparg2))
  desc, ret_features = create_descriptor(img, keypoints, brief_params);
end

function getKeypointsMatchBRIEF(  img1, img2;
                                  kparg1 = 12, kparg2=0.3, threshold=0.1,
                                  size = 256, window = 20, seed = 123 )
  #
  desc_1, ret_kps_1 = getCornersDescriptors( img1; kparg1, kparg2, size, window, seed )
  desc_2, ret_kps_2 = getCornersDescriptors( img2; kparg1, kparg2, size, window, seed )

  match_keypoints(ret_kps_1, ret_kps_2, desc_1, desc_2, threshold)
end

function getKeypointsMatchBRISK(img1, img2;
                                threshold=0.1,
                                kwargs... )
  #
  desc_1, ret_kps_1 = getFeaturesDescriptors( img1; kwargs... )
  desc_2, ret_kps_2 = getFeaturesDescriptors( img2; kwargs... )

  match_keypoints(Keypoints(ret_kps_1), Keypoints(ret_kps_2), desc_1, desc_2, threshold)
end

function calcMatchMinimize( img1::AbstractMatrix{<:Gray}, img2::AbstractMatrix{<:Gray};
                            keypoint_fnc::Function=getKeypointsMatchBRIEF,
                            kwargs...)
  #

  matches = keypoint_fnc( img1, img2; kwargs... )

  res = Optim.optimize((xyr)->cost_xyr!(xyr, matches), [0.0;0.0;0.0])
  res.minimizer, matches
end

function pushMatchWindowMedian!(iCB::AbstractVector{<:AbstractMatrix{<:Gray}},
                                trBuffer::AbstractVector{<:AbstractVector},
                                newImg::AbstractMatrix{<:Gray};
                                kwargs... )
  #
  for bimg in iCB
    trans, matches = calcMatchMinimize(newImg, bimg; kwargs...);
    push!(trBuffer, trans)
  end

  Ti = median((x->x[1]).(trBuffer))
  Tj = median((x->x[2]).(trBuffer))
  Tr = median((x->x[3]).(trBuffer))

  return [Ti; Tj; Tr]
end

Base.@kwdef struct BasicImageOnlyStabilization{N}
  imgBuf::CircularBuffer{<:AbstractMatrix{<:Gray}} = CircularBuffer{Matrix{Gray{N0f8}}}(N)
  trBuf::CircularBuffer{<:AbstractVector{Float64}} = CircularBuffer{Vector{Float64}}(N)
end

function (bi::BasicImageOnlyStabilization{N})( newImg::AbstractMatrix;
                                            offset::AbstractVector=zeros(3),
                                            roi_i = :,
                                            roi_j = :,
                                            scale_i=1,
                                            scale_j=1,
                                            scale_r=1,
                                            kwargs... ) where N
  #
  # image tracking must be in grayscale
  img_g = Gray.(newImg[roi_i,roi_j])
  # fill buffer with first image
  if length(bi.imgBuf) == 0
    for _ in 1:N
      push!(bi.imgBuf, img_g)
    end
  end
  tr = pushMatchWindowMedian!(bi.imgBuf, bi.trBuf, img_g ; kwargs... )
  # replace oldest image in buffer with new
  push!(bi.imgBuf, img_g)

  rot = ImageTransformations.recenter(_Rot.RotMatrix(scale_r*tr[3] + offset[3]), [Base.size(img_g)...] .÷ 2)  # a rotation around the center

  tform = rot ∘ Translation(scale_i*tr[1],scale_j*tr[2])
  warp(newImg, tform, axes(newImg))
end

function _stabImageSeq( fg,
                        lbls=sortDFG( ls(fg, tags=[:POSE;]) );
                        threshold=0.4, size=1024, window=50, kparg1 = 11,
                        filecount::Int=0,
                        folderpath="/tmp/caesar/philos",
                        filepath=folderpath*"/quickstab_$filecount.mp4",
                        logentry=false,
                        scale_j=0,
                        scale_r=0.8,
                        roi_i=1:900,
                        offset=[0;0;0.017*pi],
                        N::Integer=3 )
  #

  BIOS! = BasicImageOnlyStabilization{N}()

  des = listDataEntrySequence(fg, lbls[1], r"IMG_CENTER", sortDFG)
  img = fetchDataImage(fg, lbls[1], des[1])

  R_ = _Rot.RotMatrix(0.017*pi)
  img_ = warp(img, ImageTransformations.recenter(R_, center(img)));

  encoder_options = (crf=23, preset="medium")
  framerate=90
  open_video_out(filepath, img, framerate=framerate, encoder_options=encoder_options) do writer

    for lb in lbls
      println("stab ", lb)

      des = listDataEntrySequence(fg, lb, r"IMG_CENTER", sortDFG)
      imgs = fetchDataImage.(fg, lb, des)

      for (id,img) in enumerate(imgs)
        frame = BIOS!(img; roi_i, scale_j, scale_r, offset,
                      threshold, size, window, kparg1)
        #
        write(writer, frame)
        addData!(   fg, :camera, lb, Symbol(:STAB_, des[id]),
                    Caesar.toFormat(format"PNG", frame),
                    mimeType="image/png",
                    description="ImageMagick.readblob(imgBytes)" )
      end

    end

  end

  # store log
  if logentry
    open(folderpath*"/log.txt", "a") do fid
      println(fid, filepath, " : threshold=$threshold, size=$size, window=$window, kparg1=$kparg1, N=$N")
    end
  end

end

#


text:
import GreyDecision.Random
import GreyDecision.Utility
import GreyDecision.GreyNumbers: GreyNumber

@testset "Random Grey Numbers" begin

    @testset "Single Random GreyNumber" begin
        g = rand(GreyNumber)
        @test first(g) >= 0.0
        @test last(g) <= 1.0
    end

    @testset "Typed Random GreyNumber" begin
        g = rand(GreyNumber{Float64})
        @test first(g) >= 0.0
        @test last(g) <= 1.0

        h = rand(GreyNumber{UInt8})
        @test first(h) < 256
        @test last(h) >= 0
    end

    @testset "Array of random grey numbers" begin
        v = rand(GreyNumber{Float64}, 10)

        @test length(v) == 10
        @test typeof(v[1]) == GreyNumber{Float64}
        @test first(v[1]) >= 0.0
        @test last(v[1]) <= 1.0
    end

    @testset "Matrix of random grey numbers" begin
        v = rand(GreyNumber{Float64}, 10, 20)

        @test size(v) == (10, 20)
        @test typeof(v[1, 1]) == GreyNumber{Float64}
        @test first(v[1, 1]) >= 0.0
        @test last(v[1, 2]) <= 1.0
        @test first(v[2, 1]) >= 0.0
        @test last(v[2, 2]) <= 1.0
    end

end


text:
struct ZygoteAD <: ADBackend end
ADBackend(::Val{:zygote}) = ZygoteAD
function setadbackend(::Val{:zygote})
    ADBACKEND[] = :zygote
end


text:
function plotresult(result::Tuple{Dict{Symbol, Array},Array})
    observations, loglikelihoods = result

    nsamples = length(loglikelihoods)
    nvars = length(observations)
    i = 0
    for (var,samples) in observations
        subplot(nvars+1,2,2*i+1); plt[:hist](samples, color="black")
        subplot(nvars+1,2,2*i+2); scatter(1:nsamples, samples, s=1, c="black", alpha=0.1)
        i += 1
    end
    subplot(nvars+1,1,nvars+1)
    plot(1:nsamples, loglikelihoods, linewidth=1, c="black")
end


text:
"""
Generic interface for implementing inference algorithms.
An implementation of an algorithm should include the following:
1. A type specifying the algorithm and its parameters, derived from InferenceAlgorithm
2. A method of `run` function that produces results of inference. TODO: specify the format of this output
3. Methods for modelling functions specified in this file. This is where actual inference happens.

Turing translates models to thunks that call the modelling functions below at specified points.
Their dispatch is based on the value of a global variable `sampler`.
To include a new inference algorithm implement the requirements mentioned above in a separate file,
then include that file at the end of this one.
"""
type ParticleSampler{T} <: Sampler{T}
  alg ::  T
  particles :: ParticleContainer
  model :: Function
  ParticleSampler(m :: Function, a :: T) = (s = new(); s.alg = a; s.model = m; s)
end

## Concrete algorithm implementations.

include("is.jl")
include("smc.jl")
include("pgibbs.jl")
include("support/resample.jl")


text:
import LinearAlgebra: I, tr

@testset "Taylor Green Vortex Problem" begin #for q in quadratures
    q = D2Q9()

    scale = 1
    ν = 1.0 / 6.0
    problem = TaylorGreenVortex(ν, scale)
    f = LatticeBoltzmann.initialize(AnalyticalEquilibrium(), q, problem)

    # TODO: when an a thermal model is considered, the density fluctuates
    # TODO: when a thermal model is considerd we have constant density
    @testset "Hermite Pressure component" begin
        q = D2Q9()
        D = dimension(q)
        N = div(LatticeBoltzmann.order(q), 2)
        Hs = [
            [hermite(Val{n}, q.abscissae[:, i], q) for i in 1:length(q.weights)]
            for n in 1:N
        ]

        scale = 1
        ν = 1.0 / 6.0
        problem = DecayingShearFlow(ν, scale, static = false)
        problem = TaylorGreenVortex(ν, scale, static = true)

        x = 0.30
        y = 0.30
        f = LatticeBoltzmann.initial_condition(AnalyticalEquilibrium(), q, problem, x, y)

        ρ = sum(f)

        τ = problem.ν
        τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem) + 0.5

        a_f = [sum([f[idx] * Hs[n][idx] for idx in 1:length(q.weights)]) for n in 1:N]

        P = a_f[2] - (a_f[1] * a_f[1]') / ρ - ρ * I
        P = P * (1 - 1 / (2 * τ))

        σ_lb = (P - I * tr(P) / (D)) / (problem.u_max)
        σ = LatticeBoltzmann.deviatoric_tensor(q, problem, x, y, 0.0)
    end
    @testset "Pressure component" begin
        scale = 1
        ν = 2.0 / 6.0
        problem = TaylorGreenVortex(ν, scale, static = false)
        problem = DecayingShearFlow(ν, scale, static = false)

        x = 0.30
        y = 0.30
        f = LatticeBoltzmann.initial_condition(AnalyticalEquilibrium(), q, problem, x, y)

        ρ = LatticeBoltzmann.density(q, f)

        cs = 1 / q.speed_of_sound_squared
        # cs = q.speed_of_sound_squared
        d_u = problem.u_max * LatticeBoltzmann.velocity_gradient(problem, x, y, 0.0)
        τ = problem.ν
        τ = q.speed_of_sound_squared * LatticeBoltzmann.lattice_viscosity(problem) + 0.5
        for f_idx in 1:length(f)
            f[f_idx] +=
                -(q.weights[f_idx] * ρ * τ / cs) *
                sum(LatticeBoltzmann.hermite(Val{2}, q.abscissae[:, f_idx], q) .* d_u)
        end

        u = zeros(LatticeBoltzmann.dimension(q))
        LatticeBoltzmann.velocity!(q, f, ρ, u)
        T = LatticeBoltzmann.temperature(q, f, ρ, u)

        σ_lb =
            LatticeBoltzmann.momentum_flux(q, f, ρ, u) -
            I * LatticeBoltzmann.pressure(q, f, ρ, u)
        σ = LatticeBoltzmann.deviatoric_tensor(q, problem, x, y, 0.0)
    end
end


text:
# Sioux Falls network data
# http://www.bgu.ac.il/~bargera/tntp/

#Link travel time = free flow time * ( 1 + B * (flow/capacity)^Power ).
#Link generalized cost = Link travel time + toll_factor * toll + distance_factor * distance

# Traffic Assignment Data structure
type TA_Data
    network_name::String

    number_of_zones::Int64
    number_of_nodes::Int64
    first_thru_node::Int64
    number_of_links::Int64

    start_node::Array
    end_node::Array
    capacity::Array
    link_length::Array
    free_flow_time::Array
    B::Array
    power::Array
    speed_limit::Array
    toll::Array
    link_type::Array

    total_od_flow::Float64

    travel_demand::Array
    od_pairs::Array

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64
end

function load_ta_network_(out_dir, files_ID, month_w, day,  instan)
    toll_factor = 0
    distance_factor = 0
    network_data_file = files_ID * "_net_" * month_w * "_" * string(day) * '_' * instan * ".txt"
    trip_table_file = files_ID * "_trips_" * month_w * "_" * string(day) * '_' * instan * ".txt"
    best_objective = 0

    network_name = files_ID

    network_data_file =  out_dir  * "data_traffic_assignment_uni-class/" * network_data_file
    trip_table_file =   out_dir  * "data_traffic_assignment_uni-class/" * trip_table_file

    number_of_zones = 0
    number_of_links = 0
    number_of_nodes = 0
    first_thru_node = 0

    n = 0
    A = 0
    n = open(network_data_file, "r")
    A = readlines(n)
    for line in A
        #print( line)
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<NUMBER OF NODES>")
            number_of_nodes = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<FIRST THRU NODE>")
            first_thru_node = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<NUMBER OF LINKS>")
            number_of_links = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<END OF METADATA>")
       # println(line)
        end
    end

    @assert number_of_links > 0
    start_node = Array{Int}(number_of_links)
    end_node = Array{Int}(number_of_links)
    capacity = zeros(number_of_links)
    link_length = zeros(number_of_links)
    free_flow_time = zeros(number_of_links)
    B = zeros(number_of_links)
    power = zeros(number_of_links)
    speed_limit = zeros(number_of_links)
    toll = zeros(number_of_links)
    link_type = Array{Int}(number_of_links)

    idx = 1

    apa = 0
    for line in A
        if contains(line, "~")
            continue
        end

        if contains(line, ";")
            line = strip(line, '\n')
            line = strip(line, ';')
            numbers = split(line)
            start_node[idx] = parse(Int, numbers[1])
            end_node[idx] = parse(Int, numbers[2])

            capacity[idx] = parse(Float64, numbers[3])
            link_length[idx] = parse(Float64, numbers[4])

            free_flow_time[idx] = parse(Float64, numbers[5])
            B[idx] = parse(Float64, numbers[6])
            power[idx] = parse(Float64, numbers[7])
            speed_limit[idx] = parse(Float64, numbers[8])
            toll[idx] = parse(Float64, numbers[9])
            link_type[idx] = parse(Float64, numbers[10])

            idx = idx + 1
        end
    end

    ##################################################
    # Trip Table
    ##################################################

    number_of_zones_trip = 0
    total_od_flow = 0

    f = open(trip_table_file, "r")
    fe = readlines(f)

    for line in fe
        #println(line)
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones_trip = parse( Int, line[ search(line, '>')+1 : end ] )
        elseif contains(line, "<TOTAL OD FLOW>")
            total_od_flow = parse( Float64, line[ search(line, '>')+1 : end ] )
        elseif contains(line, "<END OF METADATA>")
            break
        end
    end

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = []
    origin = 0
    for line in fe
        if contains(line, "Origin")
            origin = parse( Int, split(line)[2] )
        elseif contains(line, ";")
            pairs = split(line, ";")
            for i=1:size(pairs)[1]
                if contains(pairs[i], ":")
                    pair = split(pairs[i], ":")
                    destination = parse( Int, strip(pair[1]) )
                    od_flow = parse( Float64, strip(pair[2]) )
                    travel_demand[origin, destination] = od_flow
                    push!(od_pairs, (origin, destination))
                   # println("origin=$origin, destination=$destination, flow=$od_flow")
                end
            end
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        network_name,
        number_of_zones,
        number_of_nodes,
        first_thru_node,
        number_of_links,
        start_node,
        end_node,
        capacity,
        link_length,
        free_flow_time,
        B,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective)

    return ta_data

end # end of load_network


text:
# Uses new formulation from limbdark paper.

include("s2_stable.jl")

s2diff(b,r) = s2(r,b)-convert(Float64,s2(big(r),big(b)))

# Run some tests:

using PyPlot

nb = 1024; nr = 1024
fig,axes = subplots(3,2)
# Compute r = b+-\epsilon
epsilon = 1e-8
fepsp= zeros(Float64,nb)
fepsbigp= zeros(Float64,nb)
fepsm= zeros(Float64,nb)
fepsbigm= zeros(Float64,nb)
b=linspace(epsilon,2.0,nr)
for ib=1:nb
  fepsbigp[ib] = float(s2(big(b[ib]+epsilon),big(b[ib])))
  fepsp[ib] = float(s2(b[ib]+epsilon,b[ib]))
  fepsbigm[ib] = float(s2(big(b[ib]-epsilon),big(b[ib])))
  fepsm[ib] = float(s2(b[ib]-epsilon,b[ib]))
end
ax = axes[1]
ax[:plot](b,fepsp,label=L"$s_2(b,b+10^{-8})$")
ax[:plot](b,fepsm,label=L"$s_2(b,b-10^{-8})$")
ax[:plot](b,s2.(b,b),label=L"$s_2(b,b)$")
ax[:legend](loc="upper right",fontsize=10)
ax = axes[2]
ax[:plot](b,(fepsp-s2.(b,b))/epsilon,label=L"$s_2(b,b+10^{-8})-s_2(b,b)$")
ax[:plot](b,(fepsbigp-s2.(b,b))/epsilon,linewidth=2,linestyle="--",label=L"$s_2(b,b+10^{-8})-s_2(b,b)$, big")
ax[:plot](b,(fepsm-s2.(b,b))/epsilon,label=L"$s_2(b,b-10^{-8})-s_2(b,b)$")
ax[:plot](b,(fepsbigm-s2.(b,b))/epsilon,linewidth=2,linestyle="--",label=L"$s_2(b,b-10^{-8})-s_2(b,b)$, big")
ax[:legend](loc="right",fontsize=5)
ax = axes[3]
ax[:plot](b,s2diff.(b,b+epsilon),label=L"$s_2(b,b+10^{-8})-$ big")
ax[:plot](b,s2diff.(b,b-epsilon),linewidth=2,linestyle="--",label=L"$s_2(b,b-10^{-8})-$ big")
#ax[:plot](b,fepsp-fepsbigp,label=L"$s_2(b,b+10^{-8})-$ big")
#ax[:plot](b,fepsm-fepsbigm,linewidth=2,linestyle="--",label=L"$s_2(b,b-10^{-8})-$ big")
ax[:legend](loc="lower right",fontsize=6)

# Compute b+r = 1+-\epsilon
epsilon = 1e-8
b=linspace(0.0,2.0,nr)
for ib=1:nr
  fepsbigp[ib] = float(s2(big(1.0-b[ib]+epsilon),big(b[ib])))
  fepsp[ib] = float(s2(1.0-b[ib]+epsilon,b[ib]))
  fepsbigm[ib] = float(s2(big(1.0-b[ib]-epsilon),big(b[ib])))
  fepsm[ib] = float(s2(1.0-b[ib]-epsilon,b[ib]))
end
ax = axes[4]
ax[:plot](b,fepsp,label=L"$s_2(b,1-b+10^{-8})$")
ax[:plot](b,fepsm,label=L"$s_2(b,1-b-10^{-8})$")
ax[:plot](b,s2.(1.0-b,b),label=L"$s_2(b,1-b)$")
ax[:legend](loc="lower right",fontsize=10)
ax = axes[5]
ax[:plot](b,(fepsp-s2.(1.0-b,b))/epsilon,label=L"$s_2(b,1-b+10^{-8})-s_2(b,1-b)$")
ax[:plot](b,(fepsbigp-s2.(1.0-b,b))/epsilon,linewidth=2,linestyle="--",label=L"$s_2(b,1-b+10^{-8})-s_2(b,1-b)$, big")
ax[:plot](b,(fepsm-s2.(1.0-b,b))/epsilon,label=L"$s_2(b,1-b-10^{-8})-s_2(b,1-b)$")
ax[:plot](b,(fepsbigm-s2.(1.0-b,b))/epsilon,linewidth=2,linestyle="--",linestyle="--",label=L"$s_2(b,1-b-10^{-8})-s_2(b,1-b)$, big")
ax[:legend](loc="lower right",fontsize=5)
ax = axes[6]
ax[:plot](b,s2diff.(1.0-b+epsilon,b),label=L"$s_2(b,1-b+10^{-8})-$ big")
ax[:plot](b,s2diff.(1.0-b-epsilon,b),linewidth=2,linestyle="--",label=L"$s_2(b,1-b+10^{-8})-$ big")
ax[:legend](loc="upper right",fontsize=6)
savefig("s2_machine.pdf", bbox_inches="tight")

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using BandedMatrices

# The Solve Function
function DiffEqBase.__solve(prob::BVProblem, alg::Shooting; kwargs...)
    bc = prob.bc
    u0 = deepcopy(prob.u0)
    # Form a root finding function.
    loss = function (resid, minimizer)
        uEltype = eltype(minimizer)
        tmp_prob = remake(prob,u0=minimizer)
        sol = solve(tmp_prob,alg.ode_alg;kwargs...)
        bc(resid,sol,sol.prob.p,sol.t)
        nothing
    end
    opt = alg.nlsolve(loss, u0)
    sol_prob = remake(prob,u0=opt[1])
    sol = solve(sol_prob, alg.ode_alg;kwargs...)
    if sol.retcode == opt[2]
        DiffEqBase.solution_new_retcode(sol,:Success)
    else
        DiffEqBase.solution_new_retcode(sol,:Failure)
    end
    sol
end

function DiffEqBase.__solve(prob::BVProblem, alg::Union{GeneralMIRK,MIRK}; dt=0.0, kwargs...)
    if dt<=0
        error("dt must be positive")
    end
    n = Int(cld((prob.tspan[2]-prob.tspan[1]),dt))
    x = collect(range(prob.tspan[1], stop=prob.tspan[2], length=n+1))
    S = BVPSystem(prob, x, alg_order(alg))
    if isa(prob.u0, AbstractVector{<:Number})
        copyto!.(S.y, (prob.u0,))
    elseif isa(prob.u0, AbstractVector{<:AbstractArray})
        copyto!(S.y, prob.u0)
    else
        error("u0 must be a Vector or Vector of Arrays")
    end
    tableau = constructMIRK(S)
    cache = alg_cache(alg, S)
    # Upper-level iteration
    vec_y = Array{eltype(S.y[1])}(undef, S.M*S.N)              # Vector
    reorder! = function (resid)
        # reorder the Jacobian matrix such that it is banded
        tmp_last = resid[end]
        for i in (length(resid)-1):-1:1
            resid[i+1] = resid[i]
        end
        resid[1], resid[end] = resid[end], tmp_last
    end
    loss = function (resid, minimizer)
        nest_vector!(S.y, minimizer)
        Φ!(S, tableau, cache)
        isa(prob.problem_type, TwoPointBVProblem) ? eval_bc_residual!(S) : general_eval_bc_residual!(S)
        flatten_vector!(resid, S.residual)
        reorder!(resid)
        nothing
    end

    jac_wrapper = BVPJacobianWrapper(loss)

    flatten_vector!(vec_y, S.y)
    opt = isa(prob.problem_type, TwoPointBVProblem) ? alg.nlsolve(ConstructJacobian(jac_wrapper, vec_y), vec_y) : alg.nlsolve(ConstructJacobian(jac_wrapper, S, vec_y), vec_y) # Sparse matrix is broken
    nest_vector!(S.y, opt[1])

    retcode = opt[2] ? :Success : :Failure
    DiffEqBase.build_solution(prob, alg, x, S.y, retcode = retcode)
end

#=
function solve(prob::BVProblem, alg::GeneralMIRK; dt=0.0, kwargs...)
    n = Int(cld((prob.tspan[2]-prob.tspan[1]),dt))
    x = collect(linspace(prob.tspan..., n+1))
    S = BVPSystem(prob.f, prob.bc, x, length(prob.u0), alg_order(alg))
    S.y[1] = prob.u0
    tableau = constructMIRK(S)
    cache = alg_cache(alg, S)
    # Upper-level iteration
    vec_y = Array{eltype(S.y[1])}(S.M*S.N)              # Vector
    loss = function (minimizer, resid)
        nest_vector!(S.y, minimizer)
        Φ!(S, tableau, cache)
        general_eval_bc_residual!(S)
        flatten_vector!(resid, S.residual)
        nothing
    end
    flatten_vector!(vec_y, S.y)
    opt = alg.nlsolve(loss, vec_y)
    nest_vector!(S.y, opt[1])
    retcode = opt[2] ? :Success : :Failure
    DiffEqBase.build_solution(prob, alg, x, S.y, retcode = retcode)
end
=#


text:
"""
    sliceview(A, code)
    slicecopy(A, code)

Slice array `A` according to `code`, a tuple of length `ndims(A)`,
in which `:` indicates a dimension of the slices, and `*` a dimension separating them.
For example if `code = (:,*,:)` then slices are either `view(A, :,i,:)`
or `A[:,i,:]` with `i=1:size(A,2)`.
"""
function sliceview(A::AbstractArray{T,N}, code::Tuple) where {T,N}
    N == length(code) || throw(ArgumentError("wrong code length"))
    if code == (:,*)
        collect(eachcol(A))
    elseif code == (*,:)
        collect(eachrow(A))
    elseif count(isequal(*), code) == 1
        collect(eachslice(A, dims = findfirst(isequal(*), code)))
    else
        iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(A,d) : Ref(:), Val(N))...)
        [ view(A, i...) for i in iter ]
    end
end

@doc @doc(sliceview)
function slicecopy(A::AbstractArray{T,N}, code::Tuple) where {T,N}
    N == length(code) || throw(ArgumentError("wrong code length"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(A,d) : Ref(:), Val(N))...)
    [ A[i...] for i in iter ]
end

sliceview(A::AbstractArray{T,N}) where {T,N} = sliceview(A, ntuple(i-> i==N ? (*) : (:), N))
slicecopy(A::AbstractArray{T,N}) where {T,N} = slicecopy(A, ntuple(i-> i==N ? (*) : (:), N))

"""
    glue!(B, A, code)
    copy_glue(A, code) = glue!(Array{T}(...), A, code)

Copy the contents of an array of arrays into one larger array,
un-doing `sliceview` / `slicecopy` with the same `code`.
This function can handle arbitrary codes.
(Also called `stack` or `align` elsewhere.)

    cat_glue(A, code)
    red_glue(A, code)

The same result, but calling either things like `hcat(A...)`
or things like `reduce(hcat, A)`.
The code must be sorted like `(:,:,:,*,*)`, except that `(*,:)` is allowed.

    glue(A)
    glue(A, code)

If `code` is omitted, the default is like `(:,:,*,*)`
with `ndims(first(A))` colons first, then `ndims(A)` stars.
If the inner arrays are `StaticArray`s (and the code is sorted) then it calls `static_glue`.
Otherwise it will call `red_glue`, unless code is unsuitable for that, in which case `copy_glue`.
"""
function glue(A::AbstractArray{<:AbstractArray{T,IN},ON}, code::Tuple=defaultcode(IN,ON)) where {T,IN,ON}
    if iscodesorted(code) || code == (*,:)
        red_glue(A, code)
    else
        copy_glue(A, code)
    end
end

glue(A::AbstractArray{<:Tuple,ON}, code::Tuple=defaultcode(1,ON)) where {ON} = copy_glue(A, code)

glue(A::AbstractArray{<:Number,ON}, code::Tuple=defaultcode(0,ON)) where {ON} = A

defaultcode(IN::Int, ON::Int) = ntuple(d-> d<=IN ? (:) : (*), IN+ON)

@doc @doc(glue)
@inline function red_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    if code == (:,*)
        reduce(hcat, A) # this is fast, specially optimised
    elseif code == (*,:)
        reduce(vcat, PermuteDims.(A))
    elseif code == (:,:,*)
        mat = reduce(hcat, A)
        reshape(mat, gluedsize(A, code))
    elseif iscodesorted(code)
        mat = reduce(hcat, vec(vec.(A)))
        reshape(mat, gluedsize(A, code))
    else
        throw(ArgumentError("can't glue code = $code with reduce(cat...)"))
    end
end

@doc @doc(glue)
@inline function cat_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    if code == (:,*)
        hcat(A...)
    elseif code == (*,:)
        vcat(PermuteDims.(A)...)
    elseif count(isequal(*), code) == 1 && code[end] == (*)
        cat(A...; dims = length(code))
    elseif iscodesorted(code)
        flat = cat(vec(A)...; dims = length(code)-ndims(A)+1)
        reshape(flat, gluedsize(A, code))
    else
        throw(ArgumentError("can't glue code = $code with cat(A...)"))
    end
end

@doc @doc(glue)
@inline function copy_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    B = Array{eltype(first(A))}(undef, gluedsize(A, code))
    glue!(B, A, code)
end

@doc @doc(glue)
function glue!(B::AbstractArray{T,N}, A::AbstractArray{IT,ON}, code::Tuple) where {T,N,IT,ON}
    gluecodecheck(A, code)
    N == _ndims(A) + _ndims(first(A))  || throw(DimensionMismatch("wrong size target"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(B,d) : Ref(:), Val(N))...)
    for i in iter
        copyto!(view(B, i...), A[decolonise(i)...] )
    end
    B
end

function gluecodecheck(A::AbstractArray, code::Tuple)
    colons = countcolons(code)
    inner = _ndims(first(A))
    outer = _ndims(A)
    colons == inner || throw(DimensionMismatch("wrong number of dimensions: " *
        "ndims(first(A)) == $inner not cannot be glued with code = $(pretty(code))"))
    length(code) - colons == outer || throw(DimensionMismatch("wrong number of dimensions: " *
        "ndims(A) == $outer cannot be glued with code = $(pretty(code))"))
end

@generated function decolonise(i::Tuple)
    ind = Int[]
    for k in 1:length(i.parameters)
        if i.parameters[k] != Colon
            push!(ind, k)
        end
    end
    Expr(:tuple, [Expr(:ref, :i, k) for k in ind]...)
end

@generated function countcolons(code::Tuple)
    n = 0
    for s in code.parameters
        if s == Colon
            n += 1
        end
    end
    n
end

@generated function iscodesorted(code::Tuple) # true if all colons before all stars
    colons = true
    sorted = true
    for s in code.parameters
        if s != Colon && colons # then we're at transition, change flag
            colons = false
        elseif s == Colon && !colons # then a : is following a *
            sorted = false
        end
    end
    sorted
end

@generated function gluedsize(A::AbstractArray{IT, N}, code::Tuple) where {IT, N}
    list = Any[]
    dout = 1
    din = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(first(A),$din) ))
            din += 1
        else
            push!(list, :( size(A,$dout) ))
            dout += 1
        end
    end
    :( ($(list...),) )
end

_ndims(A) = ndims(A)
_ndims(A::Tuple) = 1

using LazyStack

@inline function lazy_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    if code == (*,:)
        PermuteDims(stack(A))
    elseif iscodesorted(code)
        stack(A)
    else
        error("can't glue code = $code with LazyStack")
    end
end

using ZygoteRules # TODO add tests?

# Rules moved from SliceView.jl

@adjoint sliceview(A::AbstractArray, code::Tuple) =
    sliceview(A, code), Δ -> (glue(Δ, code), nothing)

@adjoint red_glue(A::AbstractArray, code::Tuple) =
    red_glue(A, code), Δ -> (sliceview(Δ, code), nothing)

@adjoint copy_glue(A::AbstractArray, code::Tuple) =
    copy_glue(A, code), Δ -> (sliceview(Δ, code), nothing)


text:
module Categorical_Cross_Entropy_Loss
    using LoopVectorization

    function func(output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        loss_matrix = zeros(Float32, size(output_matrix))
        @avxt for i in axes(loss_matrix, 1), j in axes(loss_matrix, 2)
            loss_matrix[i,j] = -sample_matrix[i,j]*log(max(output_matrix[i,j], 1e-8))
        end
        return loss_matrix
    end

    function prop!(δ::Array{Float32}, output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        @avxt for i in eachindex(δ)
            δ[i] = -sample_matrix[i]/max(output_matrix[i], 1e-8)
        end
    end
end


text:
module StatisticTools

export weightedmean

"""
weightedmean(values::Array{Float64,1}, errors::Array{Float64,1})

computes the weighted mean of values, using as weights w = 1/ sigma^2
"""
function weightedmean(x::Array{Float64,1}, errors::Array{Float64,1})
    @assert length(x) == length(errors)
    weights = 1./ (errors.^2)
    sigma = sqrt( 1./ sum( weights ) )
    mean = sigma^2 * sum(x .* weights)

    mean, sigma
end

end # module


text:
"""
    Shape{T}

    Abstract type; Subtypes contain geometric and visual information for a Body.
"""
abstract type Shape{T} end

"""
    EmptyShape{T} <: Shape{T}

    Contains no geometric or visual information
"""
struct EmptyShape{T} <: Shape{T}
    EmptyShape() = new{Float64}()
end

#TODO: change to MeshShape

"""
    Mesh{T} <: Shape{T}

    Contains geometric and visual information based on .obj file
"""
mutable struct Mesh{T} <: Shape{T}
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}

    path::String
    scale::SVector{3,T}
    color::RGBA

    function Mesh(path::String;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((position_offset, axis_offset))...)
        new{T}(position_offset, axis_offset, path, scale, color)
    end

    function Mesh(path::String, m::Real, J::AbstractMatrix;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            name::Symbol=Symbol("body_" * randstring(4)),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((m, J, position_offset, axis_offset))...)
        return Body(m, J; name=name, shape=new{T}(position_offset, axis_offset, path, scale, color))
    end
end

"""
    Box{T} <: Shape{T}

    Cuboid geometry

    position_offset: geometry origin offset from center of mass
    axis_offset: orientation offset from body frame
    xyz: dimensions (meters)
    scale: scaling
    color: RGBA
"""
mutable struct Box{T} <: Shape{T}
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}

    xyz::SVector{3,T}
    scale::SVector{3,T}
    color::RGBA

    function Box(x::Real, y::Real, z::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((x, y, z, position_offset, axis_offset))...)
        new{T}(position_offset, axis_offset, [x; y; z], scale, color)
    end

    function Box(x::Real, y::Real, z::Real, m::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            name::Symbol=Symbol("body_" * randstring(4)),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((x, y, z, m, position_offset, axis_offset))...)
        J = 1 / 12 * m * diagm([y^2 + z^2; x^2 + z^2; x^2 + y^2])
        return Body(m, J; name=name, shape=new{T}(position_offset, axis_offset, [x;y;z], scale, color))
    end
end

"""
    Cylinder{T} <: Shape{T}

    cylinder geometry

    position_offset: geometry origin offset from center of mass
    axis_offset: orientation offset from body frame
    rh: radius and height dimensions (meters)
    scale: scaling
    color: RGBA
"""
mutable struct Cylinder{T} <: Shape{T}
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}

    rh::SVector{2,T}
    scale::SVector{3,T}
    color::RGBA

    # Cylinder points in the z direction
    function Cylinder(r::Real, h::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((r, h, position_offset, axis_offset))...)
        new{T}(position_offset, axis_offset, [r;h], scale, color)
    end

    function Cylinder(r::Real, h::Real, m::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            name::Symbol=Symbol("body_" * randstring(4)),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((r, h, m, position_offset, axis_offset))...)
        J = 1 / 2 * m * diagm([r^2 + 1 / 6 * h^2; r^2 + 1 / 6 * h^2; r^2])
        return Body(m, J; name=name, shape=new{T}(position_offset, axis_offset, [r;h], scale, color))
    end
end

"""
    Capsule{T} <: Shape{T}

    capsule geometry

    position_offset: geometry origin offset from center of mass
    axis_offset: orientation offset from body frame
    rh: radius and height dimensions (meters)
    scale: scaling
    color: RGBA
"""
mutable struct Capsule{T} <: Shape{T}
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}

    rh::SVector{2,T}
    scale::SVector{3,T}
    color::RGBA

    # Capsule points in the z direction
    function Capsule(r::Real, h::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion= one(Quaternion),
            scale::AbstractVector=sones(3),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((r, h, position_offset, axis_offset))...)
        new{T}(position_offset, axis_offset, [r; h], scale, color)
    end

    function Capsule(r::Real, h::Real, m::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            name::Symbol=Symbol("body_" * randstring(4)),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((r, h, m, position_offset, axis_offset))...)

        mass_cylinder = π * h * r^2.0
        mass_hemisphere = π * 2.0 / 3.0 * r^3.0
        mass_total = mass_cylinder + 2 * mass_hemisphere
        Ixx_cylinder = mass_cylinder * (1.0 / 12.0 * h^2.0 + 1.0 / 4.0 * r^2.0)
        Ixx_hemisphere = 83.0 / 320 * mass_hemisphere * r^2
        d = (3.0 / 8.0 * r + 0.5 * h)
        Ixx = Ixx_cylinder + 2.0 * (Ixx_hemisphere + mass_hemisphere * d^2.0)
        Izz = 0.5 * mass_cylinder * r^2.0 + mass_hemisphere * (2.0 * r^2.0) / 5.0

        J = m * diagm([Ixx; Ixx; Izz])

        return Body(m, J; name=name, shape=new{T}(position_offset, axis_offset, [r; h], scale, color))
    end
end

"""
    Shapes{T} <: Shape{T}

    composite geometry

    shape: list of Shape objects
    position_offset: geometry origin offset from center of mass
    axis_offset: orientation offset from body frame
    xyz: dimensions (meters)
    scale: scaling
    color: RGBA
"""
mutable struct Shapes{T} <: Shape{T}
    shape::Vector
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}
    scale::SVector{3,T}
    color::RGBA

    function Shapes(shapes::Vector{Shape{T}};
        position_offset::AbstractVector=szeros(3),
        axis_offset::Quaternion=one(Quaternion),
        scale::AbstractVector=sones(3),
        name::Symbol=Symbol("body_" * randstring(4)),
        color=RGBA(0.75, 0.75, 0.75)) where T
        new{T}(shapes, position_offset, axis_offset, scale, color)
    end

    function Shapes(shapes::Vector, m::T, J;
        position_offset::AbstractVector=szeros(3),
        axis_offset::Quaternion=one(Quaternion),
        scale::AbstractVector=sones(3),
        name::Symbol=Symbol("body_" * randstring(4)),
        color=RGBA(0.75, 0.75, 0.75)) where T
        Body(m, J; name=name, shape=new{T}(shapes, position_offset, axis_offset, scale, color))
    end
end

"""
    Sphere{T} <: Shape{T}

    sphere geometry

    position_offset: geometry origin offset from center of mass
    axis_offset: orientation offset from body frame
    r: radius (meters)
    scale: scaling
    color: RGBA
"""
mutable struct Sphere{T} <: Shape{T}
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}
    r::T
    scale::SVector{3,T}
    color::RGBA

    function Sphere(r::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((r, position_offset, axis_offset))...)
        new{T}(position_offset, axis_offset, r, scale, color)
    end

    function Sphere(r::Real, m::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            name::Symbol=Symbol("body_" * randstring(4)),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((r, m, position_offset, axis_offset))...)
        J = 2 / 5 * m * diagm([r^2 for i = 1:3])
        return Body(m, J; name=name, shape=new{T}(position_offset, axis_offset, r, scale, color))
    end
end

"""
    Pyramid{T} <: Shape{T}

    pyramid geometry

    position_offset: geometry origin offset from center of mass
    axis_offset: orientation offset from body frame
    wh: width and height dimensions (meters)
    scale: scaling
    color: RGBA
"""
mutable struct Pyramid{T} <: Shape{T}
    position_offset::SVector{3,T}
    axis_offset::Quaternion{T}
    wh::SVector{2,T}
    scale::SVector{3,T}
    color::RGBA

    # Pyramid points in the z direction, Center of mass at 1/4 h
    function Pyramid(w::Real, h::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((w, h, position_offset, axis_offset))...)
        new{T}(position_offset, axis_offset, [w;h], scale, color)
    end

    function Pyramid(w::Real, h::Real, m::Real;
            position_offset::AbstractVector=szeros(3),
            axis_offset::Quaternion=one(Quaternion),
            scale::AbstractVector=sones(3),
            name::Symbol=Symbol("body_" * randstring(4)), color=RGBA(0.75, 0.75, 0.75))
        T = promote_type(quateltype.((w, h, m, position_offset, axis_offset))...)
        J = 1/80 * m * diagm([4*w^2+3*h^2;4*w^2+3*h^2;8*w^2])
        return Body(m, J; name=name, shape=new{T}(position_offset, axis_offset, [w;h], scale, color))
    end
end

function convert_shape(box::Box)
    x, y, z = Tuple(box.xyz)
    return GeometryBasics.HyperRectangle(Vec(-x / 2.0, -y / 2.0, -z / 2.0), Vec(x, y, z))
end

function convert_shape(cylinder::Cylinder)
    r, h = Tuple(cylinder.rh)
    return GeometryBasics.Cylinder(Point(0.0, 0.0, -h / 2.0),Point(0.0, 0.0, h / 2.0), r)
end

function convert_shape(sphere::Sphere)
    r = sphere.r
    return GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), r)
end

function convert_shape(pyramid::Pyramid)
    w, h = Tuple(pyramid.wh)
    return GeometryBasics.Pyramid(Point(0.0, 0.0, -h / 4.0), h, w)
end

function convert_shape(mesh::Mesh)
    return MeshFileObject(joinpath(@__DIR__, "..", "..", mesh.path))
end

function convert_shape(::EmptyShape)
    return nothing
end

function convert_shape(capsule::Capsule)
    r, h = Tuple(capsule.rh)
    p1 = Point(0.0, 0.0, -h / 2.0)
    p2 = Point(0.0, 0.0, h / 2.0)
    cyl = GeometryBasics.Cylinder(p1, p2, r)
    cap1 = GeometryBasics.Sphere(p1, r)
    cap2 = GeometryBasics.Sphere(p2, r)
    return [cyl, cap1, cap2]
end

function convert_shape(shapes::Shapes)
    geom = []
    for s in shapes.shape
        push!(geom, convert_shape(s))
    end
    return geom
end

# color
set_color!(shape::EmptyShape, color) = nothing

function set_color!(shape::Shape, color)
    shape.color = color
end

function set_color!(shapes::Shapes, color)
    for i in eachindex(shapes)
        shapes.shape[i].color = color
    end
end


text:
abstract LearnableParams

immutable NoParams <: LearnableParams end

"""
We need a wrapper for params that can hold:
- a vector view of the full parameter vector
- one or more reshaped views of that vector view

so for t::Transformation:

θ(t) --> SubArray
t.w, t.b --> view(θ)

Problem: every time we change the "master parameter vector"
we need to update the views (but... do we?? gotta check if
    just appending to the master params will invalidate the views)

Solution?  Store anonymous functions which return views of θ(t),
    and update w/b by calling those functions.

By default we can do `view(zeros(T, params_length(t)), :)`, then swap
out with a master param vector on update.
"""
type Params{T <: AbstractVector, VIEWS <: Tuple, S <: Tuple} <: LearnableParams
    θ::T
    ∇::T
    views::VIEWS
    ∇_views::VIEWS
    sizes::S
end

function Params(θ::AbstractVector, ∇::AbstractVector, sizes = ())
    views = splitview(θ, sizes)[1]
    ∇_views = splitview(∇, sizes)[1]

    # note: we pass views in so that if we need to swap these out with other views later
    #   we don't have to make a copy (type T is stable)
    n = length(θ)
    Params(view(θ,1:n), view(∇,1:n), views, ∇_views, sizes)
end

function Params{T}(::Type{T}, n::Int, sizes = ())
    θ = zeros(T, n)
    ∇ = zeros(T, n)
    Params(θ, ∇, sizes)
end

# most Learnables will just update their reference
#   but others (Chain, Graphs) will need to re-consolidate params
function reset_params!(t::Learnable, θ::AbstractVector, ∇::AbstractVector)
    p = t.params
    p.θ = θ
    p.∇ = ∇
    p.views = splitview(p.θ, p.sizes)[1]
    p.∇_views = splitview(p.∇, p.sizes)[1]
    return
end

Base.getindex(p::Params, i::Int) = p.views[i]

value(p::Params) = p.θ
grad(p::Params) = p.∇

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

function consolidate_params{T}(::Type{T}, transforms::Learnable...)
    consolidate_params(T, collect(transforms))
end

function consolidate_params{T}(::Type{T}, transforms::AbstractVector;
                                θ = nothing, ∇ = nothing)
    lens = ntuple(i -> params_length(transforms[i]), length(transforms))
    nparams = sum(lens)
    if θ == nothing
        θ = zeros(T, nparams)
        ∇ = zeros(T, nparams)
    end
    sizes = map(l -> (l,), lens)
    θs = splitview(θ, sizes)[1]
    ∇s = splitview(∇, sizes)[1]
    for (i,t) in enumerate(transforms)
        if params_length(t) > 0
            # first update the values of the new array, then reset the params with this reference
            θs[i][:] = t.params.θ
            ∇s[i][:] = t.params.∇
            reset_params!(t, θs[i], ∇s[i])
            @assert t.params.θ === θs[i]
        end
    end
    Params(θ, ∇)
end