text:
module NeuralNetwork

typealias AA AbstractArray

include("utils.jl")
include("mlp.jl")

export ConvolutionLayer, LogisticRegression, fit!, predict
export onehot, maxindices

type ConvolutionLayer
    filterw::Int
    filterh::Int
    stepw::Int
    steph::Int
    weights::Matrix
    bias::Matrix

    function ConvolutionLayer(filter_size::Int, step_size::Int,
                              weights = zeros(filter_size, filter_size))
        filterw, filterh = to_tuple_if_int(filter_size)
        stepw, steph = to_tuple_if_int(step_size)
        ConvolutionLayer(filterw, filterh, stepw, steph, weights)
    end

    function ConvolutionLayer(filterw::Int, filterh::Int,
                              stepw::Int, steph::Int,
                              weights = zeros(filterw, filterh))
        new(filterw, filterh, stepw, steph, weights)
    end
end

function forward(layer::ConvolutionLayer, image::Matrix)
    """
    Forward computation for an image
    """
    filterw, filterh = layer.filterw, layer.filterh
    startpointsx = 1:layer.stepw:(size(image, 1)-filterw+1)
    startpointsy = 1:layer.steph:(size(image, 2)-filterh+1)

    output = zeros(length(startpointsx), length(startpointsy))
    for (ox, x) in enumerate(startpointsx)
        for (oy, y) in enumerate(startpointsy)
            subimage = image[x:x+filterw-1, y:y+filterh-1]
            output[ox, oy] = sum(subimage .* layer.weights)
        end
    end
    output += layer.bias
    return output
end

function backward(layer::ConvolutionLayer)
end

end # module

max_stars_repo_name: IshitaTakeshi/NeuralNetwork

max_issues_repo_name: IshitaTakeshi/NeuralNetwork

max_forks_repo_name: IshitaTakeshi/NeuralNetwork


text:
using CUDA
using Oceananigans.Solvers: solve!, set_source_term!
using Oceananigans.Solvers: poisson_eigenvalues
using Oceananigans.Models.NonhydrostaticModels: solve_for_pressure!
using Oceananigans.Models.HydrostaticFreeSurfaceModels: _compute_w_from_continuity!
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions

function poisson_solver_instantiates(grid, planner_flag)
    solver = FFTBasedPoissonSolver(grid, planner_flag)
    return true  # Just making sure the FFTBasedPoissonSolver does not error/crash.
end

function random_divergent_source_term(grid)
    arch = architecture(grid)
    default_bcs = FieldBoundaryConditions()
    u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
    v_bcs = regularize_field_boundary_conditions(default_bcs, grid, :v)
    w_bcs = regularize_field_boundary_conditions(default_bcs, grid, :w)

    Ru = CenterField(grid, boundary_conditions=u_bcs)
    Rv = CenterField(grid, boundary_conditions=v_bcs)
    Rw = CenterField(grid, boundary_conditions=w_bcs)
    U = (u=Ru, v=Rv, w=Rw)

    Nx, Ny, Nz = size(grid)
    set!(Ru, rand(Nx, Ny, Nz))
    set!(Rv, rand(Nx, Ny, Nz))
    set!(Rw, rand(Nx, Ny, Nz))

    fill_halo_regions!(Ru, arch)
    fill_halo_regions!(Rv, arch)
    fill_halo_regions!(Rw, arch)

    # Compute the right hand side R = ∇⋅U
    ArrayType = array_type(arch)
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, U.u.data, U.v.data, U.w.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    return R, U
end

