text:
__precompile__(true)

"The module for [PCG Family](@ref)."
module PCG

# PCG
export PCGStateMCG, PCGStateOneseq, PCGStateUnique, PCGStateSetseq

# PCG Methods
export PCG_XSH_RS, PCG_XSH_RR, PCG_RXS_M_XS, PCG_XSL_RR, PCG_XSL_RR_RR

export PCGUInt, PCGMethod, PCG_LIST

export bounded_rand, advance!

const pcg_uints = (UInt8, UInt16, UInt32, UInt64, UInt128)
const PCGUInt = Union{pcg_uints...}

"""
One of PCG output method: high xorshift, followed by a random shift.

Fast.
"""
const PCG_XSH_RS = Val{:XSH_RS}
"""
One of PCG output method: high xorshift, followed by a random rotate.

Fast. Slightly better statistically than `PCG_XSH_RS`.
"""
const PCG_XSH_RR = Val{:XSH_RR}
"""
One of PCG output method: random xorshift, mcg multiply, fixed xorshift.

The most statistically powerful generator, but slower than some of the others.
"""
const PCG_RXS_M_XS = Val{:RXS_M_XS}
"""
One of PCG output method: fixed xorshift (to low bits), random rotate.

Useful for 128-bit types that are split across two CPU registers.
"""
const PCG_XSL_RR = Val{:XSL_RR}
"""
One of PCG output method: fixed xorshift (to low bits), random rotate (both parts).

Useful for 128-bit types that are split across two CPU registers. Use this in need of an invertable 128-bit
RNG.
"""
const PCG_XSL_RR_RR = Val{:XSL_RR_RR}

const pcg_methods = (PCG_XSH_RS, PCG_XSH_RR, PCG_RXS_M_XS, PCG_XSL_RR, PCG_XSL_RR_RR)

"""
The `Union` of all the PCG method types: `PCG_XSH_RS`, `PCG_XSH_RR`, `PCG_RXS_M_XS`, `PCG_XSL_RR`, and `PCG_XSL_RR_RR`.
"""
const PCGMethod = Union{pcg_methods...}

include("pcg_list.jl")

include("bases.jl")

include("main.jl")

end


text:
abstract type OuterControl <: InverterComponent end

abstract type ActivePowerControl end
abstract type ReativePowerControl end

@def outercontrol_ports begin
    state_input = [:vpll_d, :vpll_q, :ε_pll, :vd_cap, :vq_cap, :id_o, :iq_o]
    inner_input = [Vdo_var, Vdo_var, ω_freq_estimator_var]
end

"""
Parameters of a Outer-Loop controller using a virtual inertia with VSM for active power controller
and a reactive power droop controller.
# Conmutable structor
```julia
VirtualInertiaQDroop(A, R)
```

# Arguments
*  `A`::Float64 : Active power controller using virtual inertia with VSM
*  `R`::Float64 : Reactive power controller using reactive power droop
"""
mutable struct VirtualInertiaQdroop{A <: ActivePowerControl,
                                   R <: ReativePowerControl} <: OuterControl
    active_power::A
    reactive_power::R
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

        function VirtualInertiaQdroop(active_power::A,
                                     reactive_power::R) where {A <: ActivePowerControl,
                                                               R <: ReativePowerControl}

            total_states = active_power.n_states + reactive_power.n_states
            @outercontrol_ports

            new{A,R}(active_power,
                     reactive_power,
                     total_states,
                     vcat(active_power.states, reactive_power.states),
                     Ports(state_input, inner_input))

        end

end

"""
Parameters of a Virtual Inertia with SRF using VSM for active power controller
# Conmutable structor
```julia
VirtualInertia(Ta, kd, kw, ωb)
```

# Arguments
*  `Ta`::Float64 : VSM inertia constant
*  `kd`::Float64 : VSM damping constant
*  `kw`::Float64 : frequency droop gain
*  `ωb`::Float64 : rated angular frequency
"""
mutable struct VirtualInertia <: ActivePowerControl
    Ta::Float64
    kd::Float64
    kω::Float64
    ωb::Float64
    n_states::Int64
    states::Vector{Symbol}

        function VirtualInertia(Ta::Float64,
                                kd::Float64,
                                kω::Float64,
                                ωb::Float64)

            new(Ta,
                kd,
                kω,
                ωb,
                2,
                [:δω_vsm, :δθ_vsm])

        end

end

"""
Parameters of a Reactive Power droop controller
# Conmutable structor
```julia
ReactivePowerDroop(kq, ωf)
```

# Arguments
*  `kq`::Float64 : reactive power droop gain
*  `ωf`::Float64 : reactive power filter cutoff frequency (rad/sec)
"""
mutable struct ReactivePowerDroop <: ReativePowerControl
    kq::Float64
    ωf::Float64
    n_states::Int64
    states::Vector{Symbol}

        function ReactivePowerDroop(kq::Float64,
                                   ωf::Float64)

            new(kq, ωf, 1, [:qm])
        end
end

ext: jl


text:
# helper functions used by PhyloPlot to export (then plot)
# a network as a phylo R object.

"""
    makemissing!(x::AbstractVector)

Turn to `missing` any element of `x` exactly equal to -1.0.
Used for branch lengths and γs. `x` needs to accept missing values.
If not, this can be done with `allowmissing(x)`.
"""
@inline function makemissing!(x::AbstractVector)
    for i in 1:length(x)
        if x[i] == -1.0
            x[i] = missing
        end
    end
end

"""
    majoredgematrix(net::HybridNetwork)

Matrix of major edges from `net` where edge[i,1] is the number of the
parent node of edge i and edge[i,2] is the number of the child node of edge i.
Assume `nodes_changed` was updated, to list nodes in pre-order.

# Examples

```jldoctest
julia> net = readTopology("(A,(B,(C,D)));");

julia> PhyloNetworks.resetNodeNumbers!(net);

julia> PhyloNetworks.majoredgematrix(net)
6×2 Matrix{Int64}:
 5  1
 5  6
 6  2
 6  7
 7  3
 7  4
```
"""
function majoredgematrix(net::HybridNetwork)
    edge = Matrix{Int}(undef, length(net.edge)-length(net.hybrid), 2) # major edges
    i = 1 #row index for edge matrix
    for n in net.nodes_changed # topological pre-order
        !n.leaf || continue # skip leaves: associate node with children edges
        for e in n.edge
            if e.node[e.isChild1 ? 2 : 1] == n && e.isMajor
                #exclude parent edge and minor hybrid edges
                edge[i,1] = n.number # parent
                edge[i,2] = e.node[e.isChild1 ? 1 : 2].number # child
                i += 1
            end
        end
    end
    return edge
end

"""
    majoredgelength(net::HybridNetwork)

Generate vector of edge lengths of major `net` edges organized in the same order
as the `edge` matrix created via `majoredgematrix`. Considers values of `-1.0` as
missing values, recognized as NA in `R`.
Output: vector allowing for missing values.

Assume `nodes_changed` was updated, to list nodes in pre-order.

# Examples

```jldoctest
julia> net = readTopology("(((A:3.1,(B:0.2)#H1:0.3::0.9),(C,#H1:0.3::0.1):1.1),D:0.7);");

julia> directEdges!(net); preorder!(net);

julia> PhyloNetworks.majoredgelength(net)
8-element Vector{Union{Missing, Float64}}:
  missing
 0.7
  missing
 1.1
  missing
 3.1
 0.3
 0.2
```
""" #"
function majoredgelength(net::HybridNetwork)
    edgeLength = Array{Union{Float64,Missing}}(undef, length(net.edge)-length(net.hybrid))
    i=1
    for n in net.nodes_changed # topological pre-order
        if !n.leaf
            for e in n.edge # for parent major edge below
                if e.node[e.isChild1 ? 2 : 1] == n && e.isMajor
                    edgeLength[i] = e.length
                    i=i+1
                end
            end
        end
    end
    #-1.0 interpreted as missing
    makemissing!(edgeLength)
    return edgeLength
end

"""
    minorreticulationmatrix(net::HybridNetwork)

Matrix of integers, representing the minor hybrid edges in `net`.
edge[i,1] is the number of the parent node of the ith minor hybrid edge,
and edge[i,2] is the number of its child node.
Node numbers may be negative, unless they were modified by `resetNodeNumbers!`.
Assumes correct `isChild1` fields.

# Examples

```julia-repl
julia> net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);");
julia> PhyloNetworks.minorreticulationmatrix(net)
1×2 Matrix{Int64}:
 -6  3
```
""" #"
function minorreticulationmatrix(net::HybridNetwork)
    reticulation = Matrix{Int}(undef, length(net.hybrid), 2) # initialize
    j = 1 # row index, row = reticulate edge
    for e in net.edge
        if !e.isMajor # minor (hybrid) edges only
            reticulation[j,1] = getParent(e).number
            reticulation[j,2] = getChild(e).number
            j += 1
        end
    end
    return reticulation
end

"""
    minorreticulationlength(net::HybridNetwork)

Vector of lengths for the minor hybrid edges, organized in the same order
as in the matrix created via `minorreticulationmatrix`.
Replace values of `-1.0` with missing values recognized by `R`.
Output: vector allowing for missing values.

# Examples

```jldoctest
julia> net = readTopology("(((A:3.1,(B:0.2)#H1:0.4::0.9),(C,#H1:0.3::0.1):1.1),D:0.7);");

julia> PhyloNetworks.minorreticulationlength(net)
1-element Vector{Union{Missing, Float64}}:
 0.3
```
""" #"
function minorreticulationlength(net::HybridNetwork)
    reticulationLength = Vector{Union{Float64,Missing}}(undef, 0) # initialize
    for e in net.edge
        if !e.isMajor #find minor hybrid edge
            push!(reticulationLength, e.length)
        end
    end
    makemissing!(reticulationLength)
    return reticulationLength
end

"""
    minorreticulationgamma(net::HybridNetwork)

Vector of minor edge gammas (inheritance probabilities) organized in the
same order as in the matrix created via `minorreticulationmatrix`.
Considers values of `-1.0` as missing values, recognized as NA in `R`.
Output: vector allowing for missing values.

# Examples

```julia-repl
julia> net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);");

julia> PhyloNetworks.minorreticulationgamma(net)
1-element Vector{Union{Float64, Missings.Missing}}:
 0.1
```
 """ #"
