text:
module SpikingNeuralNets
using LightGraphs, SimpleWeightedGraphs, Statistics, SparseArrays, Lazy
import Base.size, Base.length, Base.Broadcast.broadcastable

export AbstractSNN,  # data-types
    size, neurons, memory, bias, spikes,
    spike, spiketimes, spiketime, spiketrain,
    voltages, voltage, configuration, input, output, inputs, outputs,
    is_input, is_output, excitors, inhibitors,  # accessor functions
    potential, # neuron potential function
    sigmoid, thresh, ramp, fire, # transfer functions
    step!  # used to iterate the SNN

"""
    AbstractSNN{VT<:Real, ST<:Integer, G<:AbstractGraph}

An abstract SNN type where VT is the voltage type,
ST is the spike train index type, and G is the graph type.

An `AbstractSNN{VT, ST, G}` must have the following elements:
    - `graph::G`
    - `I::Vector{Bool}`
    - `O::Vector{Bool}`
    - `S::Vector{Dict{ST, ST}}`
    - `V::Vector{Vector{VT}}`
"""
abstract type AbstractSNN{VT<:Real, ST<:Integer, G<:AbstractGraph{ST}} end
Broadcast.broadcastable(snn::AbstractSNN) = Ref(snn)

"""
    size(snn::AbstractSNN)::Integer

Return the number of neurons in the network `snn`.
"""
size(snn::AbstractSNN) = nv(snn.graph)

"""
    length(snn::AbstractSNN)::Integer

Return the number of neurons in the network `snn`.
"""
length(snn::AbstractSNN) = nv(snn.graph)

"""
    neurons(snn::AbstractSNN)::UnitRange{Integer}

Return the a range of indices for all neurons in the network `snn`.
"""
neurons(snn::AbstractSNN) = vertices(snn.graph)

"""
    memory(snn::AbstractSNN)::Integer

Return the size of the network `snn`'s voltage memory.
"""
memory(snn::AbstractSNN) = length(snn.V)

"""
    bias(snn::AbstractSNN)::Vector{<:Real}

Return the bias for all neurons in the network `snn`.
"""
function bias(snn::AbstractSNN{VT})::Vector{VT} where {VT<:Real}
    return zeros(VT, neurons(snn))
end

"""
    bias(snn::AbstractSNN, n::Integer)::Real

Return the bias of neuron `n` in the network `snn`.
"""
@views function bias(snn::AbstractSNN{VT}, n::Integer)::VT where {VT<:Real}
    return bias(snn)[n]
end

"""
    voltages(snn::AbstractSNN{VT})::Vector{VT} where {VT<:Real}

Get the current voltage of all neurons in the network `snn`.
"""
@views function voltages(snn::AbstractSNN{VT})::Vector{VT} where {VT<:Real}
    return snn.V[end]
end

"""
    voltage(snn::AbstractSNN{VT}, n::Integer)::VT where {VT<:Real}

Get the current voltage of neuron `n` in the network `snn`.
"""
@views function voltage(snn::AbstractSNN{VT}, n::Integer)::VT where {VT<:Real}
    return snn.V[end][n]
end

"""
    spike(snn::AbstractSNN{VT,ST}, n::Integer; t::Integer=zero(ST))::ST

Calculate whether neuron `n` in the network `snn` was spiking `t` timesteps ago.
By default, calculates whether `n` is currently spiking.
"""
function spike(snn::AbstractSNN{VT,ST}, n::Integer; t::Integer=zero(ST))::ST where {VT<:Real, ST<:Integer}
    return get(snn.S[n], ST(t), zero(ST))
end

"""
    spikes(snn::AbstractSNN, neurons::AbstractVector{<:Integer}=neurons(snn); t::Integer=0)::Vector{ST}

Calculate whether `neurons` in the network `snn` are currently spiking.
"""
function spikes(snn::AbstractSNN{VT,ST}, neurons::AbstractVector{<:Integer}=neurons(snn); t::Integer=0)::Vector{ST} where {VT<:Real, ST<:Integer}
    return spike.(snn, neurons; t=t)
end

"""
    spiketimes(snn::AbstractSNN{VT, ST}, n::Integer)::Vector{ST} where {VT<:Real, ST<:Integer}

Return the number of timesteps into the past of all spikes for neuron `n` in the network `snn`.
"""
@views function spiketimes(snn::AbstractSNN{VT, ST}, n::Integer)::Base.KeySet{ST, Dict{ST, ST}} where {VT<:Real, ST<:Integer}
    return keys(snn.S[n])
end

"""
    spiketimes(snn::AbstractSNN{VT, ST})::Vector{Vector{ST}} where {VT<:Real, ST<:Integer}

Return the number of timesteps into the past of all spikes for all neurons in the network `snn`.
"""
function spiketimes(snn::AbstractSNN{VT,ST})::Vector{Vector{ST}} where {VT<:Real, ST<:Integer}
    return @. collect(spiketimes(snn, neurons(snn)))
end

"""
    spiketime(snn::AbstractSNN{VT, ST}, n::Integer)::ST where {VT<:Real, ST<:Integer}

Return the number of timesteps into the past of the most recent spike for neuron
    `n` in the network `snn`, or typemax(ST) if `n` has not yet spiked.
"""
@views function spiketime(snn::AbstractSNN{VT,ST}, n::Integer)::ST where {VT<:Real, ST<:Integer}
    if isempty(snn.S[n])
        return typemax(ST)
    end
    return minimum(keys(snn.S[n]))
end

"""
    spiketrain(snn::AbstractSNN, t::T)::SparseMatrixCSC{Bool, T}

Return the full matrix S[neuron, time] (with t columns) of binary spike data
for all neurons in the network `snn` for the last `t` time-steps.
"""
function spiketrain(snn::AbstractSNN, t::T)::SparseMatrixCSC{Bool, T} where T<:Integer
    if t <= 0
        return spzeros(Bool, size(snn), 0)
    end
    nspikes = sum(map(length, snn.S))
    I = Vector{T}(undef, nspikes)
    J = Vector{T}(undef, nspikes)
    S = trues(nspikes)
    s = 1
    @inbounds @views for n in eachindex(snn.S)
        for spiketime in keys(snn.S[n])
            I[s] = n
            J[s] = t - spiketime
            s += 1
        end
    end

    i = J .>= 1
    return sparse(I[i], J[i], S[i], size(snn), t)
end

"""
    inputs(snn::AbstractSNN)

Return the indices of all input neurons in the network `snn`.
"""
inputs(snn::AbstractSNN) = findall(snn.I)

"""
    outputs(snn::AbstractSNN)

Return the indices of all output neurons in the network `snn`.
"""
outputs(snn::AbstractSNN) = findall(snn.O)

"""
    input(snn::AbstractSNN{VT, ST})::Vector{ST} where {VT<:Real, ST<:Integer}

Get the current spike activity of `snn`'s input neurons.
"""
function input(snn::AbstractSNN{VT, ST})::Vector{ST} where {VT<:Real, ST<:Integer}
    return spikes(snn, findall(snn.I))
end

"""
    inputs(snn::AbstractSNN{VT, ST}, n::Integer)::Vector{ST} where {VT<:Real, ST<:Integer}

Get the indices of all input neurons with projections to neuron `n`
"""
function inputs(snn::AbstractSNN{VT, ST}, n::Integer)::Vector{ST} where {VT<:Real, ST<:Integer}
    i = inneighbors(snn.graph, n)
    return i[snn.I[i]]
end

"""
    outputs(snn::AbstractSNN{VT, ST}, n::Integer)::Vector{ST} where {VT<:Real, ST<:Integer}

Get the indices of all output neurons with projections from neuron `n`
"""
function outputs(snn::AbstractSNN{VT, ST}, n::Integer)::Vector{ST} where {VT<:Real, ST<:Integer}
    o = outneighbors(snn.graph, n)
    return o[snn.O[o]]
end

"""
    output(snn::AbstractSNN{VT, ST})::Vector{ST} where {VT<:Real, ST<:Integer}

Get the current spike activity of `snn`'s output neurons.
"""
function output(snn::AbstractSNN{VT, ST})::Vector{ST} where {VT<:Real, ST<:Integer}
    return spikes(snn, findall(snn.O))
end

"""
    is_input(snn::AbstractSNN, n::Integer)::Bool

Return true if neuron `n` is an input neuron in `snn`.
"""
is_input(snn::AbstractSNN, n::Integer)::Bool = snn.I[n]

"""
    is_output(snn::AbstractSNN, n::Integer)::Bool

Return true if neuron `n` is an output neuron in `snn`.
"""
is_output(snn::AbstractSNN, n::Integer)::Bool = snn.O[n]

"""
    excitors(snn::AbstractSNN, n::integer)::Vector{<:Integer}

Return the indices of all neurons in the network `snn` with
excitatory connections to neuron `n`. If `inputs=false` (default),
exclude input neurons from the resulting vector.
"""
function excitors(snn::AbstractSNN, n::Integer; inputs=false)::Vector{<:Integer}
    i = inneighbors(snn.graph, n)
    w = view(weights(snn.graph), i, n)
    if inputs
        return i[w .> 0]
    end
    return i[@. w > 0 & !snn.I[i]]
end

"""
    inhibitors(snn::AbstractSNN, n::integer)::Vector{<:Integer}

Return the indices of all neurons in the network `snn` with
inhibitory connections to neuron `n`. If `inputs=false` (default),
exclude input neurons from the resulting vector.
"""
function inhibitors(snn::AbstractSNN, n::Integer; inputs=false)::Vector{<:Integer}
    i = inneighbors(snn.graph, n)
    w = view(weights(snn.graph), i, n)
    if inputs
        return i[w .< 0]
    end
    return i[@. w < 0 & !snn.I[i]]
end

"""
    potential(snn::AbstractSNN, n::Integer)::Real

Calculate the mean of a neuron `n`'s inputs in the
current configuration of the network `snn`.
"""
function potential(snn::AbstractSNN, n::Integer)::Real
    # ensure the vertex is valid
    if snn.I[n]
        throw("Neuron $n is an input neuron in the SNN, so its potential must be provided.")
    end

    return mean(spikes(snn, inneighbors(snn.graph, n))) + bias(snn, n)
end

"""
    potential(snn::AbstractSNN{VT, ST ,G}, n::Integer)::Real where {VT<:Real, ST<:Integer, G<:AbstractSimpleWeightedGraph}

Calculate the weighted average of a neuron `n`'s inputs in the
current configuration of the weighted network `snn`.
"""
function potential(snn::AbstractSNN{VT, ST ,G}, n::Integer)::Real where {VT<:Real, ST<:Integer, G<:AbstractSimpleWeightedGraph}
    # ensure the vertex is valid
    if  snn.I[n]
        throw("Neuron $n is an input neuron in the SNN, so its potential must be provided.")
    end

    neighbors = inneighbors(snn.graph, n)
    return sum(spikes(snn, neighbors) .* weights(snn.graph)[neighbors, n]) + bias(snn, n)
end

"""
    sigmoid(snn::AbstractSNN, n::Integer, x::Real)::Real

The sigmoid spike transfer function.
"""
sigmoid(snn::AbstractSNN, n::Integer, x::Real)::Real = 1 / (1 + exp(-x))

"""
    thresh(x::Real, θ::Real)::Real

The threshold/step transfer function.
"""
thresh(x::Real, θ::Real)::Real = x >= θ ? 1 : 0

"""
    thresh(snn::AbstractSNN, n::Integer, x::Real)::Real

The step transfer function with a threshold of 0.
"""
thresh(snn::AbstractSNN, n::Integer, x::Real)::Real = thresh(x, 0)

"""
    thresh(snn::AbstractSNN, n::Integer, x::Real, θ::Real)::Real

The threshold/step transfer function.
"""
thresh(snn::AbstractSNN, n::Integer, x::Real, θ::Real)::Real = thresh(x, θ)

"""
    ramp(x::Real)::Real

The ramp activation function (a linear rectifier).
"""
ramp(x::Real)::Real = max(0, x)

"""
    ramp(snn::AbstractSNN, n::Integer, x::Real)::Real

The ramp activation function (a linear rectifier).
"""
ramp(snn::AbstractSNN, n::Integer, x::Real)::Real = ramp(x)

"""
    fire(prob::Real, on=1, off=0)::Real

Return `true` with probability `prob`, otherwise return `false`.
"""
fire(prob::Real)::Bool = rand() <= prob

"""
    increment_keys(d::Dict{K,V})::Dict{K,V} where {K<:Real, V}

Increment the listed keys in a `Dict` `d` by one(K).
"""
function increment_keys(d::Dict{K,V})::Dict{K,V} where {K<:Real, V}
    return Dict(zip(keys(d) .+ one(K), values(d)))
end