function random_divergence_free_source_term(grid)
    default_bcs = FieldBoundaryConditions()
    u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
    v_bcs = regularize_field_boundary_conditions(default_bcs, grid, :v)
    w_bcs = regularize_field_boundary_conditions(default_bcs, grid, :w)

    # Random right hand side
    Ru = CenterField(grid, boundary_conditions=u_bcs)
    Rv = CenterField(grid, boundary_conditions=v_bcs)
    Rw = CenterField(grid, boundary_conditions=w_bcs)
    U = (u=Ru, v=Rv, w=Rw)

    Nx, Ny, Nz = size(grid)
    set!(Ru, rand(Nx, Ny, Nz))
    set!(Rv, rand(Nx, Ny, Nz))
    set!(Rw, zeros(Nx, Ny, Nz))

    arch = architecture(grid)
    fill_halo_regions!(Ru, arch, nothing, nothing)
    fill_halo_regions!(Rv, arch, nothing, nothing)
    fill_halo_regions!(Rw, arch, nothing, nothing)

    event = launch!(arch, grid, :xy, _compute_w_from_continuity!, U, grid,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    fill_halo_regions!(Rw, arch, nothing, nothing)

    # Compute the right hand side R = ∇⋅U
    ArrayType = array_type(arch)
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, Ru.data, Rv.data, Rw.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    return R
end

#####
##### Regular rectilinear grid Poisson solver
#####

function divergence_free_poisson_solution(grid, planner_flag=FFTW.MEASURE)
    arch = architecture(grid)
    ArrayType = array_type(arch)
    FT = eltype(grid)

    solver = FFTBasedPoissonSolver(grid, planner_flag)
    R, U = random_divergent_source_term(grid)

    p_bcs = FieldBoundaryConditions(grid, (Center, Center, Center))
    ϕ   = CenterField(grid, boundary_conditions=p_bcs)  # "kinematic pressure"
    ∇²ϕ = CenterField(grid, boundary_conditions=p_bcs)

    # Using Δt = 1 but it doesn't matter since velocities = 0.
    solve_for_pressure!(ϕ.data, solver, 1, U)

    compute_∇²!(∇²ϕ, ϕ, arch, grid)

    return CUDA.@allowscalar interior(∇²ϕ) ≈ R
end

#####
##### Test that Poisson solver error converges as error ~ N⁻²
#####

ψ(::Type{Bounded}, n, x) = cos(n*x/2)
ψ(::Type{Periodic}, n, x) = cos(n*x)

k²(::Type{Bounded}, n) = (n/2)^2
k²(::Type{Periodic}, n) = n^2

function analytical_poisson_solver_test(arch, N, topo; FT=Float64, mode=1)
    grid = RectilinearGrid(arch, FT, topology=topo, size=(N, N, N), x=(0, 2π), y=(0, 2π), z=(0, 2π))
    solver = FFTBasedPoissonSolver(grid)

    xC, yC, zC = nodes((Center, Center, Center), grid, reshape=true)

    TX, TY, TZ = topology(grid)
    Ψ(x, y, z) = ψ(TX, mode, x) * ψ(TY, mode, y) * ψ(TZ, mode, z)
    f(x, y, z) = -(k²(TX, mode) + k²(TY, mode) + k²(TZ, mode)) * Ψ(x, y, z)

    solver.storage .= convert(array_type(arch), f.(xC, yC, zC))

    ϕc = rhs = solver.storage
    solve!(ϕc, solver, rhs)

    ϕ = real(Array(solver.storage))

    L¹_error = mean(abs, ϕ - Ψ.(xC, yC, zC))

    return L¹_error
end

function poisson_solver_convergence(arch, topo, N¹, N²; FT=Float64, mode=1)
    error¹ = analytical_poisson_solver_test(arch, N¹, topo; FT, mode)
    error² = analytical_poisson_solver_test(arch, N², topo; FT, mode)

    rate = log(error¹ / error²) / log(N² / N¹)

    TX, TY, TZ = topo
    @info "Convergence of L¹-normed error, $(typeof(arch)), $FT, ($(N¹)³ -> $(N²)³), topology=($TX, $TY, $TZ): $rate"

    return isapprox(rate, 2, rtol=5e-3)
end

#####
##### Vertically stretched Poisson solver
#####

get_grid_size(TX, TY, TZ, Nx, Ny, Nz) = (Nx, Ny, Nz)
get_grid_size(::Type{Flat}, TY, TZ, Nx, Ny, Nz) = (Ny, Nz)
get_grid_size(TX, ::Type{Flat}, TZ, Nx, Ny, Nz) = (Nx, Nz)

get_xy_interval_kwargs(TX, TY, TZ) = (x=(0, 1), y=(0, 1))
get_xy_interval_kwargs(TX, ::Type{Flat}, TZ) = (x=(0, 1),)
get_xy_interval_kwargs(::Type{Flat}, TY, TZ) = (y=(0, 1),)

function vertically_stretched_poisson_solver_correct_answer(FT, arch, topo, Nx, Ny, zF)
    Nz = length(zF) - 1
    sz = get_grid_size(topo..., Nx, Ny, Nz)
    xy_intervals = get_xy_interval_kwargs(topo...)
    vs_grid = RectilinearGrid(arch, FT; topology=topo, size=sz, z=zF, xy_intervals...)
    solver = FourierTridiagonalPoissonSolver(vs_grid)

    p_bcs = FieldBoundaryConditions(vs_grid, (Center, Center, Center))
    ϕ   = CenterField(vs_grid, boundary_conditions=p_bcs)  # "kinematic pressure"
    ∇²ϕ = CenterField(vs_grid, boundary_conditions=p_bcs)

    R = random_divergence_free_source_term(vs_grid)

    set_source_term!(solver, R)
    ϕc = solver.storage
    solve!(ϕc, solver)

    # interior(ϕ) = solution(solver) or solution!(interior(ϕ), solver)
    CUDA.@allowscalar interior(ϕ) .= real.(solver.storage)
    compute_∇²!(∇²ϕ, ϕ, arch, vs_grid)

    return CUDA.@allowscalar interior(∇²ϕ) ≈ R
end

ext: jl


text:
# This file is a part of RadiationDetectorDSP.jl, licensed under the MIT License (MIT).

__precompile__(true)

module RadiationDetectorDSP

using ArgCheck
using ArraysOfArrays
using Distributions
using DSP
using ElasticArrays
using LinearAlgebra
using ParallelProcessingTools
using Parameters
using RadiationDetectorSignals
using Random
using RecipesBase
using SIMD
using StaticArrays
using Statistics
using StatsBase
using TypedTables

include("util.jl")
include("samples.jl")
include("filters.jl")
include("trapezoidal_filter.jl")
include("generators.jl")

end # module


text:
# This file was generated, do not modify it. # hide
PyPlot.close_figs() # hide


text:
"""
    prlyap(A, C) -> X

Solve the reverse-time periodic discrete-time Lyapunov equation

    A'σXA + C = X

where `σ` is the forward shift operator `σX(i) = X(i+1)`.

For the periodic matrices `A` and `C` with the same sampling period and commensurate periods,
the periodic matrix `X` is determined, whose period is automatically set to the
the least common period of `A` and `C`.
"""
function prlyap(A::PeriodicArray, C::PeriodicArray)
   A.Ts ≈ C.Ts || error("A and C must have the same sampling time")
   period = promote_period(A, C)
   nta = numerator(rationalize(period/A.period))
   K = nta*A.nperiod*A.dperiod
   X = pslyapd(A.M, C.M; adj = true)
   return PeriodicArray(X, period; nperiod = div(K,size(X,3)))
end
prlyap(A::PeriodicArray, C::AbstractMatrix) = prlyap(A, PeriodicArray(C, A.Ts; nperiod = 1))
"""
    pflyap(A, C) -> X

Solve the forward-time periodic discrete-time Lyapunov equation

    A'XA + C = σX

where `σ` is the forward shift operator `σX(i) = X(i+1)`.

For the periodic matrices `A` and `C` with the same sampling period and commensurate periods,
the periodic matrix `X` is determined, whose period is automatically set to the
the least common period of `A` and `C`.
"""
function pflyap(A::PeriodicArray, C::PeriodicArray)
   A.Ts ≈ C.Ts || error("A and C must have the same sampling time")
   period = promote_period(A, C)
   nta = numerator(rationalize(period/A.period))
   K = nta*A.nperiod*A.dperiod
   X = pslyapd(A.M, C.M; adj = false)
   return PeriodicArray(X, period; nperiod = div(K,size(X,3)))
end
pflyap(A::PeriodicArray, C::AbstractMatrix) = pflyap(A, PeriodicArray(C, A.Ts; nperiod = 1))
"""
    pslyapd(A, C; adj = true) -> X

Solve the periodic discrete-time Lyapunov equation.

For the square `n`-th order periodic matrices `A(i)`, `i = 1, ..., pa` and
`C(i)`, `i = 1, ..., pc`  of periods `pa` and `pc`, respectively,
the periodic solution `X(i)`, `i = 1, ..., p` of period `p = lcm(pa,pc)` of the
periodic Lyapunov equation is solved:  for `adj = true`,

    A(i)'*X(i+1)*A(i) + C(i) = X(i), i = 1, ..., p

and for `adj = false`

    A(i)*X(i)*A(i)' + C(i) = X(i+1), i = 1, ..., p.

The periodic matrices `A` and `C` are stored in the `n×n×pa` and `n×n×pc` 3-dimensional
arrays `A` and `C`, respectively, and `X` results as a `n×n×p` 3-dimensional array.

The periodic discrete analog of the Bartels-Steward method based on the periodic Schur form
of the periodic matrix `A` is employed [1].

_Reference:_

[1] A. Varga. Periodic Lyapunov equations: some applications and new algorithms.
              Int. J. Control, vol, 67, pp, 69-87, 1997.
"""
function pslyapd(A::AbstractArray{T1, 3}, C::AbstractArray{T2, 3}; adj::Bool = true) where {T1, T2}
   n = LinearAlgebra.checksquare(A[:,:,1])
   pa = size(A,3)
   pc = size(C,3)
   (LinearAlgebra.checksquare(C[:,:,1]) == n && all([issymmetric(C[:,:,i]) for i in 1:pc])) ||
      throw(DimensionMismatch("all C[:,:,i] must be $n x $n symmetric matrices"))
   p = lcm(pa,pc)

   T = promote_type(T1, T2)
   T <: BlasFloat  || (T = promote_type(Float64,T))
   A1 = T1 == T ? A : A1 = convert(Array{T,3},A)
   C1 = T2 == T ? C : C1 = convert(Array{T,3},C)

   # Reduce A to Schur form and transform C
   AS, Q, _, KSCHUR = pschur(A1)

   #X = Q'*C*Q
   X = Array{T,3}(undef, n, n, p)

   for i = 1:p
       ia = mod(i-1,pa)+1
       ic = mod(i-1,pc)+1
       ia1 = mod(i,pa)+1

       X[:,:,i] = adj ? utqu(view(C1,:,:,ic),view(Q,:,:,ia)) :
                        utqu(view(C1,:,:,ic),view(Q,:,:,ia1))
   end
   # solve A'σXA - X + C = 0
   pdlyaps!(KSCHUR, AS, X; adj)

   #X <- Q*X*Q'
   for i = 1:p
       ia = mod(i-1,pa)+1
       utqu!(view(X,:,:,i),view(Q,:,:,ia)')
   end
   return X
end
"""
     pslyapdkr(A, C; adj = true) -> X

Solve the periodic discrete-time Lyapunov matrix equation

      A'σXA + C = X, if adj = true,

or

      A*X*A' + C =  σX, if adj = false,

where `σ` is the forward shift operator `σX(i) = X(i+1)`.
The Kronecker product expansion of equations is employed. `A` and `C` are
periodic square matrices, and `A` must not have characteristic multipliers on the unit circle.
This function is not recommended for large order matrices or large periods.
"""
function pslyapdkr(A::AbstractArray{T, 3}, C::AbstractArray{T, 3}; adj = true) where {T}
    m, n, pc = size(C)
    n == LinearAlgebra.checksquare(A[:,:,1])
    m == LinearAlgebra.checksquare(C[:,:,1])
    m == n  || throw(DimensionMismatch("A and C have incompatible dimensions"))
    pa = size(A,3)
    n2 = n*n
    p = lcm(pa,pc)
    N = p*n2
    R = zeros(T, N, N)
    adj ? copyto!(view(R,N-n2+1:N,1:n2),kron(A[:,:,pa]',A[:,:,pa]')) :
          copyto!(view(R,1:n2,N-n2+1:N),kron(A[:,:,pa],A[:,:,pa]))
    (i2, j2) = adj ? (n2, n2+n2) : (n2+n2, n2)
    for i = 1:p-1
        i1 = i2-n2+1
        j1 = j2-n2+1
        ia = mod(i-1,pa)+1
        adj ? copyto!(view(R,i1:i2,j1:j2),kron(A[:,:,ia]',A[:,:,ia]')) :
              copyto!(view(R,i1:i2,j1:j2),kron(A[:,:,ia],A[:,:,ia]))
        i2 += n2
        j2 += n2
    end
    indc = mod.(0:p-1,pc).+1
    return adj ? reshape((I-R) \ (C[:,:,indc][:]), n, n, p) :
    reshape((I-R) \ (C[:,:,circshift(indc,1)][:]), n, n, p)
end
function pdlyaps!(KSCHUR::Int, A::StridedArray{T1,3}, C::StridedArray{T1,3}; adj = true) where {T1<:BlasReal}
   n = LinearAlgebra.checksquare(A[:,:,1])
   pa = size(A,3)
   pc = size(C,3)
   (LinearAlgebra.checksquare(C[:,:,1]) == n && all([issymmetric(C[:,:,i]) for i in 1:pc])) ||
      throw(DimensionMismatch("all C[:,:,i] must be $n x $n symmetric matrices"))
   rem(pc,pa) == 0 || error("the period of C must be an integer multiple of A")
   (KSCHUR <= 0 || KSCHUR > pa ) &&
         error("KSCHUR has a value $KSCHUR, which is inconsistent with A ")

   if pa == 1 && pc == 1
      lyapds!(view(A,:,:,1), view(C,:,:,1); adj)
      return C[:,:,:]
   end
   ONE = one(T1)

   # determine the dimensions of the diagonal blocks of real Schur form

   G = Array{T1,3}(undef,2,2,pc)
   W = Matrix{Float64}(undef,2,14)
   WX = Matrix{Float64}(undef,4,5)
   ba, p = MatrixEquations.sfstruct(A[:,:,KSCHUR])
   if adj
      #
      # Solve    A(j)'*X(j+1)*A(j) + C(j) = X(j) .
      #
      # The (K,L)th blocks of X(j), j = 1, ..., p are determined
      # starting from upper-left corner column by column by
      #
      #   A(j)(K,K)'*X(j+1)(K,L)*A(j)(L,L) - X(j)(K,L) = -C(j)(K,L) - R(j)(K,L)
      #
      # where
      #                K              L-1
      #   R(j)(K,L) = SUM {A(j)(I,K)'*SUM [X(j+1)(I,J)*A(j)(J,L)]}
      #               I=1             J=1
      #
      #                 K-1
      #             +  {SUM [A(j)(I,K)'*X(j+1)(I,L)]}*A(j)(L,L)
      #                 I=1
      i = 1
      @inbounds  for kk = 1:p
          dk = ba[kk]
          k = i:i+dk-1
          j = 1
          ir = 1:i-1
          for ll = 1:kk
              dl = ba[ll]
              j1 = j+dl-1
              l = j:j1
              Ckl = view(C,k,l,1:pc)
              y = view(G,1:dk,1:dl,1:pc)
              copyto!(y,Ckl)
              if kk > 1
                 # C(j+1)[l,k] = C(j+1)[l,ir]*A(j)[ir,k]
                 ic = 1:j1
                 for ii = 1:pc
                     ia = mod(ii-1,pa)+1
                     ii1 = mod(ii,pc)+1
                     mul!(view(C,l,k,ii1),view(C,l,ir,ii1),view(A,ir,k,ia))
                     #y += C(j+1)[ic,k]'*A(j)[ic,l]
                     mul!(view(y,:,:,ii),transpose(view(C,ic,k,ii1)),view(A,ic,l,ia),ONE,ONE)
                 end
              end
              Ckl[:,:,:] .= dpsylv2(adj, dk, dl, KSCHUR, view(A,k,k,1:pa), view(A,l,l,1:pa), y, W, WX)
              if ll == kk && dl == 2
                 for ii = 1:pc
                     temp = 0.5*(Ckl[1,2,ii]+Ckl[2,1,ii])
                     Ckl[1,2,ii] = temp; Ckl[2,1,ii] = temp
                 end
              end
              j += dl
              if ll < kk
                 # C(j+1)[l,k] += C(j+1)[k,l]'*A(j)[k,k]
                 for ii = 1:pc
                     ia = mod(ii-1,pa)+1
                     ii1 = mod(ii,pc)+1
                     mul!(view(C,l,k,ii1),transpose(view(C,k,l,ii1)),view(A,k,k,ia),ONE,ONE)
                 end
              end
          end
          if kk > 1
             # C(j)[ir,k] = C(j)[k,ir]'
             for ii = 1:pc
                 transpose!(view(C,ir,k,ii),view(C,k,ir,ii))
             end
          end
          i += dk
      end
   else
      #
      # Solve    A(j)*X(j)*A(j)' + C(j) = X(j+1) .
      #
      # The (K,L)th block of X(j) is determined starting from
      # bottom-right corner column by column by
      #
      #    A(j)(K,K)*X(j)(K,L)*A(j)(L,L)' - X(j+1)(K,L) = -C(j)(K,L) - R(j)(K,L)
      #
      # Where
      #
      #                 N               N
      #    R(j)(K,L) = SUM {A(j)(K,I)* SUM [X(j)(I,J)*A(j)(L,J)']} +
      #                I=K            J=L+1
      #
      #                N
      #             { SUM [A(j)(K,J)*X(j)(J,L)]}*A(j)(L,L)'
      #              J=K+1
      j = n
      for ll = p:-1:1
          dl = ba[ll]
          l = j-dl+1:j
          i = n
          ir = j+1:n
          for kk = p:-1:ll
              dk = ba[kk]
              i1 = i-dk+1
              k = i1:i
              Clk = view(C,l,k,1:pc)
              y = view(G,1:dl,1:dk,1:pc)
              copyto!(y,Clk)
              if ll < p
                 ic = i1:n
                 for ii = 1:pc
                     ia = mod(ii-1,pa)+1
                     # C(j)[k,l] = C(j)[k,ir]*A(j)[l,ir]'
                     mul!(view(C,k,l,ii),view(C,k,ir,ii),transpose(view(A,l,ir,ia)))
                     # y += (A(j)[k,ic]*C(j)[ic,l])'
                     mul!(view(y,:,:,ii),transpose(view(C,ic,l,ii)),transpose(view(A,k,ic,ia)),ONE,ONE)
                 end
              end
              Clk[:,:,:] .= dpsylv2(adj, dl, dk, KSCHUR, view(A,l,l,1:pa), view(A,k,k,1:pa), y, W, WX)
              i -= dk
              if i >= j
                 for ii = 1:pc
                     ia = mod(ii-1,pa)+1
                     # C(j)[k,l] += (A(j)[l,l]*C(j)[l,k])'
                     mul!(view(C,k,l,ii),transpose(view(C,l,k,ii)),transpose(view(A,l,l,ia)),ONE,ONE)
                 end
              else
                 break
              end
          end
          if ll < p
             ir = i+2:n
             for ii = 1:pc
                 # C(j)[ir,l] = C(j)[l,ir]'
                 transpose!(view(C,ir,l,ii),view(C,l,ir,ii))
             end
          end
          j -= dl
      end
   end
   return C[:,:,:]
end
function dpsylv2(REV::Bool, N1::Int, N2::Int, KSCHUR::Int, TL::StridedArray{T,3}, TR::StridedArray{T,3},
                 B::StridedArray{T,3}, W::AbstractMatrix{T}, WX::AbstractMatrix{T}) where {T}
#     To solve for the N1-by-N2 matrices X_j, j = 1, ..., P,
#     1 <= N1,N2 <= 2, in the P simultaneous equations:

#     if REV = true

#       TL_j'*X_(j+1)*TR_j - X_j = B_j, X_(P+1) = X_1  (1)

#     or if REV = false

#       TL_j*X_j*TR_j' - X_(j+1) = B_j, X_(P+1) = X_1  (2)

#     where TL_j is N1 by N1, TR_j is N2 by N2, B_j is N1 by N2,
#     and ISGN = 1 or -1.

#     NOTE: This routine is primarily intended to be used in conjuntion
#           with solvers for periodic Lyapunov equations. Thus, both
#           TL and TR are formed from the diagonal blocks of the same
#           matrix in periodic real Schur form. W and WX are working matrices
#           allocated only once in the caller routine.
#           TL and TR are assumed to have the same period, but
#           B may have different period than TL and TR. The period
#           of B must be an integer multiple of that of TL and TR.
#           In the interests of speed, this routine does not
#                   check the inputs for errors.

#     METHOD

#     An initial approximation X_j, j=1, ..., P, is computed by reducing
#     the system (1) or (2) to an equivalent single Lyapunov equation.
#     Then, the accuracy of the solution is iteratively improved by
#     performing forward or backward sweeps depending on the stability
#     of eigenvalues. A maximum of 30 sweeps are performed.

#     REFERENCES

#     [1] A. Varga.
#         Periodic Lyapunov equations: some applications and new algorithms.
#         Int. J. Control, vol, 67, pp, 69-87, 1997.

	IND(J,P) = mod( J-1, P ) + 1
   P = size(TL,3)
   PB = size(B,3)
   i1 = 1:N1; i2 = 1:N2;
   Xw = view(WX,:,1:4)
   Yw = view(WX,:,5)
   Z = view(W,:,7:8)
   X = Array{T,3}(undef, N1, N2, PB)

   # Quick return if possible.
   (N1 == 0 || N2 == 0) && (return X)
   if P == 1 && PB == 1
      copyto!(view(X,i1,i2,1), view(B,i1,i2,1))
      MatrixEquations.lyapdsylv2!(REV, view(X, i1, i2, 1),  N1, N2, view(TL, i1, i1, 1), view(TR, i2, i2, 1), Xw, Yw)
      return X[:,:,:]
   end
   # partition working space
   AL = view(W,:,1:2)
   AR = view(W,:,3:4)
   Q = view(W,:,5:6)
   ZOLD = view(W,:,9:10)
   ATMP = view(W,:,11:12)
   BTMP = view(W,:,13:14)

   EPSM = 2*eps(T)

#     Define

#       AL(i,j) := TL_{j+i-1}*...*TL_{j+1}*TL_{j}, AL(j,j) := I;
#       AR(i,j) := TR_{j+i-1}*...*TR_{j+1}*TR_{j}, AR(j,j) := I;

#     If REV = true, compute an initial approximation for
#     X_{KSCHUR+1} = Z by solving

#       AL(KSCHUR+P+1,KSCHUR+1)'*Z*AR(KSCHUR+P+1,KSCHUR+1)  - Z = Qr(KSCHUR+1)

#     where

#       Qr(j)  = B_{j} + TL_{j}'*Qr(j+1)*TR_{j},  Qr(KSCHUR) = B(KSCHUR),
#                for j = KSCHUR+P-1, ..., KSCHUR+1

#     If REV = false, compute an initial approximation for
#     X_{KSCHUR} = Z by solving

#       AL(KSCHUR+P,KSCHUR)*Z*AR(KSCHUR+P,KSCHUR)'  - Z = Qf(KSCHUR)

#     where

#       Qf(j)  = B_{j} + TL_{j}*Qf(j-1)*TR_{j}',  Qf(KSCHUR) = B(KSCHUR),
#                for j = KSCHUR+1, ..., KSCHUR+P-1
   x = 1
   if REV
      L1 = KSCHUR+PB-1
      L2 = KSCHUR+1
      LSTEP = -1
   else
      L1 = KSCHUR+1
      L2 = KSCHUR+PB-1
      LSTEP = 1
   end
   AL[ i1, i1 ] = TL[ i1, i1, KSCHUR ]
   AR[ i2, i2 ] = TR[ i2, i2, KSCHUR ]
   Q[i1,i2] = B[i1,i2,KSCHUR]

   if  N1 == 1 && N2 == 1
       for JJ = L1:LSTEP:L2
           J = IND( JJ, P )
           JB = IND( JJ, PB )
           X11 =  TL[ 1, 1, J ]
           Y11 = TR[ 1, 1, J ]
           AL[ 1, 1 ] = AL[ 1, 1] * X11
           AR[ 1, 1 ] = AR[ 1, 1] * Y11
           Q[ 1, 1 ] = B[ 1, 1, JB] + X11 * Q[ 1, 1 ] * Y11
       end
   elseif N1 == 1 && N2 == 2
       for JJ = L1:LSTEP:L2
           J = IND( JJ, P )
           JB = IND( JJ, PB )
           X11 =  TL[ 1, 1, J ]
           AL[ 1, 1 ] = AL[ 1, 1] * X11
           Q[1,1] = X11*Q[1,1]
           Q[1,2] = X11*Q[1,2]

           Y11 = TR[ 1, 1, J ]
           Y22 = TR[ 2, 2, J ]
           Y12 = TR[ 1, 2, J ]
           if REV
               AR[1,2] = AR[1,1]*Y12 + AR[1,2]*Y22
               AR[1,1] = AR[1,1]*Y11
               AR[2,2] = AR[2,1]*Y12 + AR[2,2]*Y22
               AR[2,1] = AR[2,1]*Y11
               Q[1,2] = Q[1,1]*Y12 + Q[1,2]*Y22 + B[ 1, 2, JB ]
               Q[1,1] = Q[1,1]*Y11 + B[ 1, 1, JB ]
           else
               AR[1,1] = Y11*AR[1,1] + Y12*AR[2,1]
               AR[2,1] = Y22*AR[2,1]
               AR[1,2] = Y11*AR[1,2] + Y12*AR[2,2]
               AR[2,2] = Y22*AR[2,2]
               Q[1,1] = Q[1,1]*Y11 + Q[1,2]*Y12 + B[ 1, 1, JB ]
               Q[1,2] = Q[1,2]*Y22 + B[ 1, 2, JB ]
           end
       end
   elseif N1 == 2 && N2 == 1
       for JJ = L1:LSTEP:L2
           J = IND( JJ, P )
           JB = IND( JJ, PB )
           X11 = TL[ 1, 1, J ]
           X12 = TL[ 1, 2, J ]
           X22 = TL[ 2, 2, J ]
           if REV
              AL[1,2] = AL[1,1]*X12 + AL[1,2]*X22
              AL[1,1] = AL[1,1]*X11
              AL[2,2] = AL[2,1]*X12 + AL[2,2]*X22
              AL[2,1] = AL[2,1]*X11
              Q[2,1] = X12*Q[1,1] + X22*Q[2,1]
              Q[1,1] = X11*Q[1,1]
           else
              AL[1,1] = X11*AL[1,1] + X12*AL[2,1]
              AL[2,1] = X22*AL[2,1]
              AL[1,2] = X11*AL[1,2] + X12*AL[2,2]
              AL[2,2] = X22*AL[2,2]
              Q[1,1] = X11*Q[1,1] + X12*Q[2,1]
              Q[2,1] = X22*Q[2,1]
            end

            Y11 = TR[ 1, 1, J ]
            AR[1,1] = AR[1,1]*Y11
            Q[1,1] = Q[1,1]*Y11 + B[ 1, 1, JB ]
            Q[2,1] = Q[2,1]*Y11 + B[ 2, 1, JB ]
       end
   elseif N1 == 2 && N2 == 2
       for JJ = L1:LSTEP:L2
          J = IND( JJ, P )
          JB = IND( JJ, PB )
          X11 = TL[ 1, 1, J ]
          X12 = TL[ 1, 2, J ]
          X22 = TL[ 2, 2, J ]
          if REV
             AL[1,2] = AL[1,1]*X12 + AL[1,2]*X22
             AL[1,1] = AL[1,1]*X11
             AL[2,2] = AL[2,1]*X12 + AL[2,2]*X22
             AL[2,1] = AL[2,1]*X11
             Q[2,1] = X12*Q[1,1] + X22*Q[2,1]
             Q[1,1] = X11*Q[1,1]
             Q[2,2] = X12*Q[1,2] + X22*Q[2,2]
             Q[1,2] = X11*Q[1,2]
          else
             AL[1,1] = X11*AL[1,1] + X12*AL[2,1]
             AL[2,1] = X22*AL[2,1]
             AL[1,2] = X11*AL[1,2] + X12*AL[2,2]
             AL[2,2] = X22*AL[2,2]
             Q[1,1] = X11*Q[1,1] + X12*Q[2,1]
             Q[2,1] = X22*Q[2,1]
             Q[1,2] = X11*Q[1,2] + X12*Q[2,2]
             Q[2,2] = X22*Q[2,2]
          end

          Y11 = TR[ 1, 1, J ]
          Y12 = TR[ 1, 2, J ]
          Y22 = TR[ 2, 2, J ]
          if REV
             AR[1,2] = AR[1,1]*Y12 + AR[1,2]*Y22
             AR[1,1] = AR[1,1]*Y11
             AR[2,2] = AR[2,1]*Y12 + AR[2,2]*Y22
             AR[2,1] = AR[2,1]*Y11
             Q[1,2] = Q[1,1]*Y12 + Q[1,2]*Y22 + B[ 1, 2, JB ]
             Q[1,1] = Q[1,1]*Y11 + B[ 1, 1, JB ]
             Q[2,2] = Q[2,1]*Y12 + Q[2,2]*Y22 + B[ 2, 2, JB ]
             Q[2,1] = Q[2,1]*Y11 + B[ 2, 1, JB ]
          else
             AR[1,1] = Y11*AR[1,1] + Y12*AR[2,1]
             AR[2,1] = Y22*AR[2,1]
             AR[1,2] = Y11*AR[1,2] + Y12*AR[2,2]
             AR[2,2] = Y22*AR[2,2]
             Q[1,1] = Q[1,1]*Y11 + Q[1,2]*Y12 + B[ 1, 1, JB ]
             Q[1,2] = Q[1,2]*Y22 + B[ 1, 2, JB ]
             Q[2,1] = Q[2,1]*Y11 + Q[2,2]*Y12 + B[ 2, 1, JB ]
             Q[2,2] = Q[2,2]*Y22 + B[ 2, 2, JB ]
          end
       end
   end
   Z[i1,i2] = -Q[i1,i2]
#     Compute X_[KSCHUR+1] (if REV=.TRUE.) or X_KSCHUR (if REV=.FALSE.).
   MatrixEquations.lyapdsylv2!(REV, view(Z,i1,i2), N1, N2, view(AL,i1,i1), view(AR,i2,i2), Xw, Yw)
   XNORM = norm(view(Z,i1,i2), Inf)
   XNORM == 0 && (XNORM = one(T))

#     Determine the type of iteration to use

   if N1 == 1
      EL = abs( AL[1,1] )
   else
      EL = sqrt( abs( AL[1,1]*AL[2,2] - AL[1,2]*AL[2,1] ) )
   end
   if N2 == 1
      ER = abs( AR[1,1] )
   else
      ER = sqrt( abs( AR[1,1]*AR[2,2] - AR[1,2]*AR[2,1] ) )
   end
   STAB = EL*ER < one(T)

#     Save initial X_[KSCHUR+1].

   ZOLD[i1,i2] = Z[i1,i2]

#     Set iteration indices.

   if REV
      if STAB
          L1 = KSCHUR + PB
          L2 = KSCHUR + 1
          LSTEP = -1
          ITIND = -1
      else
          L1 = KSCHUR + 1
          L2 = KSCHUR + PB
          LSTEP = 1
          ITIND = 0
      end
   else
      if STAB
          L1 = KSCHUR
          L2 = KSCHUR + PB - 1
          LSTEP = 1
          ITIND = 0
      else
          L1 = KSCHUR + PB - 1
          L2 = KSCHUR
          LSTEP = -1
          ITIND = -1
      end
   end
   for ITER = 1:30
       if N1 == 1 && N2 == 1
          if STAB

#             Use direct recursion.

             for JJ = L1:LSTEP:L2
               J = IND( JJ, P )
               JB = IND( JJ, PB )
               J1 = IND( JJ+ITIND+1, PB )
               Z[1,1] = TL[ 1, 1, J ] * Z[1,1] * TR[ 1, 1, J ] - B[ 1, 1, JB ]
               X[ 1, 1, J1 ] = Z[ 1, 1 ]
             end
          else

#             Use inverse recursion.

             for JJ = L1:LSTEP:L2
                  J = IND( JJ, P )
                  JB = IND( JJ, PB )
                  J1 = IND( JJ+ITIND+1, PB )
                  Z[1,1] = (Z[1,1] + B[ 1, 1, JB ] ) / TL[ 1, 1, J ] / TR[ 1, 1, J ]
                  X[ 1, 1, J1 ] = Z[ 1, 1 ]
             end
          end
       elseif N1 == 1 && N2 == 2
          if STAB

#              Use direct recursion.

             for JJ = L1:LSTEP:L2
                 J = IND( JJ, P )
                 JB = IND( JJ, PB )
                 J1 = IND( JJ+ITIND+1, PB )
                 X11 = TL[ 1, 1, J ]
                 Y11 = TR[ 1, 1, J ]
                 Y22 = TR[ 2, 2, J ]
                 if REV
                    Y12 = TR[ 1, 2, J ]
                    Y21 = TR[ 2, 1, J ]
                 else
                    Y12 = TR[ 2, 1, J ]
                    Y21 = TR[ 1, 2, J ]
                 end

                 Z[1,1] = X11*Z[1,1]
                 Z[1,2] = X11*Z[1,2]
                 TEMP = Z[1,1]*Y11 + Z[1,2]*Y21
                 Z[1,2] = Z[1,1]*Y12 + Z[1,2]*Y22 - B[ 1, 2, JB ]
                 Z[1,1] = TEMP - B[ 1, 1, JB ]
                 X[ 1, 1, J1 ] = Z[ 1, 1 ]
                 X[ 1, 2, J1 ] = Z[ 1, 2 ]
              end
          else

#              Use inverse recursion.

             for JJ = L1:LSTEP:L2
                 J = IND( JJ, P )
                 JB = IND( JJ, PB )
                 J1 = IND( JJ+ITIND+1, PB )
                 BTMP[ 1, 1 ] = ( Z[1,1] + B[1,1,JB] ) / TL[ 1, 1, J ]
                 BTMP[ 2, 1 ] = ( Z[1,2] + B[1,2,JB] ) / TL[ 1, 1, J ]
                 ATMP[ 1, 1 ] = TR[ 1, 1, J ]
                 ATMP[ 2, 2 ] = TR[ 2, 2, J ]
                 if REV
                    ATMP[ 1, 2 ] = TR[ 2, 1, J ]
                    ATMP[ 2, 1 ] = TR[ 1, 2, J ]
                 else
                    ATMP[ 1, 2 ] = TR[ 1, 2, J ]
                    ATMP[ 2, 1 ] = TR[ 2, 1, J ]
                 end

                 #CALL DGESV( 2, 1, ATMP, 2, JPIV, BTMP, 2, INFO )
                 luslv!(ATMP,view(BTMP,1:2,1:1))
                 #BTMP = ATMP\BTMP
                 Z[ 1, 1 ] = BTMP[ 1, 1 ]
                 Z[ 1, 2 ] = BTMP[ 2, 1 ]
                 X[ 1, 1, J1 ] = Z[ 1, 1 ]
                 X[ 1, 2, J1 ] = Z[ 1, 2 ]
              end
          end
       elseif N1 == 2 && N2 == 1
          if STAB

#              Use direct recursion.

             for JJ = L1:LSTEP:L2
                 J = IND( JJ, P )
                 JB = IND( JJ, PB )
                 J1 = IND( JJ+ITIND+1, PB )
                 X11 = TL[ 1, 1, J ]
                 X22 = TL[ 2, 2, J ]
                 if REV
                    X12 = TL[ 2, 1, J ]
                    X21 = TL[ 1, 2, J ]
                 else
                    X12 = TL[ 1, 2, J ]
                    X21 = TL[ 2, 1, J ]
                 end
                 Y11 = TR[ 1, 1, J ]

                 TEMP = X11*Z[1,1] + X12*Z[2,1]
                 Z[2,1] = X21*Z[1,1] + X22*Z[2,1]
                 Z[1,1] = TEMP
                 Z[1,1] = Z[1,1]*Y11 - B[ 1, 1, JB ]
                 Z[2,1] = Z[2,1]*Y11 - B[ 2, 1, JB ]
                 X[ 1, 1, J1 ] = Z[ 1, 1 ]
                 X[ 2, 1, J1 ] = Z[ 2, 1 ]
              end
          else

#              Use inverse recursion.

             for JJ = L1:LSTEP:L2
                 J = IND( JJ, P )
                 JB = IND( JJ, PB )
                 J1 = IND( JJ+ITIND+1, PB )
                 Z[ 1, 1 ] = ( Z[1,1] + B[1,1,JB] ) / TR[ 1, 1, J ]
                 Z[ 2, 1 ] = ( Z[2,1] + B[2,1,JB] ) / TR[ 1, 1, J ]
                 ATMP[ 1, 1 ] = TL[ 1, 1, J ]
                 ATMP[ 2, 2 ] = TL[ 2, 2, J ]
                 if REV
                    ATMP[ 1, 2 ] = TL[ 2, 1, J ]
                    ATMP[ 2, 1 ] = TL[ 1, 2, J ]
                 else
                    ATMP[ 1, 2 ] = TL[ 1, 2, J ]
                    ATMP[ 2, 1 ] = TL[ 2, 1, J ]
                 end

                 #CALL DGESV( 2, 1, ATMP, 2, JPIV, Z, 2, INFO )
                 luslv!(ATMP,view(Z,1:2,1:1))
                 #Z = ATMP\Z
                 X[ 1, 1, J1 ] = Z[ 1, 1 ]
                 X[ 2, 1, J1 ] = Z[ 2, 1 ]
              end
          end
       elseif N1 == 2 && N2 == 2
          if STAB

#              Use direct recursion.

             for JJ = L1:LSTEP:L2
                 J = IND( JJ, P )
                 JB = IND( JJ, PB )
                 J1 = IND( JJ+ITIND+1, PB )
                 X11 = TL[ 1, 1, J ]
                 X22 = TL[ 2, 2, J ]
                 Y11 = TR[ 1, 1, J ]
                 Y22 = TR[ 2, 2, J ]
                 if REV
                    X12 = TL[ 2, 1, J ]
                    X21 = TL[ 1, 2, J ]
                    Y12 = TR[ 1, 2, J ]
                    Y21 = TR[ 2, 1, J ]
                 else
                    X12 = TL[ 1, 2, J ]
                    X21 = TL[ 2, 1, J ]
                    Y12 = TR[ 2, 1, J ]
                    Y21 = TR[ 1, 2, J ]
                 end

                 TEMP = X11*Z[1,1] + X12*Z[2,1]
                 Z[2,1] = X21*Z[1,1] + X22*Z[2,1]
                 Z[1,1] = TEMP
                 TEMP = X11*Z[1,2] + X12*Z[2,2]
                 Z[2,2] = X21*Z[1,2] + X22*Z[2,2]
                 Z[1,2] = TEMP
                 TEMP = Z[1,1]*Y11 + Z[1,2]*Y21
                 Z[1,2] = Z[1,1]*Y12 + Z[1,2]*Y22 - B[ 1, 2, JB ]
                 Z[1,1] = TEMP - B[ 1, 1, JB ]
                 TEMP = Z[2,1]*Y11 + Z[2,2]*Y21
                 Z[2,2] = Z[2,1]*Y12 + Z[2,2]*Y22 - B[ 2, 2, JB ]
                 Z[2,1] = TEMP - B[ 2, 1, JB ]
                 X[ 1, 1, J1 ] = Z[ 1, 1 ]
                 X[ 2, 1, J1 ] = Z[ 2, 1 ]
                 X[ 1, 2, J1 ] = Z[ 1, 2 ]
                 X[ 2, 2, J1 ] = Z[ 2, 2 ]
              end
          else

#              Use inverse recursion.

             for JJ = L1:LSTEP:L2
                 J = IND( JJ, P )
                 JB = IND( JJ, PB )
                 J1 = IND( JJ+ITIND+1, PB )
                 BTMP[ 1, 1 ] = Z[ 1, 1 ] + B[ 1, 1, JB ]
                 BTMP[ 2, 1 ] = Z[ 1, 2 ] + B[ 1, 2, JB ]
                 BTMP[ 1, 2 ] = Z[ 2, 1 ] + B[ 2, 1, JB ]
                 BTMP[ 2, 2 ] = Z[ 2, 2 ] + B[ 2, 2, JB ]
                 ATMP[ 1, 1 ] = TR[ 1, 1, J ]
                 ATMP[ 2, 2 ] = TR[ 2, 2, J ]
                 if REV
                    ATMP[ 1, 2 ] = TR[ 2, 1, J ]
                    ATMP[ 2, 1 ] = TR[ 1, 2, J ]
                 else
                    ATMP[ 1, 2 ] = TR[ 1, 2, J ]
                    ATMP[ 2, 1 ] = TR[ 2, 1, J ]
                 end

                 #CALL DGESV( 2, 2, ATMP, 2, JPIV, BTMP, 2, INFO )
                 luslv!(ATMP,view(BTMP,1:2,1:2))
                 #BTMP = ATMP\BTMP
                 Z[ 1, 1 ] = BTMP[ 1, 1 ]
                 Z[ 1, 2 ] = BTMP[ 2, 1 ]
                 Z[ 2, 1 ] = BTMP[ 1, 2 ]
                 Z[ 2, 2 ] = BTMP[ 2, 2 ]

                 ATMP[ 1, 1 ] = TL[ 1, 1, J ]
                 ATMP[ 2, 2 ] = TL[ 2, 2, J ]
                 if REV
                    ATMP[ 1, 2 ] = TL[ 2, 1, J ]
                    ATMP[ 2, 1 ] = TL[ 1, 2, J ]
                 else
                    ATMP[ 1, 2 ] = TL[ 1, 2, J ]
                    ATMP[ 2, 1 ] = TL[ 2, 1, J ]
                 end

                 #CALL DGESV( 2, 2, ATMP, 2, JPIV, Z, 2, INFO )
                 luslv!(ATMP,view(Z,1:2,1:2))
                 #Z = ATMP\Z
                 X[ 1, 1, J1 ] = Z[ 1, 1 ]
                 X[ 1, 2, J1 ] = Z[ 1, 2 ]
                 X[ 2, 1, J1 ] = Z[ 2, 1 ]
                 X[ 2, 2, J1 ] = Z[ 2, 2 ]
              end
          end
       end
       DNORM = zero(T)
       for J = 1:N2
          for I = 1:N1
             DNORM = max( DNORM, abs( ZOLD[I,J] - Z[I,J] ) )
             ZOLD[ I, J ] = Z[ I, J ]
          end
       end
       #println("XNORM = $XNORM DNORM = $DNORM")
       DNORM <= EPSM*XNORM && (return -X)
   end
   @warn "iterative process not converging: solution may be inaccurate"
   return -X
#   END of DPSYLV2
end
@inline function luslv!(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T
   #
   #  fail = luslv!(A,B)
   #
   # This function is a speed-oriented implementation of a Gaussion-elimination based
   # solver of small order linear equations of the form A*X = B. The computed solution X
   # overwrites the vector B, while the resulting A contains in its upper triangular part,
   # the upper triangular factor U of its LU decomposition.
   # The diagnostic output parameter fail, of type Bool, is set to false in the case
   # of normal return or is set to true if the exact singularity of A is detected
   # or if the resulting B has non-finite components.
   #
   n, m = size(B)
   @inbounds begin
         for k = 1:n
            # find index max
            kp = k
            if k < n
                amax = abs(A[k, k])
                for i = k+1:n
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            iszero(A[kp,k]) && return true
            if k != kp
               # Interchange
               for i = 1:n
                   tmp = A[k,i]
                   A[k,i] = A[kp,i]
                   A[kp,i] = tmp
               end
               for j = 1:m
                   tmp = B[k,j]
                   B[k,j] = B[kp,j]
                   B[kp,j] = tmp
               end
            end
            # Scale first column
            Akkinv = inv(A[k,k])
            i1 = k+1:n
            Ak = view(A,i1,k)
            rmul!(Ak,Akkinv)
            # Update the rest of A and B
            for j = k+1:n
                axpy!(-A[k,j],Ak,view(A,i1,j))
            end
            for j = 1:m
                axpy!(-B[k,j],Ak,view(B,i1,j))
            end
         end
         ldiv!(UpperTriangular(A), view(B,:,:))
         return any(!isfinite, B)
   end
end


text:
# invert all the Svalbard glaciers

@everywhere debugging = false
debugging && println("\n\nDEBUG RUN !!!!!!!!\n\n")
include("svalbard-setup.jl")

@everywhere begin
    using BSON
    include("svalbard-setup.jl")
    # runs fast if true

    update_cache = true
    svalbard_nrs = debugging ? svalbard_nrs[1:20] : svalbard_nrs

    import Plots, StatPlots
    import Base.Test: @inferred
    import DataStructures: OrderedDict
    Plots.pyplot()
    #Plots.gr()
    # if strip(readstring(`hostname`))=="x230"
    #     pyplot()
    # else
    #     ENV["GKSwstype"]="100" # for headless/no display
    #     gr()
    # end
    import JLD2,JLD
    import FileIO
    import KissMCMC

    # Run it
    ########
    n_1d_theta = 3

    mcmc_niter = debugging ? 10^3 : 5*10^5
    # parameters to play with
    srand(5)

    sigma_h_model = [0.0, 10.0]
    sigma_iv_model = [0.0, 10.0]
    pm0 = BM.MPara()

    n_1d_btilde=3
    n_1d_fsl=3
    n_1d_temp=1

    fit_target = [BM.FitTarget.h, BM.FitTarget.h_iv, BM.FitTarget.length][2]
    runtyp = ["test", "testmid", "prodlow", "prod"][4]

    # used for plotting:
    theta0_with_h_iv_dict = OrderedDict(:btilde=> zeros(BM.F,n_1d_theta), # NOTE in m/y (error)
                                        :fsl=> zeros(BM.F,n_1d_theta), # error
                                        :dist_exp=> pm0.dist_exp,
                                        :h => zeros(BM.F,n_1d_theta), # for plotting
                                        :iv => zeros(BM.F,n_1d_theta), # for plotting
                                        )
    theta0_with_h_dict = deepcopy(theta0_with_h_iv_dict)
    pop!(theta0_with_h_dict, :iv)
    theta0_with_iv_dict = deepcopy(theta0_with_h_iv_dict)
    pop!(theta0_with_iv_dict, :h)

    dir_output = "output/svalbard/04dec/"*runtyp
end

!isdir(dir_output) && mkpath(dir_output)

@everywhere begin
    function var2filename(gid, sigma_h_model, sigma_iv_model)
        out = dir_output*"$(BM.getrgi(gid))_sh=$(sigma_h_model)_siv=$(sigma_iv_model)"
    end
end

#@sync @parallel for nr in svalbard_nrs[1:end÷2]
#@sync @parallel for nr in svalbard_nrs[end÷2+1:end]
@sync @parallel for nr in svalbard_nrs
    for sh in sigma_h_model, siv in sigma_iv_model
        gid = BM.RGIGlacier(nr, region_svalbard)
        fl = var2filename(gid, sh, siv)
        if debugging
            println(fl)
            tic()
        end
        try
            dt, pl = init_svalbard_gl_dt(gid)
            dt.para.opts[:pm][:sigma_h_model] = sh
            dt.para.opts[:pm][:sigma_iv_model] = siv
            gl,gb,pp,pm,pn,pl,dt = init_svalbard_gl(dt, pl)

            # Inversion
            ###########
            th0d = BM.theta0_dict_defaults(pm, n_1d_btilde, n_1d_fsl, n_1d_temp)

            (theta0, logposterior, logposterior1d, logprior, logprior_debug,
             loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
             pmcmc_defaults, fit_target) =
                BM.init_inverse(gb, pp, pm, pn,
                                runtyp=Symbol(runtyp),
                                theta0_dict=th0d,
                                fit_target=fit_target)
            pmcmc = BM.MCMCNum(;pmcmc_defaults...)

            # # Make not continuous glacier continuous:
            # if isinf(logprior_debug(theta0.th0,false)[2][2])
            #     theta0.th0[theta0.names[:btilde]] += 1
            # end

            if debugging
                @inferred logposterior(theta0.th0)
                print("Time for logposterior:  ")
                @time logposterior(theta0.th0)
                println("Value logposterior:  $(logposterior(theta0.th0)[1])")
            end

            res = if debugging
                print("Time for $(pmcmc.niter) MCMC steps:  ")
                @time BM.mcmc(logposterior, theta0, pmcmc;
                              verbose=false, use_progress_meter=false)
            else
                BM.mcmc(logposterior, theta0, pmcmc;
                        verbose=false, use_progress_meter=false)
            end
            # # h and iv averaged to where parameters are
            # h_th, iv_th = BM.h_iv_at_thetas(theta0, mc_hs1d, mc_ivs1d)

            # theta0_h_iv = BM.Theta0(theta0_with_h_iv_dict, gb, pp, pm, pn)
            # theta0_h = BM.Theta0(theta0_with_h_dict, gb, pp, pm, pn)
            # theta0_iv = BM.Theta0(theta0_with_iv_dict, gb, pp, pm, pn)

            # varnames = BM.get_varnames(theta0)
            # varnames_h_iv = BM.get_varnames(theta0_h_iv)
            # varnames_h = BM.get_varnames(theta0_h)
            # varnames_iv = BM.get_varnames(theta0_iv)

            # thetas_h_iv = vcat(thetas, h_th, iv_th)
            # thetas_h = vcat(thetas, h_th)
            # thetas_iv = vcat(thetas, iv_th)

            # #print_results(thetas, accept_ratio, names=varnames)

            # fig = BM.plottheta(thetas, theta0, toplot=:all, size=(1400,1400))
            # Plots.savefig(fig, fl*"--thetas.png")

            # fig = BM.plottheta(thetas_h, theta0_h, toplot=:all, size=(1400,1400))
            # Plots.savefig(fig, fl*"--thetas_h.png")

            # fig = BM.plottheta(thetas_iv, theta0_iv, toplot=:all, size=(1400,1400))
            # Plots.savefig(fig, fl*"--thetas_iv.png")

            # # more plots
            # # 1D results with bands
            # fig = BM.plotinv1d(gb, blobs)
            # Plots.savefig(fig, fl*"--1d.png")

            # # 2D error binned back to 1D
            # fig = BM.plotinv1d_err(gb, blobs)
            # Plots.savefig(fig, fl*"--1d_err.png")

            # # lots of 2D plots
            # fig = BM.plotinv2d(gl, blobs, size=(1400,1400))
            # Plots.savefig(fig, fl*"--2d_overview.png")

            # # error of mean IV
            # fig = BM.plotinv2d_iverr(gl, blobs, reuse=false)
            # Plots.savefig(fig, fl*"--2d_iverr.png")

            # # summaries
            # # as JLD
            # res = KissMCMC.summarize_run(thetas_h_iv, names=varnames_h_iv)
            # FileIO.save(fl*".jld2", "res", res)

            # # as text
            # fll = splitdir(fl)[2]
            # ar = hcat([[strip(string(v)),m,s] for (v,m,s) in zip(res[:var],res[:mean],res[:std])]...)
            # l1,l2,l3 = split(reprmime("text/plain", ar),'\n')[2:4]
            # write(fl*"--summary", "$fll\n    $l1\n    $l2\n    $l3\n")

            # ## forward
            # ## run forward model with most likely paras
            # res = KissMCMC.summarize_run(thetas, names=varnames_h_iv)

            # hs1d, ivs1d, hs2d, ivs2d, pm_,
            # taus1d, taus1d_l, taus2d = fwdm_fn(res[:mode])
            # fig = BM.plot_run(gb, hs1d, taus1d, ivs1d, hs2d, taus2d,
            #                   ivs2d, taus1d_l, pp, pm,
            #                   plot2d=[:h,:iv,:ivlog][1])
            # Plots.savefig(fig, fl*"--map.png")

            # FileIO.save(fl*"-map2d.jld2",
            #             "mh2d", mh2d, "sh2d", sh2d,
            #             "miv2d", miv2d, "siv2d", siv2d,
            #             "min_h2d", min_h2d, "max_h2d", max_h2d,
            #             "min_iv2d", min_iv2d, "max_iv2d", max_iv2d)

        @time bson(fl*".bson", Dict(:res => res)) # takes 2 minutes!
        catch e
            e==InterruptException() && rethrow(e)
            fll = splitdir(fl)[2]
            write(fl*"--summary", "$fll\n    error $e\n")
            println("error: $e")
        end
        if debugging
            toc()
            println("\n")
        end
    end
end


text:
module EHRAuthentication

include("Cerner/Cerner.jl")

end # module


text:
module PointDefects

   using JuLIP: deleteat!, JVecF, AbstractAtoms
   using JuLIP.ASE: rnn, extend!
   import JuLIPMaterials: cluster

   export PointDefect, Vacancy, SelfInterstitial

   abstract PointDefect
   immutable Vacancy <: PointDefect end
   immutable SelfInterstitial <: PointDefect end

   function cluster(species::AbstractString, R::Real, ::Type{Vacancy})
      at = cluster(species, R)
      deleteat!(at, 1)
      return at
   end

   interstitial_positions = Dict(
      "Si" => JVecF([-2.53312, 0.0, 0.0])
   )

   function cluster(species::AbstractString, R::Real, ::Type{SelfInterstitial})
      at = cluster(species, R)
      extend!(at, species, interstitial_positions[species])
      return at
   end

end


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
  Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, PassiveTracerFlows

makedocs(
   modules = [PassiveTracerFlows],
   doctest = false,
     clean = true,
 checkdocs = :all,
    format = :html,
   authors = "Gregory L. Wagner and Navid C. Constantinou",
  sitename = "PassiveTracerFlows.jl",
     pages = Any[
              "Home" => "index.md",
              "Modules" => Any[
                "modules/traceradvdiff.md",
              ],
              "DocStrings" => Any[
              "man/types.md",
              "man/functions.md"]
             ]
)

deploydocs( repo = "github.com/FourierFlows/PassiveTracerFlows.jl.git" )


text:
module SimulationParameters

include("physics.jl")

export Δt, nsteps, L, nx, ny, Δx, Δy, φ, K, μ, Pin, Pout,
       P₀, M, ψ, ψ₀, bc, inlet, wall, relax_steps

## CFL condition: dt ∼ dx² !!

Δt = 0.01
# Δt = 0.005
# Δt = 0.0001
# Δt = 0.0001
# nsteps = 10000
# nsteps = 50000
# nsteps = 25000

# nsteps = 500000
# nsteps = 320000
nsteps = 1000

relax_steps = 100

## Space grid [0, 2] × [0, 2]
L = 2
nx = ny = 20
# nx = ny = 30
# nx = ny = 50
# nx = ny = 150
Δx = Δy = L / nx

## Porosity
φ = 0.7

## Specific permeability
K = 1e-12

## Dynamic viscocity
μ = [1.8e-5, 2.14e-4]

## Pressures at top and bottom, and initial pressure
Pin = 1e6
Pout = 1e5
P₀ = Pout

## Molar masses
M = [MOLAR_MASS_IGAS, MOLAR_MASS_PENTANE]

## ψ = ν₁ / ν₂ = m₁ / (m₂ M₁ / M₂ + m₁)
## Molar Composition on the Inlet
ψ = 0.9

## Initial Molar Composition
ψ₀ = 0.1

## -------------------------- BC -------------------------- ##

# inlet = 1 : nx ÷ 2
# wall  = nx ÷ 2 + 1: nx
inlet = 1 : nx
wall  = nx: nx

a_u = zeros(max(nx, ny), 4, 2)
a_v = zeros(max(nx, ny), 4, 2)
a_ρ = zeros(max(nx, ny), 4, 2)
a_s = zeros(max(nx, ny), 4, 2)
a_P = zeros(max(nx, ny), 4, 2)

b_u = zeros(max(nx, ny), 4, 2)
b_v = zeros(max(nx, ny), 4, 2)
b_ρ = zeros(max(nx, ny), 4, 2)
b_s = zeros(max(nx, ny), 4, 2)
b_P = zeros(max(nx, ny), 4, 2)

## h[:][1][k] -> left
##
## h[:][2][k] -> right
##
## h[:][3][k] -> top
##
## h[:][4][k] -> bottom

k = 1 : 2

## top & bottom
@views a_v[:, :, k]     .= 1

## left & right
@views a_u[:, :, k]     .= 1

## left & right & top & bottom
@views a_ρ[:, :, k]     .= 1

## inlet at bottom [1, nx ÷ 2]
@views a_s[inlet, 4, k] .= 1

## walls & outlet
@views b_s[:, 1:3,  k]  .= 1
@views b_s[wall, 4, k]  .= 1

## outlet
@views a_P[:, 3, k]     .= 1

## inlet
@views a_P[inlet, 4, k] .= 1

## walls left & right
@views b_P[:, 1:2, k]   .= 1

## bottom wall
@views b_P[wall, 4, k]  .= 1

u_bc = BoundaryCondition(a_u, b_u)
v_bc = BoundaryCondition(a_v, b_v)
ρ_bc = BoundaryCondition(a_ρ, b_ρ)
s_bc = BoundaryCondition(a_s, b_s)
P_bc = BoundaryCondition(a_P, b_P)

bc = BoundaryConditions(u_bc, v_bc, ρ_bc, s_bc, P_bc)
end

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
module HopfieldNets
    export HopfieldNet, DiscreteHopfieldNet, ContinuousHopfieldNet
    export update!, energy, settle!, train!, associate!, storkeytrain!

    include("includes.jl")
end


text:
include("../src/pvtproperties.jl")

@testset "Gas PVT" begin

    @test LeeGasViscosity(0.7, 2000.0, 160.0, 0.75) ≈ 0.017 atol = 0.001
    @test HankinsonWithWichertPseudoCriticalTemp(0.65, 0.05, 0.08) .- (370.2, 351.4, 18.09) |> x -> all(abs.(x) .<= 1)
    @test HankinsonWithWichertPseudoCriticalPressure(0.65, 0.05, 0.08) ≈ 634.0 atol = 2
    @test PapayZFactor(634, 351.4, 600.0, 100.0) ≈ 0.921 atol = 0.01
    @test KareemEtAlZFactor(663.29, 377.59, 2000, 150) ≈ 0.8242 atol = 0.05
    @test KareemEtAlZFactor_simplified(663.29, 377.59, 2000, 150) ≈ 0.8242 atol = 0.15
    @test gasVolumeFactor(1200, 0.85, 200) ≈ 0.013 atol = 0.001
    @test gasDensity_insitu(0.916, 0.883, 346.6, 80.3) ≈ 1.79 atol = 0.01

end

@testset "Oil PVT" begin

    @test StandingSolutionGOR(30.0, 0.6, 800.0, 120.0, nothing, 2650) ≈ 121.4 atol = 1
    @test StandingSolutionGOR(41.06, 0.916, 346.6, 80.3, nothing, 2650) ≈ 109.7 atol = 1
    @test StandingBubblePoint(30, 0.75, 120, 100) ≈ 614 atol = 1
    @test StandingOilVolumeFactor(30.0, 0.6, 121.4, 800.0, 120.0) ≈ 1.07 atol = 0.05
    @test oilDensity_insitu(41.06,  0.916,  109.7,  1.05) ≈ 49.9 atol = 0.5
    @test BeggsAndRobinsonDeadOilViscosity( 37.9,  120) ≈ 4.05 atol = 0.05
    @test GlasoDeadOilViscosity(37.9, 120) ≈ 2.30 atol = 0.1
    @test ChewAndConnallySaturatedOilViscosity(2.3, 769.0) ≈ 0.638 atol = 0.1

end

@testset "Water PVT" begin

    @test waterDensity_stb(1.12) ≈ 69.9 atol = 0.1
    @test GouldWaterVolumeFactor(50.0, 120.0) ≈ 1.01 atol = 0.05

end


text:
using Pkg
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "Juice.jl tests" begin

end


text:
module CommunalHelperMusicParamTrigger

using ..Ahorn, Maple

@mapdef Trigger "CommunalHelper/MusicParamTrigger" MusicParamTrigger(
    x::Integer,
    y::Integer,
    width::Integer = Maple.defaultTriggerWidth,
    height::Integer = Maple.defaultTriggerHeight,
    param::String="",
    enterValue::Number=1.0,
    exitValue::Number=0.0,
)

const placements = Ahorn.PlacementDict(
    "Music Paramater Trigger (Communal Helper)" => Ahorn.EntityPlacement(
        MusicParamTrigger,
        "rectangle",
    ),
)

end


text:
function getgenericHamiltoniantest( Id,
    Ix,
    Iy_no_im,
    Iz,
    ω0::Vector{T},
    J_vals::Vector{T},
    J_inds::Vector{Tuple{Int,Int}}) where T <: Real
#
    N = length(ω0)
    N_couplings = length(J_vals)
    @assert N_couplings == length(J_inds)

    # pair-wise terms.
    H1 = zeros(T, 2^N, 2^N)
    for i = 1:N_couplings

        j = J_inds[i][1]
        k = J_inds[i][2]

        H1 += getmultiIjIk(Id, Ix, Iy_no_im, Iz, j, k, N, 2*π*J_vals[i])
    end

    # single terms.
    ωIzs = collect( ω0[j] .* Iz for j = 1:N )
    H2 = Kronecker.kroneckersum(ωIzs...)

    return H1+H2, H1, H2
end

function getgenericHamiltoniantest( Id,
    Ix,
    Iy_no_im,
    Iz,
    ω0::Vector{T},
    J::Matrix{T}) where T <: Real

    N = length(ω0)
    @assert size(J,1) == size(J,2) == N

    # pair-wise terms.
    H1 = zeros(T, 2^N, 2^N)
    for j = 1:N
        for k = j+1:N
            #H += (2*π*J[j,k]) .* getmultiIjIk(Id, Ix, Iy_no_im, Iz, j, k, N)
            H1 += getmultiIjIk(Id, Ix, Iy_no_im, Iz, j, k, N, 2*π*J[j,k])

        end
    end

    # single terms.
    ωIzs = collect( ω0[j] .* Iz for j = 1:N )
    H2 = Kronecker.kroneckersum(ωIzs...)

    return H1+H2, H1, H2
end

# Only the upper-right off-diagonals of J is used.
# No secular dipole-dipole.
function getgenericHamiltonian( Id,
    Ix,
    Iy_no_im,
    Iz,
    ω0::Vector{T},
    J::Matrix{T}) where T <: Real

    N = length(ω0)
    @assert size(J,1) == size(J,2) == N

    # pair-wise terms.
    H = zeros(T, 2^N, 2^N)
    for j = 1:N
        for k = j+1:N
            #H += (2*π*J[j,k]) .* getmultiIjIk(Id, Ix, Iy_no_im, Iz, j, k, N)
            H += getmultiIjIk(Id, Ix, Iy_no_im, Iz, j, k, N, 2*π*J[j,k])

        end
    end

    # single terms.
    ωIzs = collect( ω0[j] .* Iz for j = 1:N )
    H += Kronecker.kroneckersum(ωIzs...)

    return H
end

# No secular dipole-dipole.
function getgenericHamiltonian0( Id,
    Ix,
    Iy_no_im,
    Iz,
    ω0::Vector{T},
    J::Matrix{T}) where T <: Real

    N = length(ω0)
    @assert size(J,1) == size(J,2) == N

    H = zeros(T, 2^N, 2^N)

    for j = 1:N
        H += ω0[j] .* getmultiIz(Id, j, N)

        for k = j+1:N
            H += (2*π*J[j,k]) .* getmultiIjIk0(Id, Ix, Iy_no_im, Iz, j, k, N)
        end
    end

    return H
end

# No secular dipole-dipole.
function getgenericHamiltonian( Id,
    Ix,
    Iy_no_im,
    Iz,
    ω0::Vector{T},
    J_vals::Vector{T},
    J_inds::Vector{Tuple{Int,Int}}) where T <: Real

    N = length(ω0)
    N_couplings = length(J_vals)
    @assert N_couplings == length(J_inds)

    # pair-wise terms.
    H = zeros(T, 2^N, 2^N)
    for i = 1:N_couplings

        j = J_inds[i][1]
        k = J_inds[i][2]

        H += getmultiIjIk(Id, Ix, Iy_no_im, Iz, j, k, N, 2*π*J_vals[i])
    end

    # single terms.
    ωIzs = collect( ω0[j] .* Iz for j = 1:N )
    H += Kronecker.kroneckersum(ωIzs...)

    return H
end

# has secular dipole-dipole. J here is not really J-coupling, but J - d.
function getgenericHamiltonian( Id,
    Ix,
    Iy_no_im,
    Iz,
    ω0::Vector{T},
    J_vals::Vector{T},
    J_inds::Vector{Tuple{Int,Int}},
    Di_vals::Vector{T},
    Di_inds::Vector{Tuple{Int,Int}}) where T <: Real

    N = length(ω0)
    N_couplings = length(Di_vals)
    @assert N_couplings == length(Di_inds)

    # Allocate base operators.
    z_operators = Vector{Matrix{T}}(undef, N)
    fill!(z_operators, Id)

    # pair-wise terms.
    H = getgenericHamiltonian(Id, Ix, Iy_no_im, Iz, ω0, J_vals, J_inds)
    for i = 1:N_couplings

        j = Di_inds[i][1]
        k = Di_inds[i][2]

        ### z.
        z_operators[j] = (3*Di_vals[i]) .* Iz
        z_operators[k] = Iz

        # get product operator.
        H_di = Kronecker.kronecker(z_operators...)

        H += H_di
    end

    return H
end

# v must be an eigen vector of Iz. This routine does not verify this condition.
function getzangularmomentum(v::Vector{T}, Iz; tol = 1e-7)::T where T

    b = Iz*v
    m = dot(b,v)/dot(v,v)

    ## This was meant to be a sanity-check, but when N > 4, the tol becomes larger than 1e-7.
    # The exclusion of these states mean there are insufficient number of splits,
    #   thus, the area ratio does not hold, unless weak coupling is "really satisifed"?
    # Solution: just disable this check, and assume v must be an eigenvector of Iz.
    # if norm(m .* v - b) > tol
    #     return -Inf
    # end

    return m
end

function getorderofcoherence(Iz, Q::Vector{Vector{T}}) where T

    M = length(Q)
    @assert size(Iz, 1) == M == size(Iz, 2)

    M_array = collect( getzangularmomentum(Q[r], Iz) for r = 1:length(Q))

    p = Matrix{T}(undef, M, M)
    fill!(p, NaN)
    for r = 1:M
        #Mr = getzangularmomentum(Q[r], Iz)
        Mr = M_array[r]

        for s = 1:M
            #Ms = getzangularmomentum(Q[s], Iz)
            Ms = M_array[s]

            p[r,s] = Mr-Ms
            #p[r,s] = round(Mr-Ms)
        end
    end

    return p, M_array
end

function normalizetoNspins!(A::Vector{T}, N_spins::Int) where T

    factor = N_spins/sum(A)
    for l = 1:length(A)
        A[l] = A[l] * factor
    end

    return nothing
end

function getaΩ( Iz_full,
    H::Matrix{T},
    Iys_no_im_full,
    Ix_full;
    tol = 1e-7) where T <: Real

    M = size(H,1)
    @assert size(H,2) == M == size(Iz_full,1) == size(Iz_full,2)

    ### eigen states.
    λ, Q_mat = eigen(H)
    # λ[6] = -λ[6]
    # Q_mat[:,6] = -Q_mat[:,6]
    Q = collect( Q_mat[:,m] for m = 1:M )

    ### pre-compute.
    IyQs = collect( Iys_no_im_full*Q[s] for s = 1:size(Q,1)) # no im since official formula for α is div by 2im.
    IxQs = collect( Ix_full*Q[r] for r = 1:size(Q,1))

    ### get order of coherence.
    #@time p, M_array = getorderofcoherence(Iz_full, Q) # TODO speed this up.
    #println("timing coherence")
    #@time p, M_array = getorderofcoherence(Iz_full, Q)
    p, M_array = getorderofcoherence(Iz_full, Q)

    a = Vector{T}(undef, M*M)
    F = Vector{T}(undef, M*M)
    inds = Vector{Tuple{Int,Int}}(undef, M*M)

    i_out = 0

    # for r = 1:M
    #     for s = 1:M
    for (r,s) in Base.Iterators.product(1:M,1:M)

        #if Utilities.isnumericallyclose(p[r,s], -one(T), tol)
        if isnumericallyclose(round(p[r,s]), -one(T), tol)
            #println("i_out = ", i_out)
            i_out += 1

            inds[i_out] = (r,s)

            ## signal amplitude. Section A.8.2 Spin Dynamics, Eq 2.1.131 Ernst.
            t1 = -dot(conj(Q[r]), IyQs[s])
            t2 = dot(conj(Q[s]), IxQs[r])

            # factor 2 instead of 2*im since Iys_no_im_full
            #   is divided by 2 instead of 2*im.
            a[i_out] = 2*(t1 * t2)

            ## signal frequency.
            F[i_out] = λ[s] - λ[r]
        end
    end

    resize!(inds, i_out)
    resize!(a, i_out)
    resize!(F, i_out)

    return a, F, p, λ, Q, inds, M_array
end

function getcfreq(z::Vector{T}) where T
    M = length(z)

    F = Matrix{T}(undef, M, M)
    for r = 1:M
        for s = 1:M
            F[r,s] = z[s] - z[r]
        end
    end

    return F
end

# Q is operator. Evalates <r| Q |s>.
function evaltrace( Q::Kronecker.KroneckerSum{T}, r, s)::T where T
    return dot(conj(r),Q*s)
end

function evaltrace( Q::Kronecker.KroneckerSum{T},
    V)::Matrix{T} where T

    M = length(V)

    out = Matrix{T}(undef, M, M)
    for r = 1:M
    for s = 1:M
    out[r,s] = evaltrace(Q, V[r], V[s])
    end
    end

    return out
end

# tp_basis[r,:] is the elementary state vector for each of the spins.
# +1 for state 1/2, -1 for state -1/2.
function gettensorproductbasis(Iz::Matrix{T}, Id, N_spins::Int) where T
    tp_basis = Matrix{Int}(undef, 2^N_spins, N_spins)

    for n = 1:N_spins

    Is = collect( Id for j = 1:N_spins )
    Is[n] = Iz
    I_full = Kronecker.kronecker(Is...)

        for i = 1:size(I_full,1)
            m = I_full[i,i]

            if m > 0
                tp_basis[i,n] = 1
            else
                tp_basis[i,n] = -1
            end
        end
    end

    return tp_basis
end

# for single (π/2)_x pulse. # equation A.32.
function evalamplitudeA32( Iys_no_im_full,
    Ix_full,
    V::Vector{Vector{T}})::Matrix{T} where T

    M = length(V)

    out = Matrix{T}(undef, M, M)
    for r = 1:M
    for s = 1:M

    t1 = -dot(conj(V[r]), Iys_no_im_full*V[s])
    t2 = dot(conj(V[s]), Ix_full*V[r])

    # factor 2 instead of 2*im since Iys_no_im_full
    #   is divided by 2 instead of 2*im.
    out[r,s] = 2*(t1 * t2)
    end
    end

    return out
end

# no phase. common peak width λ.
function evalabsorptionspectrum(x::T, λ::T, a::Array{T,D}, Ω::Array{T,D}) where {T,D}
    L = length(a)
    @assert length(Ω) == L

    out = zero(T)
    for l = 1:L
    denominator = λ^2 + (x-Ω[l])^2
    out += a[l]*λ/denominator
    end

    return out
end

function evalabsorptionspectrumarea(lower_limit::T,
    upper_limit::T,
    λ::T,
    a::Array{T,D},
    Ω::Array{T,D}) where {T,D}

    L = length(a)
    @assert length(Ω) == L

    out = zero(T)
    for l = 1:L

    numerator = Ω[l] - lower_limit
    Ia = -atan(numerator,λ)

    numerator = Ω[l] - upper_limit
    Ib = -atan(numerator,λ)

    out += a[l]*(Ib-Ia)
    end

    return out
end


text:
"""
    cleanup(fig, exclude)

Internal function to delte all auxiliary files associated with the figure `fig` apart from those
contained in `exclude`. In particular: any data files supporting `Drawings` object and any `.gle`
script. This function is called after the `gle` engine is called (after an output image has
been generated) and only if the flag `GP_ENV["DEL_INTERM"]` is set to `true` (default).
"""
function cleanup(fig::Figure{GLE}, exclude=Vector{String}())
    # aux `.gle` folder
    rm(joinpath(GP_ENV["TMP_PATH"], fig.id * ".gle"), force=true)
    # aux `fig.id...` files
    auxfiles = Iterators.filter(f -> startswith(f, fig.id), readdir(GP_ENV["TMP_PATH"]))
    for af ∈ auxfiles
        (af ∈ exclude) || rm(joinpath(GP_ENV["TMP_PATH"], af), force=true)
    end
    return nothing
end

"""
    savefig()
    savefig(fname)
    savefig(fig; opts...)
    savefig(fig, fname; opts...)

Save the current figure (`gcf()`) or the figure `fig` to file. If the name `fname` is provided
with an extension, that name will be used and the extension will be used to determine the output
format. If no `fname` is provided, the `fig.id` will be used. Optional keyword arguments are:

* `format`: a string specifying the output format when the `fname` does not specify an extension,
possible formats are `eps|ps|pdf|svg|jpg|png`.
* `path`: a string specifying the output path where the output file should be saved. If it doesn't
exist, an attempt will be made to create it (which may fail). By default the file is saved in
the current folder.
* `res`: output file resolution (typical number is 200).
"""
function savefig(fig::Figure{GLE}, fn::String="";
                 format::String="png", path::String="", res::Int=200, _noout=false)

    GP_ENV["HAS_BACKEND"] || (@warn "No backend available to render the figure."; return)

    # by default take the figure id as name
    isempty(fn) && (fn = fig.id)
    # extract device from file name (if any)
    fn, ext = splitext(fn)
    isempty(ext) && (ext = ifelse(isempty(format), "png", format))
    ext = lowercase(ext)
    startswith(ext, ".") && (ext = ext[2:end])

    # check it's associated with a valid output type
    ext ∈ ["eps", "ps", "pdf", "svg", "jpg", "png"] ||
            throw(OptionValueError("output file type", ext))

    # read parameters from the `fig`: transparency & tex
    cairo     = ifelse(isdef(fig.transparency), `-cairo`, ``)
    texlabels = ifelse(isdef(fig.texlabels),    `-tex`,   ``)
    # path for the GLE script and the log
    fto  = joinpath(GP_ENV["TMP_PATH"], fig.id)
    fin  = fto * ".gle"
    flog = fto * ".log"
    # path for the output file
    fpo  = normpath(ifelse(isempty(path), joinpath(pwd(), fn), joinpath(path, fn)))
    # check the output dir
    dir = splitdir(fpo)[1]
    if !isdir(dir)
        @warn("The directory $dir does not exist. I will try to create it and save there.")
        mkpath(dir) # this may error
    end
    fout = fpo * ".$ext"
    # check what to do given the background
    transparent = ifelse(isnothing(fig.bgcolor) || Colors.alpha(fig.bgcolor) < 1,
                         `-transparent`, ``)
    # Assemble the figure and compile the lot note that assemble_figure writes to `fin`
    # (see also apply_gle/figure)
    assemble_figure(fig)
    glecom = pipeline(`gle -d $ext -r $res $cairo $texlabels $transparent -vb 0 -o $fout $fin`, stdout=flog, stderr=flog)
    # only for warmup: finish early (don't call gle)
    _noout && return
    # spwan the gle command, and check if it failed
    if _noout || !success(glecom)
        log = read(flog, String)
        GP_ENV["DEL_INTERM"] && cleanup(fig)
        error("GLE error: ... \n$log")
    end
    # cleanup if required and return output file name
    GP_ENV["DEL_INTERM"] && cleanup(fig, [fn * ".$ext"])
    return fout
end
savefig(fn::String=""; opts...) = savefig(gcf(), fn; opts...)

"""
    PreviewFigure

Internal type to wrap around a figure that is to be previewed. It may be previewed as an image.
"""
struct PreviewFigure
    fig::Figure
end

"""
    ImgPreview

Internal type to wrap around a figure that is to be previewed as an image. Note the distinction
with `PreviewFigure` which precedes it and ends up being shown as an image only if the continuous
preview toggle is on (`GP_ENV["CONT_PREVIEW"]`).
"""
struct ImgPreview
    fig::Figure
end

# NOTE tricks to suppress automatic opening of plot panel in atom
preview(f::Figure) = ImgPreview(f)
preview() = preview(gcf())

function Base.show(io::IO, obj::Union{PreviewFigure,DrawingHandle})
    GP_ENV["CONT_PREVIEW"] || return nothing
    fig = isa(obj, PreviewFigure) ? obj.fig : gcf()
    display(ImgPreview(fig))
    return nothing
end

function Base.show(io::IO, ::MIME"image/png", obj::ImgPreview)
    # check if can display
    (isdefined(Main, :Atom) && Main.Atom.PlotPaneEnabled.x) ||
        (isdefined(Main, :IJulia) && Main.IJulia.inited)    ||
            (@warn("Preview is only available in Juno and IJulia."); return nothing)
    # trigger a draft build
    isdir(GP_ENV["TMP_PATH"]) || mkpath(GP_ENV["TMP_PATH"])
    fname = savefig(obj.fig, "__PREVIEW__"; res=100, path=GP_ENV["TMP_PATH"])
    isnothing(fname) && return nothing
    # write to IO
    write(io, read(fname))
    GP_ENV["DEL_INTERM"] && rm(fname)
    return nothing
end


text:
using PlanktonIndividuals, Serialization

grid = RectilinearGrid(size = (1, 1, 16), x = (0,32), y = (0,32), z = (0,-32),
                              topology = (Bounded, Bounded, Bounded), halo = (2, 2, 2))

model = PlanktonModel(CPU(), grid, mode = MacroMolecularMode())

function tot_mass(nut, g)
    mass = zeros(g.Nx, g.Ny, g.Nz)
    for i in 1:g.Nx
        for j in 1:g.Ny
            for k in 1:g.Nz
                mass[i,j,k] = nut[i+g.Hx, j+g.Hy, k+g.Hz] * PlanktonIndividuals.Grids.volume(i+g.Hx, j+g.Hy, k+g.Hz, g)
            end
        end
    end
    return sum(mass)
end

TP = tot_mass(model.nutrients.PO4.data, grid) +
     tot_mass(model.nutrients.DOP.data, grid) +
     tot_mass(model.nutrients.POP.data, grid)
TP = TP + sum(model.individuals.phytos.sp1.data.PST .+
              model.individuals.phytos.sp1.data.DNA .* model.individuals.phytos.sp1.p.R_PC_DNA .+
              model.individuals.phytos.sp1.data.RNA .* model.individuals.phytos.sp1.p.R_PC_RNA)

uvel = zeros(2,1,16,11)
vvel = zeros(1,2,16,11)
wvel = zeros(1,1,17,11)

for i in 1:11
    uvel[:,:,:,i] .= randn(2,1,16) .* 1e-4
    vvel[:,:,:,i] .= randn(1,2,16) .* 1e-4
    wvel[:,:,:,i] .= randn(1,1,17) .* 1e-4
end

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10, vels=(u=uvel, v=vvel, w=wvel))

update!(sim)

TPt = tot_mass(model.nutrients.PO4.data, grid) +
      tot_mass(model.nutrients.DOP.data, grid) +
      tot_mass(model.nutrients.POP.data, grid)
TPt = TPt + sum(model.individuals.phytos.sp1.data.PST .+
               model.individuals.phytos.sp1.data.DNA .* model.individuals.phytos.sp1.p.R_PC_DNA .+
               model.individuals.phytos.sp1.data.RNA .* model.individuals.phytos.sp1.p.R_PC_RNA)

@testset "PlanktonIndividuals 1D tests:" begin
    @test isapprox(TP,TPt; atol=1e1)
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec ${JULIA:-julia} --color=yes --startup-file=no \
    -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" \
    "$@"
=#

csv_path, output = ARGS

using SnoopCompile

data = SnoopCompile.read(csv_path)
pc = SnoopCompile.parcel(reverse!(data[2]))
SnoopCompile.write(output, pc)


text:
################################################################################
# A few PerpleX calculation examples using the Julia-Perplex interface

## --- Import some useful packages
    using StatGeochem
    using Plots; gr();
    using Statistics, DelimitedFiles, SpecialFunctions

## --- Configure

    # Absolute paths to perplex resources
    perplexdir = joinpath(resourcepath,"perplex-stable")
    scratchdir = "./scratch/" # Location of directory to store output files

    # Attempt to install perplex, if not already extant
    if !isfile(joinpath(perplexdir,"vertex"))
        # Make sure resourcepath exists
        run(`mkdir -p $resourcepath`)

        # Try to compile PerpleX from source; if that fails, try to download linux binaries
        try
            # Check if there is a fortran compiler
            run(`gfortran -v`)

            # Download Perplex v6.8.7 -- known to work with interface used here
            file = download("https://storage.googleapis.com/statgeochem/perplex-6.8.7-source.zip", joinpath(resourcepath,"perplex-stable.zip"))

            # # For a more updated perplex version, you might also try
            # file = download("https://petrol.natur.cuni.cz/~ondro/perplex-sources-stable.zip", joinpath(resourcepath,"perplex-stable.zip"))

            run(`unzip -u $file -d $resourcepath`) # Extract
            system("cd $perplexdir; make") # Compile
        catch
            @warn "Failed to compile from source, trying precompiled linux binaries instead"
            run(`mkdir -p $perplexdir`)
            file = download("https://petrol.natur.cuni.cz/~ondro/Perple_X_6.8.7_Linux_64_gfortran.tar.gz","perplex-6.8.7-linux.tar.gz")
            run(`tar -xzf $file -C $perplexdir`)
        end
    end

## --- # # # # # # # # # # # # # Initial composition # # # # # # # # # # # # # #

    ## McDonough Pyrolite
    #elements =    [ "SIO2", "TIO2", "AL2O3",  "FEO",  "MNO",  "MGO",  "CAO", "NA2O",  "K2O",  "H2O",  "CO2",]
    #composition = [45.1242, 0.2005, 4.4623, 8.0723, 0.1354, 37.9043, 3.5598, 0.3610, 0.0291, 0.1511, 0.0440,]

    # # Kelemen (2014) primitive continental basalt. H2O and CO2 are guesses
    # elements =    [ "SIO2", "TIO2", "AL2O3",  "FEO",  "MNO",  "MGO",  "CAO", "NA2O",  "K2O",  "H2O",  "CO2",]
    # composition = [50.0956, 0.9564, 15.3224, 8.5103, 0.1659, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]

    # # Kelemen (2014) primitive continental basalt excluding Mn. H2O and CO2 are guesses
    # elements =    [ "SIO2", "TIO2", "AL2O3",  "FEO",  "MGO",  "CAO", "NA2O",  "K2O",  "H2O",  "CO2",]
    # composition = [50.0956, 0.9564, 15.3224, 8.5103, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]

    # Kelemen (2014) primitive continental basalt excluding Mn and Ti since most melt models can"t handle them..
    elements =    [ "SIO2", "AL2O3",  "FEO",  "MGO",  "CAO", "NA2O",  "K2O",  "H2O",  "CO2",]
    composition = [50.0956, 15.3224, 8.5103, 9.2520, 9.6912, 2.5472, 0.8588, 2.0000, 0.6000,]

    # # Average Archean basalt (EarthChem data)
    # elements =    [ "SIO2", "TIO2", "AL2O3",   "FEO",  "MNO",   "MGO",  "CAO", "NA2O",  "K2O",  "H2O",  "CO2",]
    # composition = [49.2054, 0.8401, 12.0551, 11.4018, 0.2198, 12.3997, 9.3113, 1.6549, 0.4630, 1.8935, 0.5555,]

    # # A random granite
    # elements =    [ "SIO2", "TIO2", "AL2O3",  "FEO", "MGO", "CAO", "NA2O",  "K2O",  "H2O",  "CO2",]
    # composition = [  69.16,   0.42,   14.96, 2.9064,  1.52,  1.81,    3.7,   4.95, 1.8386, 0.4407,]

## --- # # # # # # # # # # # Some solution model options # # # # # # # # # # # #
    # Emphasis on phases from Green (2016) -- developed for metabasites, includes what is probably the best (and most expensive) amphibole model. Use with hp11ver.dat
    G_solution_phases = "Augite(G)\nOpx(JH)\ncAmph(G)\noAmph(DP)\nO(JH)\nSp(JH)\nGrt(JH)\nfeldspar_B\nMica(W)\nBio(TCC)\nChl(W)\nCtd(W)\nCrd(W)\nSa(WP)\nSt(W)\nIlm(WPH)\nAtg(PN)\nT\nB\nF\nDo(HP)\nScap\nChum\nNeph(FB)\n"
    G_excludes ="ged\nfanth\ngl\nilm\nilm_nol\n"

    # Emphasis on phases from White (2014) -- developed for metapelites. Use with hp11ver.dat (though can apparenty run with hp02ver.dat without crashing)
    W_solution_phases = "Omph(HP)\nOpx(W)\ncAmph(DP)\noAmph(DP)\nO(JH)\nSp(JH)\nGt(W)\nfeldspar_B\nMica(W)\nBi(W)\nChl(W)\nCtd(W)\nCrd(W)\nSa(WP)\nSt(W) \nIlm(WPH)\nAtg(PN)\nT\nB\nF\nDo(HP)\nScap\nChum\nPu(M)\n"
    W_excludes = "andr\nts\nparg\ngl\nged\nfanth\n"

    # Emphasis on phases from Jennings and Holland (2015) -- developed for mantle melting. Use with hp11ver.dat
    JH_solution_phases = "Cpx(JH)\nOpx(JH)\ncAmph(DP)\noAmph(DP)\nO(JH)\nSp(JH)\nGrt(JH)\nfeldspar_B\nMica(W)\nBio(TCC)\nChl(W)\nCtd(W)\nCrd(W)\nSa(WP)\nSt(W)\nIlm(WPH)\nAtg(PN)\nT\nB\nF\nDo(HP)\nScap\nChum\nNeph(FB)\n"
    JH_excludes = "ts\nparg\ngl\nged\nfanth\n"

    # Emphasis on phases from Holland and Powell -- all phases can be used with hp02ver.dat.
    HP_solution_phases = "Omph(HP)\nOpx(HP)\nGlTrTsPg\nAnth\nO(HP)\nSp(HP)\nGt(HP)\nfeldspar_B\nMica(CF)\nBio(TCC)\nChl(HP)\nCtd(HP)\nSapp(HP)\nSt(HP)\nIlHm(A)\nDo(HP)\nT\nB\nF\n"
    HP_excludes = ""

## --- # # # # # # # # # # # # # Isobaric example # # # # # # # # # # # # # # # #

    # Input parameters
    P = 1000 # Pressure, bar
    T_range = [0+273.15, 1500+273.15] # Temperature range, Kelvin

    # # Configure (run build and vertex)
    # melt_model = "melt(G)"
    # @time perplex_configure_isobar(perplexdir, scratchdir, composition, elements, P, T_range,
    #     dataset="hpha11ver.dat",
    #     npoints=100,
    #     excludes=G_excludes,
    #     solution_phases=melt_model*"\n"*G_solution_phases
    # )

    melt_model = "melt(HP)"
    @time perplex_configure_isobar(perplexdir, scratchdir, composition, elements, P, T_range,
        dataset="hp02ver.dat",
        npoints=100,
        excludes=HP_excludes,
        solution_phases=melt_model*"\n"*HP_solution_phases
    )

## --- Query all properties at a single temperature -- results returned as text

    T = 850+273.15
    data_isobaric = perplex_query_point(perplexdir, scratchdir, T) |> print

## --- Query the full isobar -- results returned as dict

    bulk = perplex_query_system(perplexdir, scratchdir)             # Get system data for all temperatures. Set include_fluid = "n" to get solid+melt only
    modes = perplex_query_modes(perplexdir, scratchdir)             # || phase modes
    melt = perplex_query_phase(perplexdir, scratchdir, melt_model)  # || melt data

    # Melt wt.% seems to be slightly inaccurate; use values from modes instead
    melt["wt_pct"] = modes[melt_model]

    # Create dictionary to hold solid composition and fill it using what we know from system and melt
    solid = Dict()
    solid["wt_pct"] = 100 .- melt["wt_pct"]
    for e in ["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"]
        solid[e] = (bulk[e] - (melt[e] .* melt["wt_pct"]/100)) ./ (solid["wt_pct"]/100)
    end
    renormalize!(solid,["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"],total=100)

## --- Plot melt composition as a function of melt percent

    h = plot(xlabel="Percent melt", ylabel="Wt. % in melt", title="$melt_model + G_solution_phases, $P bar")
    i = 0
    for e in ["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"]
        plot!(h, melt["wt_pct"], melt[e], label=e, color=lines[global i += 1])
        plot!(h, melt["wt_pct"], bulk[e], label="", color=lines[i], linestyle=:dot)
    end
    plot!(h,fg_color_legend=:white, framestyle=:box)
    savefig(h,"MeltComposition.pdf")
    display(h)

## --- Plot melt composition as a function of melt SiO2

    h = plot(xlabel="Magma SIO2 (wt.%)", ylabel="Wt. % in melt", title="$melt_model + G_solution_phases, $P bar")
    i = 1
    for e in ["AL2O3","FEO","MGO","CAO","NA2O","K2O"]
        plot!(h,melt["SIO2"], melt[e], label=e, color = lines[global i +=1])
    end
    plot!(h,fg_color_legend=:white, framestyle=:box)
    savefig(h,"MeltCompositionvsSiO2.pdf")
    display(h)

## --- Plot solid composition as a function of melt percent

    h = plot(xlabel="Percent melt", ylabel="Wt. % in solid", title="$melt_model + G_solution_phases, $P bar")
    i = 0
    for e in ["SIO2","AL2O3","FEO","MGO","CAO","NA2O","K2O"]
        plot!(h, melt["wt_pct"], solid[e], label=e, color=lines[global i +=1])
    end
    plot!(h,fg_color_legend=:white, framestyle=:box, legend=:topleft)
    savefig(h,"SolidComposition.pdf")
    display(h)

## --- Plot modes of all phases as a function of temperature

    h = plot(xlabel="T (C)", ylabel="Weight percent", title="$melt_model + G_solution_phases, $P bar")
    for m in modes["elements"][3:end]
        plot!(h, modes["T(K)"] .- 273.15, modes[m], label=m)
    end
    plot!(h,fg_color_legend=:white, framestyle=:box)
    savefig(h,"PhaseModes.pdf")
    display(h)

## --- Plot modes of all phases as a function of melt percent

    h = plot(xlabel="Percent melt", ylabel="Weight percent", title="$melt_model + G_solution_phases, $P bar")
    for m in modes["elements"][3:end]
        plot!(h, modes[melt_model], modes[m], label=m)
    end
    plot!(h,fg_color_legend=:white, framestyle=:box)
    savefig(h,"PhaseModesvsF.pdf")
    display(h)

## --- Plot seismic properties

    # Query seismic properties along the whole profile
    seismic = perplex_query_seismic(perplexdir, scratchdir)
    seismic["vp/vs"][seismic["vp/vs"] .> 100] .= NaN # Exclude cases where vs drops to zero

    h = plot(xlabel="Pressure", ylabel="Property")
    plot!(h,seismic["T(K)"],seismic["vp,km/s"], label="vp,km/s")
    plot!(h,seismic["T(K)"],seismic["vp/vs"], label="vp/vs")
    plot!(h,seismic["T(K)"],seismic["rho,kg/m3"]/1000, label="rho, g/cc")
    savefig(h,"GeothermSeismicProperties.pdf")
    display(h)

## --- # # # # # # # # # # # Geothermal gradient example # # # # # # # # # # # #

    # Input parameters
    P_range = [280, 28000] # Pressure range to explore, bar (roughly 1-100 km depth)
    T_surf = 273.15 # Temperature of surface (K)
    geotherm = 0.01 # Geothermal gradient of 0.1 K/bar == about 28.4 K/km
    melt_model = ""

    # Configure (run build and vertex)
    @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
        P_range, T_surf, geotherm, dataset="hp02ver.dat", excludes=HP_excludes,
        solution_phases=HP_solution_phases, npoints=200, index=2)

    # # Alternative configuration, using hpha02ver.dat
    # @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
    #     P_range, T_surf, geotherm, dataset="hpha02ver.dat", excludes="qGL\n"*HP_excludes,
    #     solution_phases=HP_solution_phases, npoints=200, index=2)

    # # Alternative configuration, using hpha02ver.dat and new phases for metapelites
    # @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
    #     P_range, T_surf, geotherm, dataset="hpha02ver.dat", excludes="qGL\n"*W_excludes,
    #     solution_phases=W_solution_phases, npoints=200, index=2)

## --- Plot modes of all phases as a function of temperature

    modes = perplex_query_modes(perplexdir, scratchdir, index=2)             # || phase modes

    h = plot(xlabel="T (C)", ylabel="Weight percent")
    for m in modes["elements"][3:end]
        plot!(h, modes["T(K)"] .- 273.15, modes[m], label=m)
    end
    plot!(h,fg_color_legend=:white, framestyle=:box)
    savefig(h,"GeothermPhaseModes.pdf")
    display(h)

## --- Plot seismic properties

    # Query seismic properties along the whole profile
    seismic = perplex_query_seismic(perplexdir, scratchdir, index=2)
    seismic["vp/vs"][seismic["vp/vs"] .> 100] .= NaN # Exclude cases where vs drops to zero

    h = plot(xlabel="Pressure", ylabel="Property")
    plot!(h,seismic["P(bar)"],seismic["vp,km/s"], label="vp,km/s")
    plot!(h,seismic["P(bar)"],seismic["vp/vs"], label="vp/vs")
    plot!(h,seismic["P(bar)"],seismic["rho,kg/m3"]/1000, label="rho, g/cc")
    plot!(h,seismic["P(bar)"],seismic["T(K)"]/1000, label="T(K)/1000")
    savefig(h,"GeothermSeismicProperties.pdf")
    display(h)

## --- Query all properties at a single pressure

    P = 10000
    data_geotherm = perplex_query_point(perplexdir, scratchdir, P, index=2) |> print

## --- Compare seismic properties for several different geotherms, as a function of Pressure

    # Input parameters
    P_range = [280, 14000] # Pressure range to explore, bar (roughly 1-50 km depth)
    T_surf = 273.15 # Temperature of surface (K)
    dataset="hp02ver.dat"
    yelem = "vp,km/s"
    # yelem = "rho,kg/m3"

    h = plot(xlabel="Pressure (bar)", ylabel=yelem)
    for i=1:5
        geotherm = i/50
        # Configure (run build and vertex)
        @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
            P_range, T_surf, geotherm, dataset=dataset, excludes=HP_excludes,
            solution_phases=HP_solution_phases, npoints=200, index=10+i)
        # Query perplex results
        seismic = perplex_query_seismic(perplexdir, scratchdir, index=10+i)
        # Plot results
        plot!(h,seismic["P(bar)"],seismic[yelem], label="$geotherm K/bar")
        display(h)
    end

