text:
import Base.REPL: outstream
export FullScreenDialog

type InlineDialog
    w::Widget
    tty
    was_invalidated::Bool
    t
    do_redraw
    InlineDialog(w::Widget,tty) = new(w,tty,false)
end

function enable_settings(tty; query_cursor = false)
    Base.reseteof(tty)
    Base.Terminals.raw!(tty,true)
    write(tty,CSI,"?25l")
    query_cursor && write(tty,CSI,"?6n")
    write(tty,CSI,"?1002h")
    write(tty,CSI,"?1005h")
end

function disable_settings(tty)
    Base.Terminals.raw!(tty,false)
    write(tty,CSI,"?25h")
    write(tty,CSI,"?1002l")
    write(tty,CSI,"?1005l")
end

type FullScreenDialog
    w::Widget
    tty
    was_invalidated::Bool
    t
    do_redraw
    FullScreenDialog(w::Widget,tty) = new(w,tty,false)
end
invalidate(dialog::Union{InlineDialog, FullScreenDialog}) = dialog.was_invalidated = true

function create_input_loop(i,tty,s)
    @schedule begin
        try
            result = :ok
            # eof needs to go second here, because it may block
            while result != :done && !eof(tty)
                result = Base.LineEdit.match_input(i.w.ctx.focuss.keymap,nothing,tty)(s,nothing)
                i.was_invalidated && i.do_redraw()
            end
        catch e
            if !isa(e,InterruptException)
                rethrow(e)
            end
        end
    end
end

immutable IOBufferTerminal
    buf::IOBuffer
end
size(tty::IOBufferTerminal) = (80,24)
Base.displaysize(tty::IOBufferTerminal) = size(tty)
width(tty::IOBufferTerminal) = size(tty)[1]
height(tty::IOBufferTerminal) = size(tty)[2]
Base.write(tty::IOBufferTerminal,args...) = write(buf,args...)

function debug(w::Widget, emdebug = false)
    buf = IOBuffer()
    i = InlineDialog(w, IOBufferTerminal(buf))
    s = setup_screen(i)
    s.offset = (1,1)
    s.fullsize = Base.displaysize(i.tty)
    focus(i.w)
    redraw(s, i.w)
    render(buf, s)
    em = VT100.ScreenEmulator()
    em.debug = emdebug
    seekstart(buf)
    VT100.parseall!(em, buf)
    em
end

function set_widget_dialog!(w, i)
    !isdefined(w, :ctx) && return
    w.ctx.dialog = i
    for nw in children(w)
        set_widget_dialog!(nw, i)
    end
end

function setup_screen(i)
    full = isa(i,FullScreenDialog)
    TerminalUI.initialize!(i.w)
    set_widget_dialog!(i.w, i)
    local s
    if full
        s = TerminalUI.DoubleBufferedTerminalScreen(Base.displaysize(i.tty))
    else
        inline_height = min(optheight(i.w),height(i.tty)-4)
        s = TerminalUI.DoubleBufferedTerminalScreen((inline_height,width(i.tty)))
    end
    s
end

function wait(i::Union{FullScreenDialog,InlineDialog})
    full = isa(i,FullScreenDialog)
    local s = setup_screen(i)
    t = i.t = create_input_loop(i,i.tty,s)
    try
        i.do_redraw = function(args...)
            redraw(s,i.w)
            render(i.tty.out_stream,s)
        end
        enable_settings(i.tty; query_cursor = !full)
        focus(i.w)
        if !full
            println("Waiting for cursor position")
            wait(TerminalUI.curspos_condition)
            println("Done")
            pos = TerminalUI.curspos
            # Leave 1 line between the prompt and the dialog, otherwise it looks
            # squished.
            s.offset = (pos[1]+1,pos[2])
            s.fullsize = Base.displaysize(i.tty)
        end
        map(size->(resized(s,size); i.do_redraw()),monitor_resize(i.tty))
        map((args...)->render(i.tty.out_stream,s),TerminalUI.invalidated2)
        wait(t)
        afterembed(i.tty.out_stream,s)
    catch e
        bt = catch_backtrace()
        @show e
        Base.show_backtrace(STDERR, bt)
        rethrow(e)
    finally
        disable_settings(i.tty)
        close(i)
    end
end

function change_widget!(i::Union{FullScreenDialog,InlineDialog}, w)
    i.w = w
    TerminalUI.initialize!(i.w)
    set_widget_dialog!(i.w, i)
    focus(i.w)
    invalidate(i.w)
    i.do_redraw()
end

function close(d::Union{FullScreenDialog,InlineDialog})
    if !istaskdone(d.t)
        Base.throwto(d.t, InterruptException())
    end
end

function print_snapshot(i)
    full = isa(i,FullScreenDialog)
    local s = setup_screen(i)
    t = i.t = create_input_loop(i,i.tty,s)
    try
        enable_settings(i.tty; query_cursor = !full)
        if !full
            wait(TerminalUI.curspos_condition)
            pos = TerminalUI.curspos
            # Leave 1 line between the prompt and the dialog, otherwise it looks
            # squished.
            s.offset = (pos[1]+1,pos[2])
            s.fullsize = Base.displaysize(i.tty)
        end
        focus(i.w)
        redraw(s,i.w)
        render(i.tty.out_stream,s)
    finally
        disable_settings(i.tty)
        close(i)
    end
end


text:
module AlgebraCholeskyTest

using Test
using ReactiveMP
using Random

using LinearAlgebra

@testset "Helpers" begin

    @testset "cholesky related" begin
        rng = MersenneTwister(1234)

        for size in (2, 3, 4, 5, 10, 100, 1000)
            L = rand(rng, size, size)
            A = L * L'

            @test LinearAlgebra.cholesky(A).L ≈ ReactiveMP.fastcholesky(A).L
            @test inv(LinearAlgebra.cholesky(A)) ≈ inv(ReactiveMP.fastcholesky(A))
            @test inv(LinearAlgebra.cholesky(A)) ≈ ReactiveMP.cholinv(A)
            @test sqrt(A) * sqrt(A)' ≈ ReactiveMP.cholsqrt(A) * ReactiveMP.cholsqrt(A)'
            @test A ≈ ReactiveMP.cholsqrt(A) * ReactiveMP.cholsqrt(A)'
            @test logdet(A) ≈ ReactiveMP.chollogdet(A)
            @test all((inv(cholesky(A)), logdet(A)) .≈ ReactiveMP.cholinv_logdet(A))
        end

        for size in (2, 3, 4, 5, 10, 100, 1000)
            A = Diagonal(rand(rng, size))
            @test inv(A) ≈ ReactiveMP.cholinv(A)
            @test (sqrt(A) * sqrt(A)') ≈ (ReactiveMP.cholsqrt(A) * ReactiveMP.cholsqrt(A)')
            @test A ≈ (ReactiveMP.cholsqrt(A) * ReactiveMP.cholsqrt(A)')
            @test logdet(A) ≈ ReactiveMP.chollogdet(A)
            @test all((inv(A), logdet(A)) .≈ ReactiveMP.cholinv_logdet(A))
        end

        for A in rand(rng, 3)
            @test inv(A) ≈ ReactiveMP.cholinv(A)
            @test (sqrt(A) * sqrt(A)') ≈ (ReactiveMP.cholsqrt(A) * ReactiveMP.cholsqrt(A)')
            @test A ≈ (ReactiveMP.cholsqrt(A) * ReactiveMP.cholsqrt(A)')
            @test logdet(A) ≈ ReactiveMP.chollogdet(A)
            @test all((inv(A), logdet(A)) .≈ ReactiveMP.cholinv_logdet(A))
        end

    end

end

end

ext: jl


text:
using SparseArrays
using LinearAlgebra

#
import Base.exp;
import Base.sin;
import Base.cos;
function exp(A::Matrix{Complex{BigFloat}})
    Z=zero(A);
    n=size(A,1);
    Ak=Matrix{eltype(A)}(I,n,n)
    for k=0:100
        Z+= Ak/(factorial(big(k)));
        Ak[:] =Ak*A;
    end

    return Z
end
function sin(A::Matrix{BigFloat})
    return (exp(1im*A)-exp(-1im*A))/2im
end

function cos(A::Matrix{BigFloat})
    return (exp(1im*A)+exp(-1im*A))/2
end

struct HelmholtzProblem
    sigma
    Aders
    Uv
    V
    Afun
    b
    #
    dx
    xv
    s
end

function compute_Mlincomb(prob::HelmholtzProblem,sigma,Y)
    if (prob.sigma != sigma)
        @show sigma
        @show prob.sigma
        error("Wrong sigma");
    end
    z=zeros(eltype(Y),size(Y,1));
    for i=1:size(Y,2)
        z += prob.Aders[i]*Y[:,i];
    end
    return z;
end
function compute_Mder(prob::HelmholtzProblem,mu)
    return prob.Afun(mu);
end

"""

    prob=create_helmholtz(n,a,b,c,N,k_eval,beta_eval,h_eval)

Creates an object representin the discretization of

```u''(x)+(1+mu*k(x))^2*u+beta(x)*u=h(x)```

by discretizing the interval [a,b] with n discretization
points. In the domain [b,c] we have

* `beta(x)=0`
* `k(x)=k0=k(c)`.

"""
function create_helmholtz(n,a,b,c,N,k_eval,beta_eval,h_eval)

    k0=k_eval(c);

    s=2;

    # Discretize interior interval [a,b]
    xv=Vector(range(a,b,length=n+1))[2:end];
    dx=xv[2]-xv[1];

    # Second derivative
    DD=spdiagm(0 => -2*ones(n-1),
               1 => ones(n-1),
               -1=> ones(n-2))/(dx^2)

    # Create auxiliary diagonal matrices
    K_eval=spdiagm(0=>[k_eval.(xv[1:end-1]);0]);
    L=spdiagm(0=>[beta_eval.(xv[1:end-1]);0]);
    II= spdiagm(0=>[ones(n-1);0]);
    K=mu -> spdiagm(0=>[(1 .+ mu*k_eval.(xv[1:end-1])).^2;0]);

    # Boundary condition matrices
    Xn=spzeros(n,2); Xn[n,1]=1; Xn[n,2]=1;
    Yn=spzeros(n,2);
    Yn[n,1]=1;
    Yn[n,2]=3/(2*dx);
    Yn[n-1,2]=-2/dx;
    Yn[n-2,2]=1/(2*dx);

    g=mu->cos((c-b)*(I+mu*k0))
    f=mu->(I+mu*k0)\sin((c-b)*(I+mu*k0));
    F=mu->spdiagm(0=>[g(mu);f(mu)]);

    # Create the A-function

    Afun1=mu->DD+K(mu)+L;
    Afun2=mu->Xn*F(mu)*Yn';
    Afun=mu->Afun1(mu)+Afun2(mu);

    # Compute derivatives of A
    ders=Vector{typeof(DD)}(undef,N);

    ders[1]=DD+II+L;
    ders[2]=2*K_eval;
    ders[3]=2*K_eval*K_eval;

    # Compute derivatives of boundary conditions
    J=big.(diagm(0=>zeros(N), 1=>ones(N-1)));
    fders=Float64.(f(J)[1,:].*factorial.(big.(0:N-1)))
    gders=Float64.(g(J)[1,:].*factorial.(big.(0:N-1)))
    for i=1:N
        Q=Xn*spdiagm(0=>[gders[i];fders[i]])*Yn';
        if (isassigned(ders,i))
            ders[i] += Q;
        else
            ders[i]=Q;
        end
    end

    rhs=[h_eval.(xv[1:end-1]);0] # right-hand side

    #
    Uv_type=Matrix{eltype(DD)}
    Uv=Vector{Uv_type}(undef,N);

    for i=(s+2):N
        Uv[i]=Xn*spdiagm(0=>[gders[i];fders[i]]);
    end
    V=Yn;

    myproblem=HelmholtzProblem(0,ders,Uv,V,Afun,rhs,dx,xv,s);

end

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
function get_steady_state(problem::LinearOptics{Scalar}) # @memoize
    @debug "start: get steady state"

    G  = get_interaction_matrix(problem)
    Ωₙ = apply_laser_over_atoms(problem.laser, problem.atoms)
    βₛ = (im/2)*(G \ Ωₙ)

    @debug "end  : get steady state"
    return βₛ
end
function get_steady_state(problem::NonLinearOptics{MeanField}; time_max=50)
    u0 = default_evolution_initial_condition(problem)
    tspan = (0.0, time_max)

    steady_state = time_evolution(problem, u0, tspan; save_on=false)[end]
    return steady_state
end

function time_evolution(problem::LinearOptics{T}, u₀, tspan::Tuple;  kargs...) where T <: Linear
    @debug "start: time evolution - LinearOptics"
    ### parameters == constant vector and matrices
    G = get_interaction_matrix(problem)
    Ωₙ = -0.5im*apply_laser_over_atoms(problem.laser, problem.atoms)
    parameters = view(G,:,:), view(Ωₙ,:)

    ### calls for solver
    problemFunction = get_evolution_function(problem)
    prob = ODEProblem(problemFunction, u₀, tspan, parameters)
    solution = DifferentialEquations.solve(prob, VCABM(); dt=1e-10, abstol=1e-10, reltol=1e-10, kargs...)

    @debug "end  : time evolution - LinearOptics"
    return solution
end

"""
    default_evolution_initial_condition(::Scalar) = zeros(ComplexF64, N)
"""
function default_evolution_initial_condition(problem::LinearOptics{Scalar})
    zeros(ComplexF64, problem.atoms.N) # I must use "zeros" and NOT an undef array - with trash data inside
end

get_evolution_function(problem::LinearOptics{Scalar}) = Scalar!

function Scalar!(du, u, p, t)
    G, Ωₙ = p

    du[:] = G*u + Ωₙ
    return nothing
end

# function extract_solution_from_Scalar_Problem(solution)
#     nSteps = length(solution.u)
#     time_array = [solution.t[i]  for i in 1:nSteps]
#     βₜ = [solution.u[i]  for i in 1:nSteps]
#     return time_array, βₜ
# end

# ### --------------- MEAN FIELD ---------------
function time_evolution(problem::NonLinearOptics{MeanField}, u₀, tspan::Tuple;  kargs...)
    @debug "start: time evolution - NonLinearOptics"
    G = get_interaction_matrix(problem)

    #=
        I don't sum over diagonal elements during time evolution
     thus, to avoid an IF statement, I put a zero on diagonal
    =#
    G[diagind(G)] .= zero(eltype(G))

    #=
        Also, the definition for Mean Field evolution needs
     +(Γ/2), where scalar kernel return -(Γ/2).
     Just multiply by -1 fixes
    =#
    G = -G

    Ωₙ = apply_laser_over_atoms(problem.laser, problem.atoms)
    Wₙ = zeros(eltype(Ωₙ),  size(Ωₙ))
    G_βₙ = zeros(eltype(Ωₙ), size(Ωₙ))
    parameters = view(G,:,:), view(diag(G),:), view(Ωₙ,:), Wₙ, problem.laser.Δ, problem.atoms.N, G_βₙ

    ### calls for solver
    problemFunction = get_evolution_function(problem)
    prob = ODEProblem(problemFunction, u₀, tspan, parameters)
    solution = DifferentialEquations.solve(prob, VCABM(); dt=1e-10, abstol=1e-10, reltol=1e-10, kargs...)

    @debug "end  : time evolution - NonLinearOptics"
    return solution
end
get_evolution_function(problem::NonLinearOptics{MeanField}) = MeanField!

"""
    default_evolution_initial_condition(NonLinearOptics{MeanField})
β₀ = zeros(ComplexF64, atoms.N)
z₀ = 2β₀.*conj.(β₀) .- 1
"""
function default_evolution_initial_condition(problem::NonLinearOptics{MeanField})
    β₀ = zeros(ComplexF64, problem.atoms.N)
    z₀ = 2β₀.*conj.(β₀) .- 1
    u₀ = vcat(β₀, z₀)
    return u₀
end
function MeanField!(du, u, p, t)
    # parameters
    G, diagG, Ωₙ, Wₙ, Δ, N, G_βₙ = p

    βₙ = @view u[1:N]
    zₙ = @view u[N+1:end]

    #= Code below is equivalento to =#
    # Wₙ[:] .= Ωₙ/2 - im*(G*βₙ - diagG.*βₙ) # don't forget the element wise multiplication of "diagG.*βₙ"
    # @. du[1:N] = (im*Δ - Γ/2)*βₙ + im*Wₙ*zₙ
    # @. du[N+1:end] = -Γ*(1 + zₙ) - 4*imag(βₙ*conj(Wₙ))

    BLAS.gemv!('N', ComplexF64(1.0), G, βₙ, ComplexF64(0.0), G_βₙ) # == G_βₙ = G*βₙ
    @simd for i ∈ eachindex(Wₙ)
        @inbounds Wₙ[i] = Ωₙ[i]/2 - im*(G_βₙ[i] - diagG[i]*βₙ[i])
    end
    @simd for i ∈ eachindex(βₙ)
        @inbounds du[i] = (im*Δ - Γ/2)*βₙ[i] + im*Wₙ[i]*zₙ[i]
    end
    @simd for i ∈ eachindex(zₙ)
        @inbounds du[i+N] = -Γ*(1 + zₙ[i]) - 4*imag(βₙ[i]*conj(Wₙ[i]))
    end

    return nothing
end
# function extract_solution_from_MeanField_Problem(solution)
#     nTimeSteps = length(solution.u)
#     time_array = [solution.t[i]  for i in 1:nTimeSteps   ]

#     N = length(solution.u[1])÷2
#     βₜ = [solution.u[i][1:N]  for i in 1:nTimeSteps   ]
#     zₜ = [solution.u[i][(N+1):end]  for i in 1:nTimeSteps   ]
#     return time_array, βₜ, zₜ
# end


text:
"""
`extrap_prep_dim(::Type{Reflect}, d)`

First, translate x into the domain over [lbound, 2(ubound-lbound)) i.e. into twice the size of the domain.
Next, if x is now in the upper part of this ''double-domain´´, reflect over the middle to obtain a new value x' for which f(x') == f(x), but where x' is inside the domain
"""
function extrap_prep_dim(::Type{Reflect}, d)
    xs_d = Symbol("xs_", d)
    quote
        start = lbound(etp.itp, $d)
        width = ubound(etp.itp, $d) - start

        $xs_d = mod($xs_d - start, 2width) + start
        $xs_d > start + width && (xs_d = start + width - $xs_d)
    end
end

extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}) = extrap_prep_dim(Reflect, d)
function extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}, ::Val{:lo})
    xs_d = Symbol("xs_", d)
    :($xs_d < lbound(etp.itp, $d) && $(extrap_prep_dim(Reflect, d)))
