text:
using TensorKit
using statmech_tm_solver
using TensorOperations
using LinearAlgebra

io = open("result_triangular_ising_fulldiag.txt", "w")
close(io)
T = mpo_triangular_AF_ising()

T5 = ncon([T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 1]])
T5 = permute(T5, (1, 3, 5, 7, 9), (2, 4, 6, 8, 10))

T6 = ncon([T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 1]])
T6 = permute(T6, (1, 3, 5, 7, 9, 11), (2, 4, 6, 8, 10, 12))

T7 = ncon([T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 1]])
T7 = permute(T7, (1, 3, 5, 7, 9, 11, 13), (2, 4, 6, 8, 10, 12, 14))

T8 = ncon([T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 1]])
T8 = permute(T8, (1, 3, 5, 7, 9, 11, 13, 15), (2, 4, 6, 8, 10, 12, 14, 16))

T9 = ncon([T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 1]])
T9 = permute(T9, (1, 3, 5, 7, 9, 11, 13, 15, 17), (2, 4, 6, 8, 10, 12, 14, 16, 18))

T10 = ncon([T, T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 10], [10, -19, -20, 1]])
T10 = permute(T10, (1, 3, 5, 7, 9, 11, 13, 15, 17, 19), (2, 4, 6, 8, 10, 12, 14, 16, 18, 20))

T11 = ncon([T, T, T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 10], [10, -19, -20, 11], [11, -21, -22, 1]])
T11 = permute(T11, (1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21), (2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22))

T12 = ncon([T, T, T, T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 10], [10, -19, -20, 11], [11, -21, -22, 12], [12, -23, -24, 1]])
T12 = permute(T12, (1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23), (2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24))

for (n, Tn) in zip([5, 6, 7, 8, 9, 10, 11, 12], [T5, T6, T7, T8, T9, T10, T11, T12])
    w, _ = eig(Tn)
    w = diag(w.data)
    w_reals = real.(w)
    w_imags = imag.(w)

    io = open("result_triangular_ising_fulldiag.txt", "a+")
    for (w_real, w_imag) in zip(w_reals, w_imags)
        println(io, "$(n)   $(w_real)   $(w_imag) ")
    end
    close(io)
end

w, v = eig(T12);
space(v)
v0 = v.data[:, 4096]

w.data[4096, 4096]
λs=diag(w.data)

findmax(imag.(λs))

v11 = v.data[:, 11]

findmax(norm.(v0))

findall(>(0.03), norm.(v11))

v11[220]


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration
using Test

const _order = 28
const _abstol = 1.0E-20
const tT = Taylor1(_order)

@testset "Tests: dot{x}=x^2, x(0) = 1" begin
    eqs_mov(x, p, t) = x^2
    t0 = 0.0
    x0 = 1.0
    x0T = Taylor1(x0, _order)
    tT[1] = t0
    TaylorIntegration.jetcoeffs!(eqs_mov, tT, x0T, nothing)
    @test x0T.coeffs[end] == 1.0
    δt = _abstol^inv(_order-1)
    @test TaylorIntegration.stepsize(x0T, _abstol) == δt

    tv, xv = taylorinteg(eqs_mov, 1, 0.0, 1.0, _order, _abstol)
    @test length(tv) == 501
    @test length(xv) == 501
    @test xv[1] == x0
    @test tv[end] < 1.0

    tv, xv = taylorinteg(eqs_mov, x0, 0.0, 1.0, _order, _abstol, nothing)
    @test length(tv) == 501
    @test length(xv) == 501
    @test xv[1] == x0
    @test tv[end] < 1.0

    trange = 0.0:1/8:1.0
    xv = taylorinteg(eqs_mov, 1, trange, _order, _abstol)
    @test length(xv) == length(trange)
    @test typeof(xv) == Array{typeof(x0),1}
    @test xv[1] == x0
    @test isnan(xv[end])
    @test abs(xv[5] - 2.0) ≤ eps(2.0)
    tvr, xvr = taylorinteg(eqs_mov, x0, trange[1], trange[end-1], _order, _abstol)
    @test tvr[1] == trange[1]
    @test tvr[end] == trange[end-1]
    @test xvr[1] == xv[1]
    @test xvr[end] == xv[end-1]

    trange = 0.0:1/8:1.0
    xv = taylorinteg(eqs_mov, x0, trange, _order, _abstol, nothing)
    @test length(xv) == length(trange)
    @test typeof(xv) == Array{typeof(x0),1}
    @test xv[1] == x0
    @test isnan(xv[end])
    @test abs(xv[5] - 2.0) ≤ eps(2.0)
    tvr, xvr = taylorinteg(eqs_mov, x0, trange[1], trange[end-1], _order, _abstol, nothing)
    @test tvr[1] == trange[1]
    @test tvr[end] == trange[end-1]
    @test xvr[1] == xv[1]
    @test xvr[end] == xv[end-1]

    tarray = vec(trange)
    xv2 = taylorinteg(eqs_mov, x0, tarray, _order, _abstol)
    @test xv[1:end-1] == xv2[1:end-1]
    @test xv2[1:end-1] ≈ xv[1:end-1] atol=eps() rtol=0.0
    @test length(xv2) == length(tarray)
    @test typeof(xv2) == Array{typeof(x0),1}
    @test xv2[1] == x0
    @test isnan(xv2[end])
    @test abs(xv2[5] - 2.0) ≤ eps(2.0)
end

@testset "Tests: dot{x}=x^2, x(0) = 3; nsteps <= maxsteps" begin
    eqs_mov(x, p, t) = x.^2 #the ODE (i.e., the equations of motion)
    exactsol(t, x0) = x0/(1.0-x0*t) #the analytical solution
    t0 = 0.0
    tmax = 0.3
    x0 = 3.0
    x0T = Taylor1(x0, _order)
    tT[1] = t0
    TaylorIntegration.jetcoeffs!(eqs_mov, tT, x0T, nothing)
    @test x0T.coeffs[end] == 3.0^(_order+1)
    δt = (_abstol/x0T.coeffs[end-1])^inv(_order-1)
    @test TaylorIntegration.stepsize(x0T, _abstol) == δt

    tv, xv = taylorinteg(eqs_mov, x0, 0, tmax, _order, _abstol)
    @test length(tv) < 501
    @test length(xv) < 501
    @test length(tv) == 14
    @test length(xv) == 14
    @test xv[1] == x0
    @test tv[end] < 1/3
    @test tv[end] == tmax
    @test abs(xv[end]-exactsol(tv[end], xv[1])) < 2e-14

    tmax = 0.33
    tv, xv = taylorinteg(eqs_mov, x0, t0, tmax, _order, _abstol)
    @test length(tv) < 501
    @test length(xv) < 501
    @test length(tv) == 28
    @test length(xv) == 28
    @test xv[1] == x0
    @test tv[end] < 1/3
    @test tv[end] == tmax
    @test abs(xv[end]-exactsol(tv[end], xv[1])) < 5e-12
end

@testset "Test non-autonomous ODE (1): dot{x}=cos(t)" begin
    let f(x, p, t) = cos(t)
        t0 = 0//1
        tmax = 10.25*(2pi)
        abstol = 1e-20
        order = 25
        x0 = 0.0 #initial conditions such that x(t)=sin(t)
        tv, xv = taylorinteg(f, x0, t0, tmax, order, abstol)
        @test length(tv) < 501
        @test length(xv) < 501
        # @test length(xT[:,2]) < 501
        @test xv[1] == x0
        @test tv[1] == t0
        @test abs(sin(tmax)-xv[end]) < 1e-14

        tmax = 15*(2pi)
        tv, xv = taylorinteg(f, x0, t0, tmax, order, abstol)
        @test length(tv) < 501
        @test length(xv) < 501
        @test xv[1] == x0
        @test abs(sin(tmax)-xv[end]) < 1e-14
    end
end

ext: jl


text:
using StatsBase    # for sampling
using Random, Distributions    # for continuous HMM
using LinearAlgebra