## --- Compare seismic properties for several different geotherms, as a function of temperature

    # Input parameters
    P_range = [280, 16000] # Pressure range to explore, bar (roughly 1-60 km depth)
    T_surf = 273.15 # Temperature of surface (K)
    dataset="hp02ver.dat"
    yelem = "vp,km/s"
    # yelem = "rho,kg/m3"

    h = plot(xlabel="Temperature (K)", ylabel=yelem)
    for i=1:5
        geotherm = i/50
        # Configure (run build and vertex)
        @time perplex_configure_geotherm(perplexdir, scratchdir, composition, elements,
            P_range, T_surf, geotherm, dataset=dataset, excludes=HP_excludes,
            solution_phases=HP_solution_phases, npoints=200, index=10+i)
        # Query perplex results
        seismic = perplex_query_seismic(perplexdir, scratchdir, index=10+i)
        # Plot results
        plot!(h,seismic["T(K)"],seismic[yelem], label="$geotherm K/bar")
        display(h)
    end

## --- Build a pseudosection

    # Input parameters
    mingeotherm = 0.009 # Geothermal gradient, K/bar. For reference, 0.1 K/bar ≈ 28.4 K/km
    max_lith_T = 1300 # Maximum temperature of the TBL lithosphere (no point calculating beyond that)
    P_range = [280, (max_lith_T-273.15)/mingeotherm] # Pressure range to explore, bar (roughly 1-100 km depth)
    T_range = [273.15, max_lith_T] # Temperature range to explore, K
    melt_model = ""

    # Configure (run build and vertex)
    @time perplex_configure_pseudosection(perplexdir, scratchdir, composition,
        elements, P_range, T_range, dataset="hpha02ver.dat", excludes=HP_excludes,
        solution_phases=melt_model*HP_solution_phases, index=3, xnodes=200, ynodes=200)