function minorreticulationgamma(net::HybridNetwork)
    reticulationGamma = Vector{Union{Float64,Missing}}(undef, 0) #initialize
    for e in net.edge
        if !e.isMajor # minor hybrid edges only
            push!(reticulationGamma, e.gamma)
        end
    end
    makemissing!(reticulationGamma)
    return reticulationGamma
end


text:
import Distributions
import DataFrames
import RCall
import StatsBase

function simulation_conditions()
    conds = Dict{Symbol, Int}()
    conds[:n] = rand(100:2000)
    conds[:pct_over] = rand(100:1000)
    conds[:pct_under] = rand(100:1000)
    prop_majority = rand(range(0.5, stop=0.99, length=100))
    conds[:n_majority] = round(Int, prop_majority * conds[:n])
    conds[:n_minority] = conds[:n] - conds[:n_majority]
    return conds
end

function simdata(sim_conditions)
    X = rand(sim_conditions[:n], 10)
    y = vcat(zeros(sim_conditions[:n_majority]), ones(sim_conditions[:n_minority]))
    X = hcat(X, y)
    return X
end

function smote_counts_jl_array(sim_conditions)
    Xy_array = simdata(sim_conditions)
    pct_over = sim_conditions[:pct_over]
    pct_under = sim_conditions[:pct_under]
    X2_array, y2_array = ClassImbalance._smote(
        Xy_array[:, 1:10],
        Xy_array[:, 11],
        5,
        pct_over,
        pct_under,
        )
    y2_counts_array = collect(values(StatsBase.countmap(y2_array)))
    return y2_counts_array
end

function smote_counts_jl_df(sim_conditions)
    Xy_array = simdata(sim_conditions)
    pct_over = sim_conditions[:pct_over]
    pct_under = sim_conditions[:pct_under]
    X = DataFrames.DataFrame()
    p = size(Xy_array, 2)
    @assert(p >= 2)
    for j = 1:(p - 1)
        column_name = Symbol(string("x", j))
        X[column_name] = Xy_array[:, j]
    end
    y = Xy_array[:, end]
    X2, y2 = ClassImbalance._smote(
        X,
        y,
        5,
        pct_over,
        pct_under,
        )
    y2_counts_df = collect(values(StatsBase.countmap(y2)))
    return y2_counts_df
end

function smote_counts_r(sim_conditions)
    n = sim_conditions[:n]
    pct_over = sim_conditions[:pct_over]
    pct_under = sim_conditions[:pct_under]
    n_minority = sim_conditions[:n_majority]
    n_majority = sim_conditions[:n_minority]
    RCall.@rput n
    RCall.@rput pct_over
    RCall.@rput pct_under
    RCall.@rput n_minority
    RCall.@rput n_majority
    RCall.R"""
        library("DMwR")
        X <- matrix(rnorm(n*10), ncol = 10)
        X <- cbind(1, X)
        dat <- data.frame(X)

        dat$y <- c(rep(0, n_majority), rep(1, n_minority))
        dat$y <- as.factor(dat$y)

        dat2 <- SMOTE(y ~ ., dat, perc.over = pct_over, perc.under = pct_under)
        df <- data.frame(table(dat2$y))
        y2_counts_r <- array(df[, 2])
    """
    RCall.@rget dat2
    RCall.@rget y2_counts_r
    return y2_counts_r
end

function smote_comparison_array(sim_conditions)
    # Comparing to output from SMOTE() in the DMwR package
    # Xy = simdata(sim_conditions)
    y2_counts_array = sort(smote_counts_jl_array(sim_conditions))
    y2_counts_r  = sort(smote_counts_r(sim_conditions))

    @info(string("y2_counts_r: "), y2_counts_r)
    @info(string("y2_counts_array: "), y2_counts_array)

    count_1_error_absolute =
        abs(abs(y2_counts_array[1]) - abs(y2_counts_r[1]))
    count_2_error_absolute =
        abs(abs(y2_counts_array[2]) - abs(y2_counts_r[2]))
    sum_error_absolute =
        abs(abs(sum(abs.(y2_counts_array))) - abs(sum(abs.(y2_counts_r))))

    count_1_error_relative_julia_percent =
        count_1_error_absolute/abs(y2_counts_array[1])*100.0
    count_1_error_relative_r_percent =
        count_1_error_absolute/abs(y2_counts_r[1])*100.0
    count_2_error_relative_julia_percent =
        count_2_error_absolute/abs(y2_counts_array[2])*100.0
    count_2_error_relative_r_percent =
        count_2_error_absolute/abs(y2_counts_r[2])*100.0
    sum_error_relative_julia_percent =
        sum_error_absolute/abs(sum(abs.(y2_counts_array)))*100.0
    sum_error_relative_r_percent =
        sum_error_absolute/abs(sum(abs.(y2_counts_r)))*100.0

    if !all(y2_counts_array .== y2_counts_r)
        @warn(
            string(
                "The counts from our smote function did not exactly match ",
                "the counts from R ",
                "(y2_counts_array was not equal to y2_counts_r).",
                ),
            y2_counts_array,
            y2_counts_r,
            sim_conditions,
            )
        @warn(
            string("Difference between counts: "),
            count_1_error_absolute,
            count_2_error_absolute,
            sum_error_absolute,
            count_1_error_relative_julia_percent,
            count_1_error_relative_r_percent,
            count_2_error_relative_julia_percent,
            count_2_error_relative_r_percent,
            sum_error_relative_julia_percent,
            sum_error_relative_r_percent,
            )
    end

    if count_1_error_relative_julia_percent > 1.0
        @warn(
            "count_1_error_relative_julia_percent: ",
            count_1_error_relative_julia_percent,
            )
        error("count_1_error_relative_julia_percent > 1.0")
    end
    if count_1_error_relative_r_percent > 1.0
        @warn(
            "count_1_error_relative_r_percent: ",
            count_1_error_relative_r_percent,
            )
        error("count_1_error_relative_r_percent > 1.0")
    end
    if count_2_error_relative_julia_percent > 1.0
        @warn(
            "count_2_error_relative_julia_percent: ",
            count_2_error_relative_julia_percent,
            )
        error("count_2_error_relative_julia_percent > 1.0")
    end
    if count_2_error_relative_r_percent > 1.0
        @warn(
            "count_2_error_relative_r_percent: ",
            count_2_error_relative_r_percent,
            )
        error("count_2_error_relative_r_percent > 1.0")
    end
    if sum_error_relative_julia_percent > 1.0
        @warn(
            "sum_error_relative_julia_percent: ",
            sum_error_relative_julia_percent,
            )
        error("sum_error_relative_julia_percent > 1.0")
    end
    if sum_error_relative_r_percent > 1.0
        @warn(
            "sum_error_relative_r_percent: ",
            sum_error_relative_r_percent,
            )
        error("sum_error_relative_r_percent > 1.0")
    end

    status = all(y2_counts_array .== y2_counts_r)

    return status
end

function smote_comparison_df(sim_conditions)
    # Comparing to output from SMOTE() in the DMwR package
    # Xy = simdata(sim_conditions)
    y2_counts_df = sort(smote_counts_jl_df(sim_conditions))
    y2_counts_r  = sort(smote_counts_r(sim_conditions))

    @info(string("y2_counts_r: "), y2_counts_r)
    @info(string("y2_counts_df: "), y2_counts_df)

    count_1_error_absolute =
        abs(abs(y2_counts_df[1]) - abs(y2_counts_r[1]))
    count_2_error_absolute =
        abs(abs(y2_counts_df[2]) - abs(y2_counts_r[2]))
    sum_error_absolute =
        abs(abs(sum(abs.(y2_counts_df))) - abs(sum(abs.(y2_counts_r))))

    count_1_error_relative_julia_percent =
        count_1_error_absolute/abs(y2_counts_df[1])*100.0
    count_1_error_relative_r_percent =
        count_1_error_absolute/abs(y2_counts_r[1])*100.0
    count_2_error_relative_julia_percent =
        count_2_error_absolute/abs(y2_counts_df[2])*100.0
    count_2_error_relative_r_percent =
        count_2_error_absolute/abs(y2_counts_r[2])*100.0
    sum_error_relative_julia_percent =
        sum_error_absolute/abs(sum(abs.(y2_counts_df)))*100.0
    sum_error_relative_r_percent =
        sum_error_absolute/abs(sum(abs.(y2_counts_r)))*100.0

    if !all(y2_counts_df .== y2_counts_r)
        @warn(
            string(
                "The counts from our smote function did not exactly match ",
                "the counts from R ",
                "(y2_counts_df was not equal to y2_counts_r).",
                ),
            y2_counts_df,
            y2_counts_r,
            sim_conditions,
            )
        @warn(
            string("Difference between counts: "),
            count_1_error_absolute,
            count_2_error_absolute,
            sum_error_absolute,
            count_1_error_relative_julia_percent,
            count_1_error_relative_r_percent,
            count_2_error_relative_julia_percent,
            count_2_error_relative_r_percent,
            sum_error_relative_julia_percent,
            sum_error_relative_r_percent,
            )
    end

    if count_1_error_relative_julia_percent > 1.0
        @warn(
            "count_1_error_relative_julia_percent: ",
            count_1_error_relative_julia_percent,
            )
        error("count_1_error_relative_julia_percent > 1.0")
    end
    if count_1_error_relative_r_percent > 1.0
        @warn(
            "count_1_error_relative_r_percent: ",
            count_1_error_relative_r_percent,
            )
        error("count_1_error_relative_r_percent > 1.0")
    end
    if count_2_error_relative_julia_percent > 1.0
        @warn(
            "count_2_error_relative_julia_percent: ",
            count_2_error_relative_julia_percent,
            )
        error("count_2_error_relative_julia_percent > 1.0")
    end
    if count_2_error_relative_r_percent > 1.0
        @warn(
            "count_2_error_relative_r_percent: ",
            count_2_error_relative_r_percent,
            )
        error("count_2_error_relative_r_percent > 1.0")
    end
    if sum_error_relative_julia_percent > 1.0
        @warn(
            "sum_error_relative_julia_percent: ",
            sum_error_relative_julia_percent,
            )
        error("sum_error_relative_julia_percent > 1.0")
    end
    if sum_error_relative_r_percent > 1.0
        @warn(
            "sum_error_relative_r_percent: ",
            sum_error_relative_r_percent,
            )
        error("sum_error_relative_r_percent > 1.0")
    end

    status = all(y2_counts_df .== y2_counts_r)

    return status