end
function extrap_prep{N,d}(::Type{Reflect}, ::Val{N}, ::Val{d}, ::Val{:hi})
    xs_d = Symbol("xs_", d)
    :($xs_d > ubound(etp.itp, $d) && $(extrap_prep_dim(Reflect, d)))
end


text:
module CornerTransferMethods
import Base: iterate
using Parameters: @unpack
using Base.Iterators: take, enumerate, rest
using Printf: @printf, @sprintf
using LinearAlgebra: svdvals!, svd!, eigen!, Hermitian, normalize!, diag, norm
using LinearAlgebra: diagm, ishermitian, eigen, pinv
using KrylovKit: eigsolve
using TensorNetworkTensors
using TensorOperations: @tensor, scalar, checked_similar_from_indices, tensoradd!, tensorcopy, tensorcopy!

export ctm, magnetisation, isingenvironment, isingenvironmentz2
export isingctm, isingctmz2, magofβ
export ctm_kernel, ctmiterable, rotsymctmiterable, transconjctmiterable
include("ctm.jl")
export cornereoe, transfermat, clength, transferopevals
include("auxiliary-functions.jl")
export atens, asztens, atenses, atensesz2
export βofh, hofβ
export σxtensor, σztensor
include("ising.jl")

export ozofβ, oxxofβ, tfisingpropagator, tfisinghamiltonian, tfisingctm, mag, energy
include("isingmpo.jl")

end # module


text:
@testset "vector tests" begin
  @testset "vector from low level" begin
    @test (x = K(ktn(KB, 0)); eltype(x) === Bool)
    @test (x = K(ktn(KG, 0)); eltype(x) === UInt8)
    @test (x = K(ktn(KH, 0)); eltype(x) === Int16)
    @test (x = K(ktn(KI, 0)); eltype(x) === Int32)
    @test (x = K(ktn(KJ, 0)); eltype(x) === Int64)
    @test (x = K(ktn(KE, 0)); eltype(x) === Float32)
    @test (x = K(ktn(KF, 0)); eltype(x) === Float64)
    @test (x = K(ktn(KC, 0)); eltype(x) === UInt8)
    @test (x = K(ktn(KS, 0)); eltype(x) === Symbol)
    @test (x = K(ktn(KD, 0)); eltype(x) === Date)
    @test_throws ArgumentError Q.K_Chars(ktn(KB,0))
    @test String(K(kp("abc"))) == "abc"
  end  # "vector from low level"
  @testset "vector constructors" begin
    @test K_Vector([1, 2]) == K([1, 2]) == K(1:2) == [1, 2]
    @test K_Vector([:a, :b]) == K([:a, :b]) == [:a, :b]
    @test (x = K[1]; eltype(x) == Int && Array(x) == [1])
    @test (x = K[1, 2., 3]; eltype(x) == Float64 && Array(x) == [1, 2, 3])
    @test K(["", ""]) == K["", ""]
    @test eltype(Q._vector(Int, 2)) == Int
  end
  @testset "vector indexing" begin
    let x = K[1, 2]
      @test x[1] === 1
      @test_throws BoundsError x[3]
      @test_throws BoundsError x[0]
      @test (x[1] = 10; x[1] === 10)
      @test_throws BoundsError x[3] = 0
    end
  end
  @testset "vector round trip" for T in NUMBER_TYPES
    a = [typemin(T), typemax(T), zero(T)]
    x = K(a)
    @test Array(x) == a
    for n in a
      x = K(n)
      @test T(x) == n
    end
    @test String(K("αβγ")) == "αβγ"
  end
  @testset "Vector operations" begin
    x = K[1]
    @test push!(x, 2) == [1, 2]
    @test copy!(x, [10, 20]) == [10, 20]
    @test fill!(x, 0) == [0, 0]
    @test append!(x, 1:2) == [0, 0, 1, 2]
    @test empty!(x) == Int[]
  end
  @testset "string print and show" begin
    @test string(K("αβγ"), K("δ")) == "αβγδ"
    @test show_to_string(K("δ")) == """K("δ")"""
    @test show_to_string(MIME"text/plain"(), K("δ")) == """K("δ")"""
  end
  @testset "temporal" begin
    @test (x = K[Date(2002)]; x[1] = Date(2000); x.a[1] == 0)
  end
end


text:
## Show the effect of seeded PageRank on MNIST

## Load the MNIST Graph and draw it
##
using DelimitedFiles
using SparseArrays
using LinearAlgebra
ei,ej = copy.(eachcol(Int.(readdlm("mnist-train-4.edges"))))
xy = readdlm("mnist-train-4.xy")
n = maximum(ei)
A = sparse(ei,ej,1,n,n)
@assert issymmetric(A)

##
## Simple PageRank
function simplepagerank(A,α,v)
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v
  nsteps = 2*ceil(Int,log(eps(1.0))/log(α)) # upper bound on steps
  for i=1:nsteps
    x = α*(A'*(x./d)) .+ (1-α).*v
  end
  return x/sum(x) # renormalize to probability
end

## Let's show global pagerank first
pr = simplepagerank(A,0.85,ones(n)./n)
## Now show the entire graph
using GraphRecipes
using LinearAlgebra
using Plots

ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(pr)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0,alpha=0.5)
savefig("pagerank-mnist.png")
## now show Seeded PageRank
v = zeros(n)
v[1] = 1
pr = simplepagerank(A,0.85,v)

using GraphRecipes
using LinearAlgebra
using Plots

ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(pr)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0,alpha=0.5)
savefig("seeded-pagerank-mnist.png")

##
## Try the same set of nodes labeled 1 as we used for the least squares exampmle.
using MLDatasets
train_x, train_y = MNIST.traindata()
nlabels = 10
S = findall(train_y .== 1)[1:nlabels]
##
v = zeros(n)
v[S] .= 1/length(S)
pr = simplepagerank(A,0.85,v)
ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(pr)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0,alpha=0.5)
savefig("seeded-pagerank-mnist-set-1.png")

##
v = zeros(n)
v[S] .= 1/length(S)
pr = simplepagerank(A,0.5,v)
ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(pr)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0,alpha=0.5)
savefig("seeded-pagerank-mnist-set-1-alpha-half.png")

##
v = zeros(n)
v[S] .= 1/length(S)
pr = simplepagerank(A,0.99,v)
ei, ej = findnz(A)
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=0, linecolor=1, linealpha=0.2, linewidth=0.5,
  markeralpha=0.2,colorbar=false,
  axis_buffer=0.02, background=nothing)
p=sortperm(pr)
scatter!(xy[p,1],xy[p,2],marker_z=-abs.(log10.(pr[p])).^0.5,markerstrokewidth=0,alpha=0.5)
savefig("seeded-pagerank-mnist-set-1-alpha-99.png")


text:
#=
The first two consecutive numbers to have two distinct prime factors are:

14 = 2 × 7
15 = 3 × 5

The first three consecutive numbers to have three distinct prime factors are:

644 = 2² × 7 × 23
645 = 3 × 5 × 43
646 = 2 × 17 × 19.

Find the first four consecutive integers to have four distinct prime factors. What is the first of these numbers?
=#

using Primes

function calc()
  num = 4
  arr = Int64[]
  i = 7
  while length(arr) < num
    if (length(keys(factor(i))) == num) && (length(arr) == 0 || last(arr) == i-1)
      push!(arr, i)
    else
      arr = Int64[]
    end
    i += 1
  end
  arr[1]
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
module ExtremeVertexDesigns

    using Polyhedra, CDDLib, JuMP, Statistics, StatsBase, LinearAlgebra

    """
    extreme_values( model::JuMP.Model; barycenters = true)::Matrix{Float64}

        From an inequality constrained JuMP model, return a matrix containing candidate
        design points.
        Note: Only barycenters of faces are currently supported. Soon, interior barycenters will also be
        added.
    """
    function extreme_values( model::JuMP.Model; barycenters = true)::Matrix{Float64}
        #ToDo: Invite other backends?
        poly = polyhedron(model, CDDLib.Library())
        removehredundancy!(poly)
        #ToDo: Faster memcpy then vcat from pointer?
        candidates = transpose(hcat(points(poly)...))
        if barycenters
            barycenters = [ mean( incidentpoints(poly, hidx) ) for hidx in eachindex(halfspaces(poly))]
            #ToDo: Faster memcpy then vcat from pointer?
            barycenters = transpose( hcat( barycenters...) )
            grand_center = mean( barycenters, dims = 1 )
            candidates = vcat( candidates, barycenters, grand_center )
        end
        return candidates
    end
    export extreme_values

    """
    federov_D(Ndesign::Int, candidates::Matrix; k = 1e-6)

        Performs Federov exchange algorithm to select a D-optimal design from a candidate matrix.
        Follows: Cook & Nachtsheim (1980).
        Note: λ is a ridge parameter to invert a matrix. This will be swapped out for an
        invertible starting design algorithm later...
    """
    function federov_D(Ndesign::Int, candidates::Matrix; λ = 1e-6)
        Nsupport, V = size(candidates)
        @assert Ndesign < Nsupport "Design points must be < than support candidates."
        init_samples = sample(1:Nsupport, Ndesign; replace = false)
        Xdesign = candidates[init_samples,:]
        Xsupport = candidates[ setdiff(1:Nsupport, init_samples),:]
        Nsupport -= Ndesign
        #ToDo: implement invertible starting design algorithm.
        tihkinov = Diagonal(ones(V) .* λ)

        while(true)
            XdtXd       = inv( Xdesign' * Xdesign + tihkinov )
            σ_design    = argmax( [ r' * XdtXd * r for r in eachrow( Xdesign ) ] )
            σ_support   = argmax( [ r' * XdtXd * r for r in eachrow( Xsupport ) ] )
            biggest_ds  = -Inf
            pair        = (0,0)
            for nd in 1:Ndesign, ns in 1:Nsupport
                σ_ds    = Xdesign[nd,:]' * XdtXd * Xsupport[ns,:]
                ds      = σ_support - ((σ_design * σ_support) - (σ_ds ^ 2) ) - σ_design
                if ds > biggest_ds
                    biggest_ds = ds
                    pair    = (nd, ns)
                end
            end
            if biggest_ds <= 0
                break
            end
            tmp = Xdesign[first(pair),:]
            Xdesign[first(pair),:] = Xsupport[last(pair),:]
            Xsupport[last(pair),:] = tmp
        end
        return (    score = det(inv( Xdesign' * Xdesign + tihkinov )),
                    design = Xdesign )
    end
    export federov_D

    """
    K_exchange_D(k::Int, Ndesign::Int, candidates::Matrix; ridge = 1e-6)

        Performs the k-exchange algorithm to select a D-optimal design from a candidate matrix.
        Follows: Johnson & Nachtsheim (1983).
        Note: λ is a ridge parameter to invert a matrix. This will be swapped out for an
        invertible starting design algorithm later...
    """
    function k_exchange_D(k::Int, Ndesign::Int, candidates::Matrix; λ = 1e-6)
        Nsupport, V = size(candidates)
        @assert Ndesign < Nsupport "Design points must be < than support candidates."
        @assert k < Ndesign "Cannot exchange more rows(k) than are present in the design."
        init_samples = sample(1:Nsupport, Ndesign; replace = false)
        Xdesign = candidates[init_samples,:]
        Xsupport = candidates[ setdiff(1:Nsupport, init_samples),:]
        Nsupport -= Ndesign
        #ToDo: implement invertible starting design algorithm.
        tihkinov = Diagonal(ones(V) .* λ)

        while(true)
            XdtXd       = inv( Xdesign' * Xdesign + tihkinov )
            σ_design    = [ r' * XdtXd * r for r in eachrow( Xdesign ) ]
            biggest_ds = -Inf
            for nd in sortperm( σ_design )[1:k]
                biggest_ds = -Inf
                pair = (0,0)
                σ_k_design  = Xdesign[nd,:]' * XdtXd * Xdesign[nd,:]
                for ns in 1:Nsupport
                    σ_support   = Xsupport[ns,:]' * XdtXd * Xsupport[ns,:]
                    σ_ds    = Xdesign[nd,:]' * XdtXd * Xsupport[ns,:]
                    ds      = σ_support - ((σ_k_design * σ_support) - (σ_ds ^ 2) ) - σ_k_design
                    if ds > biggest_ds
                        biggest_ds = ds
                        pair    = (nd, ns)
                    end
                end
                (biggest_ds <= 0) && break
                tmp = Xdesign[first(pair),:]
                Xdesign[first(pair),:] = Xsupport[last(pair),:]
                Xsupport[last(pair),:] = tmp
            end
            (biggest_ds <= 0) && break
        end
        return (    score = det(inv( Xdesign' * Xdesign + tihkinov )),
                    design = Xdesign )
    end
    export k_exchange_D

    """
        SimplexLatticeDesign( Components::Int, Spaces::Int )
    Returns an array of tuples (`Components` in length) which represent design points.
    """
    function SimplexLatticeDesign( Components::Int, Spaces::Int )
        @assert( ( Components > 0 ) && ( Spaces > 0 ), "Components and Spaces must be positive integers!" )
        DesignPoints = Inf
        try
            DesignPoints = Int( factorial( Components + Spaces - 1) / ( factorial(Spaces) * factorial(Components - 1) ))
        catch
            @warn("Cannot estimate number of design points. May execute for the lifetime of the universe.")
        end
        PossiblePoints = reverse(0:Spaces) ./ Spaces
        DesignSpace = repeat([PossiblePoints], Components)
        Possible = [ prod for prod in Base.product(DesignSpace...) if (sum(prod) == 1) ]

        if !isinf(DesignPoints)
            @assert( DesignPoints == length(Possible) )
        end
        return Possible
    end
    return SimplexLatticeDesign

    """
        SimplexCentroidDesign( Components::Int, Order::Union{ UnitRange{ Int }, Int} )
    Returns an array of tuples (`Components` in length) which represent design points.
    Please note: the Order can be an Int or a UnitRange.
    """
    function SimplexCentroidDesign( Components::Int, Order::Union{ UnitRange{ Int }, Int} )
        HighestOrder = Order[end]
        @assert( ( Components > 0 ) && ( HighestOrder > 0 ), "Components and the highest orders must be positive integers!" )
        @assert( HighestOrder <= Components, "The highest order of the design must be less than the number of components." )
        DesignPoints = Inf
        try
            DesignPoints = 2 ^ ( Components ) - 2 ^ ( Components - HighestOrder )
            if isa(Order, Int)
                DesignPoints = 2 ^ ( HighestOrder ) - 1
            end
        catch
            @warn("Cannot estimate number of design points. May execute for the lifetime of the universe.")
        end

        PossiblePoints = zeros( Components )
        DesignSpace = []
        for order in Order
            PossiblePoints[ 1:order ] .= 1.0 / order
            push!( DesignSpace, unique( collect( Combinatorics.permutations(PossiblePoints, Components) ) ) )
        end
        DesignSpace = ( isa(Order, Int) ) ? DesignSpace[1] : reduce(vcat, DesignSpace)
        if !isinf(DesignPoints)
            println( DesignPoints )
            println(length(DesignSpace))
            @assert( DesignPoints == length(DesignSpace) )
        end
        return DesignSpace
    end
    return SimplexCentroidDesign

end


text:
##############################################################################
##
## @eachrow
##
##############################################################################

# Recursive function that traverses the syntax tree of e, replaces instances of
# ":(:(x))" with ":x[row]".
function eachrow_replace(e::Expr)
    # Traverse the syntax tree of e
    if onearg(e, :cols)
        # cols(:x) becomes cols(:x)[row]
        return Expr(:ref, Expr(:call, :cols, e.args[2]), :row)
    end

    if e.head == :.
        if e.args[1] isa QuoteNode
            e.args[1] = Expr(:ref, e.args[1], :row)
            return e
        else
            return e
        end
    end

    Expr(e.head, (isempty(e.args) ? e.args : map(eachrow_replace, e.args))...)
end

eachrow_replace(e::QuoteNode) = Expr(:ref, e, :row)

# Set the base case for helper, i.e. for when expand hits an object of type
# other than Expr (generally a Symbol or a literal).
eachrow_replace(x) = x

function eachrow_find_newcols(e::Expr, newcol_decl)
    if e.head == :macrocall && e.args[1] == Symbol("@newcol")
        ea = e.args[3]
        # expression to assign a new column to df
        return (nothing, Any[Expr(:(=), ea.args[1], Expr(:call, ea.args[2], :undef, :_N))])
    else
        if isempty(e.args)
            return (e.args, Any[])
        end
        newargs = Any[]
        for ea in e.args
            (nea, newcol) = eachrow_find_newcols(ea, newcol_decl)
            nea != nothing && push!(newargs, nea)
            nea == nothing && length(newcol) > 0 && append!(newcol_decl, newcol)
        end
        return (Expr(e.head, newargs...), newcol_decl)
    end
end

eachrow_find_newcols(x, newcol_decl) = (x, Any[])

function eachrow_helper(df, body, deprecation_warning)
    # @deprecate cannot be used because eachrow is a macro, and the @warn should not be in
    # eachrow itself because then it will be displayed when the macro is evaluated.
    deprecation_warning && @warn "`@byrow!` and `@byrow` are deprecated, use `@eachrow` instead."
    e_body, e_newcols = eachrow_find_newcols(body, Any[])
    _df = gensym()
    quote
        let $_df = $df
            local _N = nrow($_df)
            local _DF = @transform($_df, $(e_newcols...))
            $(with_helper(:_DF, :(for row = 1:_N
                $(eachrow_replace(e_body))
            end)))
            _DF
        end
    end
end

"""
    @byrow!(d, expr)

Deprecated version of `@eachrow`, see: [`@eachrow`](@ref)

Acts the exact same way. It does not change the input argument `d` in-place.
"""
macro byrow!(df, body)
    esc(eachrow_helper(df, body, true))
end

"""
    @byrow(d, expr)

Deprecated version of `@eachrow`, see: [`@eachrow`](@ref)

Acts the exact same way.
"""
macro byrow(d, body)
    esc(eachrow_helper(d, body, true))
end

"""
    @eachrow(df, body)

Act on each row of a data frame, producing a new dataframe.
Similar to

```
for row in eachrow(copy(df))
    ...
end
```

Includes support for control flow and `begin end` blocks. Since the
"environment" induced by `@eachrow df` is implicitly a single row of `df`,
use regular operators and comparisons instead of their elementwise counterparts
as in `@with`. Note that the scope within `@eachrow` is a hard scope.

`eachrow` also supports special syntax for allocating new columns. The syntax
`@newcol x::Vector{Int}` allocates a new uninitialized column `:x` with an `Vector` container
with eltype `Int`.This feature makes it easier to use `eachrow` for data
transformations. `_N` is introduced to represent the number of rows in the data frame,
`_DF` represents the `DataFrame` including added columns, and `row` represents
the index of the current row.

Changes to the rows do not affect `df` but instead a freshly allocated data frame is returned
by `@eachrow`. Also note that the returned data frame does not share columns
with `df`. See [`@eachrow!`](@ref) which employs the same syntax but modifies
the data frame in-place.

Like with `@transform`, `@eachrow` supports the use of `cols` to work with column names
stored as variables. Using `cols` with a multi-column selector, such as a `Vector` of
`Symbol`s, is currently unsupported.

### Arguments

* `df` : an `AbstractDataFrame`
* `expr` : expression operated on row by row

### Returns

The modified `AbstractDataFrame`.

### Examples

```julia
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> let x = 0
            @eachrow df begin
                if :A + :B == 3
                    x += 1
                end
            end  #  This doesn't work without the let
            x
        end
2

julia> @eachrow df begin
            if :A > :B
                :A = 0
            end
        end
3×2 DataFrame
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      2
   2 │     0      1
   3 │     0      2

julia> df2 = @eachrow df begin
           @newcol colX::Vector{Float64}
           :colX = :B == 2 ? pi * :A : :B
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> varA = :A; varB = :B;

julia> df2 = @eachrow df begin
           @newcol colX::Vector{Float64}
           :colX = cols(varB) == 2 ? pi * cols(varA) : cols(varB)
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> x = [1, 1, 1];

julia> @eachrow df begin
           x[row] = :A
       end;

julia> x
3-element Vector{Int64}:
 1
 2
 3

julia> @eachrow df begin
           @newcol m::Vector{Float64}
           :m = mean(_DF[:, row])
       end
3×3 DataFrame
 Row │ A      B      m
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  2.0
   2 │     2      1  1.66667
   3 │     3      2  1.22222

```
"""
macro eachrow(df, body)
    esc(eachrow_helper(df, body, false))
end

function eachrow!_helper(df, body)
    e_body, e_newcols = eachrow_find_newcols(body, Any[])
    _df = gensym()
    quote
        let $_df = $df
            local _N = nrow($_df)
            local _DF = @transform!($_df, $(e_newcols...))
            $(with_helper(:_DF, :(for row = 1:_N
                $(eachrow_replace(e_body))
            end)))
            _DF
        end
    end
end

"""
    @eachrow!(df, body)

Act on each row of a data frame in-place, similar to

```
for row in eachrow(df)
    ... # Actions that modify `df`.
end
```

Includes support for control flow and `begin end` blocks. Since the
"environment" induced by `@eachrow! df` is implicitly a single row of `df`,
use regular operators and comparisons instead of their elementwise counterparts
as in `@with`. Note that the scope within `@eachrow!` is a hard scope.

`eachrow!` also supports special syntax for allocating new columns. The syntax
`@newcol x::Vector{Int}` allocates a new uninitialized column `:x` with an `Vector` container
with eltype `Int`.This feature makes it easier to use `eachrow` for data
transformations. `_N` is introduced to represent the number of rows in the data frame,
`_DF` represents the `dataframe` including added columns, and `row` represents
the index of the current row.

Changes to the rows directly affect `df`. The operation will modify the
data frame in place. See [`@eachrow`](@ref) which employs the same syntax but allocates
a fresh data frame.

Like with `@transform!`, `@eachrow!` supports the use of `cols` to work with column names
stored as variables. Using `cols` with a multi-column selector, such as a `Vector` of
`Symbol`s, is currently unsupported.

### Arguments

* `df` : an `AbstractDataFrame`
* `expr` : expression operated on row by row

### Returns

The modified `AbstractDataFrame`.

### Examples

```julia
julia> using DataFramesMeta

julia> df = DataFrame(A = 1:3, B = [2, 1, 2]);

julia> let x = 0
            @eachrow! df begin
                if :A + :B == 3
                    x += 1
                end
            end  #  This doesn't work without the let
            x
        end
2

julia> df2 = copy(df);

julia> @eachrow! df2 begin
           if :A > :B
               :A = 0
           end
       end;

julia> df2
3×2 DataFrame
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      2
   2 │     0      1
   3 │     0      2

julia> df2 = copy(df);

julia> @eachrow! df2 begin
           @newcol colX::Vector{Float64}
           :colX = :B == 2 ? pi * :A : :B
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> varA = :A; varB = :B;

julia> df2 = copy(df);

julia> @eachrow! df2 begin
           @newcol colX::Vector{Float64}
           :colX = cols(varB) == 2 ? pi * cols(varA) : cols(varB)
       end
3×3 DataFrame
 Row │ A      B      colX
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  3.14159
   2 │     2      1  1.0
   3 │     3      2  9.42478

julia> x = [1, 1, 1];

julia> @eachrow! df begin
           x[row] = :A
       end;

julia> x
3-element Vector{Int64}:
 1
 2
 3

julia> @eachrow! df begin
           @newcol m::Vector{Float64}
           :m = mean(_DF[:, row])
       end
3×3 DataFrame
 Row │ A      B      m
     │ Int64  Int64  Float64
─────┼───────────────────────
   1 │     1      2  2.0
   2 │     2      1  1.66667
   3 │     3      2  1.22222
```
"""
macro eachrow!(df, body)
    esc(eachrow!_helper(df, body))
end


text:
## A script to generate a tab-separated-values file
## containing questions for a pool or test to upload
## in batch mode to BlackBoard
## Questions can be formatted using MarkDown
## or LaTeX and can be templated so that randomization
## (for pools) can be performed
## Several BlackBoard question types can be used.

using JuliaBlackBoard

# used to create .txt file with same basename and same directory as script
OPEN = JuliaBlackBoard.OPEN(@__FILE__)

# Packages used in the following examples
using Plots
using StatsBase # for mean
using GLM       # for lm
using RCall     # for calling out to R for an example

OPEN() do io

    ## ---- Question ----
    ## Simple Markdown question
    ## mt"..." is exported from Mustache
    q = mt"""
# The mean

Find the mean of

```
2 3 4
```
"""

    # use a NUM type with answer computed
    question(io, NUM,
             q,                  # question
             mean((2,3,4))       # answer
             )

    ## Note: the use of mt"""...""" makes the multi-line string into a Mustache
    ## template. For randomization problems this will be more performant.
    ## This is not necessary for the above, but is useful in general.
    ## For string interpolation, use just """...""".

    ## ---- Question ----
    ## Same as above, but we use a Mustache template
    ## so that we can do randomization
    ##
    ## The bit below: {{{:xs}}}
    ## will be substituted in by the keyword argument to `q` (defined, as `q = mt...`)
    q = mt"""
# The mean

Find the mean of

```
{{{:xs}}}
```
"""

    for xs in ((2,3,4),
               (3,4,5))
        # use a NUM type with answer computed
        question(io, NUM,
                 q(xs = join(xs, ",")), # question
                 mean(xs))              # answer
    end

    ## For more detail, when `mt` is used to create a question, then the
    ## substitution of template variables (e.g. `{{{:xs}}}`) is done by
    ## calling the question like a function with `xs` as a named argument.
    ## This allows 0 or more such subsitutions. (Variables are specified as symbol in the
    ## template). The distinction between triple and double brackets is related to HTML
    ## replacement of special entitities. Using triple braces is recommended.

    ## ---- Question ----
    ## Same question with some LaTeX code. This uses HTML to represent the LaTeX
    ## BlackBoard does not allow us to upload ŷ (unicode) or MathML so we
    ## resort to HTML as generated by `tth`.
    ##
    ## This requires `tth` (https://sourceforge.net/projects/tth/) to be available as `tth`.
    ##
    ## tth does a poor job with inline math that has height, e.g. $\hat{y}$, $\bar{x}$. The -t
    ## switch is used, as the alternative doesn't seem to work at all. It is suggested to avoid these
    ## annotations, as done here:
     q = mt"""
# The mean

Find the mean,

$$
xbar = \frac{x_1 + x_2 + \cdots + x_n}{n}.
$$

of

```
2 3 4
```
"""
    question(io, NUM, q, mean((2,3,4)))

    ## ---- Question ----
    ## Same question with some LaTeX code.
    ENV["USE_MATHJAX"] = true
    ## By setting this environment variable, questions can be rendered within
    ## BlackBoard with MathML markup **PROVIDED** one manually edits them
    ## through the web editor and saves them. (Just edit then submit, no change
    ## required
    q = mt"""
# The mean

Find the mean,

$$
\bar{x} = \frac{x_1 + x_2 + \cdots + x_n}{n}.
$$

of

```
2 3 4
```
"""
    question(io, NUM, q, mean((2,3,4)))
    ENV["USE_MATHJAX"] = false

    ## ---- Question ----
    ## Same question with LaTeX generating an image for substition. We use
    ## a different template for the question -- one
    ## that only loads an image (png) of the file generated by calling `LaTeX` on
    ## the question.
    q = mt"""
\subsection{The mean}

Find the mean,

$$
\bar{x} = \frac{x_1 + x_2 + \cdots + x_n}{n}.
$$

of

\begin{verbatim}
2 3 4
\end{verbatim}
"""

    question(io, NUM, LaTeX(q), mean((2,3,4)))

    ## ---- Question ----
    ## Same question with LaTeX generating an image for substition.
    ## However, this uses Markdown for non math-mode markup
    q = mt"""
# The mean

Find the mean,

$$
\bar{x} = \frac{x_1 + x_2 + \cdots + x_n}{n}.
$$

of

```
2 3 4
```
"""
    question(io, NUM, LaTeX′(q), mean((2,3,4)))

    ## ---- Question ----
    ## Same question with randomization.
    q = mt"""
Find the mean, $\bar{x}$, of

\begin{verbatim}
{{:xs}}
\end{verbatim}
"""

    ## Note: When viewing the above question on blackboard, there will
    ## be no identifying text in the short preview mode of the pool
    ## viewer. If might be a good idea to have a slightly more general
    ## question template, something like this where the "section" is
    ## moved into the question template
    question_tpl = mt"""
## The mean

![]({{{:latex}}})
"""
    for xs in ((1,2,3), (4,5,6))
        str = q(xs=join(xs, ", "))

        question(io, NUM,
                 question_tpl(latex=File(JuliaBlackBoard.latex_to_png(str))),
                 mean((2,3,4)))
    end

    ## ---- Question ----
    ## Using Plots to generate a plot and show inline
    ## We use `Plot(p)` to generate the plot for rendering (using Base64 encoding)

    q = mt"""
# Regression line

Estimate the slope of the regression line from the following graph

![]({{{:plot}}})
"""

    xs = rand(25)
    ys = 1 .+ 2*xs + 0.25 * rand(25)
    res = lm(@formula(y ~ x), (y=ys, x=xs))
    a, b = coef(res)
    p = scatter(xs, ys, legend=false)
    plot!(p, x -> a + b*x, extrema(xs)...)

    question(io, NUM, q(plot=Plot(p)), b, 1/10) # last number is tolerance

    ## ---- Question ----
    ## Using R to generate a plot
    ## This is almost the same
    ## Here we create a png file on the file system, so we use `File` to produce the
    ## the graph for rendering. File could be used for any png file on a local drive
    q = mt"""
# Regression line

Estimate the slope of the regression line from the following graph

![]({{{:plot}}})
"""

    pngfile = tempname() * ".png"

    @rput pngfile
    R"""
xs = rnorm(20)
ys = rnorm(20)
f = ys ~ xs
res = lm(f)

png(pngfile)
plot(f)
abline(res)
dev.off()
"""
    res = @rget res
    a, b = round.(res[:coefficients], digits=1)

    question(io, NUM,
             q(plot=File(pngfile)),   # question
             b, 0.5)                  # answer and tolerance

    ## ---- Question ----
    ## Different question types:
    ## There a some of https://help.blackboard.com/Learn/Instructor/Tests_Pools_Surveys/Reuse_Questions/Upload_Questions
    ## supported.
    ## The answers are formatted differently for each. For the above we
    ## had NUM, an answer is just a number and an optional absolute tolerance

    ## Multiple Choice
    ## answer is a  tuple of pairs ("answer 1"=>false, "answer 2"=>true)

    q = mt"""
## Multiple choice

Which is true
"""

    # tuple of pairs "possible answer" => boolean; only 1 is true
    # answers can have markup in Markdown or LaTeX (if suitably escaped)
    answers = ("The **earth** is flat"  => false,
               "The **sky** is falling" => false,
               "climate change is real" => true)

    question(io, MC, q, answers)

    ## ---- Question ----

    ## Multiple answer
    ## As with MC, only 1 or more true answers
    q = mt"""
## Multiple answer

Which is **false**
"""

    answers = ("The earth is flat"      => true,
               "The sky is falling"     => true,
               "climate change is real" => false)

    question(io, MA, q, answers)

    ## ---- Question ----
    ## True False has a boolean answer

    q = mt"""
## True/false

Is $2 > 3$?
"""
    question(io, TF, q, false)

    ## ---- Question ----
    ## Fill in the blank.
    ## Answer is a tuple of values for matching, case insensitively
    ## Would need to hand edit for other matches

    q = mt"""
## Fill in the blanks

The ____ brown fox jumped over the lazy dog.
"""

    answers = ("quick",) # only one option for an answer
    question(io, FIB, q, answers)

    ## ---- Question ----
    ## Fill in blank plus
    ##
    ## The question puts variables in between [], as in [var1]
    ## The answers are pairs ("var1"=>("possible", "exact", "answers"))
    q = mt"""
## Fill in blank -- plus!

Use a letter or name

The first three non-negative integers are [one], [two], and [three]
"""
    # Collection of pairs: "variable"=>("different", "exact", "matches", "allowed")
    answers = ("one"   => ("one", "1"),
               "two"   => ("two", "2"),
               "three" => ("three", "3")
               )
    question(io, :FIB_PLUS, q, answers)

    ## ---- Question ----
    ## ORD:
    ## Text within ( ) may be repeated for
    ## each of the answers that are part of the Ordering question. The
    ## maximum number of answers is 100. The order typed in the file
    ## is the correct order. The system will randomly order the
    ## answers.
    q = mt"""

## Ordered

The numbers in order are:
"""
    # tuple of answers
    answers = ("one", "two", "three", "four", "five")
    question(io, ORD, q, answers)

    ## ---- Question ----
    ## MAT
    ## Text within ( ) may be repeated for each of the answers that
    ## are part of the Matching question. The maximum number of
    ## answers is 100. The system will randomly order the answers and
    ## their question. When uploading a matching question, you must
    ## have a one-to-one relationship between questions and
    ## answers. If not, correct answers may be marked incorrect if
    ## more than one answer has the same value.

    q = mt"""
## Matching

Match phonetically:
"""
    # collection of pairs
    answers = ("a"   => "Eh",
               "one" => "won")
    question(io, MAT, q, answers)

    ## ---- Question ----
    ## Likert scale
    q = mt"""
# Likert scale (for surveys)

Do you like BlackBoard?
"""
    question(io, OP, q)

    ## ---- Question ----
    ## JUMBLED_SENTENCE
    ##
    ## The question marks blanks by [varname]
    ## The answer is a collection of pairs with
    ## "answer" => ("varname",) where an answer can matche more than one variable or *no*
    ## variable.
    q = mt"""

## Jumbled sentence

Question goes here [var1] (choice 1) [var2] (choice 2)
"""

    # collection of pairs "variable" => "choice in popdown menu"
    # and possible nothing => (distractors...)
    as = ("var1"  => "choice 1",
          "var2"  => "choice 2",
          nothing => ("choice 3","choice 4")) # 2 distractors

    question(io, JUMBLED_SENTENCE, q, as)

    ## Short response
    ## answer is a prompt

    q = mt"""
## Short response

Did you like this test?
"""
    question(io, SR, q, "A quick answer")

    ## ---- Question ----
    ## Essay. No context or answer is necessary
    q = mt"""
## Essay

Please write more...
"""
    question(io, ESS, q)

    ## ---- Question ----
    ## File upload. No context or answer is necessary

    q = mt"""
## File upload

Upload a screen shot. Please use a scanning feature on your phone (e.g., as provided  in the iPhone's
Notes app)
"""
    question(io, FIL, q)

end


text:
module Projections

### Imports
using Reexport
using LinearAlgebra: LinearAlgebra, UniformScaling
@reexport using StaticArrays

### Exports
export DProjections, LDProjections, MDProjections, TR, TvR, TRv,
       anyisless, anyiszero, checksizes, eachisless, getr, gett,
       llist, lptype, product, project, reduce_dof, ⋅

const ⋅ = LinearAlgebra.dot

### Implemetation

# Structs
abstract type AbstractProjections{T} end
abstract type AbstractTR end

# This is a struct that stores each projection (in spherical harmonics space)
# of a correlation function (`F` and `Fs`, for instance) for the case when the
# projections correlation matrix is diagonal in both `l` and `m`, i.e.,
# `F_{l,m;l',m'} = F_{l,m} * δ_{l,l'} * δ_{m,m'}`
struct DProjections{N, T} <: AbstractProjections{T}
    t::T
    r::SVector{N, T}
end

# This is a struct that carries information that depends only on each the
# l-values.  An important example is the projections for `Fs`, which are
# independent of `m`.  Other example is the set `l * (l + 1) * D₀` for all `l`
struct LDProjections{N, T} <: AbstractProjections{T}
    t::T
    r::SVector{N, T}
end

struct MDProjections{N, T} <: AbstractProjections{T}
    t::T
    r::SVector{N, T}
end

# Object that carries and separates
# translational and rotational information
# of a physical quantity
struct TR{T} <: AbstractTR
    t::T
    r::T
end

# This is two-field struct were the
# first component stores `K.^2 .* D₀.t`,
# and the second stores `[j * (j + 1) * D₀.r for j in 0:l]`
struct TvR{T, N} <: AbstractTR
    t::Vector{T}
    r::SVector{N, T}
end

struct TRv{T, N} <: AbstractTR
    t::T
    r::SVector{N, T}
end

# Constructors
DProjections(t::Number ) = DProjections(t, SVector{0, typeof(t)}())
LDProjections(t::Number) = LDProjections(t, SVector{0, typeof(t)}())
MDProjections(t::Number) = MDProjections(t, SVector{0, typeof(t)}())

LDProjections(v::TRv) = LDProjections(v.t, v.r)

constructorname(::Type{P}) where {P <: DProjections } = DProjections
constructorname(::Type{P}) where {P <: LDProjections} = LDProjections
constructorname(::Type{P}) where {P <: MDProjections} = MDProjections
constructorname(p::AbstractProjections) = constructorname(typeof(p))

project(x) = x
project(t::Tuple{Number}) = DProjections(first(t))
project(t::Tuple) = DProjections(first(t), SVector(Base.tail(t)))

Base.getindex(v::TvR, i) = TRv(v.t[i], v.r)
Base.length(v::TvR) = length(v.t)
Base.size(v::TvR) = size(v.t)
@inline Base.iterate(v::TvR, i = 1) =
    (i % UInt) - 1 < length(v) ? (@inbounds v[i], i + 1) : nothing

Base.eltype(::Type{TvR{T, N}}) where {T, N} = TRv{T, N}
Base.eltype(::Type{P}) where {T, P <: AbstractProjections{T}} = T

gett(v) = v.t
gett(v::Number) = v
getr(v) = v.r
getr(v::Union{Number,SVector}) = v

function checksizes(::Type{<:LDProjections{L}},
                    ::Type{<:DProjections{N}}) where {L, N}
    @assert div((L + 1) * (L + 2), 2) == N + 1
end
function checksizes(::Type{<:MDProjections{M}},
                    ::Type{<:DProjections{N}}) where {M, N}
    l = div(M + 2, 2)
    @assert div(l * (l + 1), 2) == N + 1
end
function checksizes(::Type{<:MDProjections{M}},
                    ::Type{<:LDProjections{L}}) where {L, M}
    @assert 2L == M
end
checksizes(S::Type{<:DProjections}, T::Type{<:LDProjections}) = checksizes(T, S)
checksizes(S::Type{<:DProjections}, T::Type{<:MDProjections}) = checksizes(T, S)
checksizes(S::Type{<:LDProjections}, T::Type{<:MDProjections}) = checksizes(T, S)

Base.:+(u::TR, v::TR) = TR(u.t + v.t, u.r + v.r)

Base.:-(u::TR, v::TR) = TR(u.t - v.t, u.r - v.r)

Base.:*(x::Number, p::AbstractProjections) = constructorname(p)(x * p.t, x * p.r)
Base.:*(p::AbstractProjections, x::Number) = x * p
#
Base.:*(x::Number, v::TR) = TR(x * v.t, x * v.r)
Base.:*(v::TR, x::Number) = x * v
Base.:*(u::TR, v::TR) = TR(u.t * v.t, u.r * v.r)
#
Base.:*(v::TR, p::MDProjections{2}) = @inbounds TR(v.t * (p.t + p.r[1]), v.r * p.r[2])
Base.:*(p::AbstractProjections, v::TR) = v * p

Base.:/(p::AbstractProjections, x::Number) = constructorname(p)(p.t / x, p.r / x)
#
Base.:/(v::TR, x::Number) = TR(v.t / x, v.r / x)
#
⋅(u::TR, v::TR) = u.t * v.t + u.r * v.t

for op in (:+, :-, :*)
    for P in (:DProjections, :LDProjections, :MDProjections)
        @eval begin
            Base.$op(p::($P){0}, q::($P){0}) = $P( $op(p.t, q.t) )
            Base.$op(p::($P){N}, q::($P){N}) where {N} = $P( $op(p.t, q.t), broadcast($op, p.r, q.r) )
        end
    end
end

Base.:*(v::TRv, p::DProjections) = TR(v.t * p, DProjections(zero(v.t), v.r .* p.r))

Base.:*(p::LDProjections{0}, q::DProjections{0}) = DProjections(p.t * q.t)
function Base.:*(p::LDProjections{1}, q::DProjections{2})
    t = p.t * q.t
    @inbounds begin
        r₁ = p.r[1] * q.r[1]
        r₂ = p.r[1] * q.r[2]
    end
    return DProjections(t, SVector(r₁, r₂))
end
function Base.:*(p::LDProjections{L}, q::DProjections{N}) where {L, N}
    checksizes(typeof(p), typeof(q))
    t = p.t * q.t
    T = typeof(t)
    r = zero(MVector{N, T})
    R = 1:L
    @inbounds for i in R
        r[i] = p.r.v[i] * q.r.v[i]
    end
    j = 0
    for i = (L + 1):N
        j += 1
        r[i] = p.r.v[R[j]] * q.r.v[i]
        if j == length(R)
            j = 0
            R = R[2:end]
        end
    end
    return DProjections(t, SVector(r))
end
Base.:*(p::DProjections, q::LDProjections) = q * p

product(x, y, z) = x * y * z
product(p::MDProjections{0}, q::DProjections{0}, r::LDProjections{0}) =
    MDProjections(p.t * q.t * r.t)
#
function product(p::MDProjections{2}, q::DProjections{2}, r::LDProjections{1})
    t = p.t * q.t * r.t
    @inbounds begin
        r₁ = p.r[1] * q.r[1] * r.r[1]
        r₂ = p.r[2] * q.r[2] * r.r[1]
    end
    return MDProjections(t, SVector(r₁, r₂))
end
#
function product(p::MDProjections{M}, q::DProjections{N}, r::LDProjections{L}) where {M, N, L}
    checksizes(typeof(p), typeof(q))
    checksizes(typeof(q), typeof(r))
    t = p.t * q.t * r.t
    T = typeof(t)
    s = zero(MVector{M, T})
    @inbounds for l = 1:L
        i = div(l * (l + 1), 2)
        m = l + L
        s[l] = p.r[l] * q.r[i] * r.r[l]
        s[m] = p.r[m] * q.r[i + 1] * r.r[l]
    end
    return MDProjections(t, SVector(s))
end

for f in (:exp, :inv)
    @eval begin
        Base.$f(p::AbstractProjections) = constructorname(p)($f(p.t), broadcast($f, p.r))
    end
end

Base.:+(J::UniformScaling, v::TR) = TR(J.λ + v.t, J.λ + v.r)
Base.:+(v::TR, J::UniformScaling) = J + v

Base.:+(v::TR, p::DProjections{0})  = DProjections(v.t + p.t)
Base.:+(v::TR, p::DProjections)     = DProjections(v.t + p.t, v.r .+ p.r)
Base.:+(v::TR, p::LDProjections{0}) = LDProjections(v.t + p.t)
Base.:+(v::TR, p::LDProjections)    = LDProjections(v.t + p.t, v.r .+ p.r)

Base.:-(J::UniformScaling, v::TR) = TR(J.λ - v.t, J.λ - v.r)

Base.:^(v::TR, n::Integer) = TR(v.t^n, v.r^n)
Base.inv(v::TR) = TR(inv(v.t), inv(v.r))

Base.zero(::Type{TR{T}}) where {T} = (o = zero(T); TR(o, o))
Base.zero(v::TR{T}) where {T} = (o = zero(T); TR(o, o))

Base.zero(::Type{ DProjections{0, T}}) where {T} = DProjections(zero(T))
Base.zero(::Type{LDProjections{0, T}}) where {T} = LDProjections(zero(T))
Base.zero(::Type{ DProjections{N, T}}) where {N, T} = DProjections(zero(T), zero(SVector{N, T}))
Base.zero(::Type{LDProjections{N, T}}) where {N, T} = LDProjections(zero(T), zero(SVector{N, T}))
Base.zero(::Type{MDProjections{N, T}}) where {N, T} = MDProjections(zero(T), zero(SVector{N, T}))
Base.zero(p::AbstractProjections) = constructorname(p)(zero(p.t), zero(p.r))

Base.one(v::TR) = TR(one(v.t), one(v.r))
#
Base.one(::Type{LDProjections{0, T}}) where {T} = LDProjections(one(T))
Base.one(::Type{LDProjections{L, T}}) where {L, T} = LDProjections(one(T), ones(SVector{L, T}))
Base.one(p::AbstractProjections) = constructorname(p)(one(p.t), ones(p.r))

anyiszero(v) = iszero(v)
anyiszero(v::TR{<:Number}) = iszero(v.t) || iszero(v.r)

anyisless(u, v) = u < v
anyisless(u, v::TR{<:Number}) = u < v.t || u < v.r
anyisless(u::TR{<:Number}, v) = u.t < v || u.r < v
anyisless(u::TR{<:Number}, y::TR{<:Number}) = u.t < v.t || u.t < v.r

eachisless(u, v) = u < v
eachisless(u, v::TR{<:Number}) = u < v.t && u < v.r
eachisless(u::TR{<:Number}, v) = u.t < v && u.r < v
eachisless(u::TR{<:Number}, v::TR{<:Number}) = u.t < v.t && u.r < v.r

function Base.isapprox(u::TR{T}, v::TR{T}; atol::Real = 0, rtol::Real = Base.rtoldefault(T),
                       nans::Bool = false) where {T <: Number}
    return isapprox(u.t, v.t; atol = atol, rtol = rtol, nans = nans) &&
           isapprox(u.r, v.r; atol = atol, rtol = rtol, nans = nans)
end

Base.muladd(x::Number, p::DProjections{0}, q::DProjections{0}) = DProjections(muladd(x, p.t, q.t))
function Base.muladd(x::Number, p::DProjections{2}, q::DProjections{2})
    t = muladd(x, p.t, q.t)
    @inbounds begin
        r₁ = muladd(x, p.r[1], q.r[1])
        r₂ = muladd(x, p.r[2], q.r[2])
    end
    return DProjections(t, SVector(r₁, r₂))
end
Base.muladd(x::Number, p::LDProjections{0}, q::LDProjections{0}) = LDProjections(muladd(x, p.t, q.t))
function Base.muladd(x::Number, p::LDProjections{1}, q::LDProjections{1})
    t = muladd(x, p.t, q.t)
    @inbounds begin
        r₁ = muladd(x, p.r[1], q.r[1])
    end
    return LDProjections(t, SVector(r₁))
end

reduce_dof(v) = v
reduce_dof(v::TR) = v.t + v.r
reduce_dof(v::TRv) = LDProjections(v.t, v.t .+ v.r)

lptype(::Type{DProjections{0, T}}) where {T} = LDProjections{0, T}
lptype(::Type{DProjections{2, T}}) where {T} = LDProjections{1, T}
@generated function lptype(::Type{DProjections{N, T}}) where {N, T}
    L = div(isqrt(9 + 8N) - 3, 2)
    P = LDProjections{L, T}
    return :($P)
end

@generated function mptype(::Type{DProjections{N, T}}) where {N, T}
    l = div(isqrt(9 + 8N) - 3, 2)
    V = SVector{2l, T}
    return :($V)
end

llist(::Type{DProjections{0, T}}) where {T} = SVector{0, T}()
llist(::Type{DProjections{2, T}}) where {T} = SVector(T(2))

#lprod(l) = l * (l + 1)
#@generated function lplist(::Type{DProjections{N, T}}) where {N, T}
#    L = div(isqrt(9 + 8N) - 3, 2)
#    P = LDProjections{L, T}
#    quote
#        SVector(map(lprod, 1:$L)...)
#    end
#end
#llist(::T) where {T <: DProjections} = llist(T)
#llist(x) = zero(x)

end # module


text:
#Tests of proper type generation
using PyCall

@rosimport geometry_msgs.msg: PoseStamped, Vector3
@rosimport visualization_msgs.msg: Marker
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan
@rosimport std_msgs.msg: Empty
@rosimport std_msgs.msg: Float64, String

@test_throws ErrorException @rosimport fake_msgs.msg.FakeMsg
@test_throws ErrorException @rosimport std_msgs.msg.FakeMsg
@test_throws ErrorException @rosimport nav_msgs.srv.GetPlanRequest
rostypegen()

@test isdefined(Main, :geometry_msgs)
@test isdefined(Main, :std_msgs)
@test isdefined(Main, :nav_msgs)
@test isdefined(geometry_msgs.msg, :Point)
@test isdefined(geometry_msgs.msg, :Quaternion)
@test isdefined(geometry_msgs.msg, :Pose)
@test isdefined(geometry_msgs.msg, :PoseStamped)
@test isdefined(geometry_msgs.msg, :Vector3)
@test isdefined(std_msgs.msg, :Header)
@test isdefined(std_msgs.msg, :Empty)
@test isdefined(nav_msgs.msg, :Path)
@test isdefined(nav_msgs.srv, :GetPlan)
@test isdefined(nav_msgs.srv, :GetPlanRequest)
@test isdefined(nav_msgs.srv, :GetPlanResponse)

#type generation in a non-Main module
module TestModule
    using RobotOS
    @rosimport std_msgs.msg: Float32
    rostypegen(@__MODULE__)
end
@test !isdefined(std_msgs.msg, :Float32Msg)
@test isdefined(TestModule, :std_msgs)
@test isdefined(TestModule.std_msgs.msg, :Float32Msg)

#message creation
posestamp = geometry_msgs.msg.PoseStamped()
@test typeof(posestamp.pose) == geometry_msgs.msg.Pose
@test typeof(posestamp.pose.position) == geometry_msgs.msg.Point

#service creation
boolreq = std_srvs.srv.SetBoolRequest()
boolresp = std_srvs.srv.SetBoolResponse(true, "message")
planreq = nav_msgs.srv.GetPlanRequest()
planresp = nav_msgs.srv.GetPlanResponse()
@test typeof(planreq) == nav_msgs.srv.GetPlanRequest
@test typeof(planresp) == nav_msgs.srv.GetPlanResponse

#convert to/from PyObject
posestamp.pose.position = geometry_msgs.msg.Point(1,2,3)
pypose = convert(PyObject, posestamp)
@test pypose.pose.position.x == 1.
@test pypose.pose.position.y == 2.
@test pypose.pose.position.z == 3.
pypose2 = PyObject(posestamp)
@test pypose2.pose.position.x == 1.
@test pypose2.pose.position.y == 2.
@test pypose2.pose.position.z == 3.
pose2 = convert(geometry_msgs.msg.PoseStamped, pypose)
@test pose2.pose.position.x == 1.
@test pose2.pose.position.y == 2.
@test pose2.pose.position.z == 3.
@test_throws InexactError convert(geometry_msgs.msg.Pose, pypose)

#access message enum
@test visualization_msgs.msg.Marker.CUBE == 1

#Proper array handling
path = nav_msgs.msg.Path()
@test typeof(path.poses) == Array{geometry_msgs.msg.PoseStamped,1}
push!(path.poses, posestamp)

pypath = convert(PyObject, path)
path2 = convert(nav_msgs.msg.Path, pypath)
@test typeof(path.poses) == Array{geometry_msgs.msg.PoseStamped,1}
@test path2.poses[1].pose.position.x == 1.
@test path2.poses[1].pose.position.y == 2.
@test path2.poses[1].pose.position.z == 3.

#Issue #6 - Empty message
emptymsg = std_msgs.msg.Empty()
@test length(fieldnames(typeof(emptymsg))) == 0

#Issue #7/8 - Renaming conflicting message types
@test isdefined(std_msgs.msg, :Float64Msg)
@test isdefined(std_msgs.msg, :StringMsg)
@test Publisher{std_msgs.msg.Float64Msg}("x", queue_size=10) != nothing
@test Subscriber{std_msgs.msg.Float64Msg}("x", x->x, queue_size=10) != nothing


text:
module Data_Utils
    using Statistics
    using Plots
    using FFTW
    using Statistics
    using Optim
    using Images, FileIO, ImageIO
    using Distributions
    using Printf
    using SparseArrays
    using LinearAlgebra
    using FITSIO
    # put the cwd on the path so Julia can find the module
    push!(LOAD_PATH, pwd()*"/main")
    using DHC_2DUtils

    export readdust
    export S2_uniweights
    export S2_whitenoiseweights
    export invert_covmat
    export whitenoiseweights_forlog
    export readsfd
    export readsfd_fromsrc
    export dbn_coeffs_calc
    export calc_1dps
    export get_dbn_coeffs

    function readdust(Nx)
        RGBA_img = load(pwd()*"/scratch_DF/t115_clean_small.png")
        img = reinterpret(UInt8,rawview(channelview(RGBA_img)))
        return imresize(Float64.(img[1,:,:])[1:256, 1:256], (Nx, Nx))
    end

    #Functions that add simulated noise to the true image and return covariance and std^2 for the simulations
    function S2_uniweights(im, fhash; high=25, Nsam=1000, iso=false, norm=true, smooth=false, smoothval=0.8, coeff_mask=nothing)
        #=
        Noise model: uniform[-high, high]
        =#

        (Nx, Ny)  = size(im)
        if Nx != Ny error("Input image must be square") end
        (N1iso, Nf)    = size(fhash["S1_iso_mat"])

        # fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)

        S2   = DHC_compute(im, fhash, doS2=true, doS20=false, norm=norm, iso=iso)
        Ns    = length(S2)
        S2arr = zeros(Float64, Ns, Nsam)
        println("Ns", Ns)
        for j=1:Nsam
            noise = rand(Nx,Nx).*(2*high) .- high
            init = im+noise
            if smooth init = imfilter(init, Kernel.gaussian(smoothval)) end
            S2arr[:,j] = DHC_2DUtils.DHC_compute(init, fhash, doS2=true, doS20=false, norm=norm, iso=iso)
        end
        wt = zeros(Float64, Ns)
        for i=1:Ns
            wt[i] = std(S2arr[i,:])
        end
        msub = S2arr .- mean(S2arr, dims=2)
        cov = (msub * msub')./(Nsam-1)
        if coeff_mask!=nothing
            cov = cov[coeff_mask, coeff_mask]
            wt = wt[coeff_mask]
        end
        println("Condition number of diag cov", cond(Diagonal(wt.^(2))))
        println("Output of S2weights: Condition Number of covariance", cond(cov))
        return wt.^(2), cov
    end

    function S2_whitenoiseweights(im, fhash, dhc_args; loc=0.0, sig=1.0, Nsam=1000, norm=false, smooth=false, smoothval=0.8, coeff_mask=nothing)
        #=
        Noise model: N(0, std(I))
        Output: std^2 vector, full covariance matrix
        dhc_args is now a dictionary of args to pass to DHC_compute, incl which coeffs and apd or not
        =#
        if coeff_mask==nothing error("Running without coeff_mask") end

        (Nx, Ny)  = size(im)
        if Nx != Ny error("Input image must be square") end
        (N1iso, Nf)    = size(fhash["S1_iso_mat"])

        # fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
        S2   = DHC_compute_wrapper(im, fhash, norm=norm; dhc_args...)

        Ns    = length(S2)
        S2arr = zeros(Float64, Ns, Nsam)
        println("Ns", Ns)
        for j=1:Nsam
            noise = reshape(rand(Normal(loc, sig),Nx^2), (Nx, Nx))
            init = im+noise
            if smooth init = imfilter(init, Kernel.gaussian(smoothval)) end
            S2arr[:,j] = DHC_2DUtils.DHC_compute_wrapper(init, fhash, norm=norm; dhc_args...)
        end
        wt = zeros(Float64, Ns)
        for i=1:Ns
            wt[i] = std(S2arr[i,:])
        end
        msub = S2arr .- mean(S2arr, dims=2)
        cov = (msub * msub')./(Nsam-1)
        if coeff_mask!=nothing
            cov = cov[coeff_mask, coeff_mask]
            wt = wt[coeff_mask]
        end
        println("Output of S20weights: Condition number of diag cov", cond(Diagonal(wt.^(2))))
        println("Condition Number of covariance", cond(cov))
        return wt.^(2), cov
    end

    function whitenoiseweights_forlog(oriim, fhash, dhc_args; loc=0.0, sig=1.0, Nsam=10, norm=false, smooth=false, smoothval=0.8, coeff_mask=nothing)
        (Nx, Ny)  = size(oriim)
        if norm error("Not supposed to be run with norm true") end
        if Nx != Ny error("Input image must be square") end
        (N1iso, Nf)    = size(fhash["S1_iso_mat"])
        logim = log.(oriim)

        S2 = DHC_compute_wrapper(logim, fhash; dhc_args..., norm=norm) #Iso

        Ns    = length(S2)
        S2arr = zeros(Float64, Ns, Nsam)
        println("Ns", Ns)
        for j=1:Nsam #Order: Log(Smooth(Im+Noise))
            noise = reshape(rand(Normal(loc, sig),Nx^2), (Nx, Nx))
            init = oriim+noise
            if smooth init = imfilter(init, Kernel.gaussian(smoothval)) end
            loginit = log.(init)
            S2arr[:, j]   = DHC_compute_wrapper(loginit, fhash, norm=norm; dhc_args...) #Iso
        end
        wt = zeros(Float64, Ns)
        for i=1:Ns
            wt[i] = std(S2arr[i,:])
        end
        msub = S2arr .- mean(S2arr, dims=2)
        cov = (msub * msub')./(Nsam-1)
        if coeff_mask!=nothing
            cov = cov[coeff_mask, coeff_mask]
            wt = wt[coeff_mask]
        end
        println("Output of S2 weights: Condition number of diag cov", cond(Diagonal(wt.^(2))))
        println("Condition Number of covariance", cond(cov))
        return wt.^(2), cov
    end

    function invert_covmat(cov, epsilon=nothing)
        #=
        cov: Vector if Diagonal matrix | Covariance
        Computes inverse after adding epsilon for stabilization if requested.
        =#
        isdiag = ndims(cov)==1
        dstb = zeros(size(cov)[1])
        if epsilon!=nothing
            dstb = fill(epsilon, size(cov)[1])
        end

        if isdiag #Vector
            cov = cov + dstb
            println("Cond No before inversion", cond(Diagonal(cov)))
            icov = Diagonal(cov.^(-1))
            println("Numerical error wrt Id", mean(abs.((Diagonal(cov) * icov) - I)))
        else #Matrix
            cov = Symmetric(cov + Diagonal(dstb))
            println("Cond No before inversion", cond(cov))
            icov = inv(cov)
            println("Numerical error wrt Id", mean(abs.((cov * icov) - I)))
        end

        return icov
    end

    function readsfd(nx; logbool=false)
        # read FITS file with images
        # file in /n/fink2/dfink/mldust/dust10000.fits
        #     OR  /n/fink2/dfink/mldust/dust100000.fits
        #println("fn check")
        #fname = "~/Dust-Project/WST/DHC/scratch_NM/data/dust10000.fits"
        fname = "scratch_NM/data/dust10000.fits"
        f = FITS(fname, "r")
        big = read(f[1])

        (_,__,Nslice) = size(big)
        println(Nslice, " slices")
        if logbool #BUG: Is this consistent with apodization order
            println("LOG")
            im = log.(imresize(Float64.(big), nx, nx, Nslice)) #Log of SFD images
        else
            im = imresize(Float64.(big), nx, nx, Nslice)
        end
        return im
    end

    function readsfd_fromsrc(fname, nx; logbool=false)
        # read FITS file with images
        # file in /n/fink2/dfink/mldust/dust10000.fits
        #     OR  /n/fink2/dfink/mldust/dust100000.fits
        #println("fn check")
        f = FITS(fname, "r")
        big = read(f[1])

        (_,__,Nslice) = size(big)
        println(Nslice, " slices")
        if logbool #BUG: Is this consistent with apodization order
            println("LOG")
            im = log.(imresize(Float64.(big), nx, nx, Nslice)) #Log of SFD images
        else
            im = imresize(Float64.(big), nx, nx, Nslice)
        end
        return im
    end

    #=
    function mldust(nx; logbool=true)

        # read FITS file with images
        # file in /n/fink2/dfink/mldust/dust10000.fits
        #     OR  /n/fink2/dfink/mldust/dust100000.fits
        println("fn check")
        fname = "scratch_NM/data/dust10000.fits"
        f = FITS(fname, "r")
        big = read(f[1])

        (_,__,Nslice) = size(big)
        println(Nslice, " slices")

        #nx = 96
        #im = Float64.(big[11:11+nx-1, 11:11+nx-1, :])
        #nx = 48

        # rebin it to something smaller and take the log
        if logbool #BUG: Is this consistent with apodization order
            println("LOG")
            im = log.(imresize(Float64.(big), nx, nx, Nslice)) #Log of SFD images
        else
            im = imresize(Float64.(big), nx, nx, Nslice)
        end

        covar = dust_covar_matrix(im)

        #writefits, 'covar48_new.fits', covar

        # Cholesky factorization; chol.U and chol.L will be upper,lower triangular matrix
        println("Cholesky")
        chol = cholesky(covar)

        # generate some mock maps
        Nmock = 800
        ran = rand(Normal(), nx*nx, Nmock)

        recon = reshape(chol.L*ran, nx, nx, Nmock)

        return im, covar
    end
    =#
    function sfd_pixwise_covar_matrix(im)
        # derive covariance matrix from [Nx,Nx,Nimage] array
        println("Warning: will mean-normalize the input SFD distribution")
        (nx,__,Nslice) = size(im)
        eps=1E-6

        # make it mean zero with eps noise, otherwise covar is singular
        println("Mean zero")
        for i=1:Nslice  im[:, :, i] .-= (mean(im[:, :, i])+(rand()-0.5)*eps) end

        #println("Unit variance")
        #for i=1:Nslice  im[:, :, i] ./= std(im[:, :, i]) end

        dat = reshape(im, nx*nx, Nslice)

        # covariance matrix
        println("covariance")
        covar = (dat * dat') ./ Nslice

        # Check condition number
        println("Condition number  ", cond(covar))
        return covar

    end

    function dbn_coeffs_calc(dbnimg, fhash, dhc_args, coeff_mask)
        #dbnimg: Either a set of aimges from a dbn or their log
        (N1iso,Nf) = size(fhash["S1_iso_mat"])
        (N2iso,Nfsq) = size(fhash["S2_iso_mat"])
        Ncovsamp = size(dbnimg)[3]
        s20_dbn = zeros(Float64, Ncovsamp, 2+ (dhc_args[:iso] ? N1iso+N2iso : Nf+Nf^2))
        for idx=1:Ncovsamp
            s20_dbn[idx, :] = DHC_compute_wrapper(dbnimg[:, :, idx], fhash, norm=false; dhc_args...)
        end

        s_targ_mean = mean(s20_dbn, dims=1)
        scov  = (s20_dbn .- s_targ_mean)' * (s20_dbn .- s_targ_mean) ./(Ncovsamp-1)
        return s_targ_mean[coeff_mask], diag(scov)[coeff_mask], scov[coeff_mask, coeff_mask]
    end

    function get_dbn_coeffs(dbnimg, fhash, dhc_args; coeff_mask=nothing)
        #dbnimg: Either a set of aimges from a dbn or their log
        (N1iso,Nf) = size(fhash["S1_iso_mat"])
        (N2iso,Nfsq) = size(fhash["S2_iso_mat"])
        Ncovsamp = size(dbnimg)[3]
        s20_dbn = zeros(Float64, Ncovsamp, 2+ (dhc_args[:iso] ? N1iso+N2iso : Nf+Nf^2))
        for idx=1:Ncovsamp
            s20_dbn[idx, :] = DHC_compute_wrapper(dbnimg[:, :, idx], fhash, norm=false; dhc_args...)
        end
        if coeff_mask!=nothing
            return s20_dbn[:, coeff_mask]
        else
            return s20_dbn
        end
    end

    function calc_1dps(image, kbins::Array{Float64,1})
        #Assumes uniformly spaced kbins
        Nx = size(image)[1]
        fft_zerocenter = fftshift(fft(image))
        impf = abs2.(fft_zerocenter)
        x = (collect(1:Nx) * ones(Nx)') .- (Nx/2.0)
        y = (ones(Nx) * collect(1:Nx)') .- (Nx/2.0)
        krad  = (x.^2 + y.^2).^0.5
        meanpk = zeros(size(kbins))
        kdel = kbins[2] - kbins[1]
        #println(size(meanpk), " ", kdel)
        for k=1:size(meanpk)[1]
            filt = findall((krad .>= (kbins[k] - kdel./2.0)) .& (krad .<= (kbins[k] + kdel./2.0)))
            meanpk[k] = mean(impf[filt])
        end
        return meanpk
    end

    #Transforms
    function fnlog(s::Array{Any,1})
        return Base.log.(s)
    end

    function fndlog(s::Array{Any,1})
        return 1.0 ./s
    end

end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
struct Forcing{T<:Real}
    Fx::Array{T,2}
    Fy::Array{T,2}
    H::Array{T,2}
    η_ref::Array{T,2}
    Fη::Array{T,2}
    #sst_ref::Array{T,2}
    #sst_γ::Array{T,2}
end

function Forcing{T}(P::Parameter,G::Grid) where {T<:Real}

    @unpack wind_forcing_x,wind_forcing_y,topography = P

    if wind_forcing_x == "channel"
        Fx,_ = ChannelWind(T,P,G)
    elseif wind_forcing_x == "shear"
        Fx,_ = ShearWind(T,P,G)
    elseif wind_forcing_x == "double_gyre"
        Fx,_ = DoubleGyreWind(T,P,G)
    elseif wind_forcing_x == "constant"
        Fx,_ = ConstantWind(T,P,G)
    elseif wind_forcing_x == "none"
        Fx,_ = NoWind(T,P,G)
    end

    if wind_forcing_y == "channel"
        _,Fy = ChannelWind(T,P,G)
    elseif wind_forcing_y == "shear"
        _,Fy = ShearWind(T,P,G)
    elseif wind_forcing_y == "double_gyre"
        _,Fy = DoubleGyreWind(T,P,G)
    elseif wind_forcing_y == "constant"
        _,Fy = ConstantWind(T,P,G)
    elseif wind_forcing_x == "none"
        _,Fy = NoWind(T,P,G)
    end

    if topography == "ridge"
        H,_ = Ridge(T,P,G)
    elseif topography == "ridges"
        H = Ridges(T,P,G)
    elseif topography == "seamount"
        H = Seamount(T,P,G)
    elseif topography == "flat"
        H = FlatBottom(T,P,G)
    end

    η_ref = InterfaceRelaxation(T,P,G)
    Fη = KelvinPump(T,P,G)

    return Forcing{T}(Fx,Fy,H,η_ref,Fη)
end

"""Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally
as a cosine with strongest forcing in the middle and vanishing forcing at boundaries."""
function ChannelWind(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack Δ,x_u,y_u,Lx,Ly = G
    @unpack Fx0,Fy0,H,ρ,scale = P

    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    xx_u,yy_u = meshgrid(x_u,y_u)
    Fx = (scale*Δ*Fx0/ρ/H)*cos.(π*(yy_u/Ly .- 1/2)).^2
    Fy = (scale*Δ*Fy0/ρ/H)*cos.(π*(xx_u/Lx .- 1/2)).^2

    return T.(Fx),T.(Fy)
end

"""Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally
as a hyperbolic tangent with strongest shear in the middle."""
function ShearWind(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack Δ,x_u,y_u,Lx,Ly = G
    @unpack Fx0,Fy0,H,ρ,scale = P

    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    xx_u,yy_u = meshgrid(x_u,y_u)
    Fx = (scale*Δ*Fx0/ρ/H)*tanh.(2π*(yy_u/Ly .- 1/2))
    Fy = (scale*Δ*Fy0/ρ/H)*tanh.(2π*(xx_u/Lx .- 1/2))

    return T.(Fx),T.(Fy)
end

"""Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally
with a superposition of sin & cos for a double gyre circulation.
See Cooper&Zanna 2015 or Kloewer et al 2018."""
function DoubleGyreWind(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack Δ,x_u,y_u,Lx,Ly = G
    @unpack Fx0,Fy0,H,ρ,scale = P

    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    xx_u,yy_u = meshgrid(x_u,y_u)
    Fx = (scale*Δ*Fx0/ρ/H)*(cos.(2π*(yy_u/Ly .- 1/2)) + 2*sin.(π*(yy_u/Ly .- 1/2)))
    Fy = (scale*Δ*Fy0/ρ/H)*(cos.(2π*(xx_u/Lx .- 1/2)) + 2*sin.(π*(xx_u/Lx .- 1/2)))
    return T.(Fx),T.(Fy)
end

"""Returns constant in in space forcing matrices Fx,Fy."""
function ConstantWind(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack Δ,nux,nuy,nvx,nvy = G
    @unpack Fx0,Fy0,H,ρ,scale = P

    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    Fx = T.(scale*Δ*Fx0/ρ/H)*ones(T,nux,nuy)
    Fy = T.(scale*Δ*Fy0/ρ/H)*ones(T,nvx,nvy)

    return Fx,Fy
end

"""Returns constant in in space forcing matrices Fx,Fy."""
function NoWind(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack nux,nuy,nvx,nvy = G

    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    Fx = zeros(T,nux,nuy)
    Fy = zeros(T,nvx,nvy)

    return Fx,Fy
end

"""Returns a reference state for Newtonian cooling/surface relaxation shaped as a
hyperbolic tangent to force the continuity equation."""
function InterfaceRelaxation(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack x_T,y_T,Ly,Lx = G
    @unpack η_refh,η_refw = P

    # width of a tangent is defined as the distance between 65% of its minimum value and 65% of the max
    xx_T,yy_T = meshgrid(x_T,y_T)
    η_ref = -(η_refh/2)*tanh.(2π*(Ly/(4*η_refw))*(yy_T/Ly .- 1/2))
    return T.(η_ref)
end

"""Returns a matrix of water depth for the whole domain that contains a
Gaussian seamount in the middle. Water depth, heigth and width of the
seamount are adjusted with the constants H, topofeat_height and topofeat_width."""
function Seamount(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack x_T_halo,y_T_halo,Lx,Ly = G
    @unpack topo_width,topo_height,H = P

    xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)
    bumpx = exp.(-((xx_T .- Lx/2).^2)/(2*topo_width^2))
    bumpy = exp.(-((yy_T .- Ly/2).^2)/(2*topo_width^2))

    return T.(H .- topo_height*bumpx.*bumpy)
end

"""Returns a matrix of water depth for the whole domain that contains a
meridional Gaussian ridge in the middle. Water depth, heigth and width of the
ridge are adjusted with the constants water_depth, topofeat_height and topofeat_width."""
function Ridge(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack x_T_halo,y_T_halo,Lx,Ly = G
    @unpack topo_width,topo_height,H = P

    xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)
    bumpx = exp.(-((xx_T .- Lx/2).^2)/(2*topo_width^2))
    bumpy = exp.(-((yy_T .- Ly/2).^2)/(2*topo_width^2))

    Hx = H .- topo_height*bumpx
    Hy = H .- topo_height*bumpy
    return T.(Hx),T.(Hy)
end

"""Same as Ridge() but for n ridges at various x positions."""
function Ridges(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack x_T_halo,y_T_halo,Lx,Ly = G
    @unpack topo_width,topo_height,H = P
    @unpack topo_ridges_positions = P
    n_ridges = length(topo_ridges_positions)

    xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)

    # loop over bumps in x direction
    R = zero(xx_T) .+ H

    for i in 1:n_ridges
        R .-= topo_height*exp.(-((xx_T .- topo_ridges_positions[i]*Lx).^2)/(2*topo_width^2))
    end

    return T.(R)
end

"""Returns a matrix of constant water depth H."""
function FlatBottom(::Type{T},P::Parameter,G::Grid) where {T<:Real}
    @unpack nx,ny,haloη = G
    @unpack H = P
    return fill(T(H),(nx+2*haloη,ny+2*haloη))
end

"""Returns Kelvin wave pumping forcing of the continuity equation."""
function KelvinPump(::Type{T},P::Parameter,G::Grid) where {T<:Real}

    @unpack x_T,y_T,Lx,Ly = G
    @unpack R,ϕ,Δ = G
    @unpack A,ϕk,wk = P

    xx_T,yy_T = meshgrid(x_T,y_T)

    mϕ = 2π*R/360.          # meters per degree latitude
    y0 = Ly/2 - (ϕ-ϕk)*mϕ   # y[m] for central latitude of pumping

    Fη = A*Δ*exp.(-(yy_T.-y0).^2/(2*wk^2))

    return T.(Fη)
end

"""Time evolution of forcing."""
function Ftime(::Type{T},t::Int,ω::Real) where {T<:Real}
    return convert(T,sin(ω*t))
end


text:
"""
    normalize!(A, hasIntercept)

Centers and normalizes the columns of A in place

# Arguments

- A = 2d array of floats
- hasIntercept = boolean flag indicating whether the first column of A is the
  intercept

# Value

Centers and normalizes A in place and returns 2d arrays of the column means and L2
norms of A before standardization.

"""
function normalize!(A::AbstractArray{Float64,2}, hasIntercept::Bool)

    # If including intercept, subtract the column means from all other columns
    if hasIntercept == true
        means = mean(A, dims=1)
        A[:,2:end] = A[:,2:end].-transpose(means[2:end])
    else # Otherwise, subtract the column means from all columns
        means = Array{Float64}(undef, 1, size(A,2))
    end

    # Obtain and divide by the column L2 norms
    norms = Array{Float64}(undef, 1, size(A,2))
    for j in 1:size(A,2)
        norms[j] = norm(A[:,j])
        A[:,j] = A[:,j]/norms[j]
    end

    # Return the column means and norms
    return means, norms
end

"""
    backtransform!(B, meansX, meansZ, normsX, normsZ, Y, Xold, Zold)

Back-transform coefficient estimates B in place if X and Z were standardized
prior to the estimation-- when both X and Z include intercept columns.

# Arguments

- B = 2d array of coefficient estimates B
- meansX = 2d array of column means of X, obtained prior to standardizing X
- meansZ = 2d array of column means of Z, obtained prior to standardizing Z
- normsX = 2d array of column norms of X, obtained prior to standardizing X
- normsZ = 2d array of column norms of Z, obtained prior to standardizing Z
- Y = 2d array of response matrix Y
- Xold = 2d array row covariates X prior to standardization
- Zold = 2d array column covariates Z prior to standardization

# Value

None; back-transforms B in place

"""
function backtransform!(B::AbstractArray{Float64,2},
                        meansX::AbstractArray{Float64,2},
                        meansZ::AbstractArray{Float64,2},
                        normsX::AbstractArray{Float64,2},
                        normsZ::AbstractArray{Float64,2},
                        Y::AbstractArray{Float64,2},
                        Xold::AbstractArray{Float64,2},
                        Zold::AbstractArray{Float64,2})

    # Back transform the X intercepts (row main effects)
    prodX = (meansX[:,2:end]./normsX[:,2:end])*B[2:end, 2:end]
    B[1,2:end] = (B[1,2:end]-vec(prodX))./vec(normsZ[:,2:end])/normsX[1,1]

    # Back transform the Z intercepts (column main effects)
    prodZ = B[2:end, 2:end]*transpose(meansZ[:,2:end]./normsZ[:,2:end])
    B[2:end,1] = (B[2:end,1]-prodZ)./transpose(normsX[:,2:end])/normsZ[1,1]

    # Back transform the interactions
    B[2:end, 2:end] = B[2:end, 2:end]./transpose(normsX[:,2:end])./
                                       normsZ[:,2:end]

    # Re-estimate intercept
    B[1,1] = 0
    B[1,1] = mean(Y-Xold*B*transpose(Zold))
end

"""
    backtransform!(B, hasXIntercept, hasZIntercept,
	               meansX, meansZ, normsX, normsZ)

Back-transform coefficient estimates B in place if X and Z were standardized
prior to the estimation-- when not including intercept columns for either X
or Z.

# Arguments

- B = 2d array of coefficient estimates B
- hasXIntercept = boolean flag indicating whether or not to X has an
  intercept column
- hasZIntercept = boolean flag indicating whether or not to Z has an
  intercept column
- meansX = 2d array of column means of X, obtained prior to standardizing X
- meansZ = 2d array of column means of Z, obtained prior to standardizing Z
- normsX = 2d array of column norms of X, obtained prior to standardizing X
- normsZ = 2d array of column norms of Z, obtained prior to standardizing Z

# Value

None; back-transforms B in place

"""
function backtransform!(B::AbstractArray{Float64,2},
                        hasXIntercept::Bool, hasZIntercept::Bool,
                        meansX::AbstractArray{Float64,2},
                        meansZ::AbstractArray{Float64,2},
                        normsX::AbstractArray{Float64,2},
                        normsZ::AbstractArray{Float64,2})

    # Back transform the X intercepts (row main effects), if necessary
    if hasXIntercept == true
        prodX = (meansX[:,2:end]./normsX[:,2:end])*B[2:end, 2:end]
        B[1,2:end] = (B[1,2:end]-vec(prodX))./vec(normsZ[:,2:end])/normsX[1,1]
    end

    # Back transform the Z intercepts (column main effects), if necessary
    if hasZIntercept == true
        prodZ = B[2:end, 2:end]*transpose(meansZ[:,2:end]./normsZ[:,2:end])
        B[2:end,1] = (B[2:end,1]-prodZ)./transpose(normsX[:,2:end])/
                                         normsZ[1,1]
    end

    # Back transform the interactions, if necessary
    if (hasXIntercept == true) || (hasZIntercept == true)
        B[2:end, 2:end] = B[2:end, 2:end]./transpose(normsX[:,2:end])./
                                           normsZ[:,2:end]
    end

    # Back transform the interactions if not including any main effects
    if (hasXIntercept == false) && (hasZIntercept == false)
        B = B./transpose(normsX)./normsZ
    end
end

"""
    backtransform!(B, meansX, meansZ, normsX, normsZ, Y, Xold, Zold)

Back-transform coefficient estimates B in place if X and Z were standardized
prior to the estimation-- when both X and Z include intercept columns.

# Arguments

- B = 4d array of coefficient estimates B
- meansX = 2d array of column means of X, obtained prior to standardizing X
- meansZ = 2d array of column means of Z, obtained prior to standardizing Z
- normsX = 2d array of column norms of X, obtained prior to standardizing X
- normsZ = 2d array of column norms of Z, obtained prior to standardizing Z
- Y = 2d array of response matrix Y
- Xold = 2d array row covariates X prior to standardization
- Zold = 2d array column covariates Z prior to standardization

# Value

None; back-transforms B in place

# Some notes

B is a 4d array in which each coefficient matrix is stored along the third and fourth
dimension.

"""
function backtransform!(B::AbstractArray{Float64,4},
                        meansX::AbstractArray{Float64,2},
                        meansZ::AbstractArray{Float64,2},
                        normsX::AbstractArray{Float64,2},
                        normsZ::AbstractArray{Float64,2},
                        Y::AbstractArray{Float64,2},
                        Xold::AbstractArray{Float64,2},
                        Zold::AbstractArray{Float64,2})

    # Iterate through the first dimension of B to back-transform each
	# coefficient matrix.
    for j in 1:size(B,4)
        for i in 1:size(B,3)
            # Back transform the X intercepts (row main effects)
            prodX = (meansX[:,2:end]./normsX[:,2:end])*B[2:end, 2:end,i,j]
            B[1,2:end,i,j] = (B[1,2:end,i,j]-vec(prodX))./vec(normsZ[:,2:end])/
                                                    normsX[1,1]

            # Back transform the Z intercepts (column main effects)
            prodZ = B[2:end,2:end,i,j]*transpose(meansZ[:,2:end]./normsZ[:,2:end])
            B[2:end,1,i,j] = (B[2:end,1,i,j]-prodZ)./transpose(normsX[:,2:end])/
                                                normsZ[1,1]

            # Back transform the interactions
            B[2:end,2:end,i,j] = B[2:end,2:end,i,j]./transpose(normsX[:,2:end])./
                                                normsZ[:,2:end]

            # Re-estimate intercept
            B[1,1,i,j] = 0
            B[1,1,i,j] = mean(Y-Xold*B[:,:,i,j]*transpose(Zold))
        end
    end
end

"""
    backtransform!(B, hasXIntercept, hasZIntercept,
                   meansX, meansZ, normsX, normsZ)

Back-transform coefficient estimates B in place if X and Z were standardized
prior to the estimation-- when not including intercept columns for either X
or Z.

# Arguments

- B = 4d array of coefficient estimates
- hasXIntercept = boolean flag indicating whether or not to X has an
  intercept column
- hasZIntercept = boolean flag indicating whether or not to Z has an
  intercept column
- meansX = 2d array of column means of X, obtained prior to standardizing X
- meansZ = 2d array of column means of Z, obtained prior to standardizing Z
- normsX = 2d array of column norms of X, obtained prior to standardizing X
- normsZ = 2d array of column norms of Z, obtained prior to standardizing Z

# Value

None; back-transforms B in place.

# Some notes

B is a 4d array in which each coefficient matrix is stored along the third and fourth
dimension.

"""
function backtransform!(B::AbstractArray{Float64,4},
                        hasXIntercept::Bool, hasZIntercept::Bool,
                        meansX::AbstractArray{Float64,2},
                        meansZ::AbstractArray{Float64,2},
                        normsX::AbstractArray{Float64,2},
                        normsZ::AbstractArray{Float64,2})

    # Iterate through the first dimension of B to back-transform each
	# coefficient matrix:
    for j in 1:size(B,4)
        for i in 1:size(B,3)
            # Back transform the X intercepts (row main effects), if necessary
            if hasXIntercept == true
                prodX = (meansX[:,2:end]./normsX[:,2:end])*B[2:end,2:end,i,j]
                B[1,2:end,i,j] = (B[1,2:end,i,j]-vec(prodX))./vec(normsZ[:,2:end])/
                                                        normsX[1,1]
            end

            # Back transform the Z intercepts (column main effects), if necessary
            if hasZIntercept == true
                prodZ = B[2:end,2:end,i,j]*transpose(meansZ[:,2:end]./
                    normsZ[:,2:end])
                B[2:end,1,i,j] = (B[2:end,1,i,j]-prodZ)./transpose(normsX[:,2:end])/
                                                    normsZ[1,1]
            end

            # Back transform the interactions, if necessary
            if (hasXIntercept == true) || (hasZIntercept == true)
                B[2:end,2:end,i,j] = B[2:end,2:end,i,j]./transpose(normsX[:,2:end])./
                                                    normsZ[:,2:end]
            end

            # Back transform the interactions if not including any main effects
            if (hasXIntercept == false) && (hasZIntercept == false)
                B[:,:,i,j] = B[:,:,i,j]./transpose(normsX)./normsZ
            end
        end
    end
end