## --- Query modes from a pseudosection

    geotherm = 0.015 # Geothermal gradient, K/bar. For reference, 0.1 K/bar ≈ 28.4 K/km
    @time modes = perplex_query_modes(perplexdir, scratchdir, [280, (max_lith_T-273.15)/geotherm], T_range, index=3, npoints=200)

    h = plot(xlabel="T (C)", ylabel="Weight percent")
    for m in modes["elements"][3:end]
        plot!(h, modes["T(K)"] .- 273.15, modes[m], label=m)
    end
    plot!(h,fg_color_legend=:white, framestyle=:box)
    savefig(h,"GeothermPhaseModesPseudosection.pdf")
    display(h)

## --- Query seismic properties from a pseudosection

    geotherm = 0.15 # Geothermal gradient, K/bar. For reference, 0.1 K/bar ≈ 28.4 K/km
    @time seismic = perplex_query_seismic(perplexdir, scratchdir, [280, (max_lith_T-273.15)/geotherm], T_range, index=3, npoints=200)

    h = plot(xlabel="Pressure", ylabel="Property")
    plot!(h,seismic["P(bar)"],seismic["vp,km/s"], label="vp,km/s")
    plot!(h,seismic["P(bar)"],seismic["vp/vs"], label="vp/vs")
    plot!(h,seismic["P(bar)"],seismic["rho,kg/m3"]/1000, label="rho, g/cc")
    plot!(h,seismic["P(bar)"],seismic["T(K)"]/1000, label="T(K)/1000")
    # savefig(h,"GeothermSeismicProperties.pdf")
    # display(h)