end

function run_comparisons_array(m)
    num_success = 0
    num_failure = 0
    for i = 1:m
        @info("Array simulation: $i")
        conds = simulation_conditions()
        sim_was_success = smote_comparison_array(conds)
        if sim_was_success
            num_success += 1
        else
            num_failure += 1
        end
    end

    @info(
        "Result of running array simulations: ",
        m,
        num_success,
        num_failure,
        )
    if !(num_success + num_failure == m)
        error("!(num_success + num_failure == m)")
    end
    percent_failure = num_failure/m*100.0
    @info(
        "percent_failure: ",
        percent_failure,
        )
    Test.@test(percent_failure < 1.0)
end

function run_comparisons_df(m)
    num_success = 0
    num_failure = 0
    for i = 1:m
        @info("DataFrame simulation: $i")
        conds = simulation_conditions()
        sim_was_success = smote_comparison_df(conds)
        if sim_was_success
            num_success += 1
        else
            num_failure += 1
        end
    end

    @info(
        "Result of running DataFrame simulations: ",
        m,
        num_success,
        num_failure,
        )
    if !(num_success + num_failure == m)
        error("!(num_success + num_failure == m)")
    end
    percent_failure = num_failure/m*100.0
    @info(
        "percent_failure: ",
        percent_failure,
        )
    Test.@test(percent_failure < 1.0)
end

run_comparisons_array(1000)

run_comparisons_df(1000)


text:
import YAML
config = YAML.load(open("config.yaml"))

# Points at which to create images
parameters = readdlm("parameters.dat")

parnames = config["emulator"]["parameters"]

npoints = size(parameters, 1)

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using HDF5

import PyPlot.plt
using LaTeXStrings

cmap = plt[:get_cmap]("plasma")

grd = config["grid"]
grid = Grid(grd)

# Load the parameters file from the configuration setting, and make this an object.
pars = convert_dict(config["parameters"], config["model"])
vel = pars.vel # [km/s]
npix = config["npix"] # number of pixels
species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

# Figure out how many plots we'll have.
ncols = 8
nrows = ceil(Int, npoints/ncols)

xx = 1.5 * 9
dx = 1.5
yy = (nrows + 1) * 1.5
dy = 1.5

fig, ax = plt[:subplots](nrows=nrows, ncols=ncols, figsize=(xx, yy))

for row=1:nrows
    for col=1:ncols
        iframe = col + (row - 1) * ncols

        if iframe > npoints
          # ax[row, col][:imshow](zeros((im_ny, im_nx)), cmap=cmap, vmin=0, vmax=20, extent=ext, origin="lower")
          continue
        else

          image = @sprintf("image_%3.3i.out", iframe)

          println(row, " ", col, " ", iframe)

          nu, incl = parameters[iframe,:]

          # Load image into SkyImage
          img = imread("images/" * image)
          pars.incl = incl
          skim = imToSky(img, pars.dpc)

          vvmax = maxabs(img.data)
          norm = PyPlot.matplotlib[:colors][:Normalize](0, vvmax)

          # Image needs to be flipped along RA dimension
          ext = (skim.ra[end], skim.ra[1], skim.dec[1], skim.dec[end])

          if col != 1 || row != nrows
              ax[row, col][:xaxis][:set_ticklabels]([])
              ax[row, col][:yaxis][:set_ticklabels]([])
          else
              ax[row, col][:set_xlabel](L"$\Delta \alpha$ ('')")
              ax[row, col][:set_ylabel](L"$\Delta \delta$ ('')")
          end

          #Flip the frame for Sky convention
          frame = flipdim(img.data[:,:,1], 2)

          ax[row, col][:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=cmap, norm=norm)
        end

    end
end

fig[:subplots_adjust](hspace=0.00, wspace=0.00, top=(yy - 0.5 * dy)/yy, bottom=(0.5 * dy)/yy, left=(0.5 * dx)/xx, right=(xx - 0.5 * dy)/xx)
plt[:savefig]("chmaps_all_params.png", dpi=600)

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
# # Point group 4mm (C₄ᵥ)

# ## Preamble
using LatticeTools
using Formatting: printfmt

function display_matrix(io::IO, matrix::AbstractMatrix; prefix::AbstractString="")
    width = ceil(Int, maximum(length("$item")+1 for item in matrix)/4)*4
    for row in eachrow(matrix)
        for (icol, col) in enumerate(row)
            icol == 1 && print(io, prefix)
            printfmt(io, "{:>$(width)s}", "$col")
        end
        println(io)
    end
end

# ## Irreps of 4mm
point_symmetry = project(PointSymmetryDatabase.get(13), [1 0 0; 0 1 0])
println("symmetry: ", point_symmetry.hermann_mauguin)
println("order: ", group_order(point_symmetry))
println("irreducible_representations:")
for irrep_index in 1:num_irreps(point_symmetry)
    println("- name: GM($irrep_index)")
    println("  elements:")
    for (ename, mat) in zip(element_names(point_symmetry), irrep(point_symmetry, irrep_index))
        println("  - name: \"$ename\"")
        m = Int.(real.(mat))
        if size(m, 1) == 1
            println("    matrix: $(first(m))")
        else
            println("    matrix: |-")
            display_matrix(stdout, m; prefix="      ")
        end
    end
end


text:
import JLD
using Distributed
using NetCDF

"A type for writing checkpoints."
struct Checkpointer <: OutputWriter
    dir::AbstractString
    filename_prefix::AbstractString
    output_frequency::Int
    padding::Int
end

"A type for writing NetCDF output."
mutable struct NetCDFOutputWriter <: OutputWriter
    dir::AbstractString
    folder::AbstractString
    output_frequency::Int
    padding::Int
    naming_scheme::Symbol
    compression::Int
    async::Bool
    onefile::Bool
    nctype::DataType     # NC_FLOAT (32bit), NC_DOUBLE (64bit) are denoted with integers ?!
    run_id::Int
    runfolder::AbstractString
    ncfile_u
    ncfile_v
    ncfile_w
    ncfile_T
    ncfile_S
    output_vars::Array{String,1}
end

"A type for writing Binary output."
mutable struct BinaryOutputWriter <: OutputWriter
    dir::AbstractString
    filename_prefix::AbstractString
    output_frequency::Int
    padding::Int
end

function Checkpointer(; dir=".", prefix="", frequency=1, padding=9)
    Checkpointer(dir, prefix, frequency, padding)
end

function NetCDFOutputWriter(; dir=".", folder="test", frequency=1, padding=4,
                              naming_scheme=:file_number, compression=3, async=false,
                              onefile=true, nctype=Float32, run_id=0)
    ncfile_u = nothing
    ncfile_v = nothing
    ncfile_w = nothing
    ncfile_T = nothing
    ncfile_S = nothing

    #output_vars = ["u","v","w","T","S"]
    output_vars = ["u","T"]
    runfolder = "run"*lpad(run_id,padding,"0")

    NetCDFOutputWriter(dir, folder, frequency, padding, naming_scheme,
                        compression, async, onefile, nctype, run_id, runfolder,
                        ncfile_u, ncfile_v, ncfile_w, ncfile_T, ncfile_S,
                        output_vars)
end

"Return the filename extension for the `OutputWriter` filetype."
ext(fw::OutputWriter) = throw("Not implemented.")
ext(fw::NetCDFOutputWriter) = ".nc"
ext(fw::Checkpointer) = ".jld"

filename(fw::Checkpointer, iteration) = fw.filename_prefix * "model_checkpoint_" * lpad(iteration, fw.padding, "0") * ext(fw)

function filename(fw, name, iteration)
    if fw.onefile
        if iteration == 0
            find_runs() = try
                    runlist = filter(x->startswith(x,"run"),readdir(joinpath(fw.dir,fw.folder)))
            catch
                    mkdir(joinpath(fw.dir,fw.folder))
                    runlist = filter(x->startswith(x,"run"),readdir(joinpath(fw.dir,fw.folder)))
            end
            runlist = find_runs()

            if length(runlist) > 0 # else use default value 0 for fw.run_id
                fw.run_id = maximum([parse(Int,id[4:end]) for id in runlist])+1
            end
            fw.runfolder = "run"*lpad(fw.run_id,fw.padding,"0")
            mkdir(joinpath(fw.dir,fw.folder,fw.runfolder))
        end
        name * ext(fw)
    else
        if fw.naming_scheme == :iteration
            fw.filename_prefix * name * lpad(iteration, fw.padding, "0") * ext(fw)
        elseif fw.naming_scheme == :file_number
            file_num = Int(iteration / fw.output_frequency)
            fw.filename_prefix * name * lpad(file_num, fw.padding, "0") * ext(fw)
        else
            throw(ArgumentError("Invalid naming scheme: $(fw.naming_scheme)"))
        end
    end
end

#
# Checkpointing functions
#

function write_output(model::Model{arch}, chk::Checkpointer) where arch <: Architecture
    filepath = joinpath(chk.dir, filename(chk, model.clock.iteration))

    forcing_functions = model.forcing

    # Do not include forcing functions and FFT plans. We want to avoid serializing
    # FFTW and CuFFT plans as serializing functions is not supported by JLD, and
    # seems like a tricky business in general.
    model.forcing = nothing
    model.poisson_solver = nothing

    println("WARNING: Forcing functions will not be serialized!")

    println("[Checkpointer] Serializing model to disk: $filepath")
    f = JLD.jldopen(filepath, "w", compress=true)
    JLD.@write f model
    close(f)

    println("[Checkpointer] Reconstructing FFT plans...")
    model.poisson_solver = init_poisson_solver(arch(), model.grid, model.stepper_tmp.fCC1)

    # Putting back in the forcing functions.
    model.forcing = forcing_functions

    return nothing
end