# Definition of HMM model type
struct HMM
    # the states can be string, int, etc.
    hidden_state_space::Array{<:Any, 1}
    emission_space::Array{<:Any, 1}
    initial_distribution::Array{Float64, 1}
    transition_matrix::Array{Float64, 2}
    emission_matrix::Array{Float64, 2}
    name::String
    # check parameter dimensionalities
    function HMM(S, V, μ⁰, A, B, name = "HMM_model")
        J = length(S)    # number of hidden states
        K = length(V)    # number of emitted states
        if length(μ⁰) != J
            error("Size of hidden state space: $J and
            size of initial distribution: $(length(μ⁰)) mismatch")
        elseif size(A) != (J, J)
            error("Size of hidden state space: $J disagrees with " *
            "size of transition matrix: $(size(A, 1)) × $(size(A, 2))")
        elseif size(B) != (J, K)
            error("Size of hidden state space: $J or " *
            "size of emission space: $K disagrees with " *
            "size of emission matrix: $(size(B, 1)) × $(size(B, 2))")
        elseif !all(x -> abs(x - 1.0) < 1e-10, sum(A, dims = 2))
            error("The transition matrix is not normalized")
        elseif !all(x -> abs(x - 1.0) < 1e-10, sum(B, dims = 2))
            error("The emission matrix is not normalized")
        elseif abs(sum(μ⁰) - 1.0) ≥ 1e-6
            # there will be round-off errors for μ⁰, more tolerant
            error("The initial distribution is not normalized")
        end
        new(S, V, μ⁰, A, B, name)
    end
end

# Definition of HMM with continuous emission distribution
struct Gaussian_HMM
    # Gaussian mixture model
    hidden_state_space::Array{<:Any, 1}
    initial_distribution::Array{Float64, 1}
    transition_matrix::Array{Float64, 2}
    Gaussian_center::Array{Float64, 2}
    Gaussian_deviation::Array{Float64, 2}
    mixture_coefficient::Array{Float64, 2}
    name::String
    # check parameter dimensionalities
    function Gaussian_HMM(S, μ⁰, A, M, σ², W, name = "HMM_model")
        J = length(S)    # number of hidden states
        K = size(W, 2)    # number of mixtures
        if length(μ⁰) != J
            error("Size of hidden state space: $J and " *
            "size of initial distribution: $(length(μ⁰)) mismatch")
        elseif size(A) != (J, J)
            error("Size of hidden state space: $J disagrees with " *
            "size of transition matrix: $(size(A, 1)) × $(size(A, 2))")
        elseif !all(x -> abs(x - 1.0) < 1e-10, sum(A, dims = 2))
            error("The transition matrix is not normalized")
        elseif abs(sum(μ⁰) - 1.0) ≥ 1e-6
            error("The initial distribution is not normalized")
        elseif size(M) != (J, K)
            error("The size of Gaussian center matrix: " *
            "$(size(M, 1)) × $(size(M, 2)) is incorrect " *
            "Should be $J × $K")
        elseif size(σ²) != (J, K)
            error("The size of Gaussian deviation matrix: " *
            "$(size(σ², 1)) × $(size(σ², 2)) is incorrect " *
            "Should be $J × $K")
        elseif size(W, 1) != J
            error("The size of mixture coefficient $(size(W, 1)) × $(size(W, 2)) " *
            "does not agree with the size of hidden state space $J")
        elseif !all(x -> abs(x - 1.0) < 1e-12, sum(W, dims = 2))
            error("The Gaussian mixture coefficient is not normalized")
        end
        # check Gaussian distribution parameters
        # TODO: currently not determined whether to use mixed Gaussian
        new(S, μ⁰, A, M, σ², W, name)
    end
end

# All HMM
AbstractHMM = Union{HMM, Gaussian_HMM}

function get_param(model::HMM)
    # A helper function to return the parameters

    return (
        model.initial_distribution,
        model.transition_matrix,
        model.emission_matrix
    )
end

function get_param(model::Gaussian_HMM)
    # Get the parameters of Gaussian HMM

    return (
        model.initial_distribution,
        model.transition_matrix,
        model.Gaussian_center,
        model.Gaussian_deviation,
        model.mixture_coefficient
    )
end

function get_config(model::HMM)
    # A helper function to get model configurations

    return (
        model.hidden_state_space,
        model.emission_space,
        length(model.hidden_state_space),
        length(model.emission_space)
    )
end

function get_config(model::Gaussian_HMM)
    # Get configurations for Gaussian HMM

    return (
        model.hidden_state_space,
        length(model.hidden_state_space),
        size(model.mixture_coefficient, 2)
    )
end

function init_random_HMM(J, K, type = "discrete")
    # initialize an HMM with random (but valid) parameters
    # J is the number of hidden states
    # if discrete, K is the number of observed states
    # elseif continuous, K is the number of Gaussian mixtures
    # the hiddden/observed states are represented by integers

    # the transition probability
    A = rand(Float64, (J, J))
    A ./= sum(A, dims = 2) * ones(Float64, (1, J)) # normalize
    # the initial distribution
    μ⁰ = rand(Float64, J)
    μ⁰ ./= sum(μ⁰)

    if type == "discrete"
        B = rand(Float64, (J, K))
        B ./= sum(B, dims = 2) * ones(Float64, (1, K))
        model = HMM(
            collect(0:(J - 1)),    # the hidden states start from 0
            collect(1:K),
            μ⁰,
            A,
            B,
            "random_HMM"
        )
    elseif type == "Gaussian"
        M = rand(Float64, (J, K))    # Gaussian center
        # note that the deviance is not randomized!
        σ² = ones(Float64, (J, K))    # Gaussian deviation
        W = rand(Float64, (J, K))    # mixture coefficient
        W ./= sum(W, dims = 2)
        model = Gaussian_HMM(
            collect(0:(J - 1)),
            μ⁰,
            A,
            M,
            σ²,
            W,
            "random_HMM"
        )
    end
    return model
end

function emit(model::HMM, l::Int64)
    # generate an emitted sequence and hidden state sequence
    # with given HMM model

    # parameter shorthands
    μ⁰, A, B = get_param(model)
    S, V, J, _ = get_config(model)
    # generate a sequence of emitted states and hidden states with given HMM model
    X = fill(S[1], l)    # allocate mem for hidden states
    # convert state to index
    S_idx = Dict( S[i] => i for i in 1:J )
    # generate X₀
    X[1] = sample(S, Weights(μ⁰))
    # generate the remaining hidden states
    for i = 2:l
        X[i] = sample(S, Weights(A[S_idx[X[i - 1]], :]))
    end
    # generate emitted states
    Y = [ sample(V, Weights(B[S_idx[X[i]], :])) for i in 1:l ]
    return X, Y
end

function emit(model::Gaussian_HMM, l::Int64)
    # Generate observed sequence - continuous version

    μ⁰, A, M, σ², W = get_param(model)
    S, J = get_config(model)
    # Hidden state sequence
    X_idx = zeros(Int64, l)    # allocate mem
    X_idx[1] = sample(1:J, Weights(μ⁰))
    for i = 2:l
        X_idx[i] = sample(1:J, Weights(A[X_idx[i - 1], :]))
    end
    X = S[X_idx]    # convert to hidden state space
    # Observed sequence
    d = [ MixtureModel(Normal, collect(zip(M[i, :], sqrt.(σ²[i, :]))), W[i, :])
    for i in 1:J ]
    Y = [ rand(d[i]) for i in X_idx ]
    return X, Y
end

function b(model::HMM, j::Int64, v, ϵ = 0.01)
    # The probability that v is observed given hidden state jl
    # bⱼ(v) = P(Yₙ = v|Xₙ = sⱼ)

    B = model.emission_matrix
    V = model.emission_space
    v_idx = findfirst(V .== v)    # get the index of v
    return B[j, v_idx]
end

function b(model::Gaussian_HMM, j::Int64, v, ϵ = 0.01)
    # The continuous version

    M = model.Gaussian_center
    σ² = model.Gaussian_deviation
    W = model.mixture_coefficient
    d = MixtureModel(Normal, collect(zip(M[j, :], sqrt.(σ²[j, :]))), W[j, :])
    # take the integral around a small vicinity
    return cdf(d, v + ϵ) - cdf(d, v - ϵ)
end

function forward(model::AbstractHMM, Y::AbstractArray)
    # the forward algorithm

    # variable shorthands
    μ⁰ = model.initial_distribution
    A = model.transition_matrix
    S = model.hidden_state_space
    J = length(S)    # size of hidden state space
    N = length(Y)    # length of observed seq
    # Allocate mem
    # cₙ = P(Yₙ|Y₀^{n-1}), c₀ = P(Y₀)
    # α̂ₙ(i) = P(Xₙ=sᵢ|Y₀ⁿ), α̂₀(i) = P(X₀=sᵢ|Y₀)
    C = zeros(Float64, N)
    α̂ = zeros(Float64, (J, N))
    # Initialization
    α̂[:, 1] = [ b(model, i, Y[1]) for i in 1:J ] .* μ⁰
    C[1] = sum(α̂[:, 1])
    α̂[:, 1] /= C[1]
    # Recursion
    for n = 2:N
        # notice A': summing the multiplication over columns
        α̂[:, n] = [ b(model, i, Y[n]) for i in 1:J ] .* (A' * α̂[:, n - 1])
        C[n] = sum(α̂[:, n])
        α̂[:, n] /= C[n]
    end
    return(α̂, C)
end

function log_likelihood(model::AbstractHMM, Y::AbstractArray)
    # Calculate the log-likelihood of observed (emitted) sequence
    # with given HMM model

    _, C = forward(model, Y)
    sum([ log(c) for c in C ])
end

function backward(model::AbstractHMM, Y::AbstractArray)
    # The backward algorithm

    # variable shorthands
    μ⁰ = model.initial_distribution
    A = model.transition_matrix
    S = model.hidden_state_space
    J = length(S)    # size of hidden state space
    N = length(Y)    # length of observed seq
    # Allocate mem
    # βₙ(i) = (Y_{n+1}^N | Xₙ = sᵢ), β_N = 1
    β̂ = zeros(Float64, (J, N))
    # get the C terms
    α̂, C = forward(model, Y)
    # Initialization
    β̂[:, N] .= 1/C[N]
    # Recursion
    for n = (N - 1):-1:1
        β̂[:, n] = A * ([ b(model, i, Y[n + 1]) for i in 1:J ] .* β̂[:, n + 1]) ./ C[n]
    end
    # also returns α̂ since it's free
    return β̂, α̂, C
end

function baum_welch(initial_model::HMM, Y::AbstractArray, threshold = 1e-2)
    # The Baum-Welch algorithm to infer the parameters of HMM

    # shorthands of variables
    S, V, J, K = get_config(initial_model)
    N = length(Y)    # size of observed sequence
    # convert state to index
    V_idx = Dict( V[i] => i for i in 1:K )

    # the parameters are already initialized inside initial_model
    model = initial_model   # copy the initial model
    Δ = Inf    # difference between the current and previous log-likelihood
    𝓛 = log_likelihood(model, Y)    # initial log likelihood

    while(Δ > threshold)
        # variable shorthands
        μ⁰, A, B = get_param(model)
        # the forward and backward matrix
        β̂, α̂ , C = backward(model, Y)

        # γₙ(i) is the possibility of the n_th hidden value being sᵢ
        # γₙ(i) = P(Xₙ = sᵢ|Y;θ)
        γ = C' .* α̂ .* β̂
        # ξₙ(i, j) is the possibility that the n_th hidden state is sᵢ
        # and transit to n+1_th hidden state of sⱼ
        # ξₙ(i, j) = P(Xₙ = sᵢ, X_{n + 1} = sⱼ|Y; θ)
        ξ = zeros(Float64, (J, J, N - 1))    # allocate mem
        for n  = 1:(N - 1)
            ξ[:, :, n] = α̂[:, n] * (B[:, V_idx[Y[n + 1]]] .* β̂[:, n + 1])' .* A
        end

        # update the transition matrix
        A_new = sum(ξ, dims = 3) ./ sum(γ[:, 1:(N - 1)], dims = 2)
        A_new = dropdims(A_new, dims = 3)
        # update the initial distribution
        μ⁰_new = γ[:, 1]
        # update the emission matrix
        B_new = zeros(Float64, (J, K))    # allocate mem
        for i in 1:K
            # manually add a false: the Nth col will not be accessed
            #B_new[:, i] = sum(γ[:, [Y[1:(N-1)] .== V[i]; false]], dims = 2)
            B_new[:, i] = sum(γ[:, Y .== V[i]], dims = 2)
        end
        B_new ./= sum(B_new, dims = 2)

        # update the HMM model
        model = HMM(S, V, μ⁰_new, A_new, B_new, initial_model.name)

        # update Δln(L)
        new_𝓛 = log_likelihood(model, Y)
        Δ = abs(𝓛 - new_𝓛)
        𝓛 = new_𝓛
        println("Δ:", Δ)
    end
    return model
end

function baum_welch(initial_model::Gaussian_HMM, Y::AbstractArray, threshold = 1e-2)
    # Continuous version of Baum-Welch

    # shorthands of variables
    S, J, K = get_config(initial_model)
    N = length(Y)    # size of observed sequence
    # the parameters are already initialized inside initial_model
    model = initial_model   # copy the initial model
    Δ = Inf    # difference between the current and previous log-likelihood
    𝓛 = log_likelihood(model, Y)    # initial log likelihood
    while(Δ > threshold)
        # variable shorthands
        μ⁰, A, M, σ², W = get_param(model)
        # the forward and backward matrix
        β̂, α̂ , C = backward(model, Y)

        # γₙ(i) is the possibility of the n_th hidden value being sᵢ
        # ξₙ(i, j) is the possibility that the n_th hidden state is sᵢ
        # and transit to n+1_th hidden state of sⱼ
        γ = zeros(Float64, (J, K, N))    # allocate mem
        for n = 1:N
            p = C[n] * α̂[:, n] .* β̂[:, n]
            w = zeros(Float64, (J, K))
            for j = 1:J
                for k = 1:K
                    w[j, k] = W[j, k] * pdf(Normal(M[j, k], sqrt(σ²[j, k])), Y[n])
                end
            end
            w ./= sum(w, dims = 2)
            γ[:, :, n] = p .* w
        end
        ξ = zeros(Float64, (J, J, N - 1))    # allocate mem
        for n  = 1:(N - 1)
            ξ[:, :, n] = α̂[:, n] * ([ b(model, j, Y[n + 1]) for j in 1:J ] .* β̂[:, n + 1])' .* A
        end

        # update the transition matrix
        A_new = sum(ξ, dims = 3) ./ sum(γ[:, :, 1:(N - 1)], dims = (2, 3))
        A_new = dropdims(A_new, dims = 3)    # remove singleton dimension

        # update the initial distribution
        μ⁰_new = sum(γ[:, :, 1], dims = 2)
        μ⁰_new = vec(μ⁰_new)

        # update the mixture coefficient
        W_new = sum(γ, dims = 3)
        W_new = dropdims(W_new, dims = 3)
        W_new ./= sum(W_new, dims = 2)

        # update the mean of Gaussian distribution
        M_new = zeros(Float64, (J, K))
        for k in 1:K
            M_new[:, k] = (γ[:, k, :] * Y) ./ sum(γ[:, k, :], dims = 2)
        end

        # update the deviation of Gaussian
        σ²_new = zeros(Float64, (J, K))
        for i = 1:J
            for k in 1:K
                # σ²_new[i, k] = sum(((Y .- M[i, k]).^2 .* γ[i, k, :])[1:(N - 1)]) /
                # sum(γ[i, k, 1:(N - 1)])
                σ²_new[i, k] = sum((Y .- M[i, k]).^2 .* γ[i, k, :]) / sum(γ[i, k, :])
            end
        end

        # update the HMM model
        model = Gaussian_HMM(S, μ⁰_new, A_new, M_new, σ²_new, W_new, initial_model.name)

        # update Δln(L)
        new_𝓛 = log_likelihood(model, Y)
        Δ = abs(𝓛 - new_𝓛)
        𝓛 = new_𝓛
        println("Δ:", Δ)
    end
    return model
end

function viterbi(model::AbstractHMM, Y::AbstractArray)
    # The Viterbi algorithm to infer the hidden states

    # variable shorthands
    μ⁰ = model.initial_distribution
    A = model.transition_matrix
    S = model.hidden_state_space
    J = length(S)    # size of hidden state space
    N = length(Y)    # length of observed seq
    # Φₙ(i) is the maximum log-likelihood of the hidden states
    # from 0 to n-1, ending with Xₙ = sᵢ
    # Φₙ(i) = max{i₀ to i_{n-1}} lnP(Y₀^n, X₀^{n-1}, Xₙ = sᵢ)
    Φ = zeros(Float64, (J, N))
    Ψ = zeros(Float64, (J, N))    # for traceback
    # Initialization
    Φ[:, 1] = log.(μ⁰) + log.([ b(model, i, Y[1]) for i in 1:J ])
    # Recursion
    for n = 2:N
        for i = 1:J
            idx = argmax(log.(A[:, i]) + Φ[:, n - 1])
            Ψ[i, n] = idx
            Φ[i, n] = log(b(model, i, Y[n])) + log(A[idx, i]) + Φ[idx, n - 1]
        end
    end
    # traceback
    i_max = zeros(Int64, N)
    i_max[N] = argmax(Φ[:, N])
    for n = (N - 1):-1:1
        i_max[n] = Ψ[i_max[n + 1], n + 1]
    end
    # get Sᵢ
    X_max = S[i_max]
    return(X_max)
end


text:
"""
    PkgTPL

[Repository](https://github.com/Nosferican/PkgTPL.jl)
"""
module PkgTPL
    for (root, dirs, files) ∈ walkdir(joinpath(@__DIR__))
        foreach(file -> include(joinpath(root, file)), filter!(file -> occursin(r"^\d{2}_\w+\.jl$", file), files))
    end
end


text:
@testset "DataStructures" begin
@testset "BasicBinaryTree" begin
    @testset "BasicBinaryTree: BinaryTree" begin
        tree=BinaryTree{Int}(10,0)
        insert!(tree,1,10)
        insert!(tree,1,20)
        insert!(tree,2,50)
        #=
        [1] 0
         |       \
        [2] 10  [3] 20
         |
        [4] 50
        =#
        @test tree.n == 4
        @test depth(tree,4) == 2
        @test isleaf(tree,2) == false
        @test ch(tree,1,true) == 3
        @test tree.val[left(tree)] == 50
        @test height(tree) == 2
    end
    @testset "DisjointSet: DisjointSet" begin
        set=DisjointSet(10)
        set.par[4]=1
        set.par[3]=1
        set.par[2]=3
        set.par[5]=8
        @test find(set,2) == 1
        merge!(set,5,4)
        @test find(set,5) == find(set,2)
    end
end
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
module ICARTT

include("icartt_exceptions.jl")
include("units_icartt.jl")
include("read_icartt.jl")
include("user_utils.jl")

import .ReadICARTT: read_icartt_file;
import .ICARTTUtils: get_merge_data, search_icartt_variables;
export read_icartt_file, get_merge_data, search_icartt_variables;

end # module

max_stars_repo_name: joshua-laughner/Julia-ICARTT

max_issues_repo_name: joshua-laughner/Julia-ICARTT

max_forks_repo_name: joshua-laughner/Julia-ICARTT


text:
# Matern 5/2 ARD covariance Function

"""
    Mat52Ard <: MaternARD

ARD Matern 5/2 kernel (covariance)
```math
k(x,x') = σ²(1 + √5|x-x'|/L + 5|x-x'|²/(3L²))\\exp(- √5|x-x'|/L)
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)`` and signal standard deviation ``σ`` where
``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct Mat52Ard{T<:Real} <: MaternARD
    "Inverse squared length scale"
    iℓ2::Vector{T}
    "Signal variance"
    σ2::T
    "Priors for kernel parameters"
    priors::Array
end

"""
Matern 5/2 ARD covariance Function

    Mat52Ard(ll::Vector{Real}, lσ::Real)

# Arguments
  - `ll::Vector{Real}`: vector of length scales (given on log scale)
  - `lσ::Real`: signal standard deviation (given on log scale)
"""
Mat52Ard(ll::Vector{T}, lσ::T) where T = Mat52Ard{T}(exp.(-2 .* ll), exp(2 * lσ), [])

function set_params!(mat::Mat52Ard, hyp::AbstractVector)
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat52 kernel has $(num_params(mat)) parameters, received $(length(hyp))."))
    @views @. mat.iℓ2 = exp(-2 * hyp[1:(end-1)])
    mat.σ2 = exp(2 * hyp[end])
end

get_params(mat::Mat52Ard) = [-log.(mat.iℓ2) / 2; log(mat.σ2) / 2]
get_param_names(mat::Mat52Ard) = [get_param_names(mat.iℓ2, :ll); :lσ]
num_params(mat::Mat52Ard) = length(mat.iℓ2) + 1

cov(mat::Mat52Ard, r::Number) =
    (s = √5 * r; mat.σ2 * (1 + s + s^2 / 3) * exp(-s))

dk_dll(mat::Mat52Ard, r::Real, wdiffp::Real) =
    (s = √5 * r; 5 / 3 * mat.σ2 * wdiffp *(1 + s) * exp(-s))


text:
module SoleViz

# Write your package code here.

end

max_stars_repo_name: aclai-lab/SoleVisualizations

max_issues_repo_name: aclai-lab/SoleVisualizations

max_forks_repo_name: aclai-lab/SoleVisualizations


text:
# Base.zero(pm::PermMatrix) = PermMatrix(pm.perm, zero(pm.vals))

# TODO
# to make a mat block differentiable
#YaoBlocks.niparams(x::GeneralMatrixBlock{N,N}) where N = 1<<2N
#YaoBlocks.getiparams(x::GeneralMatrixBlock) where N = (vec(x.mat)...,)

# to make a scale block differentiable
#YaoBlocks.niparams(x::Scale{<:Number}) = 1
#YaoBlocks.getiparams(x::Scale{<:Number}) = factor(x)


text:
# Scheduling isn't really possible directly with the solver yet :/
# Most of this file is based on http://www.hakank.org/julia/constraints/
# which is a website by Håkan Kjellerstrand
# please check it out :)

# by  Håkan Kjellerstrand
function cumulative(model, start, duration, resource, limit; times_max = nothing)
    tasks = [i for i in 1:length(start) if resource[i] > 0 && duration[i] > 0]
    num_tasks = length(tasks)

    times_min = minimum(round.(Int,[JuMP.lower_bound(start[i]) for i in tasks]))
    if times_max === nothing
        times_max = maximum(round.(Int,[JuMP.upper_bound(start[i])+duration[i] for i in tasks]))
    end
    for t in times_min:times_max
        b  = @variable(model, [1:num_tasks], Bin)
        for i in tasks
            # The following don't work since ConstraintSolver don't
            # support nonlinear constraints
            # @constraint(model,sum([(start[i] <= t) * (t <= start[i] + duration[i])*resource[i] for i in tasks])  <= b)

            # is this task active during this time t?
            @constraint(model, b[i] := { start[i] <= t && t < start[i]+duration[i] }) # is this task active in time t ?
        end
        # Check that there's no conflicts in time t
        @constraint(model,sum([b[i]*resource[i] for i in tasks]) <= limit)
  end
end

#=
  Furniture moving (scheduling) in Julia ConstraintSolver.jl
  From Marriott & Stuckey: "Programming with constraints", page  112f
  Changed for this test case to be easier ;)

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also his Julia page: http://www.hakank.org/julia/
=#
function furniture_moving()
    model = Model(optimizer_with_attributes(CS.Optimizer,
                                                "logging"=>[],
                                                "branch_split"=>:InHalf,
                                                "time_limit"=>50
                                        ))

    # Furniture moving
    n = 4
    # [piano, chair, bed, table]
    durations = [30,10,15,15]
    # resource needed per task
    resources = [1,1,1,1] # <- changed for this test case to solve it faster
    max_end_time = 45
    @variable(model, 0 <= start_times[1:n] <= max_end_time, Int)
    @variable(model, minimum(durations) <= end_times[1:n]  <= max_end_time, Int)
    @variable(model, 1 <= limit <= 3, Int)
    @variable(model, 0 <= max_time <= max_end_time,Int)

    for i in 1:n
        @constraint(model,end_times[i] == start_times[i] + durations[i])
    end
    @constraint(model, end_times .<= max_time)
    cumulative(model, start_times, durations, resources, limit; times_max = max_end_time)

    # Solve the problem
    @objective(model, Min, max_time)
    optimize!(model)

    status = JuMP.termination_status(model)
    @test status == MOI.OPTIMAL
    @test JuMP.value(limit) ≈ 3
    @test JuMP.value(max_time) ≈ 30
    @test JuMP.value(start_times[1]) ≈ 0
    for i = 1:n
        @test JuMP.value(start_times[i])+durations[i] ≈ JuMP.value(end_times[i])
    end
end

@testset "Scheduling" begin
@testset "Furniture" begin
    furniture_moving()
end
end


text:
using WignerSymbols
using ThreadSafeDicts

function test_serial(n::Int)
    new_cache = WignerCache()  # create a new cache
    sum = zero(BigFloat)
    for i in 1:n
        for j in 1:n
            for k in 1:n
                sum += wigner3j(new_cache, BigFloat, i, j, k, 0, 0)
            end
        end
    end
    sum
end

function test_bounded(n::Int)
    WignerSymbols.__init__()
    sum = zero(BigFloat)
    for i in 1:n
        for j in 1:n
            for k in 1:n
                sum += wigner3j(BigFloat, i, j, k, 0, 0)
            end
        end
    end
    sum
end

function test_threaded(n::Int)
    WignerSymbols.__init__()
    summands = zeros(BigFloat, Threads.nthreads())
    Threads.@threads for i in 1:n
        tid = Threads.threadid()
        for j in 1:n
            for k in 1:n
                summands[tid] += wigner3j(BigFloat, i, j, k, 0, 0)
            end
        end
    end
    sum(summands)
end

function test_threaded_specialized(n::Int)
    caches = WignerSymbols.get_thread_caches(BigFloat, n)
    summands = zeros(BigFloat, Threads.nthreads())
    Threads.@threads for i in 1:n
        tid = Threads.threadid()
        for j in 1:n
            for k in 1:n
                summands[tid] += wigner3j(caches[tid], BigFloat, i, j, k, 0, 0)
            end
        end
    end
    sum(summands)
end

maxj = 150

## serial benchmarks

@time begin
    println("default serial:  ", test_serial(maxj))
    GC.gc()
end

@time begin
    println("bounded serial:  ", test_bounded(maxj))
    GC.gc()
end

## threaded benchmarks

@time begin
    println("bounded thread:  ", test_threaded(maxj))
    GC.gc()
end

@time begin
    println("specialized th:  ", test_threaded_specialized(maxj))
    GC.gc()
end

##


text:
# DEBUG I'm trying to set this env before any possible using of CuArrays
#
# integer in bytes
# 7G = 7 * 1000 * 1000 * 1000
# FIXME does not work
#
# For 1070
g=7.0
# For 2080 Ti
# g=9.0
ENV["CUARRAYS_MEMORY_LIMIT"] = convert(Int, round(g * 1024 * 1024 * 1024))

# this manually set works
# CuArrays.usage_limit[] = parse(Int, ENV["CUARRAYS_MEMORY_LIMIT"])
# Or better call the init function
# CuArrays.__init_memory__()
# Or
# CuArrays.__init__()
# check status
# CuArrays.memory_status()

using CuArrays

using Flux
# Zygote
using Flux.Data.MNIST
using Flux: @epochs
using Flux: throttle
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated
using Random
using Statistics

using Images: channelview
using Base.Iterators: partition

# for padarray and Fill
using Images

export load_MNIST, load_CIFAR10

using MLDatasets

import Base.show
import Base.display
using FileIO
using FileIO: @format_str
using Images: channelview, colorview, RGB

struct MyImage
    arr::AbstractArray
end

"""Write a 3-channel array to PNG stream

"""
function Base.show(io::IO, ::MIME"image/png", x::MyImage)
    img = cpu(x.arr)
    size(img)[3] in [1, 3] || error("Unsupported channel size: ", size(img)[3])
    if size(img)[3] == 3
        getRGB(X) = colorview(RGB, permutedims(X, (3,1,2)))
        img = getRGB(img)
    end
    FileIO.save(FileIO.Stream(format"PNG", io), img)
    # FileIO.save(File(format"PNG", path * ".png"), img)
end

struct EmacsDisplay <: AbstractDisplay end

"""Display a PNG by writing it to tmp file and show the filename. The filename
would be replaced by an Emacs plugin.

"""
function Base.display(d::EmacsDisplay, mime::MIME"image/png", x)
    # path, io = Base.Filesystem.mktemp()
    # path * ".png"
    path = tempname() * ".png"
    open(path, "w") do io
        show(io, mime, x)
    end

    println("$(path)")
    println("#<Image: $(path)>")
end
function display(d::EmacsDisplay, x)
    display(d, "image/png", x)
end

function register_EmacsDisplay_backend()
    for d in Base.Multimedia.displays
        if typeof(d) <: EmacsDisplay
            return
        end
    end
    # register as backend
    pushdisplay(EmacsDisplay())
    # popdisplay()
    nothing
end

register_EmacsDisplay_backend()

function test_display()
    ds, test_ds = load_CIFAR10_ds(batch_size=128);
    ds, test_ds = load_MNIST_ds(batch_size=128);
    x, y = next_batch!(ds) |> gpu;
    img = x[:,:,:,1];

    # display
    displayable("image/png")
    display(MyImage(img))

    # explicit on a (new) EmacsDisplay
    displayable(EmacsDisplay(), "image/png")
    display(EmacsDisplay(), MyImage(img));

    # showable
    Base.showable(MIME("image/png"), MyImage(img))
    Base.showable(MIME("image/png"), img)
end

function test()
    train_x, train_y = MLDatasets.MNIST.traindata();
    test_x,  test_y  = MLDatasets.MNIST.testdata();
    train_x[:,:,1]
    train_y[1]
    train_x, train_y = MLDatasets.CIFAR10.traindata();
    test_x,  test_y  = MLDatasets.CIFAR10.testdata();
    train_x[:,:,:,1]
    train_y[1]
end

# There is a data loader PR: https://github.com/FluxML/Flux.jl/pull/450
mutable struct DataSetIterator
    raw_x::AbstractArray
    raw_y::AbstractArray
    index::Int
    batch_size::Int
    nbatch::Int
    DataSetIterator(x,y,batch_size) = new(x,y,0,batch_size, convert(Int, floor(size(x)[end]/batch_size)))
end

Base.show(io::IO, x::DataSetIterator) = begin
    println(io, "DataSetIterator:")
    println(io, "  batch size: ", x.batch_size)
    println(io, "  number of batches: ", x.nbatch)
    println(io, "  x data shape: ", size(x.raw_x))
    println(io, "  y data shape: ", size(x.raw_y))
    print(io, "  current index: ", x.index)
end

function next_batch!(ds::DataSetIterator)
    if (ds.index + 1) * ds.batch_size > size(ds.raw_x)[end]
        ds.index = 0
    end
    if ds.index == 0
        indices = shuffle(1:size(ds.raw_x)[end])
        ds.raw_x = ds.raw_x[:,:,:,indices]
        ds.raw_y = ds.raw_y[:,indices]
    end
    a = ds.index * ds.batch_size + 1
    b = (ds.index + 1) * ds.batch_size

    ds.index += 1
    # FIXME use view instead of copy
    # FIXME properly slice with channel-last format
    return ds.raw_x[:,:,:,a:b], ds.raw_y[:,a:b]
end

function my_convert2image(x)
    # DEPRECATED
    #
    # FIXME the digits are horizontal-major instead of vertical-major, which is
    # very weird, create problem for my visualization. I have to use
    # MLDatasets.MNIST.convert2image, which turns floating array into white
    # background digits. WTF. I would want to just rotate it.
    reshape(MLDatasets.MNIST.convert2image(reshape(x, 28,28,:)), 28,28,1,:)
end
function load_MNIST_ds(;batch_size)
    # Ideally, it takes a batch size, and should return an iterator that repeats
    # infinitely, and shuffle before each epoch. The data should be moved to GPU
    # on demand. I prefer to just move online during training.
    train_x, train_y = MLDatasets.MNIST.traindata();
    test_x,  test_y  = MLDatasets.MNIST.testdata();
    # reshape to add channel, and onehot y encoding
    #
    # I'll just permute dims to make it column major
    train_ds = DataSetIterator(reshape(permutedims(train_x, [2,1,3]), 28,28,1,:),
                               onehotbatch(train_y, 0:9), batch_size);
    # FIXME test_ds should repeat only once
    # TODO add iterator interface, for x,y in ds end
    test_ds = DataSetIterator(reshape(permutedims(test_x, [2,1,3]), 28,28,1,:),
                              onehotbatch(test_y, 0:9), batch_size);
    # x, y = next_batch!(train_ds);
    return train_ds, test_ds
end

function load_CIFAR10_ds(;batch_size)
    train_x, train_y = MLDatasets.CIFAR10.traindata();
    test_x,  test_y  = MLDatasets.CIFAR10.testdata();
    train_ds = DataSetIterator(permutedims(train_x, [2,1,3,4]),
                               onehotbatch(train_y, 0:9), batch_size);
    test_ds = DataSetIterator(permutedims(test_x, [2,1,3,4]),
                              onehotbatch(test_y, 0:9), batch_size);
    return train_ds, test_ds
end

function compute_mean_var()
    ds, test_ds = load_CIFAR10_ds(batch_size=128)
    # FIXME do I need to consider test_ds?
    train_x, train_y = MLDatasets.CIFAR10.traindata();
    test_x,  test_y  = MLDatasets.CIFAR10.testdata();
    # dim=3 is the color channel, it will return a value for each channel

    # FIXME
    # 1. the mean var of MLDatasets and Metalhead.CIFAR10 are different
    # 2. when subtracting mean, half of values will be negative?
    mean(train_x)
    mean(train_x, dims=(1,2,4))
    std(train_x, dims=(1,2,4))
    minimum(train_x)

    # NOTE: This is in back.jl
    (trainX, trainY), (valX, valY), (testX, testY) = load_CIFAR10();
    trainX = cpu.(trainX);
    size(trainX)
    size(trainX[1])
    typeof(trainX)
    mean(cat(trainX..., dims=4), dims=(1,2,4))
    std(cat(trainX..., dims=4), dims=(1,2,4))
    mean(trainX[1])
end

function flipx(p)
    function f(img)
        if rand() < p
            reverse(img, dims=2)
        else
            img
        end
    end
end

function pad_and_crop(img)
    w = size(img)[1]
    h = size(img)[2]
    # FIXME this returns OffsetArrays, with possibly negative or 0 index
    padded = padarray(img, Fill(0,(4,4,0), (4,4,0)))
    a = rand(-3:4)
    b = rand(-3:4)
    # @show a
    # @show b
    # convert OffsetArrays to normal array with parent()
    parent(padded[a:a+w-1, b:b+h-1, :])
end

function augment_one(img)
    pad_and_crop(flipx(0.5)(img))
end

struct Augment end
Flux.@treelike Augment

# - will this apply different random op to images in a batch? Yes.
# - can this apply to batch automatically? NO!
# FIXME performance, and also on CPU/GPU?
function (a::Augment)(x::AbstractArray)
    length(size(x)) == 4 || error("Dim of input must be 4, got $(length(size(x)))")
    res = map(collect(1:size(x)[end])) do i
        augment_one(x[:,:,:,i])
    end
    cat(res..., dims=4)
end

function Base.show(io::IO, l::Augment)
  print(io, "Augment()")
end

function test_augment()
    # padarray

    ds, test_ds = load_CIFAR10_ds(batch_size=16);
    x,y = next_batch!(ds);

    # some testing
    img = x[:,:,:,1];
    sample_and_view(img)
    sample_and_view(reverse(img, dims=2))
    sample_and_view(pad_and_crop(flipx(0.5)(img)))

    # testing augment layer
    aug = Augment()
    sample_and_view(x)
    sample_and_view(aug(x))
    size(aug(x))
end

function test()
    ds, test_ds = load_MNIST_ds(batch_size=128);
    ds.nbatch
    test_ds.nbatch
    x, y = next_batch!(ds);
    size(x)
    size(y)
    sample_and_view(x)
    # size(x)
    train_ds.raw_y[1:10]
    for i=1:1000
        next_batch!(train_ds)
        @show train_ds.index
    end
    60000 / 32
end

"""
Sample up to 10 images and show the image and label. If less than 10, use all.

"""
function sample_and_view(x, y=nothing, model=nothing)
    if length(size(x)) < 4
        x = x[:,:,:,:]
    end
    if typeof(x) <: TrackedArray
        x = x.data
    end
    size(x)[1] in [28,32,56] ||
        error("Image size $(size(x)[1]) not correct size. Currently support 28 or 32.")
    num = min(size(x)[4], 10)
    @info "Showing $num images .."
    if num == 0 return end
    imgs = cpu(hcat([x[:,:,:,i] for i in 1:num]...))
    # viewrepl(imgs)
    display(MyImage(imgs))
    if y != nothing
        labels = onecold(cpu(y[:,1:num]), 0:9)
        @show labels
    end
    if model != nothing
        preds = onecold(cpu(model(gpu(x))[:,1:num]), 0:9)
        @show preds
    end
    nothing
end

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
struct Constant <: Degree{0} end

"""
Constant b-splines are *nearest-neighbor* interpolations, and effectively
return `A[round(Int,x)]` when interpolating.
"""
Constant

function positions(::Constant, ax, x)  # discontinuity occurs at half-integer locations
    xm = roundbounds(x, ax)
    δx = x - xm
    fast_trunc(Int, xm), δx
end

value_weights(::Constant, δx) = (1,)
gradient_weights(::Constant, δx) = (0,)
hessian_weights(::Constant, δx) = (0,)

padded_axis(ax::AbstractUnitRange, ::BSpline{Constant}) = ax


text:
"""
`module JAC.InteractionStrength`
    ... a submodel of JAC that contains all methods for evaluating the interaction strength (reduced matrix elements)
        for various atomic interactions.
"""
module InteractionStrength

    using  GSL, ..AngularMomentum, ..Basics, ..Defaults, ..ManyElectron, ..Nuclear, ..Radial, ..RadialIntegrals
    ##x global JAC_counter = 0

    """
    `struct  InteractionStrength.XLCoefficient`  ... defines a type for coefficients of the two-electron (Breit) interaction

        + mu        ::Int64      ... Kind of integral, presently fixed to mu = 5 to keep similarity with RATIP during development.
        + nu        ::Int64      ... Rank of the integral.
        + a         ::Orbital    ... Orbitals a, b, c, d.
        + b         ::Orbital
        + c         ::Orbital
        + d         ::Orbital
        + coeff     ::Float64    ... corresponding coefficient.
    """
    struct XLCoefficient
        mu          ::Int64
        nu          ::Int64
        a           ::Orbital
        b           ::Orbital
        c           ::Orbital
        d           ::Orbital
        coeff       ::Float64
    end

    """
    `InteractionStrength.dipole(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the  <a|| d ||b>  reduced matrix element of the dipole operator for orbital functions a, b.
            A value::Float64 is returned.
    """
    function dipole(a::Orbital, b::Orbital, grid::Radial.Grid)
        wa = AngularMomentum.CL_reduced_me(a.subshell, 1, b.subshell) * RadialIntegrals.rkDiagonal(1, a, b, grid)
        return( wa )
    end

    """
    `InteractionStrength.hamiltonian_nms(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the  <a|| h_nms ||b>  reduced matrix element of the normal-mass-shift Hamiltonian for orbital
            functions a, b. A value::Float64 is returned.  For details, see Naze et al., CPC 184 (2013) 2187, Eq. (37).
    """
    function hamiltonian_nms(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        if  a.subshell.kappa != b.subshell.kappa   return( 0. )   end
        wa = RadialIntegrals.isotope_nms(a, b, nm.Z, grid) / (2 * nm.mass)
        println("**  <$(a.subshell) || h^nms || $(b.subshell)>  = $wa" )
        return( wa )
    end

    """
    `InteractionStrength.hfs_t1(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| t^(1) ||b> reduced matrix element for the HFS coupling to the magnetic-dipole moment
            of the nucleus for orbital functions a, b. A value::Float64 is returned.
    """
    function hfs_t1(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC, Eq. (49) ... test for the proper definition of the C^L tensors.
        minusa = Subshell(1, -a.subshell.kappa)
        ##x wa = - Defaults.getDefaults("alpha") * (a.subshell.kappa + b.subshell.kappa) * AngularMomentum.CL_reduced_me_rb(minusa, 1, b.subshell) *
        wb =   - (a.subshell.kappa + b.subshell.kappa) * AngularMomentum.CL_reduced_me_rb(minusa, 1, b.subshell)
        wc =   RadialIntegrals.rkNonDiagonal(-2, a, b, grid)
        wa =   wb * wc
        #
        ##x println("**  <$(a.subshell) || t1 || $(b.subshell)>  = $wa   = $wb * $wc" )
        return( wa )
    end

    """
    `InteractionStrength.hfs_t2(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| t^(2) ||b> reduced matrix element for the HFS coupling to the electric-quadrupole moment of
            the nucleus for orbital functions a, b. A value::Float64 is returned.
     """
    function hfs_t2(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC, Eq. (49) ... test for the proper definition of the C^L tensors.
        wb = - AngularMomentum.CL_reduced_me_rb(a.subshell, 2, b.subshell)
        ## wc =   RadialIntegrals.rkDiagonal(-3, a, b, grid)
        wc =   RadialIntegrals.rkDiagonal(-3, a, b, grid)
        wa =   wb * wc
        #
        ##x println("**  <$(a.subshell) || t2 || $(b.subshell)>  = $wa   = $wb * $wc" )
        return( wa )
    end

    """
    `InteractionStrength.MbaAbsorptionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp, absorption) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega, and
            within the given gauge. A value::Float64 is returned.
    """
    function MbaAbsorptionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        wa = MbaEmissionCheng(mp, gauge, omega, b, a, grid)
        wa = conj(wa)
        return( wa )
    end

    """
    `InteractionStrength.MbaEmissionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp,emission) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega, and
            within the given gauge. A value::Float64 is returned.
    """
    function MbaEmissionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        #
        if       gauge == Basics.Magnetic
            ChengI = AngularMomentum.ChengI(-kapa, kapb, AngularJ64(mp.L))
            wa     = -1.0im / sqrt(mp.L*(mp.L+1)) * ChengI * (kapb + kapa) * RadialIntegrals.GrantILplus(mp.L, q, a, b, grid::Radial.Grid)
        #
        elseif   gauge == Basics.Babushkin
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L))
            wr     = (mp.L+1) * RadialIntegrals.GrantJL(mp.L, q, a, b, grid::Radial.Grid)
            wr     = wr +  (kapb-kapa-mp.L-1) * RadialIntegrals.GrantIL0(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr +  (kapb-kapa+mp.L+1) * RadialIntegrals.GrantIL0(mp.L+1, q, b, a, grid::Radial.Grid)
            wa     = 1.0im / sqrt(mp.L*(mp.L+1)) * ChengI * wr
        #
        elseif   gauge == Basics.Coulomb
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L))
            wr     = (mp.L+1) * (kapb-kapa+mp.L) / (2mp.L+1) * RadialIntegrals.GrantIL0(mp.L-1, q, a, b, grid::Radial.Grid)
            wr     = wr +  ((mp.L+1) * (kapb-kapa+mp.L) - (2mp.L+1) * (kapb-kapa)) / (2mp.L+1) *
                                                               RadialIntegrals.GrantIL0(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr + (mp.L+1) * (kapb-kapa-mp.L) / (2mp.L+1) * RadialIntegrals.GrantIL0(mp.L-1, q, b, a, grid::Radial.Grid)
            wr     = wr +  ((mp.L+1) * (kapb-kapa-mp.L) - (2mp.L+1) * (kapb-kapa)) / (2mp.L+1) *
                                                               RadialIntegrals.GrantIL0(mp.L+1, q, b, a, grid::Radial.Grid)
            wa     = 1.0im / sqrt(mp.L*(mp.L+1)) * ChengI * wr
        else     error("stop a")
        end
        wa = conj(wa)

        return( wa )
    end

    """
    `InteractionStrength.MbaEmissionAndrey(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp, emission) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega,
            and within the given gauge. A value::Float64 is returned. At present, only the magnetic matrix elements are
            implemented.
    """
    function MbaEmissionAndrey(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        ja   = Basics.subshell_j(a.subshell);   jb   = Basics.subshell_j(b.subshell);
        #
        if       gauge == Basics.Magnetic
            wa = -1.0im * sqrt( (2mp.L+1) * (Basics.subshell_2j(a.subshell)+1) ) / sqrt(4pi) / sqrt(mp.L * (mp.L+1)) * (-1)^mp.L * (kapb + kapa)
            wa = wa * AngularMomentum.ClebschGordan(ja, AngularM64(1//2), AngularJ64(mp.L), AngularM64(0), jb, AngularM64(1//2))
            wa = wa * RadialIntegrals.GrantILplus(mp.L, q, a, b, grid::Radial.Grid)
            wa = conj(wa)
        #
        elseif   gauge == Basics.Babushkin
            wa = 0.
        #
        elseif   gauge == Basics.Coulomb
            wa = 0.
        else     error("stop a")
        end

        return( wa )
    end

    """
    `InteractionStrength.multipoleTransition(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) multipole-transition interaction strength
            <b || T^(Mp, absorption) || a> due to Johnson (2007) for the interaction with the Mp multipole component of the
            radiation field and the transition frequency omega, and within the given gauge. A value::Float64 is returned.
    """
    function multipoleTransition(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        function besselPrime_jl(L::Int64, x::Float64)    return( GSL.sf_bessel_jl(L-1, x) - (L+1)/x * GSL.sf_bessel_jl(L, x) )       end

        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        mtp  = min(size(a.P, 1), size(b.P, 1))
        !(grid.mesh == MeshGL)  &&  error("Only for MeshGL implemented so far.")
        #
        if       gauge == Basics.Magnetic
            ChengI = AngularMomentum.ChengI(-kapa, kapb, AngularJ64(mp.L));   if  abs(ChengI) < 1.0e-10  return( 0. )   end
            wa = Complex(0.)
            for  i = 2:mtp
                wa = wa + (kapa+kapb) / (mp.L+1) * GSL.sf_bessel_jl(mp.L, q * grid.r[i]) * (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * grid.wr[i]
            end
            wa = ChengI * wa
            #
        elseif   gauge == Basics.Velocity
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L));    if  abs(ChengI) < 1.0e-10  return( 0. )   end
            wa = Complex(0.)
            for  i = 2:mtp
                wa = wa - (kapa-kapb) / (mp.L+1) *
                          ( besselPrime_jl(mp.L, q * grid.r[i]) + GSL.sf_bessel_jl(mp.L, q * grid.r[i]) / (q * grid.r[i]) ) *
                          (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * grid.wr[i]
                wa = wa + mp.L * GSL.sf_bessel_jl(mp.L, q * grid.r[i]) / (q * grid.r[i]) * (a.P[i] * b.Q[i] - a.Q[i] * b.P[i]) * grid.wr[i]
            end
            wa = ChengI * wa
            #
        elseif   gauge == Basics.Length
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L));    if  abs(ChengI) < 1.0e-10  return( 0. )   end
            wa = Complex(0.)
            for  i = 2:mtp
                wa = wa + GSL.sf_bessel_jl(mp.L, q * grid.r[i]) * (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) * grid.wr[i]
                        + GSL.sf_bessel_jl(mp.L+1, q * grid.r[i]) *
                          ( (kapa-kapb) / (mp.L+1) * (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) +
                            (a.P[i] * b.Q[i] - a.Q[i] * b.P[i]) ) * grid.wr[i]
            end
            wa = ChengI * wa
            #
        else     error("stop a")
        end

        return( wa )
    end

    """
    `InteractionStrength.schiffMoment(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the  <a|| h^(Schiff-moment) ||b>  reduced matrix element of the Schiff-moment Hamiltonian for orbital
            functions a, b and for the nuclear density as given by the nuclear model. A value::Float64 is returned.
    """
    function schiffMoment(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        printstyled("\nWarning -- InteractionStrength.schiffMoment():: Not yet implemented.", color=:cyan)
        wb = 1.0 + 2.0im
        return( wb )
    end

    """
    `InteractionStrength.weakCharge(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the  <a|| h^(weak-charge) ||b>  reduced matrix element of the weak-charge Hamiltonian for orbital functions
            a, b and for the nuclear density as given by the nuclear model. A value::Float64 is returned.
    """
    function weakCharge(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        printstyled("\nWarning -- InteractionStrength.weakCharge():: Not yet implemented.", color=:cyan)
        wb = 1.0 + 2.0im
        return( wb )
    end

    """
    `InteractionStrength.XL_Breit(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the the effective Breit interaction strengths X^L_Breit (abcd) for given rank L and orbital functions
            a, b, c and d  at the given grid. A value::Float64 is returned. At present, only the zero-frequency Breit
            interaction is taken into account.
    """
    function XL_Breit(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)
        ja2 = Basics.subshell_2j(a.subshell)
        jb2 = Basics.subshell_2j(b.subshell)
        jc2 = Basics.subshell_2j(c.subshell)
        jd2 = Basics.subshell_2j(d.subshell)
        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||  L == 0
            return( 0. )
        end
        #
        xcList   = XL_Breit0_coefficients(L,a,b,c,d)
        ##x println("\nBreit strength: L=$L, nx = $(length(xcList)) ")
        XL_Breit = XL_Breit0_densities(xcList, grid)
        #
        ##x XL_Breit = 0.
        ##x for  xc in xcList
        ##x     wb = RadialIntegrals.W5_Integral(xc.mu, xc.nu, xc.a ,xc.c ,xc.b ,xc.d, grid)
        ##x     println("Breit0: nu=$(xc.nu), a=$(xc.a.subshell), b=$(xc.b.subshell), c=$(xc.c.subshell), d=$(xc.d.subshell), " *
        ##x             "coeff=$(xc.coeff), wx=$wb, wges= $(xc.coeff * wb)")
        ##x     XL_Breit = XL_Breit + xc.coeff * wb
        ##x end
        ##x #
        return( XL_Breit )
    end

    """
    `InteractionStrength.XL_Breit0_coefficients(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital)`
        ... evaluates the combinations and pre-coefficients for the zero-frequency Breit interaction
            X^L_Breit (omega=0.; abcd) for given rank L and orbital functions a, b, c and d. A list of coefficients
            xcList::Array{XLCoefficient,1} is returned.
    """
    function XL_Breit0_coefficients(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital)
        xcList = XLCoefficient[]

        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc                end
        if   abs(xc)  <  1.0e-10    return( xcList )   end

        # Consider the individual contributions from sum_nu and sum_mu. First, take T^(nu,L)_mu = R^(nu,L)_mu
        nu = L - 1

        if  rem(la+lc+nu,2) == 1   &&   rem(lb+ld+nu,2) == 1   &&   L != 0
            wa = (L+1) / ( L*(L+L-1)*(L+L+1) )
            # mu = 1
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa+L) * (d.subshell.kappa-b.subshell.kappa+L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, b, c, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, a, d, c, xcc) )   end
            # mu = 2
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa-L) * (d.subshell.kappa-b.subshell.kappa-L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, d, a, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, c, b, a, xcc) )   end
            # mu = 3
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa+L) * (d.subshell.kappa-b.subshell.kappa-L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, d, c, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, a, b, c, xcc) )   end
            # mu = 4
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa-L) * (d.subshell.kappa-b.subshell.kappa+L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, b, a, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, c, d, a, xcc) )   end
        end
        #
        nu = L

        if  rem(la+lc+nu,2) == 1   &&   rem(lb+ld+nu,2) == 1   &&   L != 0
            wa = - (a.subshell.kappa + c.subshell.kappa) * (b.subshell.kappa + d.subshell.kappa) / (L*(L+1))
            # mu = 1
            xcc = xc * wa
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, b, c, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, a, d, c, xcc) )   end
            # mu = 2
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, d, a, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, c, b, a, xcc) )   end
            # mu = 3
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, d, c, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, a, b, c, xcc) )   end
            # mu = 4
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, b, a, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, c, d, a, xcc) )   end
        end
        #
        nu = L + 1

        if  rem(la+lc+nu,2) == 1   &&   rem(lb+ld+nu,2) == 1   &&   L != 0
            wa =  L / ( (L+1)*(L+L+1)*(L+L+3) )
            # mu = 1
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa - L - 1) * (d.subshell.kappa - b.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, b, c, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, a, d, c, xcc) )   end
            # mu = 2
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa + L + 1) * (d.subshell.kappa - b.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, d, a, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, c, b, a, xcc) )   end
            # mu = 3
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa - L - 1) * (d.subshell.kappa - b.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, d, c, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, a, b, c, xcc) )   end
            # mu = 4
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa + L + 1) * (d.subshell.kappa - b.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, b, a, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, c, d, a, xcc) )   end
        end

        # Add contributions of the S^k_mu integrals
        if  rem(la+lc+L-1,2) == 1   &&   rem(lb+ld+L+1,2) == 1
            # mu = 1
            wb =  1 / ( (L+L+1)*(L+L+1) )
            xcc = xc * wb * (c.subshell.kappa - a.subshell.kappa + L) * (d.subshell.kappa - b.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, b, a, d, c,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, b, a, d, c,  -(L+L+1) / 2 * xcc) )   end
            # mu = 2
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L) * (c.subshell.kappa - a.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, a, b, c, d,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, a, b, c, d,  -(L+L+1) / 2 * xcc) )   end
            # mu = 3
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L + 1) * (c.subshell.kappa - a.subshell.kappa - L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, d, c, b, a,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, d, c, b, a,  -(L+L+1) / 2 * xcc) )   end
            # mu = 4
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa - L) * (c.subshell.kappa - a.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, c, d, a, b,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, c, d, a, b,  -(L+L+1) / 2 * xcc) )   end
            # mu = 5
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L + 1) * (c.subshell.kappa - a.subshell.kappa + L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, d, a, b, c,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, d, a, b, c,  -(L+L+1) / 2 * xcc) )   end
            # mu = 6
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa - L) * (c.subshell.kappa - a.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, a, d, c, b,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, a, d, c, b,  -(L+L+1) / 2 * xcc) )   end
            # mu = 7
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa - L - 1) * (c.subshell.kappa - a.subshell.kappa - L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, b, c, d, a,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, b, c, d, a,  -(L+L+1) / 2 * xcc) )   end
            # mu = 8
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L) * (c.subshell.kappa - a.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, c, b, a, d,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, c, b, a, d,  -(L+L+1) / 2 * xcc) )   end
        end

        return( xcList )
    end

    """
    `InteractionStrength.XL_Breit0_densities(xcList::Array{XLCoefficient,1}, grid::Radial.Grid)`
        ... computes the the effective Breit interaction strengths X^L,0_Breit (abcd) for given rank L and a list of
            orbital functions a, b, c, d and angular coefficients at the given grid. A value::Float64 is returned.
            At present, only the zero-frequency Breit interaction is taken into account.
    """
    function XL_Breit0_densities(xcList::Array{XLCoefficient,1}, grid::Radial.Grid)

        if  grid.mesh == MeshGL
            wa = 0.
            for  xc  in  xcList  ## [end:end]
                # Use the minimal extent of any involved orbitals; this need to be improved
                mtp_ac = min(size(xc.a.P, 1), size(xc.c.P, 1));   mtp_bd = min(size(xc.b.P, 1), size(xc.d.P, 1))
                for  r = 2:mtp_ac
                    for  s = 2:mtp_bd
                        wc = xc.coeff * grid.wr[r] * grid.wr[s]
                        if      s > r   continue
                        elseif  s == r
                            wa = wa + wc * (xc.a.P[r] * xc.c.Q[r]) * (grid.r[s]^xc.nu) / (grid.r[r]^(xc.nu+1)) * (xc.b.P[s] * xc.d.Q[s]) / 2.0
                        else
                            wa = wa + wc * (xc.a.P[r] * xc.c.Q[r]) * (grid.r[s]^xc.nu) / (grid.r[r]^(xc.nu+1)) * (xc.b.P[s] * xc.d.Q[s])
                        end
                    end
                end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `InteractionStrength.XL_Coulomb(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the the effective Coulomb interaction strengths X^L_Coulomb (abcd) for given rank L and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function XL_Coulomb(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( 0. )
        end
        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc    end

        ##x za = RadialIntegrals.SlaterRk_2dim(L, a, b, c, d, grid);   zb = RadialIntegrals.SlaterRk_new(L, a, b, c, d, grid)
        ##x if abs( (za-zb)/za ) > 1.0e-12    println("XL_Coulomb: Slater  za = $za   zb = $zb  ")   end

        XL_Coulomb = xc * RadialIntegrals.SlaterRk_2dim(L, a, b, c, d, grid)
        ##x XL_Coulomb = xc * RadialIntegrals.SlaterRk_new(L, a, b, c, d, grid)
        return( XL_Coulomb )
    end

    """
    `InteractionStrength.XL_Coulomb_DH(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid, lambda::Float64)`
        ... computes the the effective Coulomb-Debye-Hückel interaction strengths X^L_Coulomb_DH (abcd) for given rank L and
            orbital functions a, b, c and d at the given grid and for the given screening parameter lambda. A value::Float64 is
            returned.
    """
    function XL_Coulomb_DH(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid, lambda::Float64)
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( 0. )
        end
        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc    end

        XL_Coulomb_DH = xc * RadialIntegrals.SlaterRk_DebyeHueckel_2dim(L, a, b, c, d, grid, lambda)
        return( XL_Coulomb_DH )
    end

    """
    `InteractionStrength.XL_plasma_ionSphere(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, lambda::Float64)`
        ... computes the effective interaction strengths X^L_ion-sphere (abcd) for given rank L and orbital functions
            a, b, c and d and for the plasma parameter lambda. A value::Float64 is returned.  **Not yet implemented !**
    """
    function XL_plasma_ionSphere(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, lambda::Float64)
        error("Not yet implemented")
    end

    """
    `InteractionStrength.X1_smsA(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the the effective interaction strengths X^1_sms,A (abcd) for fixed rank 1 and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function X1_smsA(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        wa = AngularMomentum.CL_reduced_me(a.subshell, 1, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, 1, d.subshell) *
             RadialIntegrals.Vinti(a, c, grid) * RadialIntegrals.Vinti(b, d, grid) / (2 * nm.mass)
        ## println("**  <$(a.subshell) || Vinti || $(c.subshell)>  = $(RadialIntegrals.Vinti(a, c, grid)) " )
        ## println("**  <$(b.subshell) || Vinti || $(d.subshell)>  = $(RadialIntegrals.Vinti(b, d, grid)) " )
        return( wa )
    end

    """
    `InteractionStrength.X1_smsB(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the the effective interaction strengths X^1_sms,B (abcd) for fixed rank 1 and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function X1_smsB(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        wa = - AngularMomentum.CL_reduced_me(b.subshell, 1, d.subshell) * RadialIntegrals.Vinti(b, d, grid) *
               RadialIntegrals.isotope_smsB(a, c, nm.Z, grid) / (2 * nm.mass)
        return( wa )
    end

    """
    `InteractionStrength.X1_smsC(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the the effective interaction strengths X^1_sms,C (abcd) for fixed rank 1 and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function X1_smsC(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        wa = - AngularMomentum.CL_reduced_me(b.subshell, 1, d.subshell) * AngularMomentum.CL_reduced_me(a.subshell, 1, c.subshell) *
               RadialIntegrals.Vinti(b, d, grid) * RadialIntegrals.isotope_smsC(a, c, nm.Z, grid) / (2 * nm.mass)
        return( wa )
    end

    """
    `InteractionStrength.zeeman_Delta_n1(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| Delta n^(1) ||b> reduced matrix element for the Zeeman-Schwinger contribution to the coupling
            to an external magnetic field for orbital functions a, b. A value::Float64 is returned.
    """
    function zeeman_Delta_n1(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC ... test for the proper definition of the C^L tensors.
        minusa = Subshell(1, -a.subshell.kappa)
        wa = (Defaults.getDefaults("electron g-factor") - 2) / 2. * (a.subshell.kappa + b.subshell.kappa + 1) *
             AngularMomentum.CL_reduced_me(minusa, 1, b.subshell) * RadialIntegrals.rkDiagonal(0, a, b, grid)
        return( wa )
    end

    """
    `InteractionStrength.zeeman_n1(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| n^(1) ||b> reduced matrix element for the Zeeman coupling to an external magnetic field for
            orbital functions a, b. A value::Float64 is returned.
    """
    function zeeman_n1(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC ... test for the proper definition of the C^L tensors.
        minusa = Subshell(1, -a.subshell.kappa)
        wa = - AngularMomentum.CL_reduced_me(minusa, 1, b.subshell) / (2 * Defaults.getDefaults("alpha")) *
               RadialIntegrals.rkNonDiagonal(1, a, b, grid)
        return( wa )
    end

end # module


text:
include("definitions.jl")

function readSparseData(sparseDataFilename::AbstractString, dictionaryFilename::AbstractString)
#read a dataset in sparse format. One line per document, consisting of pairs of token index and token count, space separated.
    f = open(dictionaryFilename);
        dictionary = readlines(f); #an array of strings
    close(f);

    f = open(sparseDataFilename);
    documents = Array{SparseDocument}(0);
    for ln in eachline(f)
        splitLn = split(ln);
        wordCounts = Array{Int64}(length(splitLn));
        docLength = 0;
        numTokensFound = convert(Int64,floor(length(splitLn) / 2));
        tokenIndex = zeros(Int64, numTokensFound);
        tokenCount = zeros(Int64, numTokensFound);
        ind = 1;
        for i = 1:2:length(splitLn)
           tokenIndex[ind] = parse(Int64, splitLn[i]);
           tokenCount[ind] = parse(Int64, splitLn[i + 1]);
           docLength = docLength + tokenCount[ind];
           ind = ind + 1;
        end
        push!(documents, SparseDocument([], tokenIndex, tokenCount, docLength));
    end
    close(f);

    documents, dictionary
end


text:
export sumt, prodt, xlogx, randin

"safe typed sum. Returns correct zero of correct type when iterator is empty."
function sumt(::Type{T}, itr)::T where T
    isempty(itr) ? zero(T) : sum(x::T for x in itr)::T
end

"safe typed prod. Returns correct one of correct type when iterator is empty."
function prodt(::Type{T}, itr)::T where T
    isempty(itr) ? one(T) : prod(x::T for x in itr)::T
end

"x*log(x) that gives zero when x = 0. In the future use StatsFuns.xlogx instead"
xlogx(x::Float64) = (warn("Use StatsFuns.xlogy instead"); xlogy(x,x))

"x*log(y) that gives zero when x = y = 0. In the future use StatsFuns.xlogy instead"
function xlogy(x::Float64,y::Float64)
    y ≥ 0 || DomainError()
    x == 0 ? 0. : x * log(y)
end

"Volume of a simplex with vertices in V"
function simplexvol(V::Vector{Vector{Float64}})
    n = length(V) - 1
    @assert all(length.(V) .== n)
    abs(det(hcat((V[i] - V[end] for i = 1:n)...))) / factorial(n)
end

"log of binomial coefficient"
function lbinomial(i::Real, j::Real)
    0 ≤ i < Inf && 0 ≤ j < Inf || throw(DomainError())
    lgamma(i + 1) - lgamma(j + 1) - lgamma(i - j + 1)
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9bc349e8-efb0-11ea-1392-cd735828af33
# packages to load etc
begin
	using DelimitedFiles
	using NIfTI
	using Images
	using ImageView
	using LinearAlgebra
	using StatsBase
	using GLM
	using DataFrames
	using Plots
	using Distributions
	using StatsPlots
	using PlutoUI
	using Pipe
	using Printf
	md"""## loading packages...

	The notebook makes use of some existing functionality from other code which is stored in "packages" that you can install and use.

	"""
end

# ╔═╡ a04710a0-281b-11eb-3c11-5b53a085e59e
html"<button onclick=present()>Present</button>"

# ╔═╡ 8029070e-efb0-11ea-1d0a-7f402f2e9989
md"""# What are `spm`s

Statistical parametric maps... are maps / images that contain statistical values in the same x,y,z as the original image data...

They can be images of the β weights, or other statistics that are derived from them, such as

- the coefficient of determintation, $r^{2}$
- a t-statistic,  $t$, (like from a $t$-test you already know
- statistical significance, $p$-values, etc.

"""

# ╔═╡ 44032c62-2522-11eb-32d2-af427f6ce946
md"""## Doing linear regression

The linear regression analysis can be done very neatly using *linear algebra*. The details of this are beyond this 2nd year lab class, but just to give you an idea of what things look like:

$$\large \text{{\color{teal} \textbf{data}}} = \text{\bf \color{magenta} explanatory variables} \cdot \mathbf{w} $$

$$\large \mathbf{\color{teal} y} = \mathbf{\color{magenta} X}  \cdot w$$

In this equation $\mathbf{w}$, the weights, are what the want to determine. We have measured data -- and by designing the experiment in a particular way (what stimuli when, etc), we have determined the explanatory variables.

To find the weights we need to solve this equation for $\mathbf{w}$:

$$\begin{align}
y &= X \color{green}\mathbf{w}
\end{align}$$

"""

# ╔═╡ cf9fbcc4-2534-11eb-1a82-ad480057d2b0
d, h= readdlm("data/stimulus_timing.csv", ',', Float64, header = true );

# ╔═╡ 3df12514-2535-11eb-2cf3-fdfc727a993b
X = hcat(d[:,3:4], ones(size(d,1),1)) ;

# ╔═╡ b7f96162-25af-11eb-274a-7d368d09906b
md"""
## Using the `lm()` linear model function

*linear models* and regression are some of the most fundamental techniques used in stats, so they are implemented in many software tools, including `SPSS`, and packages for `R`, `Matlab`, `python`, `julia` and many others (as we have seen even in `Excel`)

Across all those tools, the models are often presented in terms of **formulas**, e.g.

  `y ~ face_response + object_response`

which stands for a model that tries to predict the data `y` as a linear combination of two explanatory variables, `face_response` and `object_response`

$$\begin{align}
\color{green}\mathbf{w} \color{black} &= X \backslash y
\end{align}$$

or in code
```julia
w = \(X,ȳ)
# or
w = pinv(X)*ȳ
```
"""

# ╔═╡ a66ed758-25bc-11eb-29ad-01f9ba1bd5be

# ╔═╡ dcb3ad64-281b-11eb-2a9c-5d9cfce9e2e2
md"""## Model:

The measured eesponse is combination of `F` and `O` response"""

# ╔═╡ ba55a268-281b-11eb-3bcd-5b948d03fdbf
md"""## Data, model, residuals..."""

# ╔═╡ dfb188e0-2836-11eb-0081-514c967cace5
md"""## How significant is the response

Turn the beta weight into a $t$ statistic and from there into a $p$ value

$$\begin{align}
t  &= \frac{\text{mean}}{\text{standard error}}
\end{align}$$

Then we need to figure out whether that  $t$-value we calculated is far enough away from 0 to be "significant".

"""

# ╔═╡ 1a0bfe62-2837-11eb-02e3-45eff0ed23f6
t_value = 5.29832/0.229458

# ╔═╡ 39ad1ca2-284a-11eb-3c64-e5e1317eb3d4
md"""$(@bind deg_freedom Slider(2:1:20; default=5, show_value=true)) """

# ╔═╡ e86ecc6e-2849-11eb-34da-0189163c776d
begin
	xx = -5.0:0.2:5.0
	plot(Distributions.TDist(deg_freedom), fill=(0, .5,:teal),
	size=(600,200), title="t Distribution (ν=$(deg_freedom))", label=false)
end

# ╔═╡ 148b4e22-25ba-11eb-28ac-4f90284e7e41
md""" ## Visualising many weights """

# ╔═╡ 6b22738e-281f-11eb-0be4-3195a4953529
md"""## And the `t` values

"""

# ╔═╡ c886d94e-25bc-11eb-0b39-8144659ec58f
md"""## Other versions of the model"""

# ╔═╡ f543c9da-281c-11eb-05bd-ad64b4005507
md"""## Only modelling the `F` response"""

# ╔═╡ 0a4f3af8-281d-11eb-3337-8fbb9998237c
md"""## Only modelling the `O` response"""

# ╔═╡ 41d2d3cc-281d-11eb-176a-c9230c75163c
md"""## Some timing / performance measurements"""

# ╔═╡ 2b2f6be4-f07f-11ea-1459-775749ca949f
md"""## Some helper functions """

# ╔═╡ 46fbc1de-13d1-11eb-2553-71dfd823a5a0
md"""...loading in the contents of a file into a **variable** called `data` works as follows:"""

# ╔═╡ cdf28f06-1449-11eb-34b1-65a1fb83f330
md"""
## Displaying images

To look at the images, we can stich them together into one big *mosaic* and turn them into *grayscale* image that the browswer can display. To do this, we need to find out what the largest (brightest pixel value) and smallest number in the data we want to display are. Because the `max` value might be an outlier, we'll use the 1st and 99th *percentile* instead -- this is often called a **robust range**.

The function `scaleFunc` then squeezes the numbers into a range $[0..1]$ which works well for displaying.
"""

# ╔═╡ 5aed662e-2525-11eb-12a3-c9c9118f6fec
demean = function(y)
	100 .* (y./mean(y) .- 1.0)
end

# ╔═╡ 88253376-1484-11eb-0967-27875c6da199
md"""because this is something we might want to do again and again, you can turn it into a **function**, which is something that takes an  *input*, changes it and spits it out again as an *output*
"""

# ╔═╡ 141dde62-efb1-11ea-34df-cd5e838d8129
begin
	# read data in from file
	data = NIfTI.niread("data/filtered_func_data.nii");
	# small technical point here: images are often stored in y,x,z
	# order (rather thatn x,y,z... so could permute data here (mutable!)

	md"""
	Succesfully loaded a dataset that has **$(ndims(data)) dimensions** of size $(size(data))

	If we look at the table of numbers for one slice (slice #12) for one timepoint (volume = 3 in the below), it looks a bit like the following:
	"""
end

# ╔═╡ f0ef494e-f80c-11ea-1163-b954c1739162
view(data, 1:10:64, 1:10:64, 12, 3);

# ╔═╡ 14f85bc0-2533-11eb-065b-93a0cccb036a
begin
	theSz = [ size(data)[1:3]... ; 3];
	β̂ = zeros(theSz...);
	tmap = zeros(theSz...);

	@time for ix = 1:theSz[1], iy = 1:theSz[2], iz = 1:theSz[3]
		# if mod(ii, 1000) == 1
		# 	println("tada")
		# end
		ols = lm(X, Float64.(demean(data[ix,iy,iz,:])))
		β̂[ix,iy,iz,:] = coef(ols)
		tmap[ix, iy, iz, :] = 	coeftable(ols).cols[3] # t values
	end

end

# ╔═╡ a8ab8ee6-f083-11ea-1316-7fbe4e125914
begin
md"""

### Calculate the beta **weights** at a voxel

We can pick the timeseries data at one particular location and see how long it takes to do that....

[plot residuals] $(@bind plotResiduals CheckBox()) – [plot model] $(@bind plotModel CheckBox())

x coord $(@bind xLoc Slider(1:size(data,1), default = 18, show_value=true))

y coord $(@bind yLoc Slider(1:size(data,2), default = 14, show_value=true))

z coord $(@bind zLoc Slider(1:size(data,3), default = 4, show_value=true))

"""
end

# ╔═╡ 04256b16-2525-11eb-30b3-eb750b5af1a2
y = data[xLoc,yLoc,zLoc,:];

# ╔═╡ 3c5138dc-25b9-11eb-058d-71a10f974857
ȳ = demean(y);

# ╔═╡ 42ace0ae-2538-11eb-3143-31d86905a747
md"""

slice: $(@bind whichSlice Slider( 1:size(data,3) ; default=1, show_value=true))

threshold: $(@bind whichThreshold Slider( 0:0.5:15 ; default=2.0, show_value=true))

"""

# ╔═╡ d0b91c60-2537-11eb-2ff2-eda55e302239
begin
	b1 = β̂[:,:,whichSlice,1]
	b2 = β̂[:,:,whichSlice,2]

	b1[abs.(b1) .<= whichThreshold] .= NaN
	b2[abs.(b2) .<= whichThreshold] .= NaN

	b3 = b1 .- b2

	heatmap(b1, color=:reds, clim = (0,5), layout = (1,2),
			title = "weights map 1",
			subplot=1,
			aspect_ratio=1, showaxis=false)
	heatmap!(b2, color=:blues, clim = (0,5), layout = (1,2),
			title = "weights map 2",
			subplot=2,
			aspect_ratio=1, showaxis=false)

		# heatmap!(b3, clim = (-15,15), layout = (1,3),
		# 	title = "w1 - w2 (f > o)",
		# 	subplot=3,
		# 	aspect_ratio=1, showaxis=false)
end

# ╔═╡ 86dbd432-281f-11eb-1aab-b53d4a226dd0
md"""

slice (t): $(@bind whichSliceT Slider( 1:size(data,3) ; default=1, show_value=true))

threshold (for t): $(@bind whichThresholdT Slider( 0:0.5:15 ; default=2.0, show_value=true))

"""

# ╔═╡ 757b4a0e-281f-11eb-32ad-cfff43e2ebe3
begin
	t1 = tmap[:,:,whichSliceT,1]
	t2 = tmap[:,:,whichSliceT,2]

	t1[abs.(t1) .<= whichThresholdT] .= NaN
	t2[abs.(t2) .<= whichThresholdT] .= NaN

	t3 = t1 .- t2
	t4 = t2 .- t1

	# heatmap(t1, color=:reds, clim = (0,5), layout = (1,1),
	# 		title = "t map 1",
	# 		subplot=1,
	# 		aspect_ratio=1, showaxis=false)
	# heatmap!(t2, color=:blues, clim = (0,5), layout = (1,1),
	# 		title = "t map 2",
	# 		subplot=2,
	# 		aspect_ratio=1, showaxis=false)
	heatmap(t3, color=cgrad([:blue, :white, :red]), clim = (-5,5), layout = (1,1),
			title = "t map (F - O)",
			subplot=1,
			aspect_ratio=1, showaxis=false)

end

# ╔═╡ 55f99666-25ba-11eb-1fd4-a15c6751527e
df = DataFrame(y = Float64.(demean(data[xLoc,yLoc, zLoc,:])),
				face_r = X[:,1], object_r = X[:,2]);

# ╔═╡ 3bf4e998-25a4-11eb-252b-6fdb6a29a92e
ols = lm(@formula(y ~ face_r + object_r), df)

# ╔═╡ 9af3adba-25ba-11eb-34a7-59079afd5231
md"""## How well does the model fit the data?

The coefficient of determination, $r^2$, tells us how much of the variance in the data is accounted for by our model. $r^2$ = $(round( r2(ols); digits=2)) means $( round( 100*r2(ols); digits=0)) % of the variance.

It can be calculated as

$$r^{2}= 1 - \frac{\text{var(residuals)}}{\text{var(data)}}$$

A value of $r^{2} = 1.0$  means that our model is perfect.

A value of $r^{2} \approx 0$ means that our model can't explain our data at all.

"""

# ╔═╡ ab4cc236-25b1-11eb-3f63-233af5cfea96
r² = r2(ols)

# ╔═╡ 10a0bd36-2837-11eb-2fe7-69a17334a3d9
ols

# ╔═╡ 821ffcee-2849-11eb-340d-2dc707fd0a8f
md""" ## But what about a $p$-value

The p-value corresponds to the area under the curve beyond the calculated t-value... i our example for `face_r` it was t = $(round(coeftable(ols).cols[3][2]; digits=2)) and a p-value of $(@sprintf "%.2g" coeftable(ols).cols[4][2])

"""

# ╔═╡ 765e7432-25a4-11eb-10b5-d3bdab410d87
olsFvO = lm(@formula(y ~ face_r - object_r), df)

# ╔═╡ 570626cc-25a5-11eb-123c-63c5b9aad411
olsF = lm(@formula(y ~ face_r), df)

# ╔═╡ 76cdcb54-25bb-11eb-0e29-9dfeae3f803d
r2(olsF)

# ╔═╡ 9fa22482-25a5-11eb-2b76-1949341e2630
olsO = lm(@formula(y ~ object_r), df)

# ╔═╡ 894cdc18-25bb-11eb-130d-93b864ddcc42
r2(olsO)

# ╔═╡ 1ef8659c-2539-11eb-184a-d95f169892a4
result = @timed @pipe data[xLoc,yLoc,zLoc,:] |>
	 					demean( _ ) |>
	 					vec( _ ) |>
	 					X \ _

# ╔═╡ 3483ee2c-2525-11eb-2967-d343eb3d16f0
begin

	p_ = plot(ȳ, lw=2, c=:black, label="data $([xLoc, yLoc, zLoc])", size=(600,200),
			xlabel="Time (volumes)", ylabel="fMRI response",
			title="r² = $(round(r²; digits=2))")

	modelFit = modelFit = X*result.value
	residuals = ȳ .- modelFit

	if plotModel
		plot!(modelFit, c=:red, lw=2.0, label="model fit")
	end

	if plotResiduals
		plot!(residuals, line=:stem, c=:blue, label="redisuals")
	end

	p_
end

# ╔═╡ 15ad8bf8-2835-11eb-1315-f5d8ea0d2bf9
md"""The *weights* for reconstructing are: **$(@sprintf "[%.2f. %.2f, %.2f]" result.value...)**"""

# ╔═╡ a48038d4-2539-11eb-3374-dbbeaed2a3a5
result.value

# ╔═╡ b3a74a5a-2539-11eb-017a-f9b3ffe38663
r_microseconds = result.time * 1000 * 1000; # in microseconds

# ╔═╡ 3e13431a-253a-11eb-265e-f11984b2b2fb
md"""The following expression feeds the data from [$(xLoc); $(yLoc); $(zLoc)] to a function called `demean()` and then to a function called `vec()` and finally does the maths using the `\` (backslash) operator:

Time taken: ≈ $(round(r_microseconds; digits=2)) microseconds

weights: $(@sprintf("[%.2f, %.2f, %.2f]", round.(result.value; digits=2)...))

"""

# ╔═╡ 13a8fe64-efb4-11ea-1ac6-dbd49975dddf
minMax = percentile(data[:], [1; 99])

# ╔═╡ 71606d26-efb4-11ea-12e2-75f85b9ea7d8
# helper function that knows how to scale images to a range that we can display
scaleFunc = scaleminmax(minMax[1], minMax[2])

# ╔═╡ c1d831e0-1448-11eb-0b18-a555cbcb18e1
@pipe data[:,:,:,3] |>
	mosaicview( _ , nrow=4; rowmajor=true) |>
	scaleFunc.( _ ) |>
	Gray.( _ )

# ╔═╡ b69a53ee-1484-11eb-0953-2b3ff81cf56c
function simpleMosaic(d, t=3, nr = 4)
	Gray.(scaleFunc.(mosaicview(d[:,:,:,t], nrow=nr; rowmajor=true)))
end

# ╔═╡ e2bea646-1484-11eb-3a26-ffc655b139c6
simpleMosaic(data) # now we can use the code this way

# ╔═╡ 072ce484-efb2-11ea-36ee-29a0147cfc04
"""
displaySlice(data, slicenum, orientation)

return and display a slice in a particular orientation. there are smarter
ways to implement the slicing through different orientations, but the code
using if/elseif is easy to folllow.

	displaySlice(data, 2, 1)

"""
displaySlice = function(e, s = 1, orientation = 3, scaling = scaleFunc)
	if  (orientation == 1)
		displayData = view(e, s, :,:,1)
	elseif (orientation == 2)
		displayData = view(e, :,s,:,1)
	elseif( orientation == 3)
		displayData = view(e, :, :, s, 1)
	else
		error("uhoh!")
	end
	rotl90(Gray.(scaling.(displayData)))
end

# ╔═╡ 8cbe6b98-144b-11eb-3d0a-35700773e22b
"""
showMosaic - arrange 4d data (at chosen timepoint) into a mosaic

    showMosaic(data,2)  # permute the dimensions (circshift twice)
						# and then arrange into roughly nxn mosaic

"""
function showMosaic(d, dim, t=1)
	@pipe d[:,:,:,t] |>
		permutedims(_, circshift([1,2,3],dim)) |>
		mosaicview(_,
			nrow=Int64(round(sqrt(size(_,3)...))),
			rowmajor=true,
			npad = 2,
			fillvalue = maximum(_)) |>
		scaleFunc.(_) |>
		Gray.(_)
end

# ╔═╡ Cell order:
# ╠═a04710a0-281b-11eb-3c11-5b53a085e59e
# ╟─8029070e-efb0-11ea-1d0a-7f402f2e9989
# ╟─f0ef494e-f80c-11ea-1163-b954c1739162
# ╟─44032c62-2522-11eb-32d2-af427f6ce946
# ╟─04256b16-2525-11eb-30b3-eb750b5af1a2
# ╟─3c5138dc-25b9-11eb-058d-71a10f974857
# ╟─cf9fbcc4-2534-11eb-1a82-ad480057d2b0
# ╟─3df12514-2535-11eb-2cf3-fdfc727a993b
# ╟─14f85bc0-2533-11eb-065b-93a0cccb036a
# ╟─b7f96162-25af-11eb-274a-7d368d09906b
# ╟─a66ed758-25bc-11eb-29ad-01f9ba1bd5be
# ╟─dcb3ad64-281b-11eb-2a9c-5d9cfce9e2e2
# ╠═3bf4e998-25a4-11eb-252b-6fdb6a29a92e
# ╟─ba55a268-281b-11eb-3bcd-5b948d03fdbf
# ╟─3483ee2c-2525-11eb-2967-d343eb3d16f0
# ╟─15ad8bf8-2835-11eb-1315-f5d8ea0d2bf9
# ╟─a8ab8ee6-f083-11ea-1316-7fbe4e125914
# ╟─9af3adba-25ba-11eb-34a7-59079afd5231
# ╠═ab4cc236-25b1-11eb-3f63-233af5cfea96
# ╟─dfb188e0-2836-11eb-0081-514c967cace5
# ╠═1a0bfe62-2837-11eb-02e3-45eff0ed23f6
# ╠═10a0bd36-2837-11eb-2fe7-69a17334a3d9
# ╟─821ffcee-2849-11eb-340d-2dc707fd0a8f
# ╠═e86ecc6e-2849-11eb-34da-0189163c776d
# ╟─39ad1ca2-284a-11eb-3c64-e5e1317eb3d4
# ╟─148b4e22-25ba-11eb-28ac-4f90284e7e41
# ╟─d0b91c60-2537-11eb-2ff2-eda55e302239
# ╟─42ace0ae-2538-11eb-3143-31d86905a747
# ╟─6b22738e-281f-11eb-0be4-3195a4953529
# ╟─757b4a0e-281f-11eb-32ad-cfff43e2ebe3
# ╟─86dbd432-281f-11eb-1aab-b53d4a226dd0
# ╟─c886d94e-25bc-11eb-0b39-8144659ec58f
# ╠═765e7432-25a4-11eb-10b5-d3bdab410d87
# ╟─f543c9da-281c-11eb-05bd-ad64b4005507
# ╠═570626cc-25a5-11eb-123c-63c5b9aad411
# ╠═76cdcb54-25bb-11eb-0e29-9dfeae3f803d
# ╟─0a4f3af8-281d-11eb-3337-8fbb9998237c
# ╠═9fa22482-25a5-11eb-2b76-1949341e2630
# ╠═894cdc18-25bb-11eb-130d-93b864ddcc42
# ╟─55f99666-25ba-11eb-1fd4-a15c6751527e
# ╟─41d2d3cc-281d-11eb-176a-c9230c75163c
# ╟─3e13431a-253a-11eb-265e-f11984b2b2fb
# ╠═1ef8659c-2539-11eb-184a-d95f169892a4
# ╠═a48038d4-2539-11eb-3374-dbbeaed2a3a5
# ╠═b3a74a5a-2539-11eb-017a-f9b3ffe38663
# ╟─2b2f6be4-f07f-11ea-1459-775749ca949f
# ╟─46fbc1de-13d1-11eb-2553-71dfd823a5a0
# ╟─cdf28f06-1449-11eb-34b1-65a1fb83f330
# ╟─c1d831e0-1448-11eb-0b18-a555cbcb18e1
# ╟─e2bea646-1484-11eb-3a26-ffc655b139c6
# ╠═9bc349e8-efb0-11ea-1392-cd735828af33
# ╠═13a8fe64-efb4-11ea-1ac6-dbd49975dddf
# ╟─71606d26-efb4-11ea-12e2-75f85b9ea7d8
# ╟─b69a53ee-1484-11eb-0953-2b3ff81cf56c
# ╟─5aed662e-2525-11eb-12a3-c9c9118f6fec
# ╟─88253376-1484-11eb-0967-27875c6da199
# ╟─141dde62-efb1-11ea-34df-cd5e838d8129
# ╟─072ce484-efb2-11ea-36ee-29a0147cfc04
# ╟─8cbe6b98-144b-11eb-3d0a-35700773e22b

max_stars_repo_name: schluppeck/hands-on-brain-data

max_issues_repo_name: schluppeck/hands-on-brain-data

max_forks_repo_name: schluppeck/hands-on-brain-data


text:
# Interface to define custom state types

"""
    posf(state)
returns the coordinates of the state in the Frenet frame.
The return type is expected to be `Frenet`.
"""
function posf end

"""
    posg(state)
returns the coordinates of the state in the global (world) frame.
The return type is expected to be a VecSE2.
"""
function posg end

"""
    vel(state)
returns the norm of the longitudinal velocity.
"""
function vel end

"""
    velf(state)
returns the velocity of the state in the Frenet frame.
"""
function velf end

"""
    velg(state)
returns the velocity of the state in the global (world) frame.
"""
function velg end


text:
module HNSW
    using LinearAlgebra
    using Reexport
    @reexport using Distances

    export HierarchicalNSW
    export add_to_graph!, set_ef!
    export knn_search

    include("neighborset.jl")
    include("visited_lists.jl")
    include("layered_graphs.jl")
    include("interface.jl")
    include("algorithm.jl")
end # module