## ---

    geotherm = 0.15 # Geothermal gradient, K/bar. For reference, 0.1 K/bar ≈ 28.4 K/km
    @time bulk = perplex_query_system(perplexdir, scratchdir, [280, (max_lith_T-273.15)/geotherm], T_range, index=3, npoints=200, include_fluid="n")

    h = plot(xlabel="Pressure", ylabel="Property")
    plot!(h,bulk["P(bar)"],bulk["vp,km/s"], label="vp,km/s")
    plot!(h,bulk["P(bar)"],bulk["vp/vs"], label="vp/vs")
    plot!(h,bulk["P(bar)"],bulk["rho,kg/m3"]/1000, label="rho, g/cc")
    plot!(h,bulk["P(bar)"],bulk["T(K)"]/1000, label="T(K)/1000")
    # savefig(h,"GeothermSeismicProperties.pdf")
    # display(h)

## --- End of File


text:
module CommonStructures

using Statistics

export ModelFit,
       BoundariesFit,
       PlaneWave,
       LayerTMMO,
       _get_beam_parameters,
       FitProcedure

## Defines the abstract type for all the fitting.
abstract type FitProcedure end

## General structure to build the indices of refraction models for otimization.
struct ModelFit{T1,T2,T3,T4,T5} <: FitProcedure where{T1<:Symbol, T2<:Real, T3<:Real, T4<:Any, T5<:Any}
    model::T1
    df::T2
    dfm::T1
    c::T3
    N::NamedTuple{(:ninc, :nhost),Tuple{Array{T4,1},Array{T5,1}}} # effective medium equations, EMA
    function ModelFit(
        model::T1;
        df::T2=1.0/3.0,
        dfm::T1=:spheres,
        c::T3=0.0,
        N::NamedTuple{(:ninc, :nhost),Tuple{Array{T4,1},Array{T5,1}}}=(ninc=[],nhost=[]),
    ) where {T1<:Symbol, T2<:Real, T3<:Real, T4<:Any, T5<:Any}
        length(N) == 2 || throw("length(N) == 2, one vector for each index of refraction:
        N=(n1, n2), where n1 and n2 are vectors.")
        return new{T1,T2,T3,T4,T5}(model, df, dfm, c, N)
    end
end

## Definition of the boundaries for the solution space grids.
struct BoundariesFit{T1,T2,T3,T4,T5} <: FitProcedure where {T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Int64}
    odlo::T1
    odup::T2
    plo::T3
    pup::T4
    Nod::T5
    Np::T5
    function BoundariesFit(
        odlo::T1, odup::T2, plo::T3, pup::T4;
        Nod::T5=20, Np::T5=20,
    ) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Int64}
        return new{T1,T2,T3,T4,T5}(odlo, odup, plo, pup, Nod, Np)
    end