function restore_from_checkpoint(filepath)
    println("Deserializing model from disk: $filepath")
    f = JLD.jldopen(filepath, "r")
    model = read(f, "model");
    close(f)

    println("Reconstructing FFT plans...")
    model.poisson_solver = init_poisson_solver(arch(model)(), model.grid, model.stepper_tmp.fCC1)

    model.forcing = Forcing(nothing, nothing, nothing, nothing, nothing)
    println("WARNING: Forcing functions have been set to nothing!")

    return model
end

#
# Binary output function
#

function write_output(model::Model, fw::BinaryOutputWriter)
    for (field, field_name) in zip(fw.fields, fw.field_names)
        filepath = joinpath(fw.dir, filename(fw, field_name, model.clock.iteration))

        println("[BinaryOutputWriter] Writing $field_name to disk: $filepath")
        if model.metadata == :CPU
            write(filepath, field.data)
        elseif model.metadata == :GPU
            write(filepath, Array(field.data))
        end
    end
end

function read_output(model::Model, fw::BinaryOutputWriter, field_name, time)
    filepath = joinpath(fw.dir, filename(fw, field_name, time_step))
    arr = zeros(model.metadata.float_type, model.grid.Nx, model.grid.Ny, model.grid.Nz)

    open(filepath, "r") do fio
        read!(fio, arr)
    end

    return arr
end

#
# NetCDF output functions
#
# Eventually, we want to permit the user to flexibly define what is outputted.
# The user-defined function that produces output may involve computations, launching kernels,
# etc; so this API needs to be designed. For now, we simply save u, v, w, and T.

function write_output(model::Model, fw::NetCDFOutputWriter)
    fields = Dict(
        "xC" => collect(model.grid.xC),
        "yC" => collect(model.grid.yC),
        "zC" => collect(model.grid.zC),
        "xF" => collect(model.grid.xF),
        "yF" => collect(model.grid.yF),
        "zF" => collect(model.grid.zF),
        "t" => collect(model.clock.timevec[1:fw.output_frequency:end]),
        "u" => Array(model.velocities.u.data),
        "v" => Array(model.velocities.v.data),
        "w" => Array(model.velocities.w.data),
        "T" => Array(model.tracers.T.data),
        "S" => Array(model.tracers.S.data)
    )

    if fw.async
        # Execute asynchronously on worker 2.
        @async remotecall(write_output_netcdf, 2, fw, fields, model.clock.iteration)
    else
        write_output_netcdf(fw, fields, model.clock.iteration)
    end
    return nothing
end

function write_output_netcdf(fw::NetCDFOutputWriter, fields, iteration)

    xC, yC, zC = fields["xC"], fields["yC"], fields["zC"]
    xF, yF, zF = fields["xF"], fields["yF"], fields["zF"]
    t = fields["t"]

    u, v, w = fields["u"], fields["v"], fields["w"]
    T, S    = fields["T"], fields["S"]

    filename(fw, "", iteration)
    filepath(name::String,fw::NetCDFOutputWriter) = joinpath(fw.dir,fw.folder,fw.runfolder,name*ext(fw))

    all_var = [u,T]

    if iteration == 0 || ~fw.onefile    # initialise for onefile or output each timestep as single file

        xC_attr = Dict("longname" => "Locations of the cell centers in the x-direction.", "units" => "m")
        yC_attr = Dict("longname" => "Locations of the cell centers in the y-direction.", "units" => "m")
        zC_attr = Dict("longname" => "Locations of the cell centers in the z-direction.", "units" => "m")

        xF_attr = Dict("longname" => "Locations of the cell faces in the x-direction.", "units" => "m")
        yF_attr = Dict("longname" => "Locations of the cell faces in the y-direction.", "units" => "m")
        zF_attr = Dict("longname" => "Locations of the cell faces in the z-direction.", "units" => "m")

        t_attr = Dict("longname" => "time", "units" => "s")

        u_attr = Dict("longname" => "Velocity in the x-direction", "units" => "m/s")
        v_attr = Dict("longname" => "Velocity in the y-direction", "units" => "m/s")
        w_attr = Dict("longname" => "Velocity in the z-direction", "units" => "m/s")
        T_attr = Dict("longname" => "Temperature", "units" => "K")
        S_attr = Dict("longname" => "Salinity", "units" => "g/kg")

        # isfile(filepath) && rm(filepath)

        if fw.async
            println("[Worker $(Distributed.myid()): NetCDFOutputWriter] Writing fields to disk:"*joinpath(fw.folder,"run"*lpad(fw.run_id,4,"0")))
        else
            println("[NetCDFOutputWriter] Writing fields to disk: "*joinpath(fw.folder,"run"*lpad(fw.run_id,4,"0")))
        end

        xC_dim = NcDim("xC",length(xC),values=xC)
        yC_dim = NcDim("yC",length(yC),values=yC)
        zC_dim = NcDim("zC",length(zC),values=zC)

        xF_dim = NcDim("xF",length(xF),values=xF)
        yF_dim = NcDim("yF",length(yF),values=yF)
        zF_dim = NcDim("zF",length(zF),values=zF)

        t_dim = NcDim("t",0,unlimited=true)

        uvar = NcVar("u",[xF_dim,yC_dim,zC_dim,t_dim], t=fw.nctype, compress=fw.compression)
        vvar = NcVar("v",[xC_dim,yF_dim,zC_dim,t_dim], t=fw.nctype, compress=fw.compression)
        wvar = NcVar("w",[xC_dim,yC_dim,zF_dim,t_dim], t=fw.nctype, compress=fw.compression)
        Tvar = NcVar("T",[xC_dim,yC_dim,zC_dim,t_dim], t=fw.nctype, compress=fw.compression)
        Svar = NcVar("S",[xC_dim,yC_dim,zC_dim,t_dim], t=fw.nctype, compress=fw.compression)
        tvar = NcVar("t",t_dim,t=Int32)   # for some reason (bug?) this can't be fw.nctype

        if "u" in fw.output_vars fw.ncfile_u = NetCDF.create(filepath("u",fw),[uvar,tvar],mode=NC_NETCDF4) end
        if "v" in fw.output_vars fw.ncfile_v = NetCDF.create(filepath("v",fw),[vvar,tvar],mode=NC_NETCDF4) end
        if "w" in fw.output_vars fw.ncfile_w = NetCDF.create(filepath("w",fw),[wvar,tvar],mode=NC_NETCDF4) end
        if "T" in fw.output_vars fw.ncfile_T = NetCDF.create(filepath("T",fw),[Tvar,tvar],mode=NC_NETCDF4) end
        if "S" in fw.output_vars fw.ncfile_S = NetCDF.create(filepath("S",fw),[Svar,tvar],mode=NC_NETCDF4) end

        if "u" in fw.output_vars NetCDF.putatt(fw.ncfile_u,"u",u_attr) end
        if "v" in fw.output_vars NetCDF.putatt(fw.ncfile_v,"v",v_attr) end
        if "w" in fw.output_vars NetCDF.putatt(fw.ncfile_w,"w",w_attr) end
        if "T" in fw.output_vars NetCDF.putatt(fw.ncfile_T,"T",T_attr) end
        if "S" in fw.output_vars NetCDF.putatt(fw.ncfile_S,"S",S_attr) end

        if "u" in fw.output_vars NetCDF.putvar(fw.ncfile_u,"t",t) end
        if "v" in fw.output_vars NetCDF.putvar(fw.ncfile_v,"t",t) end
        if "w" in fw.output_vars NetCDF.putvar(fw.ncfile_w,"t",t) end
        if "T" in fw.output_vars NetCDF.putvar(fw.ncfile_T,"t",t) end
        if "S" in fw.output_vars NetCDF.putvar(fw.ncfile_S,"t",t) end

        if "u" in fw.output_vars NetCDF.putvar(fw.ncfile_u,"u",fw.nctype.(u),start=[1,1,1,1],count=[-1,-1,-1,1]) end
        if "v" in fw.output_vars NetCDF.putvar(fw.ncfile_v,"v",fw.nctype.(v),start=[1,1,1,1],count=[-1,-1,-1,1]) end
        if "w" in fw.output_vars NetCDF.putvar(fw.ncfile_w,"w",fw.nctype.(w),start=[1,1,1,1],count=[-1,-1,-1,1]) end
        if "T" in fw.output_vars NetCDF.putvar(fw.ncfile_T,"T",fw.nctype.(T),start=[1,1,1,1],count=[-1,-1,-1,1]) end
        if "S" in fw.output_vars NetCDF.putvar(fw.ncfile_S,"S",fw.nctype.(S),start=[1,1,1,1],count=[-1,-1,-1,1]) end

        if fw.onefile
            #TODO
            # store ncfiles in fw and keep open
        else
            NetCDF.close(fw.ncfile)
            #ncclose(filepath)
        end

    else # onefile output: append to existing file
        i = iteration
        freq = fw.output_frequency

        if "u" in fw.output_vars NetCDF.putvar(fw.ncfile_u,"u",fw.nctype.(u),start=[1,1,1,Int(i/freq)+1],count=[-1,-1,-1,1]) end
        if "v" in fw.output_vars NetCDF.putvar(fw.ncfile_v,"v",fw.nctype.(v),start=[1,1,1,Int(i/freq)+1],count=[-1,-1,-1,1]) end
        if "w" in fw.output_vars NetCDF.putvar(fw.ncfile_w,"w",fw.nctype.(w),start=[1,1,1,Int(i/freq)+1],count=[-1,-1,-1,1]) end
        if "T" in fw.output_vars NetCDF.putvar(fw.ncfile_T,"T",fw.nctype.(T),start=[1,1,1,Int(i/freq)+1],count=[-1,-1,-1,1]) end
        if "S" in fw.output_vars NetCDF.putvar(fw.ncfile_S,"S",fw.nctype.(S),start=[1,1,1,Int(i/freq)+1],count=[-1,-1,-1,1]) end

        # flush everything to file to view nc file while model is running
        #TODO this is currently not working
        if "u" in fw.output_vars NetCDF.sync(fw.ncfile_u) end
        if "v" in fw.output_vars NetCDF.sync(fw.ncfile_v) end
        if "w" in fw.output_vars NetCDF.sync(fw.ncfile_w) end
        if "T" in fw.output_vars NetCDF.sync(fw.ncfile_T) end
        if "S" in fw.output_vars NetCDF.sync(fw.ncfile_S) end
    end

    return nothing