"""
    cycle!(snn::AbstractSNN{VT, ST, G}, s::Vector{<:ST}, v::Vector{<:VT})::AbstractSNN{VT,ST,G} where {VT<:Real, ST<:Integer, G<:AbstractGraph}

    Push `v` to the front of `snn`'s voltage memory, deleting any excess
    voltages at the end of `snn`'s memory.
    Add current spike times (zeros) to `snn`'s spike train memory if `s[n] = true`,
    and increment previous spike times by 1.
"""
function cycle!(snn::AbstractSNN{VT, ST, G}, s::AbstractVector{<:ST}, v::AbstractVector{<:VT}; threaded=false)::AbstractSNN{VT,ST,G} where {VT<:Real, ST<:Integer, G<:AbstractGraph}
    # add the new voltages
    if memory(snn) == 1
        @inbounds snn.V[1] = v
    elseif memory(snn) > 1
        push!(snn.V, v)
        while length(snn.V) > memory(snn)
            popfirst!(snn.V)
        end
    end

    # increment spike times and add the new spikes
    if threaded
        @inbounds Threads.@threads for n in neurons(snn)
            snn.S[n] = increment_keys(snn.S[n])
            if !iszero(s[n])
                snn.S[n][zero(ST)] = s[n]
            end
        end
    else
        @inbounds for n in neurons(snn)
            snn.S[n] = increment_keys(snn.S[n])
            if !iszero(s[n])
                snn.S[n][zero(ST)] = s[n]
            end
        end
    end
    return snn
end

"""
    step!(snn::AbstractSNN{VT,ST,G}, input::AbstractVector{<:ST}=zeros(ST, sum(snn.I));
          pot::Function=potential, transfer::Function=thresh, f=Bool, threaded=false)::AbstractSNN{VT,ST,G} where {VT<:Real, ST<:Integer, G<:AbstractGraph}

Calculate `snn`'s new firing configuration as `f(p(pot(snn, input)))`
and cycle `snn`'s memory accordingly.
"""
function step!(snn::AbstractSNN{VT,ST,G}, input::AbstractVector{<:ST}=zeros(ST, sum(snn.I)); pot::Function=potential, transfer::Function=thresh, f=Bool, threaded=false)::AbstractSNN{VT,ST,G} where {VT<:Real, ST<:Integer, G<:AbstractGraph}
    # Allocate vectors to store the new spikes/voltages
    s = Vector{ST}(undef, size(snn))
    v = Vector{VT}(undef, size(snn))

    # retrieve the new input configuration
    I = findall(snn.I)
    @inbounds v[I] .= zero(VT)
    @inbounds s[I] = input

    # compute the new configuration
    if threaded
        @inbounds Threads.@threads for n in findall(.!snn.I)
            v[n] = pot(snn, n)
            s[n] = f(transfer(snn, n, v[n]))
        end
    else
        @inbounds for n in findall(.!snn.I)
            v[n] = pot(snn, n)
            s[n] = f(transfer(snn, n, v[n]))
        end
    end

    # add the new configuration and remove excess configuration memory
    return cycle!(snn, s, v; threaded=threaded)
end

"""
    step!(snn::AbstractSNN{VT,ST,G}, input::AbstractMatrix{<:ST}, stateFns::Vararg{Function, N};
          pot::Function=potential, transfer::Function=thresh, f=Bool, threaded=false)::Tuple{AbstractSNN{VT,ST,G}, Vararg{Vector{Any}, N}} where {VT<:Real, ST<:Integer, G<:AbstractGraph, N}

For each column in `input[neuron, time]`, step! through `snn`. Return the updated
SNN and the history of the value of each function in `stateFns` through each step.
"""
function step!(snn::AbstractSNN{VT,ST,G}, input::AbstractMatrix{<:ST}, stateFns::Vararg{Function, N}; pot::Function=potential, transfer::Function=thresh, f=Bool, threaded=false)::Tuple{AbstractSNN{VT,ST,G}, Vararg{Vector{Any}, N}} where {VT<:Real, ST<:Integer, G<:AbstractGraph, N}
    T = size(input, 2)
    state = Tuple{Vararg{Vector{Any}, N}}(Vector{Any}(undef,T) for i in 1:length(stateFns))
    @inbounds for t in 1:T
        step!(snn, input[:,t], pot=pot, transfer=transfer, f=f, threaded=threaded)
        for i in eachindex(state)
            state[i][t] = stateFns[i](snn)
        end
    end

    return snn, state...
end

"""
    step!(snn::AbstractSNN{VT,ST,G}, inputIter, stateFns::Vararg{Function, N};
          pot::Function=potential, transfer::Function=sigmoid, f=fire, threaded=false)::Tuple{AbstractSNN{VT,ST,G}, Vararg{Vector{Any}, N}} where {VT<:Real, ST<:Integer, G<:AbstractGraph, N}

For each value in an iterable `inputIter`, step! through `snn`. Return the updated
SNN and the history of the value of each function in `stateFns` through each step.
"""
function step!(snn::AbstractSNN{VT,ST,G}, inputIter, stateFns::Vararg{Function, N}; pot::Function=potential, transfer::Function=thresh, f=Bool, threaded=false)::Tuple{AbstractSNN{VT,ST,G}, Vararg{Vector{Any}, N}} where {VT<:Real, ST<:Integer, G<:AbstractGraph, N}
    T = length(inputIter)
    state = Tuple{Vararg{Vector{Any}, N}}(Vector{Any}(undef,T) for i in 1:length(stateFns))
    @inbounds for (t, input) in zip(1:T, inputIter)
        step!(snn, input, pot=pot, transfer=transfer, f=f, threaded=threaded)
        for i in eachindex(state)
            state[i][t] = stateFns[i](snn)
        end
    end

    return snn, state...
end

"""
    step!(snn::AbstractSNN{VT,ST,G}, iterations::Integer, stateFns::Vararg{Function, N};
          pot::Function=potential, transfer::Function=sigmoid, f=fire, threaded=false)::Tuple{AbstractSNN{VT,ST,G}, Vararg{Vector{Any}, N}} where {VT<:Real, ST<:Integer, G<:AbstractGraph, N}

For each `iteration` iterations, step! through `snn` with constant input. Return
the updated SNN and the history of the value of each function in `stateFns` through each step.
"""
function step!(snn::AbstractSNN{VT,ST,G}, iterations::Integer, stateFns::Vararg{Function, N};
               input::Vector{<:ST}=zeros(ST, sum(snn.I)),
               pot::Function=potential, transfer::Function=thresh, f=Bool, threaded=false)::Tuple{AbstractSNN{VT,ST,G}, Vararg{Vector{Any}, N}} where {VT<:Real, ST<:Integer, G<:AbstractGraph, N}
    return step!(snn, constantly(iterations, input), stateFns...; pot=pot, transfer=transfer, f=f, threaded=threaded)
end

include("SNNs.jl")
include("LIFs.jl")
include("utils/GraphUtils.jl")
include("utils/SNNPlot.jl")
include("inputs/Inputs.jl")
end


text:
# reminder, to type H*, do H\^+
#cd("/Users/Frank/.julia/dev/Bridge")
outdir="output/bolus/"

using Bridge, StaticArrays, Distributions
using Bridge:logpdfnormal
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models
using DelimitedFiles
using DataFrames
using CSV
using RCall

sk = 0 # skipped in evaluating loglikelihoodν

νHparam = true
simlongpath = true
obs_scheme =["full","firstcomponent"][2]

# settings in case of νH - parametrisation
ϵ = 10^(-3)
Σdiagel = 10^(-5)

# settings sampler
iterations = 500
skip_it = 50# 1000
subsamples = 0:skip_it:iterations

ρ = 0.0#95

L = @SMatrix [.5 .5 ]
#L = @SMatrix [1. 0. ; 0.0 1.0]

m, d = size(L)
Σ = SMatrix{m,m}(Σdiagel*I)
dp = 2 # dprime

# specify target process
struct Diffusion <: ContinuousTimeProcess{ℝ{2}}
    α::Float64
    β::Float64
    λ::Float64
    k::Float64
    σ1::Float64
    σ2::Float64
end

# pdfemg(x, μ, σ, λ) = λ/2*exp(λ/2*(2μ + λ*σ^2 - 2x)).*erfc((μ + λ*σ^2 - x)/(sqrt(2)*σ))
# dose(t, c) = pdfemg(t, c...)*c[3]
# dose(t, c) = 1. *(t > c)

Bridge.b(t, x, P::Diffusion) = ℝ{2}(P.α*dose(t) -(P.λ + P.β)*x[1] + (P.k-P.λ)*x[2],  P.λ*x[1] - (P.k-P.λ)*x[2])
#Bridge.b(t, x, P::Diffusion) = ℝ{2}(0.0 , 0.0)#P.α*dose(t)
Bridge.σ(t, x, P::Diffusion) = @SMatrix [P.σ1 0.0 ;0.0  P.σ1]
Bridge.constdiff(::Diffusion) = true

# specify auxiliary process
struct DiffusionAux <: ContinuousTimeProcess{ℝ{2}}
    α::Float64
    β::Float64
    λ::Float64
    k::Float64
    #c::Float64
    σ1::Float64
    σ2::Float64
end

Random.seed!(42)
Bridge.B(t, P::DiffusionAux) = @SMatrix [ -P.λ - P.β P.k-P.λ ;  P.λ  P.λ-P.k]
Bridge.β(t, P::DiffusionAux) = ℝ{2}(P.α*dose(t),0.0)

#Bridge.B(t, P::DiffusionAux) = @SMatrix [ 0.0  0.0 ;  0.0 0.0]
#Bridge.β(t, P::DiffusionAux) = ℝ{2}(0.0,0.0)#P.α*dose(t)
Bridge.σ(t, P::DiffusionAux) = @SMatrix [P.σ1 0.0 ;0.0   P.σ2]
Bridge.constdiff(::DiffusionAux) = true
Bridge.b(t, x, P::DiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P)
Bridge.a(t, P::DiffusionAux) = Bridge.outer(Bridge.σ(t, P))

#dose(t) = 2*(t/20)/(1+(t/20)^2)
dose(t) = 2*(t/2)/(1+(t/2)^2)
FT = 70.; VB = 20.; PS = VE = 15.; HE = 0.4; DT = 2.4 # Favetto-Samson
DT = 1. # out choice

P = Diffusion(FT/(1-HE), FT/(VB*(1-HE)), PS/(VB*(1-HE)), PS/(VB*(1-HE)) + PS/VE, sqrt(2),0.2)
Pt = DiffusionAux(P.α, P.β,P.λ,P.k,P.σ1,P.σ2)