end

## Definition of the beam type.
abstract type LightSource end
struct PlaneWave{T1,T2,T3} <: LightSource where {T1<:Real,T2<:Real,T3<:Real}
    λ::AbstractVector{T1}
    θ::AbstractVector{T2}
    p::T3
    function PlaneWave(
        λ::AbstractVector{T1}, θ::AbstractVector{T2};
        p::T3=0.5,
    ) where {T1<:Real,T2<:Real,T3<:Real}
        0.0<=p<=1.0 || throw("The polarisation, p argument, must lie between 0 and 1: 0.0<=p<=1.0")
        return new{T1,T2,T3}(λ, θ, p)
    end
end

## Material type definition with geometrical (physical) and optical thickness.
abstract type Material end
struct LayerTMMO{T1,T2,T3,T4} <: Material where {T1<:ComplexF64, T2<:Symbol, T3<:Real, T4<:Number}
    n::Array{T1,1}
    type::T2
    d::T3
    nλ0::Array{T4}
    function LayerTMMO(
        n::AbstractArray{T1,1};
        type::T2=:GT,
        d::T3=1.0/4.0,
        nλ0::Array{T4}=[-eps()],
    ) where {T1<:ComplexF64, T2<:Symbol, T3<:Real, T4<:Number}
        return new{T1,T2,T3,T4}(n, type, d, nλ0)
    end