end

function read_output(fw::NetCDFOutputWriter, field_name, iter)
    filepath = joinpath(fw.dir, filename(fw, "", iter))
    println("[NetCDFOutputWriter] Reading fields from disk: $filepath")

    if fw.onefile
        #TODO read only timestep iter
        field_data = ncread(filepath, field_name)
    else
        field_data = ncread(filepath, field_name)
    end

    ncclose(filepath)
    return field_data
end


text:
using Documenter, StrayCopses

makedocs(;
    modules=[StrayCopses],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/Arkoniak/StrayCopses.jl/blob/{commit}{path}#L{line}",
    sitename="StrayCopses.jl",
    authors="Andrey Oskin"
)

deploydocs(;
    repo="github.com/Arkoniak/StrayCopses.jl",
)


text:
# == PRIVATE METHODS BELOW HERE ======================================================================================= #
function _extract_species_from_phrase(phrase::String)::Array{String,1}

    # cut around the +'s
    species_array = Array{String,1}()
    species_substring_array = split(phrase, "+")

    # process each substring -
    for species_substring in species_substring_array

        # get the last value after we split for * -
        species_value = string(last(split(species_substring, "*")))

        # grab -
        push!(species_array, species_value)
    end

    # return -
    return species_array
end

function _extract_stochiometric_coefficient_from_phrase(phrase::String, species::String)::Float64

    # check: does this species occur in this phrase?
    if (occursin(species, phrase) == false)
        return 0.0
    end

    # ok: if we get here, then we have this species -
    # split around the +
    plus_split_product_array = split(phrase, "+")

    # process each component of the phrase -
    for species_substring in plus_split_product_array

        # get the values after we split for * -
        tmp_value_array = string.(split(species_substring, "*"))

        # ok: the last element will always be the species. If we have two elements the first is the st coeff -
        if (last(tmp_value_array) == species)
            if (length(tmp_value_array) == 2)
                return parse(Float64, first(tmp_value_array))
            else
                return 1.0
            end
        end
    end

    # default returns 0.0 -
    return 0.0
end

function _process_left_phrase(phrase::String, species::String)::Float64

    # get coeff -
    stm_coeff = _extract_stochiometric_coefficient_from_phrase(phrase, species)

    # ok: if non-zero, then multiply by -1
    if (stm_coeff > 0)
        return -1.0 * stm_coeff
    end

    # default -
    return stm_coeff
end

function _process_right_phrase(phrase::String, species::String)::Float64
    return _extract_stochiometric_coefficient_from_phrase(phrase, species)
end

# == PRIVATE METHODS ABOVE HERE ======================================================================================= #