if simlongpath
    # Simulate one long path
    # Random.seed!(2)
    x0 = ℝ{2}(0.0, 0.0)
    #x0 = ℝ{2}(-8.0, 1.0)
    T_long = 2.0
    dt = 0.001
    tt_long = 0.:dt:T_long
    W_long = sample(tt_long, Wiener{ℝ{dp}}())
    X_long = solve(Euler(), x0, W_long, P)
    # Extract partial observations
    lt = length(tt_long)
    # obsnum = 10
    # if obsnum > 2
    #     obsind = sort(sample(2:lt-1,obsnum-2,replace=false))
    #     pushfirst!(obsind,1)
    #     push!(obsind,lt)
    # elseif obsnum==2
    #     obsind = [1, lt]
    # else
    #     error("provide valid number of observations ")
    # end
    obsnum = 10
    obsind = 1:(lt÷obsnum):lt
    obsnum = length(obsind)

    _pairs(collection) = Base.Generator(=>, keys(collection), values(collection))
    V_ = SamplePath(collect(_pairs(X_long))[obsind])
    V = SamplePath(V_.tt, map(y -> (L*y)[1] .+ (cholesky(Σ)).U' * randn(m), V_.yy))
end

obsnum = length(V)
segnum = obsnum-1

longpath = [Any[tt_long[j], d, X_long.yy[j][d]] for d in 1:2, j in 1:5:length(X_long) ][:]
obs = [Any[V.tt[j], dind, V.yy[j][dind]] for dind in 1:1, j in 1:length(V) ][:]
obsDf = DataFrame(time=map(x->x[1],obs), component = map(x->x[2],obs),value=map(x->x[3],obs) )
longpathDf = DataFrame(time=map(x->x[1],longpath), component = map(x->x[2],longpath),value=map(x->x[3],longpath) )

function gpupdate(ν::SVector, H⁺::SMatrix, Σ, L, v)
    if all(diag(H⁺) .== Inf)
        H⁺_ = SMatrix(inv(L' * inv(Σ) * L))
        V_ = (L' * inv(Σ) * L)\(L' * inv(Σ) *  v)
        return V_, H⁺_
    else
        Z = I - H⁺*L'*inv(Σ + L*H⁺*L')*L
        return SVector(Z*H⁺*L'*inv(Σ)*v + Z*ν), Z*H⁺
    end
end

tX = ℝ{d}
tW = ℝ{dp} #Float64 #ℝ{dp}# #
typeX = SamplePath{tX}
typeW = SamplePath{tW}
typePo = Bridge.PartialBridgeνH{ℝ{d},Diffusion,DiffusionAux,ℝ{d},SArray{Tuple{d,d},Float64,2,d^2}}
XX = Vector{typeX}(undef, segnum)
XX2 = Vector{typeX}(undef, segnum)
WW = Vector{typeW}(undef, segnum)
Po = Vector{typePo}(undef, segnum)
Poᵒ = Vector{typePo}(undef, segnum) # needed when parameter estimation is done

# initialise
νend = SVector{d}(zeros(d))
Hend⁺ = SMatrix{d,d}(I/ϵ)
H⁺i = Vector{typeof(Hend⁺)}(undef, segnum)

νend, Hend⁺ = gpupdate(νend, Hend⁺ , Σ,L,V.yy[end])

dtimp = 0.0001  # mesh width for imputed paths

τ(t, T0, Tend) = T0 +(t-T0)*(2-(t-T0)/(Tend-T0))

# solve backward recursion on [0,T]

for i in segnum:-1:1
    # update on interval (t[i-1],t[i])
    tt_ = τ.(V.tt[i]:dtimp:V.tt[i+1],V.tt[i],V.tt[i+1])
    XX[i] = Bridge.samplepath(tt_, zero(tX)) # initialise
    XX2[i] = Bridge.samplepath(tt_, zero(tX)) # initialise
    WW[i] = Bridge.samplepath(tt_, zero(tW)) # initialise
    H⁺i[i] = Hend⁺
    Po[i], νend, Hend⁺ = Bridge.partialbridgeνH(tt_, P, Pt, νend, Hend⁺)
    Poᵒ[i] = Po[i]
    global νend, Hend⁺ = gpupdate(νend, Hend⁺, Σ, L, V.yy[i])
end
H⁺i[1] = Hend⁺

#elapsed_time= @elapsed begin

####################### MH algorithm ###################

# Initialisation, forward simulate on [0,T] a bridge
xstart = νend #+ √(Hend⁺) * randn(d)  # note that this is really ν(0)
for i in 1:segnum
    tt = Po[i].tt
    WW[i] = sample(tt,Wiener{tW}())
    Bridge.solve!(Euler(), XX[i], xstart, WW[i], Po[i])
    xstart = XX[i].yy[end] # starting point for next segment
end

XXinit = vcat(XX...)
longpath = [Any[tt_long[j], d, X_long.yy[j][d]] for d in 1:2, j in 1:5:length(X_long) ][:]
initpath = [Any[XXinit.tt[j], d, XXinit.yy[j][d]] for d in 1:2, j in 1:5:length(XXinit) ][:]
obs = [Any[V.tt[j], dind, V.yy[j][dind]] for dind in 1:1, j in 1:length(V) ][:]
obsDf = DataFrame(time=map(x->x[1],obs), component = map(x->x[2],obs),value=map(x->x[3],obs) )
longpathDf = DataFrame(time=map(x->x[1],longpath), component = map(x->x[2],longpath),value=map(x->x[3],longpath) )
initpathDf = DataFrame(time=map(x->x[1],initpath), component = map(x->x[2],initpath),value=map(x->x[3],initpath) )

@rput obsDf
@rput longpathDf
@rput initpathDf

R"""
library(ggplot2)
library(tidyverse)

longpathDf$component <- as.factor(longpathDf$component)
p <- ggplot() +
  ylab("") + geom_path(aes(x=time,y=value,colour=as.factor(component)),data=initpathDf)+
   geom_point(aes(x=time,y=value),data=obsDf,colour="red")+
  geom_path(aes(x=time,y=value,colour=component),data=longpathDf)+theme_minimal()+
  theme(legend.position="bottom")
  #facet_wrap(~component,ncol=1,scales='free_y') +

"""

XXo = deepcopy(XX)
WWo = deepcopy(WW)
#WW2 = deepcopy(WW)

# save some of the paths
XXsave = Any[]
if 0 in subsamples
    push!(XXsave, deepcopy(XX))
end

acc = 0
accparams = 0
mhsteps = 0
mhstepsparams = 0
Hrightmost⁺ = H⁺i[segnum]
νrightmost = Po[segnum].ν[end]
Hzero⁺ = SMatrix{d,d}(0.01*I)

C = [0.0]

for iter in 1:iterations

    finished = false
    klow = 1
    ind = segnum:-1:1
    kup = obsnum

    xstart = XX[1].yy[1]
    xstarto = xstart
    updateparams =  false#
    while !finished
        # update a block

        if updateparams
            ind = (segnum==2) ? (1:1) : (segnum:-1:1)
            kup = obsnum
        else
            segnum_update = sample(1:obsnum-klow) #obsnum-1   ##obsnum-1#   # number of segments to update
            kup = klow + segnum_update  # update on interval [t(klow), t(kup)]
            ind = (segnum_update==1) ? (1:1) : ((kup-1):-1:klow) # indices of segments to update
        end
        hasbegin = ind[end]==1
        hasend = ind[1]==segnum

        νend  = hasend ? νrightmost :  XX[ind[1]].yy[end]  # initialise on rightmost segment
        Hend⁺ = hasend ? Hrightmost⁺ : Hzero⁺ # initialise on rightmost segment
        νendᵒ = νend
        Hend⁺ᵒ = Hend⁺

        α = P.α
        if updateparams
            αᵒ = α + 0.01 * randn()
        else
            αᵒ = α
        end

        Pᵒ = Diffusion(αᵒ,P.β,P.λ,P.k,P.σ1,P.σ2)
        Ptᵒ = DiffusionAux(αᵒ,P.β,P.λ,P.k,P.σ1,P.σ2)
        Pt = DiffusionAux(α, P.β,P.λ,P.k,P.σ1,P.σ2)

        # compute guiding term
        for i in ind
            tt = Po[i].tt
            Po[i], νend, Hend⁺ = Bridge.partialbridgeνH(tt, P, Pt, νend, Hend⁺)
            νend, Hend⁺ = gpupdate(νend, Hend⁺, Σ, L, V.yy[i])
            if updateparams
                Poᵒ[i], νendᵒ, Hend⁺ᵒ = Bridge.partialbridgeνH(tt, Pᵒ, Ptᵒ, νendᵒ, Hend⁺ᵒ)
                νendᵒ, Hend⁺ᵒ = gpupdate(νendᵒ, Hend⁺ᵒ, Σ, L, V.yy[i])
            else
                Poᵒ[i] = Po[i]
                νendᵒ, Hend⁺ᵒ = νend, Hend⁺
            end
        end

        # simulate guided proposal
        for i in reverse(ind)
            tt = Po[i].tt
            if !updateparams
                sample!(WWo[i], Wiener{ℝ{2}}())
                WWo[i].yy .= ρ * WW[i].yy + sqrt(1-ρ^2) * WWo[i].yy
            else
                WWo[i].yy .= WW[i].yy
            end
            if updateparams
                xstarto = xstart = XX[1].yy[1]
            elseif i==1
                xstart = XX[1].yy[1]
                u = randn()
                xstarto = xstart + 0.1 * ℝ{2}(u, -u)
            else
                xstarto = xstart = XX[i-1].yy[end]
            end

            # at this point, either WW = WWo or Po == Poᵒ (if updateparams=true)
            solve!(Euler(), XX2[i], xstart, WW[i], Po[i])
            solve!(Euler(), XXo[i], xstarto, WWo[i], Poᵒ[i])
        end
        # compute loglikelihood
        diffll = 0.0
        for i in ind
            diffll += llikelihood(LeftRule(), XXo[i],  Poᵒ[i]) - llikelihood(LeftRule(), XX2[i],  Po[i])
        end
        if hasbegin
             diffll += logpdfnormal(xstarto-νendᵒ, Bridge.symmetrize(Hend⁺ᵒ))-logpdfnormal(xstart-νend, Bridge.symmetrize(Hend⁺))               # plus possibly log q(X0|X0o) = log q(X0o|X0)
        end
        if updateparams
            diffll += -(V.tt[end]-V.tt[1]) *(trace(Bridge.B(0.0, Ptᵒ)) - trace(Bridge.B(0.0,Pt))) +
            logpdf(Gamma(2,1),cᵒ) - logpdf(Gamma(2,1),c)
        end

        print("iter  diff_ll: ",round(diffll, digits=3))
        # MH step
        if log(rand()) <= diffll
            print("✓")
            for i in ind
                 XX[i], XXo[i] = XXo[i], XX[i]
                 WW[i], WWo[i] = WWo[i], WW[i]
            end
            P = Pᵒ
            Pt = Ptᵒ
            push!(C, αᵒ)
            if updateparams
                accparams += 1
            else
                acc += 1

            end
        end
        klow = kup # adjust counter
        finished = (klow==obsnum) # all segments have been updated (this is correct when a single sweep is done for parameter update)
        mhsteps += !updateparams
        mhstepsparams += updateparams
    end
    println()
    if iter in subsamples
        push!(XXsave, deepcopy(XX))
    end

end

@info "Done."*"\x7"^6

ave_acc_perc = 100*round(acc/mhsteps, digits=2)

limp = length(vcat(XXsave[1]...))
iterates = [Any[s,  vcat(XXsave[i]...).tt[j], dind, vcat(XXsave[i]...).yy[j][dind]] for dind in 1:d, j in 1:10:limp, (i,s) in enumerate(subsamples) ][:]
iteratesaverage = [Any[s,  vcat(XXsave[i]...).tt[j], mean(vcat(XXsave[i]...).yy[j])] for j in 1:10:limp, (i,s) in enumerate(subsamples) ][:]

write2csv = false
if write2csv
    # write long path to csv file
    f = open(outdir*"longforward_pendulum.csv","w")
    headl = "time, component, value \n"
    write(f, headl)
    writedlm(f, longpath, ',')
    close(f)
    # write mcmc iterates to csv file
    fn = outdir*"iterates-"*obs_scheme*".csv"
    f = open(fn,"w")
    headl = "iteration, time, component, value \n"
    write(f, headl)
    writedlm(f,iterates,",")
    close(f)
    # write observations to csv file
    fn = outdir*"observations.csv"
    f = open(fn,"w")
    headl = "time, component, value \n"
    write(f, headl)
    writedlm(f,obs,",")
    close(f)
end

println("Average acceptance percentage: ",ave_acc_perc,"\n")

# plotting the results
iteratesDf = DataFrame(iteration = map(x->x[1],iterates), time=map(x->x[2],iterates), component = map(x->x[3],iterates),value=map(x->x[4],iterates) )
iteratesaverageDf = DataFrame(iteration = map(x->x[1],iteratesaverage), time=map(x->x[2],iteratesaverage), value = map(x->x[3],iteratesaverage) )
@rput iteratesDf
@rput iteratesaverageDf

R"""
library(ggplot2)
library(tidyverse)
iteratesDfsub <- iteratesDf #%>% filter(iteration > 5)
p <- ggplot() +
  geom_path(mapping=aes(x=time,y=value,colour=iteration,group=iteration),data=iteratesDfsub) +
  scale_colour_gradient(low='green',high='blue')+
   ylab("") + geom_point(aes(x=time,y=value),data=obsDf,colour="red")+
   geom_path(aes(x=time,y=value),data=longpathDf,colour="yellow")+
   facet_wrap(~component,ncol=1,scales='free_y') +
theme_minimal()

show(p)
"""

R"""
ggplot() +
  geom_path(mapping=aes(x=time,y=value,colour=iteration,group=iteration),data=iteratesaverageDf) +
  scale_colour_gradient(low='green',high='blue')+
   ylab("") + geom_point(aes(x=time,y=value),data=obsDf,colour="red")+
   #geom_path(aes(x=time,y=value),data=longpathDf,colour="yellow")+
   theme_minimal()
"""

writeinfo = true
if writeinfo
    # write info to txt file
    fn = outdir*"info.txt"
    f = open(fn,"w")
    #write(f, "Elapsed time: ", string(elapsed_time),"\n")
    #write(f, "Choice of observation schemes: ",obs_scheme,"\n")
    #write(f, "Easy conditioning (means going up to 1 for the rough component instead of 2): ",string(easy_conditioning),"\n")
    write(f, "Number of iterations: ",string(iterations),"\n")
    write(f, "Skip every ",string(skip_it)," iterations, when saving to csv","\n\n")
    write(f, "Starting point: ",string(x0),"\n")
    write(f, "End time T: ", string(V.tt[end]),"\n")
    #write(f, "Endpoint v: ",string(v),"\n")
    write(f, "Noise Sigma: ",string(Σ),"\n")
    write(f, "Regularisation parameter epsilon", string(ϵ),"\n")
    write(f, "L: ",string(L),"\n\n")
    write(f, "Mesh width: ",string(dtimp),"\n")
    write(f, "rho (Crank-Nicholsen parameter: ",string(ρ),"\n")
    write(f, "Average acceptance percentage: ",string(ave_acc_perc),"\n\n")
    write(f, "Backward type parametrisation in terms of nu and H? ",string(νHparam),"\n")
    close(f)
end

ext: jl


text:
module PendulumWithFixSimulation

using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

Bar = Model(
    m  = 0.1,
    Lx = 0.1,
    Ly = Par(value=:(0.2*Lx)),
    Lz = Par(value=:(0.2*Lx)),
    vmat1 = VisualMaterial(color="Teal", transparency=0.5),
    vmat2 = VisualMaterial(color="Red"),
    frame0 = Object3D(feature=Solid(shape=Beam(axis=1, length=:Lx, width=:Ly, thickness=:Lz),
                                    massProperties=MassProperties(mass=:m),
                                    visualMaterial=:(vmat1))),
	frame1 = Object3D(parent=:frame0,
                      translation=:[-Lx/2, 0.0, 0.0],
                      feature=Visual(shape=Cylinder(axis=3, diameter=:(Ly/2), length=:(1.2*Lz)),
                                     visualMaterial=:(vmat2))),
    frame2 = Object3D(parent=:frame0, translation=:[Lx/2, 0.0, 0.0])
)

Pendulum = Model(
    m = 1.0,
    Lx = 0.1,
    world = Object3D(feature=Scene() ),
    worldFrame = Object3D(parent=:world,
                          feature=Visual(shape=CoordinateSystem(length=:(Lx/2)))),
    bar0 = Bar | Map(m=:m, Lx=:Lx),
    rev = RevoluteWithFlange(obj1=:world, obj2=:(bar0.frame1)),
    bar1 = Bar | Map(m=:m, Lx=:Lx),
    fix1 = Fix(obj1=:(bar0.frame2), obj2=:(bar1.frame1), translation=:[0.0, -Lx/4, 0.0]),
    bar2 = Bar | Map(m=:m, Lx=:Lx),
    fix2 = Fix(obj1=:(bar0.frame2), obj2=:(bar2.frame1), translation=:[0.0, Lx/4, 0.0], rotation=[0.0, 0.0, 45*pi/180]),
    capsule1 = Object3D(feature=Solid(shape=Capsule(axis=3, length=:(Lx/2), diameter=:(Lx/4)), solidMaterial="DryWood", visualMaterial="BlueTransparent")),
    fix3 = Fix(obj1=:(bar1.frame2), obj2=:capsule1, translation=:[Lx/4, 0.0, 0.0], rotation=[90*pi/180, 0.0, 0.0]),
    capsule2 = Object3D(feature=Solid(shape=Capsule(axis=3, length=:(Lx/2), diameter=:(Lx/4)), solidMaterial="DryWood", visualMaterial="GreenTransparent")),
    fix4 = Fix(obj1=:capsule2, obj2=:(bar2.frame2), translation=:[-Lx/4, 0.0, 0.0], rotation=[-90*pi/180, 0.0, 0.0])
)

PendulumWithFix = Model(
    pendulum = buildModia3D(Pendulum | Map(Lx=1.0, m=2.0, rev=Map(phi=Var(init=1.0)))),

    damper = Damper | Map(d=100.0),
    fixed = Fixed,
    connect = :[(damper.flange_b, pendulum.rev.flange),
                (damper.flange_a, fixed.flange)]
)

pendulumWithFix = @instantiateModel(PendulumWithFix, unitless=true)

stopTime = 10.0
requiredFinalStates = [-1.6164397926928633, 0.4670955914743122]
simulate!(pendulumWithFix, stopTime=stopTime, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(pendulumWithFix, ["pendulum.rev.flange.phi", "pendulum.rev.variables[1]"], figure=1)

end


text:
#--------------------
#--------------------
#### Symbolic
#--------------------
abstract type Symbolic{T} end

###
### Uni-type design
###

@enum ExprType::UInt8  SYM TERM ADD MUL POW DIV

const Metadata = Union{Nothing,Base.ImmutableDict{DataType,Any}}
const NO_METADATA = nothing

sdict(kv...) = Dict{Any, Any}(kv...)

using Base: RefValue
const EMPTY_ARGS = []
const EMPTY_HASH = RefValue(UInt(0))
const NOT_SORTED = RefValue(false)
const EMPTY_DICT = sdict()
const EMPTY_DICT_T = typeof(EMPTY_DICT)

@compactify show_methods=false begin
    @abstract struct BasicSymbolic{T} <: Symbolic{T}
        metadata::Metadata     = NO_METADATA
    end
    struct Sym{T} <: BasicSymbolic{T}
        name::Symbol           = :OOF
    end
    struct Term{T} <: BasicSymbolic{T}
        f::Any                 = identity  # base/num if Pow; issorted if Add/Dict
        arguments::Vector{Any} = EMPTY_ARGS
        hash::RefValue{UInt}   = EMPTY_HASH
    end
    struct Mul{T} <: BasicSymbolic{T}
        coeff::Any             = 0         # exp/den if Pow
        dict::EMPTY_DICT_T     = EMPTY_DICT
        hash::RefValue{UInt}   = EMPTY_HASH
        arguments::Vector{Any} = EMPTY_ARGS
        issorted::RefValue{Bool} = NOT_SORTED
    end
    struct Add{T} <: BasicSymbolic{T}
        coeff::Any             = 0         # exp/den if Pow
        dict::EMPTY_DICT_T     = EMPTY_DICT
        hash::RefValue{UInt}   = EMPTY_HASH
        arguments::Vector{Any} = EMPTY_ARGS
        issorted::RefValue{Bool} = NOT_SORTED
    end
    struct Div{T} <: BasicSymbolic{T}
        num::Any               = 1
        den::Any               = 1
        simplified::Bool       = false
        arguments::Vector{Any} = EMPTY_ARGS
    end
    struct Pow{T} <: BasicSymbolic{T}
        base::Any              = 1
        exp::Any               = 1
        arguments::Vector{Any} = EMPTY_ARGS
    end
end

function exprtype(x::BasicSymbolic)
    @compactified x::BasicSymbolic begin
        Term => TERM
        Add  => ADD
        Mul  => MUL
        Div  => DIV
        Pow  => POW
        Sym  => SYM
        _    => error_on_type()
    end
end

# Same but different error messages
@noinline error_on_type() = error("Internal error: unreachable reached!")
@noinline error_sym() = error("Sym doesn't have a operation or arguments!")
@noinline error_property(E, s) = error("$E doesn't have field $s")

# We can think about bits later
# flags
const SIMPLIFIED = 0x01 << 0

#@inline is_of_type(x::BasicSymbolic, type::UInt8) = (x.bitflags & type) != 0x00
#@inline issimplified(x::BasicSymbolic) = is_of_type(x, SIMPLIFIED)

function ConstructionBase.setproperties_object(obj::BasicSymbolic{T}, patch)::BasicSymbolic{T} where T
    nt = getproperties(obj)
    nt_new = merge(nt, patch)
    Unityper.rt_constructor(obj){T}(;nt_new...)
end

###
### TermInterface
###
using TermInterface

TermInterface.exprhead(x::Symbolic) = :call
TermInterface.symtype(x::Number) = typeof(x)
@inline TermInterface.symtype(::Symbolic{T}) where T = T

# We're returning a function pointer
@inline function TermInterface.operation(x::BasicSymbolic)
    @compactified x::BasicSymbolic begin
        Term => x.f
        Add  => (+)
        Mul  => (*)
        Div  => (/)
        Pow  => (^)
        Sym  => error_sym()
        _    => error_on_type()
    end
end

function TermInterface.arguments(x::BasicSymbolic)
    args = unsorted_arguments(x)
    @compactified x::BasicSymbolic begin
        Add => @goto ADDMUL
        Mul => @goto ADDMUL
        _   => return args
    end
    @label ADDMUL
    if !x.issorted[]
        sort!(args, lt = <ₑ)
        x.issorted[] = true
    end
    return args
end
function TermInterface.unsorted_arguments(x::BasicSymbolic)
    @compactified x::BasicSymbolic begin
        Term => return x.arguments
        Add  => @goto ADDMUL
        Mul  => @goto ADDMUL
        Div  => @goto DIV
        Pow  => @goto POW
        Sym  => error_sym()
        _    => error_on_type()
    end

    @label ADDMUL
    E = exprtype(x)
    args = x.arguments
    isempty(args) || return args
    siz = length(x.dict)
    idcoeff = E === ADD ? iszero(x.coeff) : isone(x.coeff)
    sizehint!(args, idcoeff ? siz : siz + 1)
    idcoeff || push!(args, x.coeff)
    if isadd(x)
        for (k, v) in x.dict
            push!(args, k * v)
        end
    else # MUL
        for (k, v) in x.dict
            push!(args, unstable_pow(k, v))
        end
    end
    return args

    @label DIV
    args = x.arguments
    isempty(args) || return args
    sizehint!(args, 2)
    push!(args, x.num)
    push!(args, x.den)
    return args

    @label POW
    args = x.arguments
    isempty(args) || return args
    sizehint!(args, 2)
    push!(args, x.base)
    push!(args, x.exp)
    return args
end

TermInterface.istree(s::BasicSymbolic) = !issym(s)
@inline isa_SymType(T::Val{S}, x) where {S} = x isa BasicSymbolic ? Unityper.isa_type_fun(Val(SymbolicUtils.BasicSymbolic), T, x) : false
TermInterface.issym(x::BasicSymbolic) = isa_SymType(Val(:Sym), x)
isterm(x) = isa_SymType(Val(:Term), x)
ismul(x)  = isa_SymType(Val(:Mul), x)
isadd(x)  = isa_SymType(Val(:Add), x)
ispow(x)  = isa_SymType(Val(:Pow), x)
isdiv(x)  = isa_SymType(Val(:Div), x)

###
### Base interface
###

Base.isequal(::Symbolic, x) = false
Base.isequal(x, ::Symbolic) = false
Base.isequal(::Symbolic, ::Symbolic) = false

function Base.isequal(a::BasicSymbolic{T}, b::BasicSymbolic{S}) where {T,S}
    a === b && return true

    E = exprtype(a)
    E === exprtype(b) || return false

    T === S || return false

    if E === SYM
        nameof(a) === nameof(b)
    elseif E === ADD || E === MUL
        a.coeff == b.coeff && isequal(a.dict, b.dict)
    elseif E === DIV
        isequal(a.num, b.num) && isequal(a.den, b.den)
    elseif E === POW
        isequal(a.exp, b.exp) && isequal(a.base, b.base)
    elseif E === TERM
        a1 = arguments(a)
        a2 = arguments(b)
        isequal(operation(a), operation(b)) &&
            length(a1) == length(a2) &&
            all(isequal(l, r) for (l, r) in zip(a1, a2))
    else
        error_on_type()
    end
end

Base.one( s::Symbolic) = one( symtype(s))
Base.zero(s::Symbolic) = zero(symtype(s))

Base.nameof(s::BasicSymbolic) = issym(s) ? s.name : error("None Sym BasicSymbolic doesn't have a name")

## This is much faster than hash of an array of Any
hashvec(xs, z) = foldr(hash, xs, init=z)
function Base.hash(s::BasicSymbolic, salt::UInt)
    E = exprtype(s)
    if E === SYM
        hash(nameof(s), salt ⊻ 0x4de7d7c66d41da43)
    elseif E === ADD || E === MUL
        !iszero(salt) && return hash(hash(s, zero(UInt64)), salt)
        h = s.hash[]
        !iszero(h) && return h
        hashoffset = isadd(s) ? 0xaddaddaddaddadda : 0xaaaaaaaaaaaaaaaa
        h′= hash(hashoffset, hash(s.coeff, hash(s.dict, salt)))
        s.hash[] = h′
        return h′
    elseif E === DIV
        return hash(s.num, hash(s.den, salt ⊻ 0x334b218e73bbba53))
    elseif E === POW
        hash(s.exp, hash(s.base, salt ⊻ 0x2b55b97a6efb080c))
    elseif E === TERM
        !iszero(salt) && return hash(hash(s, zero(UInt)), salt)
        h = s.hash[]
        !iszero(h) && return h
        op = operation(s)
        oph = op isa Function ? nameof(op) : op
        h′ = hashvec(arguments(s), hash(oph, salt))
        s.hash[] = h′
        return h′
    else
        error_on_type()
    end
end

###
### Constructors
###

function Sym{T}(name::Symbol; kw...) where T
    Sym{T}(; name=name, kw...)
end

function Term{T}(f, args; kw...) where T
    if eltype(args) !== Any
        args = convert(Vector{Any}, args)
    end

    Term{T}(;f=f, arguments=args, hash=Ref(UInt(0)), kw...)
end

function Term(f, args; metadata=NO_METADATA)
    Term{_promote_symtype(f, args)}(f, args, metadata=metadata)
end

function Add(::Type{T}, coeff, dict; metadata=NO_METADATA, kw...) where T
    if isempty(dict)
        return coeff
    elseif _iszero(coeff) && length(dict) == 1
        k,v = first(dict)
        if _isone(v)
            return k
        else
            coeff, dict = makemul(v, k)
            return Mul(T, coeff, dict)
        end
    end

    Add{T}(; coeff, dict, hash=Ref(UInt(0)), metadata, arguments=[], issorted=RefValue(false), kw...)
end

function Mul(T, a, b; metadata=NO_METADATA, kw...)
    isempty(b) && return a
    if _isone(a) && length(b) == 1
        pair = first(b)
        if _isone(last(pair)) # first value
            return first(pair)
        else
            return unstable_pow(first(pair), last(pair))
        end
    else
        coeff = a
        dict = b
        Mul{T}(; coeff, dict, hash=Ref(UInt(0)), metadata, arguments=[], issorted=RefValue(false), kw...)
    end
end

const Rat = Union{Rational, Integer}

function ratcoeff(x)
    if ismul(x)
        ratcoeff(x.coeff)
    elseif x isa Rat
        (true, x)
    else
        (false, NaN)
    end
end
ratio(x::Integer,y::Integer) = iszero(rem(x,y)) ? div(x,y) : x//y
ratio(x::Rat,y::Rat) = x//y
function maybe_intcoeff(x)
    if ismul(x)
        if x.coeff isa Rational && isone(x.coeff.den)
            Mul{symtype(x)}(; coeff=x.coeff.num, dict=x.dict, x.metadata, arguments=[], issorted=RefValue(false))
        else
            x
        end
    elseif x isa Rational
        isone(x.den) ? x.num : x
    else
        x
    end
end

function Div{T}(n, d, simplified=false; metadata=nothing) where {T}
    if T<:Number && !(T<:SafeReal)
        n, d = quick_cancel(n, d)
    end
    _iszero(n) && return zero(typeof(n))
    _isone(d) && return n

    if isdiv(n) && isdiv(d)
        return Div{T}(n.num * d.den, n.den * d.num)
    elseif isdiv(n)
        return Div{T}(n.num, n.den * d)
    elseif isdiv(d)
        return Div{T}(n * d.den, d.num)
    end

    d isa Number && _isone(-d) && return -1 * n
    n isa Rat && d isa Rat && return n // d # maybe called by oblivious code in simplify

    # GCD coefficient upon construction
    rat, nc = ratcoeff(n)
    if rat
        rat, dc = ratcoeff(d)
        if rat
            g = gcd(nc, dc) * sign(dc) # make denominator positive
            invdc = ratio(1, g)
            n = maybe_intcoeff(invdc * n)
            d = maybe_intcoeff(invdc * d)
            if d isa Number
                _isone(d) && return n
                _isone(-d) && return -1 * n
            end
        end
    end

    Div{T}(; num=n, den=d, simplified, arguments=[], metadata)
end

function Div(n,d, simplified=false; kw...)
    Div{promote_symtype((/), symtype(n), symtype(d))}(n, d, simplified; kw...)
end

@inline function numerators(x)
    isdiv(x) && return numerators(x.num)
    istree(x) && operation(x) === (*) ? arguments(x) : Any[x]
end

@inline denominators(x) = isdiv(x) ? numerators(x.den) : Any[1]

function (::Type{<:Pow{T}})(a, b; metadata=NO_METADATA) where {T}
    _iszero(b) && return 1
    _isone(b) && return a
    Pow{T}(; base=a, exp=b, arguments=[], metadata)
end

function Pow(a, b; metadata=NO_METADATA)
    Pow{promote_symtype(^, symtype(a), symtype(b))}(makepow(a, b)..., metadata=metadata)
end

function toterm(t::BasicSymbolic{T}) where T
    E = exprtype(t)
    if E === SYM || E === TERM
        return t
    elseif E === ADD || E === MUL
        args = []
        for (k, coeff) in t.dict
            push!(args, coeff == 1 ? k : Term{T}(E === MUL ? (^) : (*), Any[coeff, k]))
        end
        Term{T}(operation(t), args)
    elseif E === DIV
        Term{T}(/, Any[t.num, t.den])
    elseif E === POW
        Term{T}(^, [t.base, t.exp])
    else
        error_on_type()
    end
end

"""
    makeadd(sign, coeff::Number, xs...)

Any Muls inside an Add should always have a coeff of 1
and the key (in Add) should instead be used to store the actual coefficient
"""
function makeadd(sign, coeff, xs...)
    d = sdict()
    for x in xs
        if isadd(x)
            coeff += x.coeff
            _merge!(+, d, x.dict, filter=_iszero)
            continue
        end
        if x isa Number
            coeff += x
            continue
        end
        if ismul(x)
            k = Mul(symtype(x), 1, x.dict)
            v = sign * x.coeff + get(d, k, 0)
        else
            k = x
            v = sign + get(d, x, 0)
        end
        if iszero(v)
            delete!(d, k)
        else
            d[k] = v
        end
    end
    coeff, d
end

function makemul(coeff, xs...; d=sdict())
    for x in xs
        if ispow(x) && x.exp isa Number
            d[x.base] = x.exp + get(d, x.base, 0)
        elseif x isa Number
            coeff *= x
        elseif ismul(x)
            coeff *= x.coeff
            _merge!(+, d, x.dict, filter=_iszero)
        else
            v = 1 + get(d, x, 0)
            if _iszero(v)
                delete!(d, x)
            else
                d[x] = v
            end
        end
    end
    (coeff, d)
end

unstable_pow(a, b) = a isa Integer && b isa Integer ? (a//1) ^ b : a ^ b

function makepow(a, b)
    base = a
    exp = b
    if ispow(a)
        base = a.base
        exp = a.exp * b
    end
    return (base, exp)
end

function term(f, args...; type = nothing)
    if type === nothing
        T = _promote_symtype(f, args)
    else
        T = type
    end
    Term{T}(f, Any[args...])
end

"""
    unflatten(t::Symbolic{T})
Binarizes `Term`s with n-ary operations
"""
function unflatten(t::Symbolic{T}) where{T}
    if istree(t)
        f = operation(t)
        if f == (+) || f == (*)   # TODO check out for other n-ary --> binary ops
            a = arguments(t)
            return foldl((x,y) -> Term{T}(f, Any[x, y]), a)
        end
    end
    return t
end

unflatten(t) = t

"""
    similarterm(t, f, args, symtype; metadata=nothing)

Create a term that is similar in type to `t`. Extending this function allows packages
using their own expression types with SymbolicUtils to define how new terms should
be created. Note that `similarterm` may return an object that has a
different type than `t`, because `f` also influences the result.

## Arguments

- `t` the reference term to use to create similar terms
- `f` is the operation of the term
- `args` is the arguments
- The `symtype` of the resulting term. Best effort will be made to set the symtype of the
  resulting similar term to this type.
"""
TermInterface.similarterm(t::Type{<:Symbolic}, f, args; metadata=nothing, exprhead=:call) =
    similarterm(t, f, args, _promote_symtype(f, args); metadata=metadata, exprhead=exprhead)

TermInterface.similarterm(t::Type{<:Symbolic}, f::Symbol, args; metadata=nothing, exprhead=:call) =
    TermInterface.similarterm(t, eval(f), args; metadata=metadata, exprhead=exprhead)

function TermInterface.similarterm(t::Type{<:BasicSymbolic{<:Number}}, f, args, symtype; metadata=nothing, exprhead=:call)
    if f isa Symbol
        return Term{_promote_symtype(eval(f), args)}(eval(f), args; metadata=metadata)
    end
    T = symtype
    if T === nothing
        T = _promote_symtype(f, args)
    end
    if f === (+)
        Add(T, makeadd(1, 0, args...)...; metadata=metadata)
    elseif f == (*)
        Mul(T, makemul(1, args...)...; metadata=metadata)
    elseif f == (/)
        @assert length(args) == 2
        Div{T}(args...; metadata=metadata)
    elseif f == (^) && length(args) == 2
        Pow{T}(makepow(args...)...; metadata=metadata)
    else
        Term{T}(f, args, metadata=metadata)
    end
end

###
### Tree print
###

import AbstractTrees

struct TreePrint
    op
    x
end

function AbstractTrees.children(x::BasicSymbolic)
    if isterm(x) || ispow(x)
        return arguments(x)
    elseif isadd(x) || ismul(x)
        children = Any[x.coeff]
        for (key, coeff) in pairs(x.dict)
            if coeff == 1
                push!(children, key)
            else
                push!(children, TreePrint(isadd(x) ? (:*) : (:^), (key, coeff)))
            end
        end
        return children
    end
end

AbstractTrees.children(x::TreePrint) = [x.x[1], x.x[2]]

print_tree(x; show_type=false, maxdepth=Inf, kw...) = print_tree(stdout, x; show_type=show_type, maxdepth=maxdepth, kw...)

function print_tree(_io::IO, x::BasicSymbolic; show_type=false, kw...)
    if isterm(x) || isadd(x) || ismul(x) || ispow(x) || isdiv(x)
        AbstractTrees.print_tree(_io, x; withinds=true, kw...) do io, y, inds
            if istree(y)
                print(io, operation(y))
            elseif y isa TreePrint
                print(io, "(", y.op, ")")
            else
                print(io, y)
            end
            if !(y isa TreePrint) && show_type
                print(io, " [", typeof(y), "]")
            end
        end
    end
end

# I think we already have this defined on 201?
#=
function TermInterface.istree(t::Type{<:Symbolic})
    TermInterface.issym(t) ? false : true
end
=#

###
### Metadata
###
TermInterface.metadata(s::Symbolic) = s.metadata
TermInterface.metadata(s::Symbolic, meta) = Setfield.@set! s.metadata = meta

function hasmetadata(s::Symbolic, ctx)
    metadata(s) isa AbstractDict && haskey(metadata(s), ctx)
end

function getmetadata(s::Symbolic, ctx)
    md = metadata(s)
    if md isa AbstractDict
        md[ctx]
    else
        throw(ArgumentError("$s does not have metadata for $ctx"))
    end
end

function getmetadata(s::Symbolic, ctx, default)
    md = metadata(s)
    md isa AbstractDict ? get(md, ctx, default) : default
end

# pirated for Setfield purposes:
using Base: ImmutableDict
Base.ImmutableDict(d::ImmutableDict{K,V}, x, y)  where {K, V} = ImmutableDict{K,V}(d, x, y)

assocmeta(d::Dict, ctx, val) = (d=copy(d); d[ctx] = val; d)
function assocmeta(d::Base.ImmutableDict, ctx, val)::ImmutableDict{DataType,Any}
    # optimizations
    # If using upto 3 contexts, things stay compact
    if isdefined(d, :parent)
        d.key === ctx && return @set d.value = val
        d1 = d.parent
        if isdefined(d1, :parent)
            d1.key === ctx && return @set d.parent.value = val
            d2 = d1.parent
            if isdefined(d2, :parent)
                d2.key === ctx && return @set d.parent.parent.value = val
            end
        end
    end
    Base.ImmutableDict{DataType, Any}(d, ctx, val)
end

function setmetadata(s::Symbolic, ctx::DataType, val)
    if s.metadata isa AbstractDict
        @set s.metadata = assocmeta(s.metadata, ctx, val)
    else
        # fresh Dict
        @set s.metadata = Base.ImmutableDict{DataType, Any}(ctx, val)
    end
end

### Metatheory.jl e-graph rewriting integration

"""
    SymtypeAnalysis

This abstract type is used to identify the EGraph analysis
that keeps track of symtype through an EGraph. This must
be added to every EGraph that is used in SymbolicUtils.
"""
abstract type SymtypeAnalysis <: AbstractAnalysis end
_getsymtype(T::Type{<:Symbolic{X}}) where X = X
_getsymtype(T::Type{X}) where {X} = X
EGraphs.make(an::Type{SymtypeAnalysis}, g::EGraph, n::ENodeLiteral) = symtype(n.value)
EGraphs.make(an::Type{SymtypeAnalysis}, g::EGraph, n::ENodeTerm{T}) where {T} = _getsymtype(T)
EGraphs.join(an::Type{SymtypeAnalysis}, A, B) = Union{A, B}

# TODO JOIN egraph analysis
TermInterface.symtype(ec::EClass) = getdata(ec, SymtypeAnalysis, Any)

function to_symbolic(x)
    Base.depwarn("`to_symbolic(x)` is deprecated, define the interface for your " *
                 "symbolic structure using `istree(x)`, `operation(x)`, `arguments(x)` " *
                 "and `similarterm(::YourType, f, args, symtype)`", :to_symbolic, force=true)

    x
end

###
###  Pretty printing
###
const show_simplified = Ref(false)

isnegative(t::Real) = t < 0
function isnegative(t)
    if istree(t) && operation(t) === (*)
        coeff = first(arguments(t))
        return isnegative(coeff)
    end
    return false
end

# Term{}
setargs(t, args) = Term{symtype(t)}(operation(t), args)
cdrargs(args) = setargs(t, cdr(args))

print_arg(io, x::Union{Complex, Rational}; paren=true) = print(io, "(", x, ")")
isbinop(f) = istree(f) && !istree(operation(f)) && Base.isbinaryoperator(nameof(operation(f)))
function print_arg(io, x; paren=false)
    if paren && isbinop(x)
        print(io, "(", x, ")")
    else
        print(io, x)
    end
end
print_arg(io, s::String; paren=true) = show(io, s)
function print_arg(io, f, x)
    f !== (*) && return print_arg(io, x)
    if Base.isbinaryoperator(nameof(f))
        print_arg(io, x, paren=true)
    else
        print_arg(io, x)
    end
end

function remove_minus(t)
    !istree(t) && return -t
    @assert operation(t) == (*)
    args = arguments(t)
    @assert args[1] < 0
    Any[-args[1], args[2:end]...]
end

function show_add(io, args)
    negs = filter(isnegative, args)
    nnegs = filter(!isnegative, args)
    for (i, t) in enumerate(nnegs)
        i != 1 && print(io, " + ")
        print_arg(io, +,  t)
    end

    for (i, t) in enumerate(negs)
        if i==1 && isempty(nnegs)
            print_arg(io, -, t)
        else
            print(io, " - ")
            show_mul(io, remove_minus(t))
        end
    end
end

function show_pow(io, args)
    base, ex = args

    if base isa Real && base < 0
        print(io, "(")
        print_arg(io, base)
        print(io, ")")
    else
        print_arg(io, base, paren=true)
    end
    print(io, "^")
    print_arg(io, ex, paren=true)
end

function show_mul(io, args)
    length(args) == 1 && return print_arg(io, *, args[1])

    minus = args[1] isa Number && args[1] == -1
    unit = args[1] isa Number && args[1] == 1

    paren_scalar = (args[1] isa Complex && !_iszero(imag(args[1]))) ||
                   args[1] isa Rational ||
                   (args[1] isa Number && !isfinite(args[1]))

    nostar = minus || unit ||
            (!paren_scalar && args[1] isa Number && !(args[2] isa Number))

    for (i, t) in enumerate(args)
        if i != 1
            if i==2 && nostar
            else
                print(io, "*")
            end
        end
        if i == 1 && minus
            print(io, "-")
        elseif i == 1 && unit
        else
            print_arg(io, *, t)
        end
    end
end

function show_ref(io, f, args)
    x = args[1]
    idx = args[2:end]

    istree(x) && print(io, "(")
    print(io, x)
    istree(x) && print(io, ")")
    print(io, "[")
    for i=1:length(idx)
        print_arg(io, idx[i])
        i != length(idx) && print(io, ", ")
    end
    print(io, "]")
end

function show_call(io, f, args)
    fname = istree(f) ? Symbol(repr(f)) : nameof(f)
    binary = Base.isbinaryoperator(fname)
    if binary
        for (i, t) in enumerate(args)
            i != 1 && print(io, " $fname ")
            print_arg(io, t, paren=true)
        end
    else
        if issym(f)
            Base.show_unquoted(io, nameof(f))
        else
            Base.show(io, f)
        end
        print(io, "(")
        for i=1:length(args)
            print(io, args[i])
            i != length(args) && print(io, ", ")
        end
        print(io, ")")
    end
end

function show_term(io::IO, t)
    if get(io, :simplify, show_simplified[])
        return print(IOContext(io, :simplify=>false), simplify(t))
    end

    f = operation(t)
    args = arguments(t)

    if f === (+)
        show_add(io, args)
    elseif f === (*)
        show_mul(io, args)
    elseif f === (^)
        show_pow(io, args)
    elseif f === (getindex)
        show_ref(io, f, args)
    else
        show_call(io, f, args)
    end

    return nothing
end

showraw(io, t) = Base.show(IOContext(io, :simplify=>false), t)
showraw(t) = showraw(stdout, t)

function Base.show(io::IO, v::BasicSymbolic)
    if TermInterface.issym(v)
        Base.show_unquoted(io, v.name)
    else
        show_term(io, v)
    end
end

###
### Symbolic function / type inference
###

"""
    promote_symtype(f, Ts...)

The result of applying `f` to arguments of [`symtype`](#symtype) `Ts...`

```julia
julia> promote_symtype(+, Real, Real)
Real

julia> promote_symtype(+, Complex, Real)
Number

julia> @syms f(x)::Complex
(f(::Number)::Complex,)

julia> promote_symtype(f, Number)
Complex
```

When constructing [`Term`](#Term)s without an explicit symtype,
`promote_symtype` is used to figure out the symtype of the Term.
"""
promote_symtype(f, Ts...) = Any

#---------------------------
#---------------------------
#### Function-like variables
#---------------------------

struct FnType{X<:Tuple,Y} end

(f::Symbolic{<:FnType})(args...) = Term{promote_symtype(f, symtype.(args)...)}(f, Any[args...])

function (f::Symbolic)(args...)
    error("Sym $f is not callable. " *
          "Use @syms $f(var1, var2,...) to create it as a callable.")
end

"""
    promote_symtype(f::FnType{X,Y}, arg_symtypes...)

The output symtype of applying variable `f` to arugments of symtype `arg_symtypes...`.
if the arguments are of the wrong type then this function will error.
"""
function promote_symtype(f::BasicSymbolic{<:FnType{X,Y}}, args...) where {X, Y}
    if X === Tuple
        return Y
    end

    # This is to handle `Tuple{T} where T`, so we cannot reliably query the type
    # parameters of the `Tuple` in `FnType`.
    t = Tuple{args...}
    if !(t <: X)
        error("$t is not a subtype of $X.")
    end
    return Y
end

function Base.show(io::IO, f::Symbolic{<:FnType{X,Y}}) where {X,Y}
    print(io, nameof(f))
    # Use `Base.unwrap_unionall` to handle `Tuple{T} where T`. This is not the
    # best printing, but it's better than erroring.
    argrepr = join(map(t->"::"*string(t), Base.unwrap_unionall(X).parameters), ", ")
    print(io, "(", argrepr, ")")
    print(io, "::", Y)
end

@inline isassociative(op) = op === (+) || op === (*)

function _promote_symtype(f, args)
    if TermInterface.issym(f)
        promote_symtype(f, map(symtype, args)...)
    else
        if length(args) == 0
            promote_symtype(f)
        elseif length(args) == 1
            promote_symtype(f, symtype(args[1]))
        elseif length(args) == 2
            promote_symtype(f, symtype(args[1]), symtype(args[2]))
        elseif isassociative(f)
            mapfoldl(symtype, (x,y) -> promote_symtype(f, x, y), args)
        else
            promote_symtype(f, map(symtype, args)...)
        end
    end
end

###
### Macro
###

"""
    @syms <lhs_expr>[::T1] <lhs_expr>[::T2]...

For instance:

    @syms foo::Real bar baz(x, y::Real)::Complex

Create one or more variables. `<lhs_expr>` can be just a symbol in which case
it will be the name of the variable, or a function call in which case a function-like
variable which has the same name as the function being called. The Sym type, or
in the case of a function-like Sym, the output type of calling the function
can be set using the `::T` syntax.

# Examples:

- `@syms foo bar::Real baz::Int` will create
variable `foo` of symtype `Number` (the default), `bar` of symtype `Real`
and `baz` of symtype `Int`
- `@syms f(x) g(y::Real, x)::Int h(a::Int, f(b))` creates 1-arg `f` 2-arg `g`
and 2 arg `h`. The second argument to `h` must be a one argument function-like
variable. So, `h(1, g)` will fail and `h(1, f)` will work.
"""
macro syms(xs...)
    defs = map(xs) do x
        n, t = _name_type(x)
        T = esc(t)
        nt = _name_type(x)
        n, t = nt.name, nt.type
        :($(esc(n)) = Sym{$T}($(Expr(:quote, n))))
    end
    Expr(:block, defs...,
         :(tuple($(map(x->esc(_name_type(x).name), xs)...))))
end

function syms_syntax_error()
    error("Incorrect @syms syntax. Try `@syms x::Real y::Complex g(a) f(::Real)::Real` for instance.")
end

function _name_type(x)
    if x isa Symbol
        return (name=x, type=Number)
    elseif x isa Expr && x.head === :(::)
        if length(x.args) == 1
            return (name=nothing, type=x.args[1])
        end
        lhs, rhs = x.args[1:2]
        if lhs isa Expr && lhs.head === :call
            # e.g. f(::Real)::Unreal
            type = map(x->_name_type(x).type, lhs.args[2:end])
            return (name=lhs.args[1], type=:($FnType{Tuple{$(type...)}, $rhs}))
        else
            return (name=lhs, type=rhs)
        end
    elseif x isa Expr && x.head === :ref
        ntype = _name_type(x.args[1]) # a::Number
        N = length(x.args)-1
        return (name=ntype.name,
                type=:(Array{$(ntype.type), $N}),
                array_metadata=:(Base.Slice.(($(x.args[2:end]...),))))
    elseif x isa Expr && x.head === :call
        return _name_type(:($x::Number))
    else
        syms_syntax_error()
    end
end

###
### Arithmetic
###
const SN = Symbolic{<:Number}
# TODO Reviewme this is necessary for Metatheory.jl egraph rewriting
# integration. Constructors of `Add, Mul, Pow...` from Base (+, *, ^, ...)
# Should now accepts EClasses as arguments.
const SN_EC = Union{SN, EClass}

_merge(f::F, d, others...; filter=x->false) where F = _merge!(f, Dict{Any,Any}(d), others...; filter=filter)

function _merge!(f::F, d, others...; filter=x->false) where F
    acc = d
    for other in others
        for (k, v) in other
            v = f(v)
            ak = get(acc, k, nothing)
            if ak !== nothing
                v = ak + v
            end
            if filter(v)
                delete!(acc, k)
            else
                acc[k] = v
            end
        end
    end
    acc
end

function mapvalues(f, d1::AbstractDict)
    d = copy(d1)
    for (k, v) in d
        d[k] = f(k, v)
    end
    d
end

add_t(a::Number,b::Number) = promote_symtype(+, symtype(a), symtype(b))
add_t(a,b) = promote_symtype(+, symtype(a), symtype(b))
sub_t(a,b) = promote_symtype(-, symtype(a), symtype(b))
sub_t(a) = promote_symtype(-, symtype(a))

import Base: (+), (-), (*), (//), (/), (\), (^)
function +(a::SN_EC, b::SN_EC)
    if isadd(a) && isadd(b)
        return Add(add_t(a,b),
                   a.coeff + b.coeff,
                   _merge(+, a.dict, b.dict, filter=_iszero))
    elseif isadd(a)
        coeff, dict = makeadd(1, 0, b)
        return Add(add_t(a,b), a.coeff + coeff, _merge(+, a.dict, dict, filter=_iszero))
    elseif isadd(b)
        return b + a
    end
    coeff, dict = makeadd(1, 0, a, b)
    Add(add_t(a,b), coeff, dict)
end

function +(a::Number, b::SN_EC)
    iszero(a) && return b
    if isadd(b)
        Add(add_t(a,b), a + b.coeff, b.dict)
    else
        Add(add_t(a,b), makeadd(1, a, b)...)
    end
end

+(a::SN_EC, b::Number) = b + a

+(a::SN_EC) = a

function -(a::SN_EC)
    isadd(a) ? Add(sub_t(a), -a.coeff, mapvalues((_,v) -> -v, a.dict)) :
    Add(sub_t(a), makeadd(-1, 0, a)...)
end

function -(a::SN_EC, b::SN_EC)
    isadd(a) && isadd(b) ? Add(sub_t(a,b),
                               a.coeff - b.coeff,
                               _merge(-, a.dict,
                                      b.dict,
                                      filter=_iszero)) : a + (-b)
end

-(a::Number, b::SN_EC) = a + (-b)
-(a::SN_EC, b::Number) = a + (-b)

mul_t(a,b) = promote_symtype(*, symtype(a), symtype(b))
mul_t(a) = promote_symtype(*, symtype(a))

*(a::SN_EC) = a

function *(a::SN_EC, b::SN_EC)
    # Always make sure Div wraps Mul
    if isdiv(a) && isdiv(b)
        Div(a.num * b.num, a.den * b.den)
    elseif isdiv(a)
        Div(a.num * b, a.den)
    elseif isdiv(b)
        Div(a * b.num, b.den)
    elseif ismul(a) && ismul(b)
        Mul(mul_t(a, b),
            a.coeff * b.coeff,
            _merge(+, a.dict, b.dict, filter=_iszero))
    elseif ismul(a) && ispow(b)
        if b.exp isa Number
            Mul(mul_t(a, b),
                a.coeff, _merge(+, a.dict, Base.ImmutableDict(b.base=>b.exp), filter=_iszero))
        else
            Mul(mul_t(a, b),
                a.coeff, _merge(+, a.dict, Base.ImmutableDict(b=>1), filter=_iszero))
        end
    elseif ispow(a) && ismul(b)
        b * a
    else
        Mul(mul_t(a,b), makemul(1, a, b)...)
    end
end

function *(a::Number, b::SN_EC)
    if iszero(a)
        a
    elseif isone(a)
        b
    elseif isdiv(b)
        Div(a*b.num, b.den)
    elseif isone(-a) && isadd(b)
        # -1(a+b) -> -a - b
        T = promote_symtype(+, typeof(a), symtype(b))
        Add(T, b.coeff * a, Dict{Any,Any}(k=>v*a for (k, v) in b.dict))
    else
        Mul(mul_t(a, b), makemul(a, b)...)
    end
end

###
### Div
###

/(a::Union{SN_EC,Number}, b::SN_EC) = Div(a, b)

*(a::SN_EC, b::Number) = b * a

\(a::SN_EC, b::Union{Number, SN_EC}) = b / a

\(a::Number, b::SN_EC) = b / a

/(a::SN_EC, b::Number) = (b isa Integer ? 1//b : inv(b)) * a

//(a::Union{SN_EC, Number}, b::SN_EC) = a / b

//(a::SN_EC, b::T) where {T <: Number} = (one(T) // b) * a

###
### Pow
###

function ^(a::SN_EC, b)
    if b isa Number && iszero(b)
        # fast path
        1
    elseif b isa Number && b < 0
        Div(1, a ^ (-b))
    elseif ismul(a) && b isa Number
        coeff = unstable_pow(a.coeff, b)
        Mul(promote_symtype(^, symtype(a), symtype(b)),
            coeff, mapvalues((k, v) -> b*v, a.dict))
    else
        Pow(a, b)
    end
end

^(a::Number, b::SN_EC) = Pow(a, b)


text:
using GPMLj, Test, Random, PyCall

Random.seed!(123)

@testset "Likelihoods" begin
    @testset "Bernoulli" begin
        likelihood = gpflow.likelihoods.Bernoulli()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Bernoulli
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "Beta" begin
        likelihood = gpflow.likelihoods.Beta()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Beta
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "Exponential" begin
        likelihood = gpflow.likelihoods.Exponential()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Exponential
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "Gamma" begin
        likelihood = gpflow.likelihoods.Gamma()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Gamma
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "Gaussian" begin
        likelihood = gpflow.likelihoods.Gaussian()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Gaussian
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "GaussianMC" begin
        @test_throws "GaussianMC not implemented" gpflow.likelihoods.GaussianMC()
    end

    @testset "Likelihood" begin
        likelihood = gpflow.likelihoods.Likelihood()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Likelihood
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "MonteCarloLikelihood" begin
        likelihood = gpflow.likelihoods.MonteCarloLikelihood()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.MonteCarloLikelihood
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "MultiClass" begin
        likelihood = gpflow.likelihoods.MultiClass(3)
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.MultiClass
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "Ordinal" begin
        likelihood = gpflow.likelihoods.Ordinal([1,2])
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Ordinal
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "Poisson" begin
        likelihood = gpflow.likelihoods.Poisson()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.Poisson
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "RobustMax" begin
        likelihood = gpflow.likelihoods.RobustMax(3)
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.RobustMax
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "SoftMax" begin
        likelihood = gpflow.likelihoods.SoftMax(3)
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.SoftMax
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    @testset "StudentT" begin
        likelihood = gpflow.likelihoods.StudentT()
        @test typeof(likelihood)<:gpflow.AbstractLikelihood
        @test typeof(likelihood)<:gpflow.likelihoods.StudentT
        @test typeof(likelihood.o)<:PyObject
        temp = likelihood.o
        @test instantiate!(likelihood) == temp
    end

    # TODO
    # @testset "SwitchedLikelihood" begin
    #     likelihood = gpflow.likelihoods.SwitchedLikelihood()
    #     @test typeof(likelihood)<:gpflow.AbstractLikelihood
    #     @test typeof(likelihood)<:gpflow.likelihoods.SwitchedLikelihood
    #     @test typeof(likelihood.o)<:PyObject
    #     temp = likelihood.o
    #     @test instantiate!(likelihood) == temp
    # end

end #module


text:
function solve()
    x = readline()
    s = split(x, '.')
    s[1]
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# backward-compatibility to ColorTypes < v0.9 or Colors < v0.11

struct ArrayLF{T,N} <: AbstractArray{T,N}
    A::Array{T,N}
end
Base.IndexStyle(::Type{A}) where {A<:ArrayLF} = IndexLinear()
Base.size(A::ArrayLF) = size(A.A)
Base.getindex(A::ArrayLF, i::Int) = A.A[i]
Base.setindex!(A::ArrayLF, val, i::Int) = A.A[i] = val

struct ArrayLS{T,N} <: AbstractArray{T,N}
    A::Array{T,N}
end
Base.IndexStyle(::Type{A}) where {A<:ArrayLS} = IndexCartesian()
Base.size(A::ArrayLS) = size(A.A)
Base.getindex(A::ArrayLS{T,N}, i::Vararg{Int,N}) where {T,N} = A.A[i...]
Base.setindex!(A::ArrayLS{T,N}, val, i::Vararg{Int,N}) where {T,N} = A.A[i...] = val

@testset "channel_view" begin

@testset "grayscale" begin
    a = rand(2,3)
    @test channel_view(a) === a

    a0 = [Gray(N0f8(0.2)), Gray(N0f8(0.4))]
    for (a, LI) in ((copy(a0), IndexLinear()),
                    (ArrayLF(copy(a0)), IndexLinear()),
                    (ArrayLS(copy(a0)), IndexCartesian()))
        v = @inferred(channel_view(a))
        @test @inferred(IndexStyle(v)) == LI
        @test ndims(v) == 1
        @test size(v) == (2,)
        @test eltype(v) == N0f8
        @test @inferred(color_view(Gray, v)) === a
        @test parent(parent(v)) === a
        @test v[1] == N0f8(0.2)
        @test v[2] == N0f8(0.4)
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[3]
        v[1] = 0.8
        @test a[1] === Gray(N0f8(0.8))
        @test_throws BoundsError (v[0] = 0.6)
        @test_throws BoundsError (v[3] = 0.6)
        c = similar(v)
        @test eltype(c) == N0f8 && ndims(c) == 1
        @test length(c) == 2
        c = similar(v, 3)
        @test eltype(c) == N0f8 && ndims(c) == 1
        @test length(c) == 3
        c = similar(v, Float32)
        @test eltype(c) == Float32 && ndims(c) == 1
        @test length(c) == 2
        c = similar(v, Float16, (5,5))
        @test eltype(c) == Float16 && ndims(c) == 2
        @test size(c) == (5,5)
    end
end

@testset "RGB, HSV, etc" begin
    for T in (RGB, BGR, XRGB, RGBX, HSV, Lab, XYZ)
        a0 = [T(0.1,0.2,0.3), T(0.4, 0.5, 0.6)]
        for a in (copy(a0),
                  ArrayLS(copy(a0)))
            v = @inferred(channel_view(a))
            @test ndims(v) == 2
            @test size(v) == (3,2)
            @test eltype(v) == Float64
            if T in (RGB, HSV, Lab, XYZ)
                @test @inferred(color_view(T, v)) == a && color_view(T, v) isa typeof(a)
            else
                @test @inferred(color_view(T, v)) == a
            end
            @test v[1] == v[1,1] == 0.1
            @test v[2] == v[2,1] == 0.2
            @test v[3] == v[3,1] == 0.3
            @test v[4] == v[1,2] == 0.4
            @test v[5] == v[2,2] == 0.5
            @test v[6] == v[3,2] == 0.6
            @test_throws BoundsError v[0,1]
            @test_throws BoundsError v[4,1]
            @test_throws BoundsError v[2,0]
            @test_throws BoundsError v[2,3]
            v[2] = 0.8
            @test a[1] == T(0.1,0.8,0.3)
            v[2,1] = 0.7
            @test a[1] == T(0.1,0.7,0.3)
            @test_throws BoundsError (v[0,1] = 0.7)
            @test_throws BoundsError (v[4,1] = 0.7)
            @test_throws BoundsError (v[2,0] = 0.7)
            @test_throws BoundsError (v[2,3] = 0.7)
            c = similar(v)
            @test size(c) == (3,2) && eltype(c) == Float64
            c = similar(v, (3,4))
            @test size(c) == (3,4) && eltype(c) == Float64
            c = similar(v, Float32)
            @test size(c) == (3,2) && eltype(c) == Float32
            c = similar(v, Float16, (3,5,5))
            @test size(c) == (3,5,5) && eltype(c) == Float16
        end
    end
    #= FIXME
    a = reshape([RGB(1,0,0)])  # 0-dimensional
    v = @inferred(channel_view(a))
    @test axes(v) === (Base.OneTo(3),)
    v = @inferred(channel_view(a))
    @test axes(v) === (Base.OneTo(3),)
    =#
end

@testset "Gray+Alpha" begin
    for T in (AGray, GrayA)
        a = [T(0.1f0,0.2f0), T(0.3f0,0.4f0), T(0.5f0,0.6f0)]
        v = @inferred(channel_view(a))
        @test @inferred(color_view(T, v)) == a
        @test ndims(v) == 2
        @test size(v) == (2,3)
        @test eltype(v) == Float32
        @test v[1] == v[1,1] == 0.1f0
        @test v[2] == v[2,1] == 0.2f0
        @test v[3] == v[1,2] == 0.3f0
        @test v[4] == v[2,2] == 0.4f0
        @test v[5] == v[1,3] == 0.5f0
        @test v[6] == v[2,3] == 0.6f0
        @test_throws BoundsError v[0,1]
        @test_throws BoundsError v[3,1]
        @test_throws BoundsError v[2,0]
        @test_throws BoundsError v[2,4]
        v[2] = 0.8
        @test a[1] == T(0.1f0,0.8f0)
        v[2,1] = 0.7
        @test a[1] == T(0.1f0,0.7f0)
        @test_throws BoundsError (v[0,1] = 0.7)
        @test_throws BoundsError (v[3,1] = 0.7)
        @test_throws BoundsError (v[2,0] = 0.7)
        @test_throws BoundsError (v[2,4] = 0.7)
        c = similar(v)
        @test eltype(c) == Float32
        @test size(c) == (2,3)
        c = similar(v, (2,4))
        @test eltype(c) == Float32
        @test size(c) == (2,4)
        c = similar(v, Float64)
        @test eltype(c) == Float64
        @test size(c) == (2,3)
        c = similar(v, Float16, (2,5,5))
        @test eltype(c) == Float16
        @test size(c) == (2,5,5)
    end
end

@testset "Alpha+RGB, HSV, etc" begin
    for T in (ARGB,
              ABGR,
              AHSV,
              ALab,
              AXYZ,
              RGBA,
              BGRA,
              LabA,
              XYZA)
        a = [T(0.1,0.2,0.3,0.4), T(0.5,0.6,0.7,0.8)]
        v = @inferred(channel_view(a))
        @test ndims(v) == 2
        @test size(v) == (4,2)
        @test eltype(v) == Float64
        @test v[1] == v[1,1] == 0.1
        @test v[2] == v[2,1] == 0.2
        @test v[3] == v[3,1] == 0.3
        @test v[4] == v[4,1] == 0.4
        @test v[5] == v[1,2] == 0.5
        @test v[6] == v[2,2] == 0.6
        @test v[7] == v[3,2] == 0.7
        @test v[8] == v[4,2] == 0.8
        @test_throws BoundsError v[0,1]
        @test_throws BoundsError v[5,1]
        @test_throws BoundsError v[2,0]
        @test_throws BoundsError v[2,3]
        v[2] = 0.9
        @test a[1] == T(0.1,0.9,0.3,0.4)
        v[2,1] = 0.7
        @test a[1] == T(0.1,0.7,0.3,0.4)
        @test_throws BoundsError (v[0,1] = 0.7)
        @test_throws BoundsError (v[5,1] = 0.7)
        @test_throws BoundsError (v[2,0] = 0.7)
        @test_throws BoundsError (v[2,3] = 0.7)
        c = similar(v)
        @test eltype(c) == Float64
        @test size(c) == (4,2)
        c = similar(v, (4,4))
        @test eltype(c) == Float64
        @test size(c) == (4,4)
        c = similar(v, Float32)
        @test eltype(c) == Float32
        @test size(c) == (4,2)
        c = similar(v, Float16, (4,5,5))
        @test eltype(c) == Float16
        @test size(c) == (4,5,5)
    end

    #= FIXME using AxisIndices.OffsetArray
    @testset "Non-1 indices" begin
        a = OffsetArray(rand(RGB{N0f8}, 3, 5), -1:1, -2:2)
        v = @inferred(channel_view(a))
        @test @inferred(axes(v)) == (IdentityUnitRange(1:3), IdentityUnitRange(-1:1), IdentityUnitRange(-2:2))
        @test @inferred(v[1,0,0]) === a[0,0].r
        a = OffsetArray(rand(Gray{Float32}, 3, 5), -1:1, -2:2)
        v = @inferred(channel_view(a))
        @test @inferred(axes(v)) == (IdentityUnitRange(-1:1), IdentityUnitRange(-2:2))
        @test @inferred(v[0,0]) === gray(a[0,0])
        @test @inferred(v[5]) === gray(a[5])
        v[5] = -1
        @test v[5] === -1.0f0
    end
    =#
end

end

@testset "color_view" begin

@testset "grayscale" begin
    a0 = [N0f8(0.2), N0f8(0.4)]
    for (a, LI) in ((copy(a0), IndexLinear()),
                    (ArrayLF(copy(a0)), IndexLinear()),
                    (ArrayLS(copy(a0)), IndexCartesian()))
        @test_throws MethodError color_view(a)
        v = @inferred(color_view(Gray, a))
        @test @inferred(IndexStyle(v)) == LI
        @test ndims(v) == 1
        @test size(v) == (2,)
        @test eltype(v) == Gray{N0f8}
        @test @inferred(channel_view(v)) === a
        @test parent(parent(v)) === a
        @test v[1] == Gray(N0f8(0.2))
        @test v[2] == Gray(N0f8(0.4))
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[3]
        v[1] = 0.8
        @test a[1] === N0f8(0.8)
        @test_throws BoundsError (v[0] = 0.6)
        @test_throws BoundsError (v[3] = 0.6)
        c = similar(v)
        @test eltype(c) == Gray{N0f8} && ndims(c) == 1
        @test length(c) == 2
        c = similar(v, 3)
        @test eltype(c) == Gray{N0f8} && ndims(c) == 1
        @test length(c) == 3
        c = similar(v, Gray{Float32})
        @test eltype(c) == Gray{Float32}
        @test length(c) == 2
        c = similar(v, Gray{Float16}, (5,5))
        @test eltype(c) == Gray{Float16}
        @test size(c) == (5,5)
        c = similar(v, Float32)
        @test isa(c, Array{Float32, 1})
    end
    # two dimensional images and linear indexing
    a0 = N0f8[0.2 0.4; 0.6 0.8]
    for (a, LI) in ((copy(a0), IndexLinear()),
                    (ArrayLF(copy(a0)), IndexLinear()),
                    (ArrayLS(copy(a0)), IndexCartesian()))
        @test_throws MethodError color_view(a)
        v = @inferred(color_view(Gray, a))
        @test @inferred(IndexStyle(v)) == LI
        @test ndims(v) == 2
        @test size(v) == (2,2)
        @test eltype(v) == Gray{N0f8}
        @test @inferred(channel_view(v)) === a
        @test parent(parent(v)) === a
        @test v[1] == Gray(N0f8(0.2))
        @test v[2] == Gray(N0f8(0.6))
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[5]
        v[1] = 0.9
        @test a[1] === N0f8(0.9)
        @test_throws BoundsError (v[0] = 0.6)
        @test_throws BoundsError (v[5] = 0.6)
    end
end

@testset "RGB, HSV, etc" begin
    for T in (RGB, BGR, XRGB, RGBX, HSV, Lab, XYZ)
        a0 = [0.1 0.2 0.3; 0.4 0.5 0.6]'
        for a in (copy(a0),
                  ArrayLS(copy(a0)))
            v = @inferred(color_view(T,a))
            @test @inferred(channel_view(v)) === a
            @test ndims(v) == 1
            @test size(v) == (2,)
            @test eltype(v) == T{Float64}
            @test v[1] == T(0.1,0.2,0.3)
            @test v[2] == T(0.4,0.5,0.6)
            @test_throws BoundsError v[0]
            @test_throws BoundsError v[3]
            v[2] = T(0.8, 0.7, 0.6)
            @test a == [0.1 0.2 0.3; 0.8 0.7 0.6]'
            @test_throws BoundsError (v[0] = T(0.8, 0.7, 0.6))
            @test_throws BoundsError (v[3] = T(0.8, 0.7, 0.6))
            c = similar(v)
            @test size(c) == (2,)
            c = similar(v, 4)
            @test size(c) == (4,)
            c = similar(v, T{Float32})
            @test size(c) == (2,)
            c = similar(v, T)
            @test size(c) == (2,)
            c = similar(v, T{Float16}, (5,5))
            @test size(c) == (5,5)
            c = similar(v, RGB24)
            @test eltype(c) == RGB24
            @test size(c) == size(v)
        end
    end
end

@testset "Gray+Alpha" begin
    for T in (AGray, GrayA)
        a = [0.1f0 0.2f0; 0.3f0 0.4f0; 0.5f0 0.6f0]'
        v = @inferred(color_view(T, a))
        @test ndims(v) == 1
        @test size(v) == (3,)
        @test eltype(v) == T{Float32}
        @test @inferred(channel_view(v)) === a
        @test v[1] == T(0.1f0, 0.2f0)
        @test v[2] == T(0.3f0, 0.4f0)
        @test v[3] == T(0.5f0, 0.6f0)
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[4]
        v[2] = T(0.8, 0.7)
        @test a[1,2] === 0.8f0
        @test a[2,2] === 0.7f0
        @test_throws BoundsError (v[0] = T(0.8,0.7))
        @test_throws BoundsError (v[4] = T(0.8,0.7))
        c = similar(v)
        @test eltype(c) == T{Float32}
        @test size(c) == (3,)
        c = similar(v, (4,))
        @test eltype(c) == T{Float32}
        @test size(c) == (4,)
        c = similar(v, T{Float64})
        @test eltype(c) == T{Float64}
        @test size(c) == (3,)
        c = similar(v, T{Float16}, (5,5))
        @test eltype(c) == T{Float16}
        @test size(c) == (5,5)
    end
end

@testset "Alpha+RGB, HSV, etc" begin
    a = rand(ARGB{N0f8}, 5, 5)
    vc = @inferred(channel_view(a))
    @test eltype(@inferred(color_view(ARGB, vc))) == ARGB{N0f8}
    cvc = @inferred(color_view(RGBA, vc))
    @test all(cvc .== a)

    for T in (ARGB,
              ABGR,
              AHSV,
              ALab,
              AXYZ,
              RGBA,
              BGRA,
              HSVA,
              LabA,
              XYZA)
        a = [0.1 0.2 0.3 0.4; 0.5 0.6 0.7 0.8]'
        v = @inferred(color_view(T, a))
        @test eltype(v) == T{Float64}
        @test @inferred(channel_view(v)) === a
        @test ndims(v) == 1
        @test size(v) == (2,)
        @test eltype(v) == T{Float64}
        @test v[1] == T(0.1,0.2,0.3,0.4)
        @test v[2] == T(0.5,0.6,0.7,0.8)
        @test_throws BoundsError v[0]
        @test_throws BoundsError v[3]
        v[2] = T(0.9,0.8,0.7,0.6)
        @test a[1,2] == 0.9
        @test a[2,2] == 0.8
        @test a[3,2] == 0.7
        @test a[4,2] == 0.6
        @test_throws BoundsError (v[0] = T(0.9,0.8,0.7,0.6))
        @test_throws BoundsError (v[3] = T(0.9,0.8,0.7,0.6))
        c = similar(v)
        @test eltype(c) == T{Float64}
        @test size(c) == (2,)
        c = similar(v, 4)
        @test eltype(c) == T{Float64}
        @test size(c) == (4,)
        c = similar(v, T{Float32})
        @test eltype(c) == T{Float32}
        @test size(c) == (2,)
        c = similar(v, T{Float16}, (5,5))
        @test eltype(c) == T{Float16}
        @test size(c) == (5,5)
    end

    #= FIXME using AxisIndices.OffsetArray
    @testset "Non-1 indices" begin
        a = OffsetArray(rand(3, 3, 5), 1:3, -1:1, -2:2)
        v = @inferred(color_view(RGB, a))
        @test @inferred(axes(v)) == (IdentityUnitRange(-1:1), IdentityUnitRange(-2:2))
        @test @inferred(v[0,0]) === RGB(a[1,0,0], a[2,0,0], a[3,0,0])
        a = OffsetArray(rand(3, 3, 5), 0:2, -1:1, -2:2)
        @test_throws DimensionMismatch color_view(RGB, a)
    end
    =#
end

@testset "grayscale" begin
    A = NamedAxisArray{(:y, :x)}(rand(Gray{N0f8}, 4, 5));
    # FIXME
    #@test summary(A) == "2-dimensional AxisArray{Gray{N0f8},2,...} with axes:\n    :y, Base.OneTo(4)\n    :x, Base.OneTo(5)\nAnd data, a 4×5 Array{Gray{N0f8},2} with eltype Gray{Normed{UInt8,8}}"
    cv = channel_view(A);
    @test named_axes(cv) == (y = 1:4, x = 1:5)
    @test spatial_order(cv) == (:y, :x)
    @test_throws ArgumentError channeldim(cv)
end

@testset "color" begin
    A = NamedAxisArray{(:y, :x)}(rand(RGB{N0f8}, 4, 5));
    cv = channel_view(A);
    @test @inferred(named_axes(cv)) == (color = 1:3, y = 1:4, x = 1:5)
    @test @inferred(spatial_order(cv)) == (:y, :x)
    @test @inferred(channeldim(cv)) == 1
    p = permuteddimsview(cv, (2,3,1))
    @test @inferred(named_axes(p)) == (y = 1:4, x = 1:5, color = 1:3)
    @test channeldim(p) == 3
end

@testset "nested" begin
    A = NamedAxisArray(rand(RGB{N0f8}, 4, 5), y = range(1, step=2, length=4), x=1:5);
    P = permuteddimsview(A, (2, 1));
    @test @inferred(pixel_spacing(P)) == (1, 2)
    M = mappedarray(identity, A)
    @test @inferred(pixel_spacing(M)) == (2, 1)
    s = u"s" # global const
    μm = u"μm" # global const
    tax = range(0.0s, step=0.1s, length=11)
    A = NamedAxisArray(rand(N0f16, 4, 5, 11),
                  y = range(1μm, step=2μm, length=4),
                  x = range(1μm, step=1μm, length=5),
                  time = tax)
    P = permuteddimsview(A, (3, 1, 2));
    M = mappedarray(identity, A);
    @test @inferred(pixel_spacing(P)) == @inferred(pixel_spacing(M)) == (2μm, 1μm)
    @test @inferred(time_keys(P)) == @inferred(time_keys(M)) == tax
    @test has_timedim(P)
    @test spatialdims(P) == (2, 3)
    @test spatialdims(M) == (1, 2)
    @test spatial_order(P) == spatial_order(M) == (:y, :x)
    @test @inferred(spatial_size(P)) == @inferred(spatial_size(M)) == (4, 5)
    @test_throws ErrorException assert_timedim_last(P)
    assert_timedim_last(M)
    A = NamedAxisArray(rand(N0f16, 11, 5, 4),
                       (time = tax,
                        x = range(1μm, step=1μm, length=5),
                        y = range(1μm, step=2μm, length=4)))
    P = permuteddimsview(A, (3, 2, 1))
    M = mappedarray(identity, A)
    @test @inferred(spatial_keys(P)) == ((1:2:7)μm, (1:5)μm)
    @test @inferred(pixel_spacing(P)) == (2μm, 1μm)
    @test @inferred(pixel_spacing(M)) == (1μm, 2μm)
    @test @inferred(time_keys(P)) == @inferred(time_keys(M)) == tax
    @test has_timedim(P)
    @test spatialdims(P) == (1, 2)
    @test spatialdims(M) == (2, 3)
    @test spatial_order(P) == (:y, :x)
    @test spatial_order(M) == (:x, :y)
    @test @inferred(spatial_size(P)) == (4, 5)
    @test @inferred(spatial_size(M)) == (5, 4)
    assert_timedim_last(P)
    @test_throws ErrorException assert_timedim_last(M)
end

end

nothing