end

# function meanSquaredError(
#     X::Array{T1}, Xexp::Array{T1};
#     σ::Array{T1}=ones(size(Xexp)),
# ) where {T1<:Float64}
#     # mse = sum(mean(abs2.(X .- Xexp)./σ, dims=1))
#     _absq = abs2.(X .- Xexp)./σ
#     _mse = mean(_absq, dims=1)
#     _sse = sum(_absq, dims=1)
#     mse = sqrt(sum(_mse + _sse)/size(Xexp, 1))
#     # mse = mean(abs2.((X .- Xexp)./σ))
#     # mse = sum(abs2.((X .- Xexp)./σ))
#     return mse
# end

## Returns the parameters with the correct format for the beam suitable for the calculations.
function _get_beam_parameters(beam::PlaneWave)
    beam_ = PlaneWave(float.(beam.λ), deg2rad.(beam.θ))
    x = (λLen=length(beam.λ), θLen=length(beam.θ))
    return beam_, x
end

end # module


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct BaseMachine <: Machine
        R::Float64
        Xd_p::Float64
        eq_p::Float64
        MVABase::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of an Automatic Voltage Regulator Type II -  Typical static exciter model

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: (0, nothing)
- `Xd_p::Float64`: Reactance after EMF in machine per unit, validation range: (0, nothing)
- `eq_p::Float64`: Fixed EMF behind the impedance, validation range: (0, nothing)
- `MVABase::Float64`: Nominal Capacity in MVA, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct BaseMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in machine per unit"
    Xd_p::Float64
    "Fixed EMF behind the impedance"
    eq_p::Float64
    "Nominal Capacity in MVA"
    MVABase::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function BaseMachine(R, Xd_p, eq_p, MVABase, ext=Dict{String, Any}(), )
    BaseMachine(R, Xd_p, eq_p, MVABase, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function BaseMachine(; R, Xd_p, eq_p, MVABase, ext=Dict{String, Any}(), )
    BaseMachine(R, Xd_p, eq_p, MVABase, ext, )
end

# Constructor for demo purposes; non-functional.
function BaseMachine(::Nothing)
    BaseMachine(;
        R=0,
        Xd_p=0,
        eq_p=0,
        MVABase=0,
        ext=Dict{String, Any}(),
    )
end

"""Get BaseMachine R."""
get_R(value::BaseMachine) = value.R
"""Get BaseMachine Xd_p."""
get_Xd_p(value::BaseMachine) = value.Xd_p
"""Get BaseMachine eq_p."""
get_eq_p(value::BaseMachine) = value.eq_p
"""Get BaseMachine MVABase."""
get_MVABase(value::BaseMachine) = value.MVABase
"""Get BaseMachine ext."""
get_ext(value::BaseMachine) = value.ext
"""Get BaseMachine states."""
get_states(value::BaseMachine) = value.states
"""Get BaseMachine n_states."""
get_n_states(value::BaseMachine) = value.n_states
"""Get BaseMachine internal."""
get_internal(value::BaseMachine) = value.internal


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
include(joinpath("..", "common.jl"))

name = "FFMPEG"

# The products that we will ensure are always built
products = [
    ExecutableProduct("ffmpeg", :ffmpeg),
    ExecutableProduct("ffprobe", :ffprobe),
    LibraryProduct(["libavcodec", "avcodec"], :libavcodec),
    LibraryProduct(["libavdevice", "avdevice"], :libavdevice),
    LibraryProduct(["libavfilter", "avfilter"], :libavfilter),
    LibraryProduct(["libavformat", "avformat"], :libavformat),
    LibraryProduct(["libavresample", "avresample"], :libavresample),
    LibraryProduct(["libavutil", "avutil"], :libavutil),
    LibraryProduct(["libpostproc", "postproc"], :libpostproc),
    LibraryProduct(["libswresample", "swresample"], :libswresample),
    LibraryProduct(["libswscale", "swscale"], :libswscale),
]

# Dependencies that must be installed before this package can be built.
# TODO: Theora once it's available
dependencies = [
    BuildDependency("nv_codec_headers_jll"),
    Dependency("libass_jll"; compat="0.15.1"),
    Dependency("libfdk_aac_jll"),
    Dependency("FriBidi_jll"),
    Dependency("FreeType2_jll"),
    Dependency("LAME_jll"),
    Dependency("libvorbis_jll"),
    Dependency("Ogg_jll"),
    BuildDependency("LibVPX_jll"), # We use the static archive
    Dependency("x264_jll"; compat="~2021.05.05"),
    Dependency("x265_jll"; compat="~3.5"),
    Dependency("Bzip2_jll"; compat="1.0.8"),
    Dependency("Zlib_jll"),
    Dependency("OpenSSL_jll"),
    Dependency("Opus_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script(; ffplay=false), platforms, products, dependencies; preferred_gcc_version=preferred_gcc_version, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using BandedMatrices, BlockBandedMatrices, ArrayLayouts, BlockArrays, FillArrays, CuArrays, GPUArrays, LinearAlgebra
import BlockBandedMatrices: _BandedBlockBandedMatrix
import BandedMatrices: _BandedMatrix

function finitedifferences!(data::AbstractMatrix{T}) where T
    N = nblocks(data,2)
    sc = N^2 * one(T)
    for J = 1:N
        data.blocks[1,J][1,:] .= zero(T)
        data.blocks[2,J][1,:] .= sc
        data.blocks[3,J][1,:] .= zero(T)
    end
    for J = 1:N
        data.blocks[1,J][2,:] .= sc
        data.blocks[2,J][2,:] .= -4sc
        data.blocks[3,J][2,:] .= sc
    end
    for J = 1:N
        data.blocks[1,J][3,:] .= zero(T)
        data.blocks[2,J][3,:] .= sc
        data.blocks[3,J][3,:] .= zero(T)
    end
end

function finitedifferences(::Type{Mat}, N) where Mat<:AbstractMatrix
    data = BlockMatrix{eltype(Mat),Mat}(undef, Fill(3,3), Fill(N,N))
    finitedifferences!(data)
    _BandedBlockBandedMatrix(data, (Fill(N,N), Fill(N,N)), (1,1), (1,1))
end

finitedifferences(::Type{Arr}, N) where Arr<:AbstractArray = FiniteDifference(Arr{2}, N)

T = Float32
N = 10
x = y = range(zero(T), stop=one(T), length=N)

u₀ = BlockVector{T, Vec}(undef, Fill(N,N))

b = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N))
c = BlockVector{T, CuArray{T,1}}(Zeros{T}(N^2), Fill(N,N))

function mul!(c, A::BandedBlockBandedMatrix{T,<:BlockMatrix{<:Any,<:CuArray}}, b) where T
    fill!(c, 0)

    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    for J = 1:N, K = max(1,J-l):min(N,J+u)
        B = _BandedMatrix(A.data.blocks[K-J+u+1,J],N,λ,μ)
        muladd!(one(T), B, b.blocks[J], one(T), c.blocks[K])
    end
    c
end

function mul!(c, A::BandedBlockBandedMatrix{<:Any,<:BlockMatrix{<:Any,<:CuArray}}, b)
    fill!(c, 0)

    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    for J = 1:N, K = max(1,J-l):min(N,J+u)
        B = _BandedMatrix(A.data.blocks[K-J+u+1,J],N,λ,μ)
        ArrayLayouts.materialize!(MulAdd(1f0, B, b.blocks[J], 1f0, c.blocks[K]))
    end
    c
end

function mul2!(c, A::BandedBlockBandedMatrix{<:Any,<:BlockMatrix{<:Any,<:CuArray}}, x)
    fill!(c, 0)

    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    @inbounds for b = -u:0
        for J = 1-b:M
            K = J + b
            B = _BandedMatrix(A.data.blocks[b+u+1,J],N,λ,μ)
            ArrayLayouts.materialize!(MulAdd(1f0, B, x.blocks[J], 1f0, c.blocks[K]))
        end
    end

    @inbounds for b = 1:l
        for J = 1:M-b
            K = J + b
            B = _BandedMatrix(A.data.blocks[b+u+1,J],N,λ,μ)
            ArrayLayouts.materialize!(MulAdd(1f0, B, x.blocks[J], 1f0, c.blocks[K]))
        end
    end

    c
end

function mul3!(c, A::BandedBlockBandedMatrix{T,<:BlockMatrix{<:Any}}, x) where T
    # fill!(c, 0)

    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    @inbounds for b = -u:0
        for J = 1-b:M
            K = J + b
            BLAS.gbmv!('N', N, λ, μ,  one(T), A.data.blocks[b+u+1,J], x.blocks[J], one(T), c.blocks[K])
        end
    end

    @inbounds for b = 1:l
        for J = 1:M-b
            K = J + b
            BLAS.gbmv!('N', N, λ, μ,  one(T), A.data.blocks[b+u+1,J], x.blocks[J], one(T), c.blocks[K])
        end
    end

    c
end

function mul_diag!(c, A::BandedBlockBandedMatrix{<:Any,<:BlockMatrix{<:Any,<:CuArray}}, x)
    # fill!(c, 0)

    l, u = blockbandwidths(A)
    λ, μ = subblockbandwidths(A)
    N,M = nblocks(A)

    for J = 1:M
        BLAS.gbmv!('N', N, λ, μ,  1f0, A.data.blocks[2,J], x.blocks[J], 1f0, c.blocks[J])
    end
    c
end

N = 4000;
T = Float32
A = finitedifferences(CuArray{T,2},N)
x = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N))
c = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N))