# == PUBLIC METHODS BELOW HERE ======================================================================================== #
function build_stoichiometric_matrix(reactionTable::DataFrame, speciesSymbolArray::Array{String,1})::Some

    try

        # initialize -
        reaction_id_array = build_reaction_id_array(reactionTable) |> check
        species_symbol_array = speciesSymbolArray
        number_of_reactions = length(reaction_id_array)
        number_of_species = length(species_symbol_array)
        stoichiometric_matrx = zeros(number_of_species, number_of_reactions)

        # main -
        for species_index = 1:number_of_species # rows

            # grab species symbol -
            species_symbol = species_symbol_array[species_index]

            # does this symbol appear in a reaction?
            for reaction_index = 1:number_of_reactions

                # get the L and R phrases -
                L = reactionTable[reaction_index,:forward]
                R = reactionTable[reaction_index,:reverse]

                # get the stcoeff from the L phrase -
                L_st_coeff = _process_left_phrase(L, species_symbol)

                # get the stcoeff from the R phrase -
                R_st_coeff = _process_right_phrase(R, species_symbol)

                # update the stoichiometric_matrx -
                stoichiometric_matrx[species_index,reaction_index] = (L_st_coeff + R_st_coeff)
            end
        end

        # return -
        return Some(stoichiometric_matrx)
    catch error
        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_stoichiometric_matrix(reactionTable::DataFrame)::Some

    try

        # initialize -
        species_symbol_array = build_species_symbol_array(reactionTable) |> check
        stoichiometric_matrx = build_stoichiometric_matrix(reactionTable, species_symbol_array) |> check

        # return -
        return Some(stoichiometric_matrx)
    catch error
        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_flux_bounds_array(reactionTable::DataFrame; defaultFluxBoundValue::Float64=100.0)::Some

    try

        # initialize -
        reaction_id_array = build_reaction_id_array(reactionTable) |> check
        number_of_reactions = length(reaction_id_array)
        flux_bounds_array = Array{Float64,2}(undef, number_of_reactions, 2)

        # iterate through the reaction id's and determine if the reaction is reversible. If so, then build the bounds array
        # using the default values -
        for (index, id_value) in enumerate(reaction_id_array)

            # is this reaction reversible?
            is_reversible = reactionTable[index,:reversibility]
            if (is_reversible == true)
                flux_bounds_array[index,1] = -1.0 * defaultFluxBoundValue
                flux_bounds_array[index,2] = defaultFluxBoundValue
            elseif (is_reversible == false)
                flux_bounds_array[index,1] = 0.0
                flux_bounds_array[index,2] = defaultFluxBoundValue
            end
        end

        # return -
        return Some(flux_bounds_array)
    catch error
        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_species_bounds_array(speciesTable::DataFrame)::Some

    try

        # initialize -
        (number_of_species, _) = size(speciesTable)
        species_bounds_array = Array{Float64,2}(undef, number_of_species, 2)

        # main -
        for species_index = 1:number_of_species

            # get data for bounds array -
            lower_bound = speciesTable[species_index, :lower_bound]
            upper_bound = speciesTable[species_index, :upper_bound]

            # package -
            species_bounds_array[species_index,1] = lower_bound
            species_bounds_array[species_index,2] = upper_bound
        end

        # return -
        return Some(species_bounds_array)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_species_table(reactionTable::DataFrame)::Some

    try

        # initialize -
        species_symbol_array = build_species_symbol_array(reactionTable) |> check
        df_species_table = DataFrame(symbol=String[], name=Union{String,Missing}[], lower_bound=Float64[], upper_bound=Float64[])

        # main -
        for species_symbol in species_symbol_array

            # create a tuple w/the data row -
            data_row = (species_symbol, missing, 0.0, 0.0)

            # add to the df -
            push!(df_species_table, data_row)
        end

        # return -
        return Some(df_species_table)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_species_symbol_array(reactionTable::DataFrame)::Some

    try

        # initialize -
        species_symbol_array = Array{String,1}()
        reaction_phrase_array = Array{String,1}()
        (number_of_reactions, _) = size(reactionTable)

        # populate the array of reaction phrases -
        for reaction_index = 1:number_of_reactions

            # get the L and R phrases -
            L = reactionTable[reaction_index,:forward]
            R = reactionTable[reaction_index,:reverse]

            # capture -
            push!(reaction_phrase_array, L)
            push!(reaction_phrase_array, R)
        end

        # ok, so now lets populate the tmp species set -
        for reaction_phrase in reaction_phrase_array

            # get species sub array -
            species_in_reaction_phrase = _extract_species_from_phrase(reaction_phrase)

            # push into set -
            for tmp_species_symbol in species_in_reaction_phrase
                if (in(tmp_species_symbol, species_symbol_array) == false)
                    push!(species_symbol_array, tmp_species_symbol)
                end
            end
        end

        # return -
        return Some(species_symbol_array)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_objective_coefficient_array(reactionTable::DataFrame)::Some

    try

        # initialize -
        (number_of_reactions, _) = size(reactionTable)

        # default behavior -
        obj_coeff_array = zeros(number_of_reactions)

        # return -
        return Some(obj_coeff_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_reaction_id_array(reactionTable::DataFrame)::Some

    try

        # get the id col -
        id_col = reactionTable[!,:id]

        # return -
        return Some(id_col)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_transport_reaction_table(species_pair_array::Array{Pair{String,String},1})::Some

    try

        # initailize -
        id_array = String[]
        forward_reaction_string = String[]
        reverse_reaction_string = String[]
        reversibility_array = Bool[]
        ec_number_array = Union{Missing,String}[]
        lower_bound_array = Union{Missing,Float64}[]
        upper_bound_array = Union{Missing,Float64}[]

        # ok, so we have an array of pairs - assume all transport reactions are *reversible* -
        for transport_pair in species_pair_array

            # get the a => b
            first_species = transport_pair.first
            second_species = transport_pair.second

            # write the record -
            push!(id_array, "$(first_species)_exchange")
            push!(forward_reaction_string, first_species)
            push!(reverse_reaction_string, second_species)
            push!(reversibility_array, true)
            push!(ec_number_array, missing)
            push!(lower_bound_array, -Inf)
            push!(upper_bound_array, Inf)
        end

        # package into DataFrame -
        reaction_dataframe = DataFrame(id=id_array, forward_reaction=forward_reaction_string, reverse_reaction=reverse_reaction_string, reversibility=reversibility_array,
            ec_number=ec_number_array, lower_bound=lower_bound_array, upper_bound=upper_bound_array)

        # return -
        return Some(reaction_dataframe)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function build_atom_matrix(formula_array::Array{String,1};
        elements=["C","H","N","O","P","S"])::Some

    try

        # initialize -
        number_of_elements = length(elements)
        number_of_compounds = length(formula_array)
        atom_matrix = Array{Number,2}(undef,number_of_elements,number_of_compounds)

        # ok, so let's process this list of formulas -
        for (formula_index,formula) in enumerate(formula_array)

            # parse -
            atom_dictionary = parse_molecular_formula_string(formula) |> check

            # fill in the elements -
            for (element_index,element) in enumerate(elements)

                # check: do we have this element?
                get!(atom_dictionary, element, 0)

                # go -
                atom_matrix[element_index, formula_index] = atom_dictionary[element]
            end
        end

        # return -
        return Some(atom_matrix)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function parse_molecular_formula_string(formula::String)::Some

    try

        # test -
        atom_dictionary = Dict();
        local_array = Array{Char,1}()

        # turn string into char array -
        formula_char_array = collect(formula);

        # add extra 1 if last char is a letter -
        if (isnumeric(last(formula_char_array))== false)
            push!(formula_char_array,'1')
        end

        # read from the bottom -
        reverse!(formula_char_array)

        # how many chars do we have?
        while (isempty(formula_char_array) == false)

            # clean out the array from the last pass -
            empty!(local_array)

            # grab the next value -
            next_value = pop!(formula_char_array)
            if (isnumeric(next_value) == false)

                # we have an element -> read until I hit another element -
                is_ok_to_loop = true
                while (is_ok_to_loop)

                    if (isempty(formula_char_array) == true)
                        break
                    end

                    read_one_ahead = pop!(formula_char_array)
                    if (isnumeric(read_one_ahead) == true)
                        push!(local_array,read_one_ahead)
                    else

                        # ok: so if we get here - then we read the next char, but it was a
                        # letter (element) - so we need to puch it back on the stack
                        push!(formula_char_array,read_one_ahead)

                        # reverse -
                        is_ok_to_loop = false
                    end
                end

                # we need to turn local array into a string -
                buffer = ""
                [buffer*=string(x) for x in local_array]
                atom_dictionary[string(next_value)] = parse(Int64,buffer)
            end
        end

        # return -
        return Some(atom_dictionary)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

# == PUBLIC METHODS ABOVE HERE ======================================================================================== #


text:
module Job
using Cairn.Client

using Cairn.RPC: PrefixKey, Key, KeySpan

immutable t
    # NOT TRUE; FUTURE WORK. ALL MUST BE .jl
    # if *_executable ends in .jl, then it will try to be
    # imported. Otherwise, it will be piped via stdin
    dfs::Vector{Cairn.Client.t}
    #=
    A mapfn executable should have a single function declaration
    that takes an original file name, an actual file name, a results dir
    it should read from the actual file name
    and write output files of the form <results dir>/hash(ouput key) or <results dir>/key
    depending on the key type
    =#
    map_executable::Key
    reduce_executable::Key
    files::KeySpan
    workers::Int64
    result_replication::Int64
    result_dir::PrefixKey
    job_id::ASCIIString
    function t(dfs::Cairn.Client.t, map_exe::ASCIIString,
               reduce_exe::ASCIIString, data::FileSpan,
               local_data::Vector{ASCIIString}, workers::Int64,
               exe_replication::Int64, upload_replication::Int64,
               result_replication::Int64,
               results_dir::PrefixKey = PrefixKey(string(uuid4())),
               job_id::ASCIIString = string(uuid4()))
        # upload code to the DFS
        map_exe_name = string(uuid4())
        map_exe_key = "$job_id.$map_exe_name.$map_exe"
        open(map_exe, "r") do f
            mm_data = Mmap.mmap(f, Vector{UInt8})
            Cairn.Client.upload(dfs, map_exe_key,  mm_data, exe_replication)
        end
        reduce_exe_name = string(uuid4())
        reduce_exe_key = "$job_id.$reduce_exe_name.$reduce_exe"
        open(reduce_exe, "r") do f
            mm_data = Mmap.mmap(f, Vector{UInt8})
            Cairn.Client.upload(dfs, reduce_exe_key,  mm_data, exe_replication)
        end

        # upload files to the DFS (if not already there)
        @sync for (i,key) = enumerate(local_data)
            @async open(key, "r") do f
                data = Mmap.mmap(f, Vector{UInt8})
                Cairn.Client.upload(dfs, "$job_id.uploaded.$i",  data, upload_replication)
            end

        end
        if length(local_data) != 0
            data = [data; PrefixKey("$job_id.uploaded.")]
        end

        new(Cairn.Client.getMasters(dfs), map_exe_key, reduce_exe_key, data, workers, results_replication, results_dir, job_id)
    end
    t(job::t, sub::FileSpan) = new(job.dfs, job.map_executable, job.reduce_executable, sub, 0,
    job.result_replication, job.result_dir, job.job_id)
end

end


text:
export write_extxyz
export read_extxyz

import ExtXYZ

"""
    write_extxyz(file, atoms, R, cell)

Uses ExtXYZ.jl to write a .extxyz file

R can be either a single configuration (Matrix) or many (Vector{<:Matrix}).
"""
write_extxyz(file, atoms, R::Matrix, cell) =
    ExtXYZ.write_frame(file, to_extxyz_dict(atoms, R, cell))

write_extxyz(file, atoms, R::Vector{<:Matrix}, cell) =
    ExtXYZ.write_frames(file, to_extxyz_dict.(Ref(atoms), R, Ref(cell)))

"""
    read_extxyz(file)

Uses ExtXYZ.jl to read a .extxyz file
"""
function read_extxyz(file)

    dict = ExtXYZ.read_frame(file)
    atoms, R, cell = from_extxyz_dict(dict)
    positions = [R]

    for dict in ExtXYZ.read_frames(file, Iterators.countfrom(2))
        _, R, _ = from_extxyz_dict(dict)
        push!(positions, R)
    end

    return atoms, positions, cell
end

function to_extxyz_dict(atoms::Atoms, R::Matrix, cell::PeriodicCell)

    dict = Dict{String, Any}()
    dict["N_atoms"] = length(atoms)

    dict["pbc"] = cell.periodicity

    dict["cell"] = au_to_ang.(permutedims(cell.vectors, (2,1)))

    dict["arrays"] = Dict{String,Any}()
    dict["arrays"]["species"] = String.(atoms.types) |> Array
    dict["arrays"]["pos"] = au_to_ang.(R)

    dict["info"] = Dict{String,Any}()

    return dict
end

function from_extxyz_dict(dict::Dict{String,Any})

    atoms = Atoms(Symbol.(dict["arrays"]["species"]))
    R = ang_to_au.(dict["arrays"]["pos"])
    c = permutedims(dict["cell"], (2,1))
    cell = PeriodicCell{eltype(c)}(ang_to_au.(c), dict["pbc"])

    return atoms, R, cell
end


text:
using BlockArrays, LinearAlgebra, Base64

struct PartiallyImplementedBlockVector <: AbstractBlockArray{Float64,1} end

@testset "partially implemented block array" begin
    # the error thrown before was incorrect
    A = PartiallyImplementedBlockVector()
    try
        getblock(A, 1)
    catch err
        @test err isa ErrorException && err.msg == "getblock for PartiallyImplementedBlockVector is not implemented"
    end
    try
        getblock!(zeros(5), A, Block(1))
    catch err
        @test err isa ErrorException && err.msg == "getblock! for PartiallyImplementedBlockVector is not implemented"
    end
    try
        BlockArrays.setblock!(A, zeros(5), Block(1))
    catch err
        @test err isa ErrorException && err.msg == "setblock! for PartiallyImplementedBlockVector is not implemented"
    end
    try
        BlockArrays.blocksize(A, 2)
    catch err
        @test err isa ErrorException && err.msg == "blocksizes for PartiallyImplementedBlockVector is not implemented"
    end
end

@testset "Array block interface" begin
    @test blocksizes(1) == BlockArrays.BlockSizes{0}()
    @test 1[Block()] == 1

    A = randn(5)
    @test blocksizes(A) == BlockArrays.BlockSizes([5])
    @test A[Block(1)] == A
    view(A,Block(1))[1] = 2
    @test A[1] == 2
    @test_throws BoundsError A[Block(2)]

    A = randn(5,5)
    @test A[Block(1,1)] == A
end

@testset "Triangular/Symmetric/Hermitian block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (1:4))
    A .= reshape(1:length(A), size(A))

    @test blocksizes(UpperTriangular(A)) == blocksizes(Symmetric(A)) == blocksizes(A)
    @test UpperTriangular(A)[Block(2,2)] == UpperTriangular(A[2:3,2:3])
    @test UpperTriangular(A)[Block(2,3)] == A[2:3,4:6]
    @test UpperTriangular(A)[Block(3,2)] == zeros(3,2)
    @test Symmetric(A)[Block(2,2)] == Symmetric(A[2:3,2:3])
    @test Symmetric(A)[Block(2,3)] == A[2:3,4:6]
    @test Symmetric(A)[Block(3,2)] == transpose(A[2:3,4:6])
    @test Hermitian(A)[Block(2,2)] == Hermitian(A[2:3,2:3])
    @test Hermitian(A)[Block(2,3)] == A[2:3,4:6]
    @test Hermitian(A)[Block(3,2)] == A[2:3,4:6]'

    @test stringmime("text/plain", UpperTriangular(A)) == "10×10 UpperTriangular{Complex{Float64},PseudoBlockArray{Complex{Float64},2,Array{Complex{Float64},2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}}}:\n 1.0+0.0im  │  11.0+0.0im  21.0+0.0im  │  31.0+0.0im  41.0+0.0im  51.0+0.0im  │  61.0+0.0im  71.0+0.0im  81.0+0.0im   91.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │  12.0+0.0im  22.0+0.0im  │  32.0+0.0im  42.0+0.0im  52.0+0.0im  │  62.0+0.0im  72.0+0.0im  82.0+0.0im   92.0+0.0im\n     ⋅      │       ⋅      23.0+0.0im  │  33.0+0.0im  43.0+0.0im  53.0+0.0im  │  63.0+0.0im  73.0+0.0im  83.0+0.0im   93.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │  34.0+0.0im  44.0+0.0im  54.0+0.0im  │  64.0+0.0im  74.0+0.0im  84.0+0.0im   94.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅      45.0+0.0im  55.0+0.0im  │  65.0+0.0im  75.0+0.0im  85.0+0.0im   95.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅      56.0+0.0im  │  66.0+0.0im  76.0+0.0im  86.0+0.0im   96.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │  67.0+0.0im  77.0+0.0im  87.0+0.0im   97.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅      78.0+0.0im  88.0+0.0im   98.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅      89.0+0.0im   99.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅           ⋅      100.0+0.0im"
end

@testset "Adjoint/Transpose block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (2:5))
    A .= randn.() .+ randn.().*im

    @test blocksizes(A') == BlockArrays.BlockSizes(2:5, 1:4)
    @test blocksizes(Transpose(A)) == BlockArrays.BlockSizes(2:5, 1:4)

    @test A'[Block(2,2)] == A[Block(2,2)]' == A[2:3,3:5]'
    @test transpose(A)[Block(2,2)] == transpose(A[2:3,3:5])
    @test A'[Block(2,3)] == A[Block(3,2)]'
    @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

    @test BlockArray(A') == A'
    @test BlockArray(transpose(A)) == transpose(A)
end

@testset "Diagonal BlockArray" begin
    A = mortar(Diagonal(fill([1 2],2)))
    @test A isa BlockMatrix{Int,Diagonal{Matrix{Int}, Vector{Matrix{Int}}}}
    @test A[Block(1,2)] == [0 0]
    @test_throws BlockBoundsError A[Block(1,3)]
    @test A == [1 2 0 0; 0 0 1 2]
end

@testset "non-standard BlockSizes" begin
    A = BlockArray([1 2; 3 4], BlockArrays.BlockSizes((1:3, 1:3)))
    @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},BlockArrays.BlockSizes{2,NTuple{2,UnitRange{Int}}}}
    A = BlockArray([1 2; 3 4], BlockArrays.BlockSizes((1:3, [1,2,3])))
    @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},BlockArrays.BlockSizes{2,Tuple{UnitRange{Int},Vector{Int}}}}
end


text:
##### Beginning of file

function _print_welcome_message()::Nothing
    mirrorupdater_version::VersionNumber = version()
    mirrorupdater_pkgdir::String = package_directory()
    @info(string("This is MirrorUpdater, version ",mirrorupdater_version,),)
    @info(string("For help, please visit https://github.com/UnofficialJuliaMirror/MirrorUpdater.jl",),)
    @debug(string("MirrorUpdater package directory: ",mirrorupdater_pkgdir,),)
    return nothing
end

##### End of file


text:
# Tables.jl interface

Tables.istable(::Type{<:TableRowIterator}) = true
Tables.rowaccess(::Type{<:TableRowIterator}) = true
Tables.rows(itr::TableRowIterator) = itr
Tables.schema(itr::TableRowIterator) = Tables.Schema(itr.index.column_labels, fill(Any, length(itr.index.column_labels)))
Tables.columnnames(tr::TableRow) = tr.index.column_labels
Tables.getcolumn(tr::TableRow, nm::Symbol) = getdata(tr, nm)
Tables.getcolumn(tr::TableRow, i::Integer) = getdata(tr, i)

function table_to_arrays(x)
    _as_vector(y::AbstractVector) = y
    _as_vector(y) = collect(y)
    columns = Any[_as_vector(c) for c in Tables.Columns(x)]
    names = collect(Symbol, Tables.columnnames(x))
    return columns, names
end

"""
    writetable(filename, table; [overwrite], [sheetname])

Write Tables.jl `table` to the specified filename
"""
writetable(filename::AbstractString, x; kw...) = writetable(filename, table_to_arrays(x)...; kw...)

"""
    writetable(filename::AbstractString, tables::Vector{Pair{String, T}}; overwrite::Bool=false)
    writetable(filename::AbstractString, tables::Pair{String, Any}...; overwrite::Bool=false)
"""
function writetable(filename::AbstractString, tables::Vector{<:Pair}; kw...)
    data = [(name, table_to_arrays(x)...) for (name, x) in tables]
    return writetable(filename, data; kw...)
end
writetable(filename::AbstractString, tables::Pair{<:String, <:Any}...; kw...) = writetable(filename, collect(tables); kw...)

"""
    writetable!(sheet::Worksheet, table; anchor_cell::CellRef=CellRef("A1")))

Write Tables.jl `table` to the specified sheet
"""
writetable!(sheet::Worksheet, x, kw...) = writetable!(sheet, table_to_arrays(x)...; kw...)


text:
"""Transit function to call either the rhs_linear or the rhs_nonlinear."""
function rhs!(  u::Array{T,2},
                v::Array{T,2},
                η::Array{T,2},
                Diag::DiagnosticVars{T,Tprog},
                S::ModelSetup{T,Tprog},
                t::Int) where {T,Tprog}

    @unpack dynamics = S.parameters

    if dynamics == "linear"
        rhs_linear!(u,v,η,Diag,S,t)
    else
        rhs_nonlinear!(u,v,η,Diag,S,t)
    end
end

"""Tendencies du,dv of

        ∂u/∂t = qhv - ∂(1/2*(u²+v²) + gη)/∂x + Fx
        ∂v/∂t = -qhu - ∂(1/2*(u²+v²) + gη)/∂y + Fy

the nonlinear shallow water equations."""
function rhs_nonlinear!(u::AbstractMatrix,
                        v::AbstractMatrix,
                        η::AbstractMatrix,
                        Diag::DiagnosticVars,
                        S::ModelSetup,
                        t::Int)

    @unpack h,h_u,h_v,U,V = Diag.VolumeFluxes
    @unpack H = S.forcing
    @unpack ep = S.grid

    UVfluxes!(u,v,η,Diag,S)              # U,V needed for PV advection and in the continuity equation
    if S.grid.nstep_advcor == 0              # evaluate every RK substep
        advection_coriolis!(u,v,η,Diag,S)    # PV and non-linear Bernoulli terms
    end
    PVadvection!(Diag,S)                 # advect the PV with U,V

    # Bernoulli potential - recalculate for new η, KEu,KEv are only updated in advection_coriolis
    @unpack p,KEu,KEv,dpdx,dpdy = Diag.Bernoulli
    @unpack g,scale,scale_inv = S.constants
    g_scale = g*scale
    bernoulli!(p,KEu,KEv,η,g_scale,ep,scale_inv)
    ∂x!(dpdx,p)
    ∂y!(dpdy,p)

    # adding the terms
    momentum_u!(Diag,S,t)
    momentum_v!(Diag,S,t)
end

"""Tendencies du,dv of

        ∂u/∂t = fv - g∂η/∂x + Fx
        ∂v/∂t = -fu - g∂η/∂y + Fy

the linear shallow water equations."""
function rhs_linear!(   u::AbstractMatrix,
                        v::AbstractMatrix,
                        η::AbstractMatrix,
                        Diag::DiagnosticVars,
                        S::ModelSetup,
                        t::Int)

    @unpack h,h_u,h_v,U,V,dUdx,dVdy = Diag.VolumeFluxes
    @unpack g,scale = S.constants
    @unpack ep = S.grid

    # Pressure gradient
    @unpack dpdx,dpdy = Diag.Bernoulli
    g_scale = g*scale
    ∂x!(dpdx,g_scale*η)
    ∂y!(dpdy,g_scale*η)

    # Coriolis force
    @unpack qhv,qhu,v_u,u_v = Diag.Vorticity
    @unpack f_u,f_v = S.grid
    Ixy!(v_u,v)
    Ixy!(u_v,u)
    fv!(qhv,f_u,v_u,ep)
    fu!(qhu,f_v,u_v,ep)

    # adding the terms
    momentum_u!(Diag,S,t)
    momentum_v!(Diag,S,t)
end

""" Update advective and Coriolis tendencies."""
function advection_coriolis!(   u::Array{T,2},
                                v::Array{T,2},
                                η::Array{T,2},
                                Diag::DiagnosticVars{T,Tprog},
                                S::ModelSetup{T,Tprog}) where {T,Tprog}

    @unpack h = Diag.VolumeFluxes
    @unpack h_q,dvdx,dudy = Diag.Vorticity
    @unpack u²,v²,KEu,KEv = Diag.Bernoulli
    @unpack ep,f_q = S.grid

    Ixy!(h_q,h)

    # off-diagonals of stress tensor ∇(u,v)
    ∂x!(dvdx,v)
    ∂y!(dudy,u)

    # non-linear part of the Bernoulli potential
    speed!(u²,v²,u,v)
    Ix!(KEu,u²)
    Iy!(KEv,v²)

    # Potential vorticity update
    PV!(Diag,S)

    @unpack q = Diag.Vorticity          # now load what just has been calculated

    # Linear combinations of the potential vorticity q
    if S.parameters.adv_scheme == "Sadourny"
        @unpack q_u,q_v = Diag.Vorticity
        Iy!(q_u,q)
        Ix!(q_v,q)
    elseif S.parameters.adv_scheme == "ArakawaHsu"
        @unpack qα,qβ,qγ,qδ = Diag.ArakawaHsu
        AHα!(qα,q)
        AHβ!(qβ,q)
        AHγ!(qγ,q)
        AHδ!(qδ,q)
    end
end

"""Layer thickness h obtained by adding sea surface height η to bottom height H."""
function thickness!(h::AbstractMatrix,η::AbstractMatrix,H::AbstractMatrix)
    m,n = size(h)
    @boundscheck (m,n) == size(η) || throw(BoundsError())
    @boundscheck (m,n) == size(H) || throw(BoundsError())

    @inbounds for i in eachindex(η)
        h[i] = η[i] + H[i]
    end
end

"""Squared velocities u²,v²."""
function speed!(u²::AbstractMatrix,
                v²::AbstractMatrix,
                u::AbstractMatrix,
                v::AbstractMatrix)
    m,n = size(u²)
    @boundscheck (m,n) == size(u) || throw(BoundsError())

    @inbounds for i in eachindex(u)
        u²[i] = u[i]^2
    end

    m,n = size(v²)
    @boundscheck (m,n) == size(v) || throw(BoundsError())

    @inbounds for i in eachindex(v)
        v²[i] = v[i]^2
    end
end

"""Bernoulli potential p = 1/2*(u² + v²) + gη."""
function bernoulli!(p::Array{T,2},
                    KEu::Array{T,2},
                    KEv::Array{T,2},
                    η::Array{T,2},
                    g::T,
                    ep::Int,
                    scale_inv::T) where {T<:AbstractFloat}

    m,n = size(p)
    @boundscheck (m+ep,n+2) == size(KEu) || throw(BoundsError())
    @boundscheck (m+2,n) == size(KEv) || throw(BoundsError())
    @boundscheck (m,n) == size(η) || throw(BoundsError())

    one_half_scale_inv = convert(T,0.5)*scale_inv

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            p[i,j] = one_half_scale_inv*(KEu[i+ep,j+1] + KEv[i+1,j]) + g*η[i,j]
        end
    end
end