A_cpu = finitedifferences(Array{T,2},N)
x_cpu = BlockVector{T, Array{T,1}}(undef, Fill(N,N))
c_cpu = BlockVector{T, Array{T,1}}(undef, Fill(N,N))

@time mul3!(c, A, x);
@time mul3!(c_cpu, A_cpu, x_cpu);

GPUArrays.synchronize(A)
x = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N))
c = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N))
b .= randn(T,N^2)

function finitedifference_2d(n)
    h = 1/n
    D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
    D_xx = BandedBlockBandedMatrix(Kron(D², Eye(n)))
    D_yy = BandedBlockBandedMatrix(Kron(Eye(n), D²))
    D_xx + D_yy
end

A_cpu = finitedifference_2d(N)
b_cpu = randn(N^2);
c_cpu = similar(b_cpu)
@time c_cpu .= Mul(A_cpu, b_cpu)

## large bandwidths
L = U = 20
    N = 500
    A = BandedBlockBandedMatrix{Float64}(undef, (Fill(N,N), Fill(N,N)), (L, U), (L, U))
    A.data.blocks .= randn.()

    b = randn(N^2);
    c = similar(b)
@time c .= Mul(A,b)

function cu_rand(::Type{T}, N, (L, U)) where T
    data = BlockMatrix{T,CuArray{T,2}}(undef, Fill(L+U+1,L+U+1), Fill(N,N))
    _BandedBlockBandedMatrix(data, (Fill(N,N), Fill(N,N)), (L,U), (L,U))
end

T = Float32;
N = 500; L = U = 20; A = cu_rand(T, N, (L,U)); b = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N));
    c = BlockVector{T, CuArray{T,1}}(undef, Fill(N,N));

A_s = sparse(A)

@time A_s*b;

n = 1_000_000; A = CuArray{Float32}(3,n); B = similar(A); C = similar(A);
    x = CuArray{Float32}(n); y = similar(x); z = similar(x);
    a = similar(x);  b = similar(x); c = similar(x);
    @time multgbmv!((A,B,C), (x,y,z), (a,b,c), n);

function multgbmv!((A,B,C), (x,y,z), (a,b,c), n)
   BLAS.gbmv!('N', n, 1, 1, 1f0, A, x, 0f0, a)
   BLAS.gbmv!('N', n, 1, 1, 1f0, B, y, 0f0, b)
   BLAS.gbmv!('N', n, 1, 1, 1f0, C, z, 0f0, x)
end


text:
# ---
# title: 467. Unique Substrings in Wraparound String
# id: problem467
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/unique-substrings-in-wraparound-string/description/>
# hidden: true
# ---
#
# Consider the string `s` to be the infinite wraparound string of
# "abcdefghijklmnopqrstuvwxyz", so `s` will look like this:
# "...zabcdefghijklmnopqrstuvwxyzabcdefghijklmnopqrstuvwxyzabcd....".
#
# Now we have another string `p`. Your job is to find out how many unique non-
# empty substrings of `p` are present in `s`. In particular, your input is the
# string `p` and you need to output the number of different non-empty substrings
# of `p` in the string `s`.
#
# **Note:** `p` consists of only lowercase English letters and the size of p
# might be over 10000.
#
# **Example 1:**
#
#
#
#     Input: "a"
#     Output: 1
#
#     Explanation: Only the substring "a" of string "a" is in the string s.
#
#
# **Example 2:**
#
#
#
#     Input: "cac"
#     Output: 2
#     Explanation: There are two substrings "a", "c" of string "cac" in the string s.
#
#
# **Example 3:**
#
#
#
#     Input: "zab"
#     Output: 6
#     Explanation: There are six substrings "z", "a", "b", "za", "ab", "zab" of string "zab" in the string s.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end