"""Coriolis term f*v. """
function fv!(   qhv::AbstractMatrix,
                f_u::AbstractMatrix,
                v_u::AbstractMatrix,
                ep::Int)

    m,n = size(qhv)
    @boundscheck (m,n) == size(f_u) || throw(BoundsError())
    @boundscheck (m+4-ep,n+2) == size(v_u) || throw(BoundsError())
    @boundscheck ep == 1 || ep == 0 || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            qhv[i,j] = f_u[i,j]*v_u[i+2-ep,j+1]
        end
    end
end

"""Coriolis term f*u. """
function fu!(   qhu::AbstractMatrix,
                f_v::AbstractMatrix,
                u_v::AbstractMatrix,
                ep::Int)

    m,n = size(qhu)
    @boundscheck (m,n) == size(f_v) || throw(BoundsError())
    @boundscheck (m+2+ep,n+4) == size(u_v) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            qhu[i,j] = f_v[i,j]*u_v[i+1+ep,j+2]
        end
    end
end

"""Sum up the tendencies of the non-diffusive right-hand side for the u-component."""
function momentum_u!(   Diag::DiagnosticVars{T,Tprog},
                        S::ModelSetup,
                        t::Int) where {T,Tprog}

    @unpack du = Diag.Tendencies
    @unpack qhv = Diag.Vorticity
    @unpack dpdx = Diag.Bernoulli
    @unpack Fx = S.forcing
    @unpack ep,halo = S.grid

    m,n = size(du) .- (2halo,2halo)     # cut off the halo
    @boundscheck (m,n) == size(qhv) || throw(BoundsError())
    @boundscheck (m+2-ep,n+2) == size(dpdx) || throw(BoundsError())
    @boundscheck (m,n) == size(Fx) || throw(BoundsError())

    if S.parameters.seasonal_wind_x
        @unpack ωFx = S.constants
        Fxt = Ftime(T,t,ωFx)
    else
        Fxt = one(T)
    end

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            du[i+2,j+2] = (Tprog(qhv[i,j]) - Tprog(dpdx[i+1-ep,j+1])) + Tprog(Fxt*Fx[i,j])
        end
    end
end

"""Sum up the tendencies of the non-diffusive right-hand side for the v-component."""
function momentum_v!(   Diag::DiagnosticVars{T,Tprog},
                        S::ModelSetup,
                        t::Int) where {T,Tprog}

    @unpack dv = Diag.Tendencies
    @unpack qhu = Diag.Vorticity
    @unpack dpdy = Diag.Bernoulli
    @unpack Fy = S.forcing
    @unpack halo = S.grid

    m,n = size(dv) .- (2halo,2halo)     # cut off the halo
    @boundscheck (m,n) == size(qhu) || throw(BoundsError())
    @boundscheck (m+2,n+2) == size(dpdy) || throw(BoundsError())

    if S.parameters.seasonal_wind_y
        @unpack ωFy = S.constants
        Fyt = Ftime(T,t,ωFy)
    else
        Fyt = one(T)
    end

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
             dv[i+2,j+2] = -(Tprog(qhu[i,j]) + Tprog(dpdy[i+1,j+1])) + Tprog(Fyt*Fy[i,j])
        end
    end
end

"""Zonal mass flux U = uh."""
function Uflux!(U::AbstractMatrix{T},
                u::AbstractMatrix{T},
                h_u::AbstractMatrix{T},
                ep::Int,
                scale_inv::T) where {T<:AbstractFloat}

    m,n = size(U)
    @boundscheck (m,n) == size(h_u) || throw(BoundsError())
    @boundscheck (m+2+ep,n+2) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            U[i,j] = u[1+ep+i,1+j]*h_u[i,j]*scale_inv
        end
    end
end

"""Meridional mass flux V = vh."""
function Vflux!(V::AbstractMatrix{T},
                v::AbstractMatrix{T},
                h_v::AbstractMatrix{T},
                scale_inv::T) where {T<:AbstractFloat}

    m,n = size(V)
    @boundscheck (m,n) == size(h_v) || throw(BoundsError())
    @boundscheck (m+2,n+2) == size(v) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            V[i,j] = v[i+1,j+1]*h_v[i,j]*scale_inv
        end
    end
end

"""Calculate the mass/volume fluxes U,V from u,v,η."""
function UVfluxes!( u::AbstractMatrix,
                    v::AbstractMatrix,
                    η::AbstractMatrix,
                    Diag::DiagnosticVars,
                    S::ModelSetup)

    @unpack h,h_u,h_v,U,V = Diag.VolumeFluxes
    @unpack H = S.forcing
    @unpack ep = S.grid
    @unpack scale_inv = S.constants

    thickness!(h,η,H)
    Ix!(h_u,h)
    Iy!(h_v,h)

    # mass or volume flux U,V = uh,vh
    Uflux!(U,u,h_u,ep,scale_inv)
    Vflux!(V,v,h_v,scale_inv)
end


text:
using JuliaUtils
using Test

@testset "LinearLS - two datapoints" begin
    x = [1, 2]
    y = [4, 5]
    ls = linearls(x, y)
    @test ls.m == 1
    @test ls.q == 3
end

@testset "Newton - exp" begin
    f(x) = exp(x) - 1
    x0 = 2
    x = newton(f, x0, tol=1e-20, niters=10^3)
    @test isapprox(x, 0, atol=1e-10)
end


text:
# utilities for managing Pumi lbraries

module PumiConfig

using MPI

export CppBool, CONFIG_PATHS, HAVE_SIMMETRIX

# make MeshEntity* passable by MPI
if sizeof(Ptr{Void}) == sizeof(Int32)
  T = Int32
elseif  sizeof(Ptr{Void}) == sizeof(Int64)
  T = Int64
else
  error("cannot find compatible size of Ptr{Void}")
end

if !haskey(MPI.mpitype_dict, Ptr{Void})
  MPI.mpitype_dict[Ptr{Void}] = MPI.mpitype_dict[T]
end

# make PdePumiInterface findable
push!(LOAD_PATH, dirname(@__FILE__))

"""
  Typealias for C++ Bool.  This is implementation defined, so it might need
  to be manually changed for different systems
"""
const CppBool = UInt8  # this is implementation defined

"""
  Holds paths that were determined at configure time
"""
global const CONFIG_PATHS = Dict{String, String}()

"""
  Reads the PUMI_RPATH file and updates [`CONFIG_PATHS`](@ref)
"""
function readConfigFile()

  # this is basically like RPATH
  fname = joinpath(dirname(@__FILE__), "PUMI_RPATH")

  if !isfile(fname)
    error("PUMI_RPATH file not found: please rebuild PumiInterface.jl")
  end

  for line in eachline(fname)
    substr = split(chomp(line))
    @assert length(substr) == 2
    CONFIG_PATHS[substr[1]] = substr[2]
  end

  # String manipulation ins CMake is not easy, so CMake writes whatever
  # paths are convenient and this function updates them to what we really want
  updatePaths()

  return nothing
end

function updatePaths()

  # PUMI_CMAKE_PATH -> PUMI_LIBDIR
  str = CONFIG_PATHS["PUMI_CMAKE_PATH"]
  suffix = "cmake/SCOREC/SCORECConfig.cmake"
  CONFIG_PATHS["PUMI_LIBDIR"] =  replace(str, suffix, "")

  # PUMIINTERFACE_PREFIX -> PUMIINTERFACE_LIBDIR
  str = CONFIG_PATHS["PUMIINTERFACE_PREFIX"]
  CONFIG_PATHS["PUMIINTERFACE_LIBDIR"] = joinpath(str, "lib")

  return nothing
end

readConfigFile()

global const HAVE_SIMMETRIX = CONFIG_PATHS["HAVE_SIMMETRIX"] == "ON"
global const HAVE_SIMMETRIX = CONFIG_PATHS["HAVE_SIMMETRIX"] == "ON"

end  # end module


text:
using .Settings: Config

export plot

mutable struct Plotter
    o::PyObject
end
Plotter(config::Config) = Plotter(qha_plotting.Plotter(config))

@pyinterface Plotter

function plot(cfgfile)
    config = from_yaml(cfgfile)

    plotter = Plotter(config)

    desired_pressures_gpa =
        range(config["P_MIN"]; length = config["NTV"], step = config["DELTA_P"])
    config["DESIRED_PRESSURES_GPa"] = desired_pressures_gpa

    results_folder = expanduser(config["output_directory"])
    calculation_option = Dict(
        "F" => "f_tp",
        "G" => "g_tp",
        "H" => "h_tp",
        "U" => "u_tp",
        "V" => "v_tp",
        "Cv" => "cv_tp_jmolk",
        "Cp" => "cp_tp_jmolk",
        "Bt" => "bt_tp_gpa",
        "Btp" => "btp_tp",
        "Bs" => "bs_tp_gpa",
        "alpha" => "alpha_tp",
        "gamma" => "gamma_tp",
    )

    file_ftv_fitted = joinpath(results_folder, "f_tv_fitted_ev_ang3.txt")
    config["f_tv_fitted"] = file_ftv_fitted

    file_ftv_non_fitted = joinpath(results_folder, "f_tv_nonfitted_ev_ang3.txt")
    config["f_tv_non_fitted"] = file_ftv_non_fitted

    file_ptv_gpa = joinpath(results_folder, "p_tv_gpa.txt")
    config["p_tv_gpa"] = file_ptv_gpa

    plotter.fv_pv()

    for idx in config["thermodynamic_properties"]
        if idx in ("F", "G", "H", "U")
            attr_name = calculation_option[idx] * "_" * config["energy_unit"]
            file_name = attr_name * ".txt"
            file_dir = joinpath(results_folder, file_name)
            config[idx] = file_dir
            plotter.plot2file(idx)
        elseif idx == "V"
            v_ang3 = calculation_option[idx] * "_ang3"
            file_name_ang3 = v_ang3 * ".txt"
            file_dir_ang3 = joinpath(results_folder, file_name_ang3)
            config[idx] = file_dir_ang3
            plotter.plot2file(idx)
        elseif idx in ("Cv", "Cp", "Bt", "Btp", "Bs", "alpha", "gamma")
            attr_name = calculation_option[idx]
            file_name = attr_name * ".txt"
            file_dir = joinpath(results_folder, file_name)
            config[idx] = file_dir
            plotter.plot2file(idx)
        end
    end
end