text:
# This file imports all serialization and deserialization functions
# from internal modules.
import Serialization

"""
    serialize_bin(stream::IO, model)

Serialize an mlpack model `model` in the binary boost::serialization
format to the given `stream`.  Example:

```julia
_, model, _, _ = mlpack.logistic_regression(training=x, labels=y)
mlpack.serialize_bin(open("model.bin", "w"), model)
```

The model can later be loaded with `mlpack.deserialize_bin()`, or even
from mlpack bindings in other languages.  However, the format used
here will *not* work with Julia's `Serialization` package!  Use
`Serialization.serialize()` instead.
"""
function serialize_bin(stream::IO, model) end

"""
    deserialize_bin(stream::IO, model_type::Type)
Deserialize an mlpack model type from an input stream.  Specify the
type of the model manually with the `t` argument.  Example usage:

```julia
lr_model = mlpack.deserialize_bin(stream, LogisticRegression)
```

Only use this if you have saved the model in the boost::serialization
binary format using `serialize_bin()` or an mlpack binding in another
language!  If you used `Serialization.serialize()` to serialize your
model, then use `Serialization.deserialize()` to deserialize it.

Then, the returned model can be passed to appropriate mlpack functions
for machine learning tasks.
"""
function deserialize_bin(stream::IO, t::Type) end

serialize_bin(stream::IO, model::AdaBoostModel) =
    _Internal.adaboost_internal.serializeAdaBoostModel(stream, model)
deserialize_bin(stream::IO, ::Type{AdaBoostModel}) =
    _Internal.adaboost_internal.deserializeAdaBoostModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::AdaBoostModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, AdaBoostModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{AdaBoostModel})
  deserialize_bin(s.io, AdaBoostModel)
end

serialize_bin(stream::IO, model::ApproxKFNModel) =
    _Internal.approx_kfn_internal.serializeApproxKFNModel(stream, model)
deserialize_bin(stream::IO, ::Type{ApproxKFNModel}) =
    _Internal.approx_kfn_internal.deserializeApproxKFNModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::ApproxKFNModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, ApproxKFNModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{ApproxKFNModel})
  deserialize_bin(s.io, ApproxKFNModel)
end

serialize_bin(stream::IO, model::BayesianLinearRegression) =
    _Internal.bayesian_linear_regression_internal.serializeBayesianLinearRegression(stream, model)
deserialize_bin(stream::IO, ::Type{BayesianLinearRegression}) =
    _Internal.bayesian_linear_regression_internal.deserializeBayesianLinearRegression(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::BayesianLinearRegression)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, BayesianLinearRegression)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{BayesianLinearRegression})
  deserialize_bin(s.io, BayesianLinearRegression)
end

serialize_bin(stream::IO, model::CFModel) =
    _Internal.cf_internal.serializeCFModel(stream, model)
deserialize_bin(stream::IO, ::Type{CFModel}) =
    _Internal.cf_internal.deserializeCFModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::CFModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, CFModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{CFModel})
  deserialize_bin(s.io, CFModel)
end

serialize_bin(stream::IO, model::DSModel) =
    _Internal.decision_stump_internal.serializeDSModel(stream, model)
deserialize_bin(stream::IO, ::Type{DSModel}) =
    _Internal.decision_stump_internal.deserializeDSModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::DSModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, DSModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{DSModel})
  deserialize_bin(s.io, DSModel)
end

serialize_bin(stream::IO, model::DecisionTreeModel) =
    _Internal.decision_tree_internal.serializeDecisionTreeModel(stream, model)
deserialize_bin(stream::IO, ::Type{DecisionTreeModel}) =
    _Internal.decision_tree_internal.deserializeDecisionTreeModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::DecisionTreeModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, DecisionTreeModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{DecisionTreeModel})
  deserialize_bin(s.io, DecisionTreeModel)
end

serialize_bin(stream::IO, model::DTree) =
    _Internal.det_internal.serializeDTree(stream, model)
deserialize_bin(stream::IO, ::Type{DTree}) =
    _Internal.det_internal.deserializeDTree(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::DTree)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, DTree)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{DTree})
  deserialize_bin(s.io, DTree)
end

serialize_bin(stream::IO, model::FastMKSModel) =
    _Internal.fastmks_internal.serializeFastMKSModel(stream, model)
deserialize_bin(stream::IO, ::Type{FastMKSModel}) =
    _Internal.fastmks_internal.deserializeFastMKSModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::FastMKSModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, FastMKSModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{FastMKSModel})
  deserialize_bin(s.io, FastMKSModel)
end

serialize_bin(stream::IO, model::GMM) =
    _Internal.gmm_train_internal.serializeGMM(stream, model)
deserialize_bin(stream::IO, ::Type{GMM}) =
    _Internal.gmm_train_internal.deserializeGMM(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::GMM)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, GMM)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{GMM})
  deserialize_bin(s.io, GMM)
end

serialize_bin(stream::IO, model::HMMModel) =
    _Internal.hmm_train_internal.serializeHMMModel(stream, model)
deserialize_bin(stream::IO, ::Type{HMMModel}) =
    _Internal.hmm_train_internal.deserializeHMMModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::HMMModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, HMMModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{HMMModel})
  deserialize_bin(s.io, HMMModel)
end

serialize_bin(stream::IO, model::HoeffdingTreeModel) =
    _Internal.hoeffding_tree_internal.serializeHoeffdingTreeModel(stream, model)
deserialize_bin(stream::IO, ::Type{HoeffdingTreeModel}) =
    _Internal.hoeffding_tree_internal.deserializeHoeffdingTreeModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::HoeffdingTreeModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, HoeffdingTreeModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{HoeffdingTreeModel})
  deserialize_bin(s.io, HoeffdingTreeModel)
end

serialize_bin(stream::IO, model::KDEModel) =
    _Internal.kde_internal.serializeKDEModel(stream, model)
deserialize_bin(stream::IO, ::Type{KDEModel}) =
    _Internal.kde_internal.deserializeKDEModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::KDEModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, KDEModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{KDEModel})
  deserialize_bin(s.io, KDEModel)
end

serialize_bin(stream::IO, model::LARS) =
    _Internal.lars_internal.serializeLARS(stream, model)
deserialize_bin(stream::IO, ::Type{LARS}) =
    _Internal.lars_internal.deserializeLARS(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::LARS)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, LARS)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{LARS})
  deserialize_bin(s.io, LARS)
end

serialize_bin(stream::IO, model::LinearRegression) =
    _Internal.linear_regression_internal.serializeLinearRegression(stream, model)
deserialize_bin(stream::IO, ::Type{LinearRegression}) =
    _Internal.linear_regression_internal.deserializeLinearRegression(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::LinearRegression)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, LinearRegression)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{LinearRegression})
  deserialize_bin(s.io, LinearRegression)
end

serialize_bin(stream::IO, model::LinearSVMModel) =
    _Internal.linear_svm_internal.serializeLinearSVMModel(stream, model)
deserialize_bin(stream::IO, ::Type{LinearSVMModel}) =
    _Internal.linear_svm_internal.deserializeLinearSVMModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::LinearSVMModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, LinearSVMModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{LinearSVMModel})
  deserialize_bin(s.io, LinearSVMModel)
end

serialize_bin(stream::IO, model::LocalCoordinateCoding) =
    _Internal.local_coordinate_coding_internal.serializeLocalCoordinateCoding(stream, model)
deserialize_bin(stream::IO, ::Type{LocalCoordinateCoding}) =
    _Internal.local_coordinate_coding_internal.deserializeLocalCoordinateCoding(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::LocalCoordinateCoding)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, LocalCoordinateCoding)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{LocalCoordinateCoding})
  deserialize_bin(s.io, LocalCoordinateCoding)
end

serialize_bin(stream::IO, model::LogisticRegression) =
    _Internal.logistic_regression_internal.serializeLogisticRegression(stream, model)
deserialize_bin(stream::IO, ::Type{LogisticRegression}) =
    _Internal.logistic_regression_internal.deserializeLogisticRegression(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::LogisticRegression)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, LogisticRegression)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{LogisticRegression})
  deserialize_bin(s.io, LogisticRegression)
end

serialize_bin(stream::IO, model::LSHSearch) =
    _Internal.lsh_internal.serializeLSHSearch(stream, model)
deserialize_bin(stream::IO, ::Type{LSHSearch}) =
    _Internal.lsh_internal.deserializeLSHSearch(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::LSHSearch)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, LSHSearch)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{LSHSearch})
  deserialize_bin(s.io, LSHSearch)
end

serialize_bin(stream::IO, model::NBCModel) =
    _Internal.nbc_internal.serializeNBCModel(stream, model)
deserialize_bin(stream::IO, ::Type{NBCModel}) =
    _Internal.nbc_internal.deserializeNBCModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::NBCModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, NBCModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{NBCModel})
  deserialize_bin(s.io, NBCModel)
end

serialize_bin(stream::IO, model::KNNModel) =
    _Internal.knn_internal.serializeKNNModel(stream, model)
deserialize_bin(stream::IO, ::Type{KNNModel}) =
    _Internal.knn_internal.deserializeKNNModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::KNNModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, KNNModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{KNNModel})
  deserialize_bin(s.io, KNNModel)
end

serialize_bin(stream::IO, model::KFNModel) =
    _Internal.kfn_internal.serializeKFNModel(stream, model)
deserialize_bin(stream::IO, ::Type{KFNModel}) =
    _Internal.kfn_internal.deserializeKFNModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::KFNModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, KFNModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{KFNModel})
  deserialize_bin(s.io, KFNModel)
end

serialize_bin(stream::IO, model::PerceptronModel) =
    _Internal.perceptron_internal.serializePerceptronModel(stream, model)
deserialize_bin(stream::IO, ::Type{PerceptronModel}) =
    _Internal.perceptron_internal.deserializePerceptronModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::PerceptronModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, PerceptronModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{PerceptronModel})
  deserialize_bin(s.io, PerceptronModel)
end

serialize_bin(stream::IO, model::ScalingModel) =
    _Internal.preprocess_scale_internal.serializeScalingModel(stream, model)
deserialize_bin(stream::IO, ::Type{ScalingModel}) =
    _Internal.preprocess_scale_internal.deserializeScalingModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::ScalingModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, ScalingModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{ScalingModel})
  deserialize_bin(s.io, ScalingModel)
end

serialize_bin(stream::IO, model::RandomForestModel) =
    _Internal.random_forest_internal.serializeRandomForestModel(stream, model)
deserialize_bin(stream::IO, ::Type{RandomForestModel}) =
    _Internal.random_forest_internal.deserializeRandomForestModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::RandomForestModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, RandomForestModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{RandomForestModel})
  deserialize_bin(s.io, RandomForestModel)
end

serialize_bin(stream::IO, model::RANNModel) =
    _Internal.krann_internal.serializeRANNModel(stream, model)
deserialize_bin(stream::IO, ::Type{RANNModel}) =
    _Internal.krann_internal.deserializeRANNModel(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::RANNModel)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, RANNModel)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{RANNModel})
  deserialize_bin(s.io, RANNModel)
end

serialize_bin(stream::IO, model::SoftmaxRegression) =
    _Internal.softmax_regression_internal.serializeSoftmaxRegression(stream, model)
deserialize_bin(stream::IO, ::Type{SoftmaxRegression}) =
    _Internal.softmax_regression_internal.deserializeSoftmaxRegression(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::SoftmaxRegression)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, SoftmaxRegression)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{SoftmaxRegression})
  deserialize_bin(s.io, SoftmaxRegression)
end

serialize_bin(stream::IO, model::SparseCoding) =
    _Internal.sparse_coding_internal.serializeSparseCoding(stream, model)
deserialize_bin(stream::IO, ::Type{SparseCoding}) =
    _Internal.sparse_coding_internal.deserializeSparseCoding(stream)

function Serialization.serialize(s::Serialization.AbstractSerializer,
                                 model::SparseCoding)
  Serialization.writetag(s.io, Serialization.OBJECT_TAG)
  Serialization.serialize(s, SparseCoding)
  serialize_bin(s.io, model)
end

function Serialization.deserialize(s::Serialization.AbstractSerializer,
                                   ::Type{SparseCoding})
  deserialize_bin(s.io, SparseCoding)
end

ext: jl


text:
module FastLegendre

# From
# O(1) Computation of Legendre polynomials and Gauss-Legendre nodes and weights for parallel computing
# Bogaert, Ignace and Michiels, Bart and Fostier, Jan
# SIAM JOURNAL ON SCIENTIFIC COMPUTING
# 34, 2012
# http://dx.doi.org/10.1137/110855442
#
# as mentioned https://github.com/JuliaMath/SpecialFunctions.jl/issues/124
#
# Only Float64
#
# Gives a speed up over Clenshaw, which is expected, as this is O(1), not O(n):
#
# n= 1000 Clenshaw might be faster
# n= 10_000 10x faster
# n= 100_000 10x faster
# n= 1_000_000 500x faster
# n= 10_000_000 10,000x faster
# (measurement of basis(P,n)(x) compared to Basis(P,n)(x)

using FastGaussQuadrature
using SpecialFunctions

# 4 regions of evaluation
# specialized for Float64 values
function fastlegendre(n, x)
    n == 0 && return one(x)
    n == 1 && return x
    n < 100 && return direct_32(n, x)
    x < 0 && return (iseven(n) ? 1 : -1) * fastlegendre(n, -x)
    θ = acos(x)
    (n + 1) * sin(θ) > 25 && return asy_32(n, x)
    asy_33(n, x)
end

## ---- n < 100 case

# precompute first 99
const xᵢs = [FastGaussQuadrature.gausslegendre(n)[1][1:(n ÷ 2)] for n in 1:99]
const xᵢs⁻¹ = [inv.(xᵢ) for xᵢ in xᵢs]
const xᵢs2⁻¹ = [inv.(1 .- xᵢ .^ 2) for xᵢ in xᵢs]

# precompute
const Cn = zeros(Float64, 99)
Cn[1] = 1.0
for i in 1:49
    n = 2i
    Cn[n] = sqrt(pi) / gamma(i + 1) / gamma((1 - n) / 2)
    Cn[n + 1] = (n + 1) * Cn[n]
end

# use cosine or sine
function direct_32(n, x)
    if abs(x) > 0.707
        if iseven(n)
            direct_sine(Val(true), n, x)
        else
            direct_sine(Val(false), n, x)
        end
    else
        if iseven(n)
            direct_cosine(Val(true), n, x)
        else
            direct_cosine(Val(false), n, x)
        end
    end
end

#x = cos(theta) 1 - x^2 = sintheta0
function direct_sine(::Val{true}, n, x)
    tot = one(Float64)
    xs = xᵢs2⁻¹[n]
    for aᵢ in xs
        tot *= 1 - (1 - x^2) * aᵢ
    end
    tot
end

direct_sine(::Val{false}, n, x) = x * direct_sine(Val(true), n, x)

function direct_cosine(p::Val{true}, n, x)
    tot = Cn[n]
    xs = xᵢs⁻¹[n]
    for aᵢ in xs
        tot *= (1 - (x * aᵢ)^2)
    end
    tot
end
direct_cosine(p::Val{false}, n, x) = x * direct_cosine(Val(true), n, x)

## ----
const Cl0s = [gamma(n + 1) / gamma(n + 3 / 2) for n in 1:9]
τ(x) = evalpoly(
    1 / x,
    (
        1 / 1,
        0,
        -1 / 64,
        0,
        21 / 8192,
        0,
        -671 / 524288,
        0,
        180323 / 80323134,
        0,
        -20898423 / 8589934592,
        -7426362705 / 1099511627776,
    ),
)

function Cl0(l)
    1 <= l < 10 && return Cl0s[l]
    1 / sqrt(l + 3 / 4) * τ(l + 3 / 4)
end

function asy_32(n, x)
    M = 17
    θ = acos(x)
    val = sqrt((2 / pi) * (1 / sin(θ)))
    tot = 0.0
    sθ = sin(θ)
    λ = Cl0(n)
    for m in 0:(M - 1)
        α = (n + m + 1 / 2) * θ - (m + 1 / 2) * pi / 2
        tot += λ * cos(α)
        λ *= (m + 1 / 2)^2 / 2 / (m + 1) / (n + m + 3 / 2) / sθ
    end
    val * tot
end

## ---
h(n, y) = y^n * besselj(n, y)

const cns = [
    [1 / 8, -1 / 12],
    [11 / 384, -7 / 160, 1 / 160],
    [173 / 15360, -101 / 3584, 671 / 80640, -61 / 120960],
    [22931 / 3440640, -90497 / 3870720, 217 / 20480, -1261 / 967680, 1261 / 29030400],
    [
        1319183 / 247726080,
        -10918993 / 454164480,
        +1676287 / 113541120,
        -7034857 / 2554675200,
        +1501 / 8110080,
        -79 / 20275200,
    ],
    [
        233526463 / 43599790080,
        -1396004969 / 47233105920,
        +2323237523 / 101213798400,
        -72836747 / 12651724800,
        +3135577 / 5367398400,
        -1532789 / 61993451520,
        +66643 / 185980354560,
    ],
]

function fn(i, y)
    n = 2i
    cs = cns[i]
    tot = 0.0
    Δ = i - 1
    for j in i:(2i)
        tot += cs[j - Δ] * h(j, y)
    end
    tot
end

function asy_33(n, x)
    v = n + 1 / 2
    θ = acos(x)
    y = v * θ

    # Eq (3.19)
    val = besselj(0, y)

    λ = 1 / v^2
    for i in 1:6
        val += fn(i, y) * λ
        λ /= v^2
    end
    val
end

end


text:
function set_knapsacktest_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        termination_status = MOI.OPTIMAL,
        primal_status = MOI.FEASIBLE_POINT,
        dual_status = MOI.NO_SOLUTION,
        variable_primal = [1.0, 0.0, 0.0, 1.0, 1.0],
        constraint_primal = [9.0]
    )
    return
end

function set_int1test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        termination_status = MOI.OPTIMAL,
        primal_status = MOI.FEASIBLE_POINT,
        dual_status = MOI.NO_SOLUTION,
        variable_primal = [4.0, 5.0, 1.0],
        constraint_primal = [10.0, 15.0]
    )
    return
end

function set_int3test_solutions!(solver)
    # SOLVE 1
    LQOI.set_solution!(solver,
        termination_status = MOI.OPTIMAL,
        primal_status = MOI.FEASIBLE_POINT,
        dual_status = MOI.NO_SOLUTION,
        variable_primal = [1.0, -0.0, -0.0, -0.0, -0.0, 1.0, -0.0, -0.0, -0.0, -0.0, -0.0],
        constraint_primal = [0.9875]
    )
    # SOLVE 2
    LQOI.set_solution!(solver,
        termination_status = MOI.OPTIMAL,
        primal_status = MOI.FEASIBLE_POINT,
        dual_status = MOI.NO_SOLUTION,
        variable_primal = [1.0, -0.0, -0.0, -0.0, -0.0, 1.0, -0.0, -0.0, -0.0, -0.0, -0.0],
        constraint_primal = [0.9875]
    )
    return
end


text:
module Nystrom

using StaticArrays
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using SpecialFunctions
using QuadGK
using RecipesBase
using TimerOutputs
using Base.Threads

# packages to generate the geometry
using GmshSDK
using ParametricSurfaces

# base utilities
using WavePropBase
using WavePropBase.Utils
using WavePropBase.Geometry
using WavePropBase.Interpolation
using WavePropBase.Integration
using WavePropBase.Mesh
using WavePropBase.Trees
WavePropBase.@import_interface

# Interpolated Factored Green Function Method
# import IFGF
# import IFGF: IFGFOperator, centered_factor

export
    # Gmsh related stuff
    GmshSDK,
    @gmsh,
    gmsh,
    # re-export useful stuff
    Geometry,
    ParametricSurfaces,
    Laplace,
    Helmholtz,
    Elastostatic,
    Maxwell,
    # types
    NystromMesh,
    NystromDOF,
    SingleLayerKernel,
    DoubleLayerKernel,
    AdjointDoubleLayerKernel,
    HyperSingularKernel,
    SingleLayerPotential,
    DoubleLayerPotential,
    SingleLayerOperator,
    DoubleLayerOperator,
    AdjointDoubleLayerOperator,
    HyperSingularOperator,
    Density,
    TangentialDensity,
    # functions
    qweights,
    trace,
    ncross,
    coords,
    qcoords,
    γ₀,
    γ₁
    # IFGFCompressor

WavePropBase.@export_interface

include("utils.jl")
include("nystrommesh.jl")
include("kernels.jl")
include("density.jl")
include("potential.jl")
include("integraloperator.jl")
include("lebedevpoints.jl")
include("dim.jl")
include("gausskronrod.jl")
include("maxwellwrappers.jl")
include("discreteoperator.jl")
# include("ifgf.jl")

end # module

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
#=

* This file contains derivative work from the Geophysical Fluid Dynamics
* Laboratory GitHub project for the FV3 dynamical core
* (https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere). It is a stand alone
* subset of the original specifically licensed for use under the
* Apache License (https://www.apache.org/licenses/LICENSE-2.0).
* The application of the Apache License to this standalone work does not imply
* nor shall it be construed as any change to the original source license.

=#

module SwCoreModule

include("./Interpolate.jl")

using Printf, NCDatasets, Revise, OffsetArrays, Dates, Parameters, DataStructures

using .Interpolate

export print_state, write_state, c_sw!, interpolate_state!

#------------------------------------------------------------------
# print_state
#
# Prints statistics for the kernel state variables
#------------------------------------------------------------------
function print_state(message::String, data, io::IOStream)

    println(io, "TEST")
    println(io, "TEST ", "="^115)
    @printf(io, "%-5s %29s\n", "TEST ", message)
    println(io, "TEST ", "="^115)
    @printf(io, "%-5s %14s %19s %19s %19s %19s %19s\n","TEST", "Variable", "Min", "Max", "First", "Last", "RMS")
    println(io, "TEST ", "-"^115)

    # print the variables
    print_2d_variable(io, "rarea", data.rarea)
    print_2d_variable(io, "rarea_c", data.rarea_c)
    print_3d_variable(io, "sin_sg", data.sin_sg)
    print_3d_variable(io, "cos_sg", data.cos_sg)
    print_2d_variable(io, "sina_v", data.sina_v)
    print_2d_variable(io, "cosa_v", data.cosa_v)
    print_2d_variable(io, "sina_u", data.sina_u)
    print_2d_variable(io, "cosa_u", data.cosa_u)
    print_2d_variable(io, "fC", data.fC)
    print_2d_variable(io, "rdxc", data.rdxc)
    print_2d_variable(io, "rdyc", data.rdyc)
    print_2d_variable(io, "dx", data.dx)
    print_2d_variable(io, "dy", data.dy)
    print_2d_variable(io, "dxc", data.dxc)
    print_2d_variable(io, "dyc", data.dyc)
    print_2d_variable(io, "cosa_s", data.cosa_s)
    print_2d_variable(io, "rsin_u", data.rsin_u)
    print_2d_variable(io, "rsin_v", data.rsin_v)
    print_2d_variable(io, "rsin2", data.rsin2)
    print_2d_variable(io, "dxa", data.dxa)
    print_2d_variable(io, "dya", data.dya)

    print_3d_variable(io,"delp", data.delp)
    print_3d_variable(io,"delpc", data.delpc)
    print_3d_variable(io,"pt", data.pt)
    print_3d_variable(io,"ptc", data.ptc)
    print_3d_variable(io,"u", data.u)
    print_3d_variable(io,"v", data.v)
    print_3d_variable(io,"w", data.w)
    print_3d_variable(io,"uc", data.uc)
    print_3d_variable(io,"vc", data.vc)
    print_3d_variable(io,"ua", data.ua)
    print_3d_variable(io,"va", data.va)
    print_3d_variable(io,"wc", data.wc)
    print_3d_variable(io,"ut", data.ut)
    print_3d_variable(io,"vt", data.vt)
    print_3d_variable(io,"divg_d", data.divg_d)

    println(io, "TEST ", "-"^115)
    println(io, "TEST")

    return nothing

end # function print_state

#------------------------------------------------------------------
# print_2d_variable
#
# Prints statistics for a 2d state variable
#------------------------------------------------------------------
function print_2d_variable(io::IOStream,name::String, data::OffsetArray)

    min =     minimum(data)
    max =     maximum(data)
    first =   data[begin, begin]
    last =    data[end, end]
    rms =     sqrt(sum(data.^2.) / length(data))

    @printf(io, "%-5s %14s %19.10E %19.10E %19.10E %19.10E %19.10E\n","TEST", name, min, max , first, last, rms)

    return nothing

end # function print_2d_variable

#------------------------------------------------------------------
# print_3d_variable
#
# Prints statistics for a 3d state variable
#------------------------------------------------------------------
function print_3d_variable(io::IOStream, name::String, data::OffsetArray)

    min =    minimum(data)
    max =    maximum(data)
    first =  data[begin, begin, begin]
    last =   data[end, end, end]
    rms =    sqrt(sum(data.^2.) / length(data))

    @printf(io, "%-5s %14s %19.10E %19.10E %19.10E %19.10E %19.10E\n","TEST", name, min, max , first, last, rms)

    return nothing

end # function print_3d_variable

#------------------------------------------------------------------
# write_state
#
# Write state to NetCDF file
#------------------------------------------------------------------
# create a NetCDF file from the data in netCDF-4 format
function write_state(outputfilename::String, data)

  date = Dates.now()
  dateString = Dates.format(date, "yyyy-mm-dd HH:MM:SS")

  ds = NCDataset(outputfilename, "c" ; format=:netcdf4, attrib=OrderedDict(
    "creation_date"             => dateString,
    "kernel_name"               => "c_sw",
    "isd"                       => data.isd,
    "ied"                       => data.ied,
    "jsd"                       => data.jsd,
    "jed"                       => data.jed,
    "is"                        => data.is,
    "ie"                        => data.ie,
    "js"                        => data.js,
    "je"                        => data.je,
    "nord"                      => data.nord,
    "npx"                       => data.npx,
    "npy"                       => data.npy,
    "npz"                       => data.npz,
    "dt2"                       => data.dt2,
    "sw_corner"                 => data.sw_corner,
    "se_corner"                 => data.se_corner,
    "nw_corner"                 => data.nw_corner,
    "ne_corner"                 => data.ne_corner,
  ))

  # define Dimensions
  defDim(ds, "ni", data.ni)
  defDim(ds, "nj", data.nj)
  defDim(ds, "nip1", data.nip1)
  defDim(ds, "njp1", data.njp1)
  defDim(ds, "nk", data.nk)
  defDim(ds, "nk9", data.nk9)

  # Declare variables
  ncrarea =           defVar(ds, "rarea", Float64, ("ni", "nj"))
  ncrarea_c =         defVar(ds, "rarea_c", Float64, ("nip1", "njp1"))
  ncsin_sg =          defVar(ds, "sin_sg", Float64, ("ni", "nj", "nk9"))
  nccos_sg =          defVar(ds, "cos_sg", Float64, ("ni", "nj", "nk9"))
  ncsina_v =          defVar(ds, "sina_v", Float64, ("ni", "njp1"))
  nccosa_v =          defVar(ds, "cosa_v", Float64, ("ni", "njp1"))
  ncsina_u =          defVar(ds, "sina_u", Float64, ("nip1", "nj"))
  nccosa_u =          defVar(ds, "cosa_u", Float64, ("nip1", "nj"))
  ncfC =              defVar(ds, "fC", Float64, ("nip1", "njp1"))
  ncrdxc =            defVar(ds, "rdxc", Float64, ("nip1", "nj"))
  ncrdyc =            defVar(ds, "rdyc", Float64, ("ni", "njp1"))
  ncdx =              defVar(ds, "dx", Float64, ("ni", "njp1"))
  ncdy =              defVar(ds, "dy", Float64, ("nip1", "nj"))
  ncdxc =             defVar(ds, "dxc", Float64, ("nip1", "nj"))
  ncdyc =             defVar(ds, "dyc", Float64, ("ni", "njp1"))
  nccosa_s =          defVar(ds, "cosa_s", Float64, ("ni", "nj"))
  ncrsin_u =          defVar(ds, "rsin_u", Float64, ("nip1", "nj"))
  ncrsin_v =          defVar(ds, "rsin_v", Float64, ("ni", "njp1"))
  ncrsin2 =           defVar(ds, "rsin2", Float64, ("ni", "nj"))
  ncdxa =             defVar(ds, "dxa", Float64, ("ni", "nj"))
  ncdya =             defVar(ds, "dya", Float64, ("ni", "nj"))
  ncdelp =            defVar(ds, "delp", Float64, ("ni", "nj", "nk"))
  ncdelpc =           defVar(ds, "delpc", Float64, ("ni", "nj", "nk"))
  ncpt =              defVar(ds, "pt", Float64, ("ni", "nj", "nk"))
  ncptc =             defVar(ds, "ptc", Float64, ("ni", "nj", "nk"))
  ncu =               defVar(ds, "u", Float64, ("ni", "njp1", "nk"))
  ncv =               defVar(ds, "v", Float64, ("nip1", "nj", "nk"))
  ncw =               defVar(ds, "w", Float64, ("ni", "nj", "nk"))
  ncuc =              defVar(ds, "uc", Float64, ("nip1", "nj", "nk"))
  ncvc =              defVar(ds, "vc", Float64, ("ni", "njp1", "nk"))
  ncua =              defVar(ds, "ua", Float64, ("ni", "nj", "nk"))
  ncva =              defVar(ds, "va", Float64, ("ni", "nj", "nk"))
  ncwc =              defVar(ds, "wc", Float64, ("ni", "nj", "nk"))
  ncut =              defVar(ds, "ut", Float64, ("ni", "nj", "nk"))
  ncvt =              defVar(ds, "vt", Float64, ("ni", "nj", "nk"))
  ncdivg_d =          defVar(ds, "divg_d", Float64, ("nip1", "njp1", "nk"))

  # Define variables
  ncrarea[:,:] =    data.rarea
  ncrarea_c[:,:] =  data.rarea_c
  ncsin_sg[:,:,:] = data.sin_sg
  nccos_sg[:,:,:] = data.cos_sg
  ncsina_v[:,:] =   data.sina_v
  nccosa_v[:,:] =   data.cosa_v
  ncsina_u[:,:] =   data.sina_u
  nccosa_u[:,:] =   data.cosa_u
  ncfC[:,:] =       data.fC
  ncrdxc[:,:] =     data.rdxc
  ncrdyc[:,:] =     data.rdyc
  ncdx[:,:] =       data.dx
  ncdy[:,:] =       data.dy
  ncdxc[:,:] =      data.dxc
  ncdyc[:,:] =      data.dyc
  nccosa_s[:,:] =   data.cosa_s
  ncrsin_u[:,:] =   data.rsin_u
  ncrsin_v[:,:] =   data.rsin_v
  ncrsin2[:,:] =    data.rsin2
  ncdxa[:,:] =      data.dxa
  ncdya[:,:] =      data.dya
  ncdelp[:,:,:] =   data.delp
  ncdelpc[:,:,:] =  data.delpc
  ncpt[:,:,:] =     data.pt
  ncptc[:,:,:] =    data.ptc
  ncu[:,:,:] =      data.u
  ncv[:,:,:] =      data.v
  ncw[:,:,:] =      data.w
  ncuc[:,:,:] =     data.uc
  ncvc[:,:,:] =     data.vc
  ncua[:,:,:] =     data.ua
  ncva[:,:,:] =     data.va
  ncwc[:,:,:] =     data.wc
  ncut[:,:,:] =     data.ut
  ncvt[:,:,:] =     data.vt
  ncdivg_d[:,:,:] = data.divg_d

  close(ds)

end # function write_state

#------------------------------------------------------------------
# c_sw
#
# Top-level routine for c_sw kernel
#------------------------------------------------------------------
function c_sw!(data, k::Int64)

  # Unpack arguments
  @unpack sw_corner, se_corner, nw_corner, ne_corner,
  rarea, rarea_c, sin_sg, cos_sg, sina_v, cosa_v,
  sina_u, cosa_u, fC, rdxc, rdyc, dx, dy, dxc, dyc,
  cosa_s, rsin_u, rsin_v, rsin2, dxa, dya,
  delpc, delp, ptc, pt, u, v, w, uc, vc, ua, va, wc, ut,
  vt, divg_d, dt2  = data

  # Unpack local variables
  @unpack isd, ied, jsd, jed, is, ie, js, je,
          nord, npx, npy, npz, dt2 = data

  # re-dimension variables to the first two dimensions
  delpc =        @view delpc[:,:,k]
  delp =         @view delp[:,:,k]
  ptc =          @view ptc[:,:,k]
  pt =           @view pt[:,:,k]
  u =            @view u[:,:,k]
  v =            @view v[:,:,k]
  w =            @view w[:,:,k]
  uc =           @view uc[:,:,k]
  vc =           @view vc[:,:,k]
  ua =           @view ua[:,:,k]
  va =           @view va[:,:,k]
  wc =           @view wc[:,:,k]
  ut =           @view ut[:,:,k]
  vt =           @view vt[:,:,k]
  divg_d =       @view divg_d[:,:,k]

  #helper function to index dimensions
  indexdim(a, b) = b - a + 1

  # define local variables arrays and index them
  vort = OffsetArray{Float64, 2}(undef, is-1:ie+1, js-1:je+1)
  ke   = OffsetArray{Float64, 2}(undef, is-1:ie+1, js-1:je+1)
  fx   = OffsetArray{Float64, 2}(undef, is-1:ie+2, js-1:je+1)
  fx1  = OffsetArray{Float64, 2}(undef, is-1:ie+2, js-1:je+1)
  fx2  = OffsetArray{Float64, 2}(undef, is-1:ie+2, js-1:je+1)
  fy   = OffsetArray{Float64, 2}(undef, is-1:ie+1, js-1:je+2)
  fy1  = OffsetArray{Float64, 2}(undef, is-1:ie+1, js-1:je+2)
  fy2  = OffsetArray{Float64, 2}(undef, is-1:ie+1, js-1:je+2)

  # Convert corners from Int32 to Bool for Julia conditional statements
  sw_corner = !iszero(sw_corner)
  se_corner = !iszero(se_corner)
  ne_corner = !iszero(ne_corner)
  nw_corner = !iszero(nw_corner)

  iep1 = ie + 1
  jep1 = je + 1

  d2a2c_vect!(data, k)

  if nord > 0
      divergence_corner!(data, k)
  end

  for j = js-1:jep1
      for i = is-1:iep1+1
          if ut[i, j] > 0.0
              ut[i, j] = dt2 * ut[i, j] * dy[i, j] * sin_sg[i-1, j, 3]
          else
              ut[i, j] = dt2 * ut[i, j] * dy[i, j] * sin_sg[i, j, 1]
          end
      end
  end
  for j = js-1:je+2
      for i = is-1:iep1
          if vt[i, j] > 0.0
              vt[i, j] = dt2 * vt[i, j] * dx[i, j] * sin_sg[i, j-1, 4]
          else
              vt[i, j] = dt2 * vt[i, j] * dx[i, j] * sin_sg[i, j, 2]
          end
      end
  end

  # -----------------
  #  Transport delp:
  # -----------------
  # Xdir:
  fill2_4corners!(delp, pt, 1, sw_corner, se_corner, ne_corner, nw_corner, npx, npy)
  fill_4corners!(w, 1, sw_corner, se_corner, ne_corner, nw_corner, npx, npy)

  for j = js-1:je+1
      for i = is-1:ie+2
          if ut[i, j] > 0.
              fx1[i, j] = delp[i-1, j]
              fx[i, j]  = pt[i-1, j]
              fx2[i, j] = w[i-1, j]
          else
              fx1[i, j] = delp[i, j]
              fx[i, j]  = pt[i, j]
              fx2[i, j] = w[i, j]
          end
          fx1[i, j] = ut[i, j]  * fx1[i, j]
          fx[i, j]  = fx1[i, j] * fx[i, j]
          fx2[i, j] = fx1[i, j] * fx2[i, j]
      end
  end

  # Ydir:
  fill2_4corners!(delp, pt, 2, sw_corner, se_corner, ne_corner, nw_corner, npx, npy)
  fill_4corners!(w, 2, sw_corner, se_corner, ne_corner, nw_corner, npx, npy)

  for j = js-1:je+2
      for i = is-1:ie+1
          if vt[i, j] > 0.
              fy1[i, j] = delp[i, j-1]
              fy[i, j]  = pt[i, j-1]
              fy2[i, j] = w[i, j-1]
          else
              fy1[i, j] = delp[i, j]
              fy[i, j]  = pt[i, j]
              fy2[i, j] = w[i, j]
          end
          fy1[i, j] = vt[i, j]  * fy1[i, j]
          fy[i, j]  = fy1[i, j] * fy[i, j]
          fy2[i, j] = fy1[i, j] * fy2[i, j]
      end
  end
  for j = js-1:je+1
      for i = is-1:ie+1
          delpc[i, j] = delp[i, j] + (fx1[i, j] - fx1[i+1, j] + fy1[i, j] - fy1[i, j+1]) *
              rarea[i, j]
          ptc[i, j] = (pt[i, j] * delp[i, j] + (fx[i, j] - fx[i+1, j] +
                       fy[i, j] - fy[i, j+1]) * rarea[i, j]) / delpc[i, j]
          wc[i, j] = (w[i, j] * delp[i, j] + (fx2[i, j] - fx2[i+1, j] +
                       fy2[i, j] - fy2[i, j+1]) * rarea[i, j]) / delpc[i, j]
      end
  end

  # ----------------------------------------------------------------------------
  #  Compute KE:
  #  Since uc = u*, i.e. the covariant wind perpendicular to the face edge, if
  #  we want to compute kinetic energy we will need the true coordinate-parallel
  #  covariant wind, computed through u = uc*sina + v*cosa.
  #  Use the alpha for the cell KE is being computed in.
  # ----------------------------------------------------------------------------
  for j = js-1:jep1
      for i = is-1:iep1
          if ua[i, j] > 0.
              if i == 1
                  ke[1, j] = uc[1, j] * sin_sg[1, j, 1] + v[1, j] * cos_sg[1, j, 1]
              elseif i == npx
                  ke[i, j] = uc[npx, j] * sin_sg[npx, j, 1] + v[npx, j] * cos_sg[npx, j, 1]
              else
                  ke[i, j] = uc[i, j]
              end
          else
              if i == 0
                  ke[0, j] = uc[1, j] * sin_sg[0, j, 3] + v[1, j] * cos_sg[0, j, 3]
              elseif i == npx - 1
                  ke[i, j] = uc[npx, j] * sin_sg[npx-1, j, 3] + v[npx, j] * cos_sg[npx-1, j, 3]
              else
                  ke[i, j] = uc[i+1, j]
              end
          end
      end
  end
  for j = js-1:jep1
      for i = is-1:iep1
          if va[i, j] > 0.
              if j == 1
                  vort[i, 1] = vc[i, 1] * sin_sg[i, 1, 2] + u[i, 1] * cos_sg[i, 1, 2]
              elseif j == npy
                  vort[i, j] = vc[i, npy] * sin_sg[i, npy, 2] + u[i, npy] * cos_sg[i, npy, 2]
              else
                  vort[i, j] = vc[i, j]
              end
          else
              if j == 0
                  vort[i, 0] = vc[i, 1] * sin_sg[i, 0, 4] + u[i, 1] * cos_sg[i, 0, 4]
              elseif j == npy - 1
                  vort[i, j] = vc[i, npy] * sin_sg[i, npy-1, 4] + u[i, npy] * cos_sg[i, npy-1, 4]
              else
                  vort[i, j] = vc[i, j+1]
              end
          end
      end
  end

  dt4 = 0.5 * dt2
  for j = js-1:jep1
      for i = is-1:iep1
          ke[i, j] = dt4 * (ua[i, j] * ke[i, j] + va[i, j] * vort[i, j])
      end
  end

  # ------------------------------------------------------------------
  # Compute circulation on C grid
  # To consider using true co - variant winds at face edges?
  # ------------------------------------------------------------------
  for j = js-1:je+1
      for i = is:ie+1
          fx[i, j] = uc[i, j] * dxc[i, j]
      end
  end

  for j = js:je+1
      for i = is-1:ie+1
          fy[i, j] = vc[i, j] * dyc[i, j]
      end
  end

  for j = js:je+1
      for i = is:ie+1
          vort[i, j] = fx[i, j-1] - fx[i, j] - fy[i-1, j] + fy[i, j]
      end
  end

  # Remove the extra term at the corners:
  if sw_corner
      vort[1, 1] = vort[1, 1] + fy[0, 1]
  end
  if se_corner
      vort[npx, 1] = vort[npx, 1] - fy[npx, 1]
  end
  if ne_corner
      vort[npx, npy] = vort[npx, npy] - fy[npx, npy]
  end
  if nw_corner
      vort[1, npy] = vort[1, npy] + fy[0, npy]
  end

  # ------------------------------------------------------------------
  # Compute absolute vorticity
  # ------------------------------------------------------------------
  for j = js:je+1
      for i = is:ie+1
          vort[i, j] = fC[i, j] + rarea_c[i, j] * vort[i, j]
      end
  end

  # ------------------------------------------------------------------------------
  # Transport absolute vorticity:

  # To go from v to contravariant v at the edges, we divide by sin_sg;
  # but we  must multiply by sin_sg to get the proper flux.
  # These cancel, leaving us with fy1 = dt2 * v at the edges.
  # [For the same reason we only divide by sin instead of sin^2 in the interior]
  # ------------------------------------------------------------------------------
  for j = js:je
      #DEC$ VECTOR ALWAYS
      for i = is:iep1
          if i == 1 || i == npx
              fy1[i, j] = dt2 * v[i, j]
          else
              fy1[i, j] = dt2 * (v[i, j] - uc[i, j] * cosa_u[i, j]) / sina_u[i, j]
          end
          if fy1[i, j] > 0.0
              fy[i, j] = vort[i, j]
          else
              fy[i, j] = vort[i, j+1]
          end
      end
  end
  for j = js:jep1
      if j == 1 || j == npy
          #DEC$ VECTOR ALWAYS
          for i = is:ie
              fx1[i, j] = dt2 * u[i, j]
              if fx1[i, j] > 0.0
                  fx[i, j] = vort[i, j]
              else
                  fx[i, j] = vort[i+1, j]
              end
          end
      else
          #DEC$ VECTOR ALWAYS
          for i = is:ie
              fx1[i, j] = dt2 * (u[i, j] - vc[i, j] * cosa_v[i, j]) / sina_v[i, j]
              if fx1[i, j] > 0.0
                  fx[i, j] = vort[i, j]
              else
                  fx[i, j] = vort[i+1, j]
              end
          end
      end
  end

  # Update time - centered winds on the C - Grid
  for j = js:je
      for i = is:iep1
          uc[i, j] = uc[i, j] + fy1[i, j] * fy[i, j] + rdxc[i, j] * (ke[i-1, j] - ke[i, j])
      end
  end
  for j = js:jep1
      for i = is:ie
          vc[i, j] = vc[i, j] - fx1[i, j] * fx[i, j] + rdyc[i, j] * (ke[i, j-1] - ke[i, j])
      end
  end

end # function c_sw

# ------------------------------------------------------------------
#  divergence_corner
# ------------------------------------------------------------------
function divergence_corner!(data, k::Int64)

  # unpack arguments
  @unpack sw_corner, se_corner, ne_corner, nw_corner,
  rarea_c, sin_sg, cos_sg, dxc, dyc, u, v, ua, va, divg_d,
  is, ie, js, je, npx, npy = data

  # re-dimension variables to the first two dimensions
  u =         @view u[:,:,k]
  v =         @view v[:,:,k]
  ua =        @view ua[:,:,k]
  va =        @view va[:,:,k]
  divg_d =    @view divg_d[:,:,k]

  # Convert corners from Int32 to Bool for Julia conditional statements
  sw_corner = !iszero(sw_corner)
  se_corner = !iszero(se_corner)
  ne_corner = !iszero(ne_corner)
  nw_corner = !iszero(nw_corner)

  indexdim(a, b) = b - a + 1

  # local variables
  uf = OffsetArray(Array{Float64, 2}(undef, indexdim(is-2,ie+2) , indexdim(js-1, je+2)), is-2:ie+2, js-1:je+2)
  vf = OffsetArray(Array{Float64, 2}(undef, indexdim(is-1,ie+2) , indexdim(js-2, je+2)), is-1:ie+2, js-2:je+2)

  is2 = max(2, is)
  ie1 = min(npx - 1, ie + 1)

  #    9---4---8
  #    |       |
  #    1   5   3
  #    |       |
  #    6---2---7
  for j = js:je+1
    if j == 1 || j == npy
      for i = is-1:ie+1
        uf[i, j] = u[i, j] * dyc[i, j] * 0.5 * (sin_sg[i, j-1, 4] + sin_sg[i, j, 2])
      end
    else
      for i = is-1:ie+1
        uf[i, j] = (u[i, j] - 0.25 * (va[i, j-1] + va[i, j]) *
        (cos_sg[i, j-1, 4] + cos_sg[i, j, 2])) *
        dyc[i, j] * 0.5 * (sin_sg[i, j-1, 4] + sin_sg[i, j, 2])
      end
    end
  end

  for j = js-1:je+1
    for i = is2:ie1
      vf[i, j] = (v[i, j] - 0.25 * (ua[i-1, j] + ua[i, j]) *
      (cos_sg[i-1, j, 3] + cos_sg[i, j, 1])) *
      dxc[i, j] * 0.5 * (sin_sg[i-1, j, 3] + sin_sg[i, j, 1])
    end
    if is == 1
      vf[1, j] = v[1, j] * dxc[1, j] * 0.5 * (sin_sg[0, j, 3] + sin_sg[1, j, 1])
    end
    if ie + 1 == npx
      vf[npx, j] = v[npx, j] * dxc[npx, j] * 0.5 * (sin_sg[npx-1, j, 3] + sin_sg[npx, j, 1])
    end
  end

  for j = js:je+1
    for i = is:ie+1
      divg_d[i, j] = vf[i, j-1] - vf[i, j] + uf[i-1, j] - uf[i, j]
    end
  end

  # Remove the extra term at the corners:
  if sw_corner
    divg_d[1, 1] = divg_d[1, 1] - vf[1, 0]
  end
  if se_corner
    divg_d[npx, 1] = divg_d[npx, 1] - vf[npx, 0]
  end
  if ne_corner
    divg_d[npx, npy] = divg_d[npx, npy] + vf[npx, npy]
  end
  if nw_corner
    divg_d[1, npy] = divg_d[1, npy] + vf[1, npy]
  end

  for j = js:je+1
    for i = is:ie+1
      divg_d[i, j] = rarea_c[i, j] * divg_d[i, j]
    end
  end
  return nothing
end # function divergence_corner

# ------------------------------------------------------------------
# d2a2c_vect
#
#  There is a limit to how far this routine can fill uc and vc in the
#  halo, and so either mpp_update_domains or some sort of boundary
#  routine [extrapolation, outflow, interpolation from a nested grid]
#  is needed after c_sw is completed if these variables are needed
#  in the halo
# -------------------------------------------------------------------
function d2a2c_vect!(data, k::Int64)

  @unpack sw_corner, se_corner, ne_corner, nw_corner,
  sin_sg, cosa_u, cosa_v, cosa_s, rsin_u, rsin_v,
  rsin2, dxa, dya, u, v, ua, va, uc, vc, ut, vt,
  is, js, ie, je, isd, ied, jsd, jed, npx, npy = data

  # "slice the 3d arrays to 2d arrays"
  u =       @view u[:,:,k]
  v =       @view v[:,:,k]
  uc =      @view uc[:,:,k]
  vc =      @view vc[:,:,k]
  ua =      @view ua[:,:,k]
  va =      @view va[:,:,k]
  ut =      @view ut[:,:,k]
  vt =      @view vt[:,:,k]

  # Convert corners from Int32 to Bool for Julia conditional statements
  sw_corner = !iszero(sw_corner)
  se_corner = !iszero(se_corner)
  ne_corner = !iszero(ne_corner)
  nw_corner = !iszero(nw_corner)

  # 4-pt Lagrange interpolation
  a1 =  0.5625
  a2 = -0.0625
  # volume-conserving cubic with 2nd drv=0 at end point:
  c1 = -2 / 14
  c2 = 11 / 14
  c3 =  5 / 14

  indexdim(a, b) = b - a + 1

  # Local Variables
  utmp = OffsetArray(fill(1.e8, (indexdim(isd,ied),indexdim(jsd,jed))), isd:ied, jsd:jed)
  vtmp = OffsetArray(fill(1.e8, (indexdim(isd,ied),indexdim(jsd,jed))), isd:ied, jsd:jed)

  # -----------------------------
  # Interior:
  # -----------------------------
  for j = max(4, js - 1) : min(npy-4, je+1)
    for i = max(4, isd) : min(npx-4, ied)
      utmp[i, j] = a2 * (u[i, j-1] + u[i, j+2]) + a1 * (u[i, j] + u[i, j+1])
    end
  end
  for j = max(4, jsd) : min(npy-4, jed)
    for i = max(4, is - 1) : min(npx-4, ie+1)
      vtmp[i, j] = a2 * (v[i-1, j] + v[i+2, j]) + a1 * (v[i, j] + v[i+1, j])
    end
  end

  # -----------------------------
  # edges:
  # -----------------------------
  if js == 1 || jsd < 4
    for j = jsd:3
      for i = isd:ied
        utmp[i, j] = 0.5 * (u[i, j] + u[i, j+1])
        vtmp[i, j] = 0.5 * (v[i, j] + v[i+1, j])
      end
    end
  end

  if je + 1 == npy || jed >= npy - 4
    for j = npy-3:jed
      for i = isd:ied
        utmp[i, j] = 0.5 * (u[i, j] + u[i, j+1])
        vtmp[i, j] = 0.5 * (v[i, j] + v[i+1, j])
      end
    end
  end

  if is == 1 || isd < 4
    for j = max(4, jsd): min(npy-4, jed)
      for i = isd:3
        utmp[i, j] = 0.5 * (u[i, j] + u[i, j+1])
        vtmp[i, j] = 0.5 * (v[i, j] + v[i+1, j])
      end
    end
  end

  if ie + 1 == npx || ied >= npx - 4
    for j = max(4, jsd): min(npy-4, jed)
      for i = npx-3:ied
        utmp[i, j] = 0.5 * (u[i, j] + u[i, j+1])
        vtmp[i, j] = 0.5 * (v[i, j] + v[i+1, j])
      end
    end
  end

  for j = js-2:je+2
    for i = is-2:ie+2
      ua[i, j] = (utmp[i, j] - vtmp[i, j] * cosa_s[i, j]) * rsin2[i, j]
      va[i, j] = (vtmp[i, j] - utmp[i, j] * cosa_s[i, j]) * rsin2[i, j]
    end
  end

  # A - > C
  # -----------------------------
  # Fix the edges
  # -----------------------------
  # Xdir:
  if sw_corner
    for i = -2:0
      utmp[i, 0] = -vtmp[0, 1-i]
    end
  end
  if se_corner
    for i = 0:2
      utmp[npx+i, 0] = vtmp[npx, i+1]
    end
  end
  if ne_corner
    for i = 0:2
      utmp[npx+i, npy] = -vtmp[npx, je-i]
    end
  end
  if nw_corner
    for i = -2:0
      utmp[i, npy] = vtmp[0, je+i]
    end
  end

  ifirst = max(3,     is-1)
  ilast  = min(npx-2, ie+2)

  # ----------------------------------------------------------
  # 4th order interpolation for interior points:
  # -----------------------------------------------------------
  for j = js-1:je+1
    for i = ifirst:ilast
      uc[i, j] = a1 * (utmp[i-1, j] + utmp[i, j]) + a2 * (utmp[i-2, j] + utmp[i+1, j])
      ut[i, j] = (uc[i, j] - v[i, j] * cosa_u[i, j]) * rsin_u[i, j]
    end
  end

  # Xdir:
  if sw_corner
    ua[-1, 0] = -va[0, 2]
    ua[ 0, 0] = -va[0, 1]
  end
  if se_corner
    ua[npx,   0] = va[npx, 1]
    ua[npx+1, 0] = va[npx, 2]
  end
  if ne_corner
    ua[npx,   npy] = -va[npx, npy-1]
    ua[npx+1, npy] = -va[npx, npy-2]
  end
  if nw_corner
    ua[-1, npy] = va[0, npy-2]
    ua[ 0, npy] = va[0, npy-1]
  end

  if is == 1
    for j = js-1:je+1
      uc[0, j] = c1 * utmp[-2, j] + c2 * utmp[-1, j] + c3 * utmp[0, j]
      ut[1, j] = edge_interpolate4!(ua[-1:2, j], dxa[-1:2, j])
      #Want to  the UPSTREAM value
      if ut[1, j] < 0.0
        uc[1, j] = ut[1, j] * sin_sg[0, j, 3]
      else
        uc[1, j] = ut[1, j] * sin_sg[1, j, 1]
      end
      uc[2, j] = c1 * utmp[3, j] + c2 * utmp[2, j] + c3 * utmp[1, j]
      ut[0, j] = (uc[0, j] - v[0, j] * cosa_u[0, j]) * rsin_u[0, j]
      ut[2, j] = (uc[2, j] - v[2, j] * cosa_u[2, j]) * rsin_u[2, j]
    end

    if ie + 1 == npx
      for j = js-1:je+1
        uc[npx-1, j] =
        c1 * utmp[npx-3, j] + c2 * utmp[npx-2, j] + c3 * utmp[npx-1, j]
        i = npx
        ut[i, j] = 0.25 * (-ua[i-2, j] + 3.0 * (ua[i-1, j] + ua[i, j]) - ua[i+1, j])
        ut[i, j] = edge_interpolate4!(ua[i-2:i+1, j], dxa[i-2:i+1, j])
        if ut[i, j] < 0.0
          uc[i, j] = ut[i, j] * sin_sg[i-1, j, 3]
        else
          uc[i, j] = ut[i, j] * sin_sg[i, j, 1]
        end
        uc[npx+1, j] = c3 * utmp[npx, j] + c2 * utmp[npx+1, j] + c1 * utmp[npx+2, j]
        ut[npx-1, j] = (uc[npx-1, j] - v[npx-1, j] * cosa_u[npx-1, j]) * rsin_u[npx-1, j]
        ut[npx+1, j] = (uc[npx+1, j] - v[npx+1, j] * cosa_u[npx+1, j]) * rsin_u[npx+1, j]
      end
    end

  end

  # ---------
  # Ydir:
  # ---------
  if sw_corner
    for j = -2:0
      vtmp[0, j] = -utmp[1-j, 0]
    end
  end
  if nw_corner
    for j = 0:2
      vtmp[0, npy+j] = utmp[j+1, npy]
    end
  end
  if se_corner
    for j = -2:0
      vtmp[npx, j] = utmp[ie+j, 0]
    end
  end
  if ne_corner
    for j = 0:2
      vtmp[npx, npy+j] = -utmp[ie-j, npy]
    end
  end
  if sw_corner
    va[0, -1] = -ua[2, 0]
    va[0,  0] = -ua[1, 0]
  end
  if se_corner
    va[npx,  0] = ua[npx-1, 0]
    va[npx, -1] = ua[npx-2, 0]
  end
  if ne_corner
    va[npx,   npy] = -ua[npx-1, npy]
    va[npx, npy+1] = -ua[npx-2, npy]
  end
  if nw_corner
    va[0,   npy] = ua[1, npy]
    va[0, npy+1] = ua[2, npy]
  end

  for j = js-1:je+2
    if j == 1
      for i = is-1:ie+1
        vt[i, j] = edge_interpolate4!(va[i, -1:2], dya[i, -1:2])
        if vt[i, j] < 0.0
          vc[i, j] = vt[i, j] * sin_sg[i, j-1, 4]
        else
          vc[i, j] = vt[i, j] * sin_sg[i, j, 2]
        end
      end
    elseif j == 0 || j == npy - 1
      for i = is-1:ie+1
        vc[i, j] = c1 * vtmp[i, j-2] + c2 * vtmp[i, j-1] + c3 * vtmp[i, j]
        vt[i, j] = (vc[i, j] - u[i, j] * cosa_v[i, j]) * rsin_v[i, j]
      end
    elseif j == 2 || j == npy + 1
      for i = is-1:ie+1
        vc[i, j] = c1 * vtmp[i, j+1] + c2 * vtmp[i, j] + c3 * vtmp[i, j-1]
        vt[i, j] = (vc[i, j] - u[i, j] * cosa_v[i, j]) * rsin_v[i, j]
      end
    elseif j == npy
      for i = is-1:ie+1
        vt[i, j] = 0.25 * (-va[i, j-2] + 3.0 * (va[i, j-1] + va[i, j]) - va[i, j+1])
        vt[i, j] = edge_interpolate4!(va[i, j-2:j+1], dya[i, j-2:j+1])
        if vt[i, j] < 0.0
          vc[i, j] = vt[i, j] * sin_sg[i, j-1, 4]
        else
          vc[i, j] = vt[i, j] * sin_sg[i, j, 2]
        end
      end
    else
      # 4th order interpolation for interior points:
      for i = is-1:ie+1
        vc[i, j] = a2 * (vtmp[i, j-2] + vtmp[i, j+1]) + a1 * (vtmp[i, j-1] + vtmp[i, j])
        vt[i, j] = (vc[i, j] - u[i, j] * cosa_v[i, j]) * rsin_v[i, j]
      end
    end
  end
  return nothing
end # function d2a2c_vect

#------------------------------------------------------------------
# edge_interpolate4
#------------------------------------------------------------------
function edge_interpolate4!(ua, dxa)

  u0L = 0.5 * ((2.0 * dxa[2] + dxa[1]) * ua[2] - dxa[2] * ua[1]) / (dxa[1] + dxa[2])
  u0R = 0.5 * ((2.0 * dxa[3] + dxa[4]) * ua[3] - dxa[3] * ua[4]) / (dxa[3] + dxa[4])

  return u0L + u0R

  # This is the original edge - interpolation code, which makes
  # a relatively small increase in the error in unstretched case 2.
  #
  # edge_interpolate4 = 0.25 * [3 * [ua[2] + ua[3]] - [ua[1] + ua[4]]]

end #function edge_interpolate4

#----------------------------------------------------------------------------------
# fill2_4corners
#
#  This function fills the 4 corners of the scalar fields only as needed by c_core
#----------------------------------------------------------------------------------
function fill2_4corners!(q1, q2, dir, sw_corner, se_corner, ne_corner, nw_corner, npx, npy)

  if dir == 1
    if sw_corner
      q1[-1, 0] = q1[0, 2]
      q1[ 0, 0] = q1[0, 1]
      q2[-1, 0] = q2[0, 2]
      q2[ 0, 0] = q2[0, 1]
    end
    if se_corner
      q1[npx+1, 0] = q1[npx, 2]
      q1[npx,   0] = q1[npx, 1]
      q2[npx+1, 0] = q2[npx, 2]
      q2[npx,   0] = q2[npx, 1]
    end
    if nw_corner
      q1[ 0, npy] = q1[0, npy-1]
      q1[-1, npy] = q1[0, npy-2]
      q2[ 0, npy] = q2[0, npy-1]
      q2[-1, npy] = q2[0, npy-2]
    end
    if ne_corner
      q1[npx,   npy] = q1[npx, npy-1]
      q1[npx+1, npy] = q1[npx, npy-2]
      q2[npx,   npy] = q2[npx, npy-1]
      q2[npx+1, npy] = q2[npx, npy-2]
    end

  end

  if dir == 2
    if sw_corner
      q1[0,  0] = q1[1, 0]
      q1[0, -1] = q1[2, 0]
      q2[0,  0] = q2[1, 0]
      q2[0, -1] = q2[2, 0]
    end
    if se_corner
      q1[npx,  0] = q1[npx-1, 0]
      q1[npx, -1] = q1[npx-2, 0]
      q2[npx,  0] = q2[npx-1, 0]
      q2[npx, -1] = q2[npx-2, 0]
    end
    if nw_corner
      q1[0, npy]   = q1[1, npy]
      q1[0, npy+1] = q1[2, npy]
      q2[0, npy]   = q2[1, npy]
      q2[0, npy+1] = q2[2, npy]
    end
    if ne_corner
      q1[npx, npy]   = q1[npx-1, npy]
      q1[npx, npy+1] = q1[npx-2, npy]
      q2[npx, npy]   = q2[npx-1, npy]
      q2[npx, npy+1] = q2[npx-2, npy]
    end

  end

  return nothing
end # function fill2_4corners

#----------------------------------------------------------------------------------
# fill_4corners
#
#  This function fills the 4 corners of the scalar fields only as needed by c_core
#----------------------------------------------------------------------------------
function fill_4corners!(q, dir, sw_corner, se_corner, ne_corner, nw_corner, npx, npy)

  if dir == 1
    if sw_corner
      q[-1, 0] = q[0, 2]
      q[ 0, 0] = q[0, 1]
    end
    if se_corner
      q[npx+1, 0] = q[npx, 2]
      q[npx,   0] = q[npx, 1]
    end
    if nw_corner
      q[ 0, npy] = q[0, npy-1]
      q[-1, npy] = q[0, npy-2]
    end
    if ne_corner
      q[npx,   npy] = q[npx, npy-1]
      q[npx+1, npy] = q[npx, npy-2]
    end

  end

  if dir == 2
    if sw_corner
      q[0,  0] = q[1, 0]
      q[0, -1] = q[2, 0]
    end
    if se_corner
      q[npx,  0] = q[npx-1, 0]
      q[npx, -1] = q[npx-2, 0]
    end
    if nw_corner
      q[0, npy  ] = q[1, npy]
      q[0, npy+1] = q[2, npy]
    end
    if ne_corner
      q[npx, npy  ] = q[npx-1, npy]
      q[npx, npy+1] = q[npx-2, npy]
    end

  end

  return nothing
end #function fill_4corners

#----------------------------------------------------------------------------------
# interpolate_state
#
#  Increase the database by interpFactor.
#----------------------------------------------------------------------------------
function interpolate_state!(original_state, interpFactor::Int64)

  @unpack rarea, rarea_c, sin_sg, cos_sg, sina_v, cosa_v,
  sina_u, cosa_u, fC, rdxc, rdyc, dx, dy, dxc, dyc,
  cosa_s, rsin_u, rsin_v, rsin2, dxa, dya,
  delpc, delp, ptc, pt, u, v, w, uc, vc, ua, va, wc, ut,
  vt, divg_d  = original_state

  @unpack isd, ied, jsd, jed, is, ie, js, je, nord, npx, npy, npz = original_state

  odims = Array{Int64}(undef, 6)
  idims = Array{Int64}(undef, 6)

  odims[1] = isd
  odims[2] = ied
  odims[3] = jsd
  odims[4] = jed

  interpolateCalculateSpace2D!(odims, interpFactor, idims)

  new_rsin2 = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rsin2, odims, new_rsin2, idims, interpFactor)
  new_dxa = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(dxa, odims, new_dxa, idims, interpFactor)
  new_dya = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(dya, odims, new_dya, idims, interpFactor)
  new_cosa_s = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(cosa_s, odims, new_cosa_s, idims, interpFactor)
  new_rarea = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rarea, odims, new_rarea, idims, interpFactor)

  new_ied = idims[2] # Beginning subscripts are unchanged.
  new_jed = idims[4] # No changes to the 3rd dimension.

  odims[1] = isd
  odims[2] = ied
  odims[3] = jsd
  odims[4] = jed + 1

  interpolateCalculateSpace2D!(odims, interpFactor, idims)

  new_sina_v = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(sina_v, odims, new_sina_v, idims, interpFactor)
  new_cosa_v = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(cosa_v, odims, new_cosa_v, idims, interpFactor)
  new_rsin_v = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rsin_v, odims, new_rsin_v, idims, interpFactor)
  new_rdyc = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rdyc, odims, new_rdyc, idims, interpFactor)
  new_dx = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(dx, odims, new_dx, idims, interpFactor)
  new_dyc = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(dyc, odims, new_dyc, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied + 1
  odims[3] = jsd
  odims[4] = jed

  interpolateCalculateSpace2D!(odims, interpFactor, idims)

  new_sina_u = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(sina_u, odims, new_sina_u, idims, interpFactor)
  new_cosa_u = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(cosa_u, odims, new_cosa_u, idims, interpFactor)
  new_rdxc = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rdxc, odims, new_rdxc, idims, interpFactor)
  new_dy = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(dy, odims, new_dy, idims, interpFactor)
  new_dxc = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(dxc, odims, new_dxc, idims, interpFactor)
  new_rsin_u = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rsin_u, odims, new_rsin_u, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied + 1
  odims[3] = jsd
  odims[4] = jed + 1

  interpolateCalculateSpace2D!(odims, interpFactor, idims)

  new_rarea_c = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(rarea_c, odims, new_rarea_c, idims, interpFactor)
  new_fC = OffsetArray{Float64, 2}(undef, idims[1]:idims[2], idims[3]:idims[4])
  interpolateArray2D!(fC, odims, new_fC, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied
  odims[3] = jsd
  odims[4] = jed
  odims[5] = 1
  odims[6] = 9

  interpolateCalculateSpace3D!(odims, interpFactor, idims)

  new_sin_sg = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(sin_sg, odims, new_sin_sg, idims, interpFactor)
  new_cos_sg = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(cos_sg, odims, new_cos_sg, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied
  odims[3] = jsd
  odims[4] = jed
  odims[5] = 1
  odims[6] = npz

  interpolateCalculateSpace3D!(odims, interpFactor, idims)

  new_delpc = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(delpc, odims, new_delpc, idims, interpFactor)
  new_delp = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(delp, odims, new_delp, idims, interpFactor)
  new_ptc = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(ptc, odims, new_ptc, idims, interpFactor)
  new_pt = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(pt, odims, new_pt, idims, interpFactor)
  new_w = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(w, odims, new_w, idims, interpFactor)
  new_ua = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(ua, odims, new_ua, idims, interpFactor)
  new_va = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(va, odims, new_va, idims, interpFactor)
  new_wc = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(wc, odims, new_wc, idims, interpFactor)
  new_ut = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(ut, odims, new_ut, idims, interpFactor)
  new_vt = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(vt, odims, new_vt, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied + 1
  odims[3] = jsd
  odims[4] = jed
  odims[5] = 1
  odims[6] = npz

  interpolateCalculateSpace3D!(odims, interpFactor, idims)

  new_v = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(v, odims, new_v, idims, interpFactor)
  new_uc = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(uc, odims, new_uc, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied
  odims[3] = jsd
  odims[4] = jed + 1
  odims[5] = 1
  odims[6] = npz

  interpolateCalculateSpace3D!(odims, interpFactor, idims)

  new_u = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(u, odims, new_u, idims, interpFactor)
  new_vc = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(vc, odims, new_vc, idims, interpFactor)

  odims[1] = isd
  odims[2] = ied + 1
  odims[3] = jsd
  odims[4] = jed + 1
  odims[5] = 1
  odims[6] = npz

  interpolateCalculateSpace3D!(odims, interpFactor, idims)

  new_divg_d = OffsetArray{Float64, 3}(undef, idims[1]:idims[2], idims[3]:idims[4], idims[5]: idims[6])
  interpolateArray3D!(divg_d, odims, new_divg_d, idims, interpFactor)

  # Exchange the old dimensions to the new dimensions.
  ied = new_ied
  jed = new_jed
  is = 1
  ie = new_ied - 3
  js = 1
  je = new_jed - 3
  nord = 1
  npx = ie + 1
  npy = je + 1
  npz = 127

  # Redefine the dimensions in original_state Struct
  setfield!(original_state, :ied,  Int32(new_ied))
  setfield!(original_state, :jed,  Int32(new_jed))
  setfield!(original_state, :is,   Int32(1))
  setfield!(original_state, :ie,   Int32(new_ied - 3))
  setfield!(original_state, :js,   Int32(1))
  setfield!(original_state, :je,   Int32(new_jed - 3))
  setfield!(original_state, :nord, Int32(1))
  setfield!(original_state, :npx,  Int32(ie +1))
  setfield!(original_state, :npy,  Int32(je + 1))
  setfield!(original_state, :npz,  Int32(127))

  # Redefine the NetCDF dimensions in original_state Struct
  setfield!(original_state, :ni,   abs(isd - ied - 1))
  setfield!(original_state, :nj,   abs(jsd - jed - 1))
  setfield!(original_state, :nip1, abs(isd - ied - 2))
  setfield!(original_state, :njp1, abs(jsd - jed - 2))

  # Reassign the newly dimensioned and interpolated variables as OffsetArrays which slice off the elements past ied + 1 & jed +1, and are negatively indexed
  setfield!(original_state, :rarea , OffsetArray(new_rarea[  isd : ied   , jsd: jed      ],  isd : ied   , jsd: jed      ))
  setfield!(original_state, :rarea_c,OffsetArray(new_rarea_c[isd : ied+1 , jsd: jed+1    ],  isd : ied+1 , jsd: jed+1    ))
  setfield!(original_state, :sin_sg, OffsetArray(new_sin_sg[ isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :cos_sg, OffsetArray(new_cos_sg[ isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :sina_v, OffsetArray(new_sina_v[ isd : ied   , jsd: jed+1    ],  isd : ied   , jsd: jed+1    ))
  setfield!(original_state, :cosa_v, OffsetArray(new_cosa_v[ isd : ied   , jsd: jed+1    ],  isd : ied   , jsd: jed+1    ))
  setfield!(original_state, :sina_u, OffsetArray(new_sina_u[ isd : ied+1 , jsd: jed      ],  isd : ied+1 , jsd: jed      ))
  setfield!(original_state, :cosa_u, OffsetArray(new_cosa_u[ isd : ied+1 , jsd: jed      ],  isd : ied+1 , jsd: jed      ))
  setfield!(original_state, :fC    , OffsetArray(new_fC[     isd : ied+1 , jsd: jed+1    ],  isd : ied+1 , jsd: jed+1    ))
  setfield!(original_state, :rdxc  , OffsetArray(new_rdxc[   isd : ied+1 , jsd: jed      ],  isd : ied+1 , jsd: jed      ))
  setfield!(original_state, :rdyc  , OffsetArray(new_rdyc[   isd : ied   , jsd: jed+1    ],  isd : ied   , jsd: jed+1    ))
  setfield!(original_state, :dx    , OffsetArray(new_dx[     isd : ied   , jsd: jed+1    ],  isd : ied   , jsd: jed+1    ))
  setfield!(original_state, :dy    , OffsetArray(new_dy[     isd : ied+1 , jsd: jed      ],  isd : ied+1 , jsd: jed      ))
  setfield!(original_state, :dxc   , OffsetArray(new_dxc[    isd : ied+1 , jsd: jed      ],  isd : ied+1 , jsd: jed      ))
  setfield!(original_state, :dyc   , OffsetArray(new_dyc[    isd : ied   , jsd: jed+1    ],  isd : ied   , jsd: jed+1    ))
  setfield!(original_state, :cosa_s, OffsetArray(new_cosa_s[ isd : ied   , jsd: jed      ],  isd : ied   , jsd: jed      ))
  setfield!(original_state, :rsin_u, OffsetArray(new_rsin_u[ isd : ied+1 , jsd: jed      ],  isd : ied+1 , jsd: jed      ))
  setfield!(original_state, :rsin_v, OffsetArray(new_rsin_v[ isd : ied   , jsd: jed+1    ],  isd : ied   , jsd: jed+1    ))
  setfield!(original_state, :rsin2 , OffsetArray(new_rsin2[  isd : ied   , jsd: jed      ],  isd : ied   , jsd: jed      ))
  setfield!(original_state, :dxa   , OffsetArray(new_dxa[    isd : ied   , jsd: jed      ],  isd : ied   , jsd: jed      ))
  setfield!(original_state, :dya   , OffsetArray(new_dya[    isd : ied   , jsd: jed      ],  isd : ied   , jsd: jed      ))
  setfield!(original_state, :delpc , OffsetArray(new_delpc[  isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :delp  , OffsetArray(new_delp[   isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :ptc   , OffsetArray(new_ptc[    isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :pt    , OffsetArray(new_pt[     isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :u     , OffsetArray(new_u[      isd : ied   , jsd: jed+1, : ],  isd : ied   , jsd: jed+1, : ))
  setfield!(original_state, :v     , OffsetArray(new_v[      isd : ied+1 , jsd: jed  , : ],  isd : ied+1 , jsd: jed  , : ))
  setfield!(original_state, :w     , OffsetArray(new_w[      isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :uc    , OffsetArray(new_uc[     isd : ied+1 , jsd: jed  , : ],  isd : ied+1 , jsd: jed  , : ))
  setfield!(original_state, :ua    , OffsetArray(new_ua[     isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :vc    , OffsetArray(new_vc[     isd : ied   , jsd: jed+1, : ],  isd : ied   , jsd: jed+1, : ))
  setfield!(original_state, :va    , OffsetArray(new_va[     isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :wc    , OffsetArray(new_wc[     isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :ut    , OffsetArray(new_ut[     isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :vt    , OffsetArray(new_vt[     isd : ied   , jsd: jed  , : ],  isd : ied   , jsd: jed  , : ))
  setfield!(original_state, :divg_d, OffsetArray(new_divg_d[ isd : ied+1 , jsd: jed+1, : ],  isd : ied+1 , jsd: jed+1, : ))

  return nothing

end #function interpolate_state

end # module sw_core_mod


text:
function combined(op, kernel_function)
    function combined_kernel_function(u, v, k=nothing)
        ker = kernel_function(u, v)
        isnothing(k) ? ker : op.(k, ker)
    end
end

# Kernel functions

"""
    gaussiankernel(u, v)

Compute the Gaussian kernel, given by the formula
```
    κ(x,y) = exp(-‖x-y‖²/2).
```
`u` and `v` are matrices, where each column is a datapoint.
"""
function gaussiankernel(u, v)
    uu = sum(abs2, u, dims=1)
    vv = sum(abs2, v, dims=1)
    acc = u' * v
    @. acc = exp(acc - uu' / 2 - vv / 2)
    return acc
end

function rrule(::typeof(gaussiankernel), u, v)
    r = gaussiankernel(u, v)
    # r̄ is short for ∂l / ∂r
    function gaussiankernel_pullback(r̄)
        m = r .* r̄ # pullback the exponential
        ū = v * m'
        ū .-= sum(m', dims=1) .* u
        v̄ = u * m
        v̄ .-= sum(m, dims=1) .* v
        return NO_FIELDS, ū, v̄
    end
    return r, gaussiankernel_pullback
end

"""
    additivegaussiankernel(u, v, k=nothing)

Compute the radial basis function kernel (see [`gaussiankernel`](@ref)) on `u`, `v`.
If `k` is not `nothing`, add it to the result.
"""
const additivegaussiankernel = combined(+, gaussiankernel)

"""
    additivegaussiankernel(u, v, k=nothing)

Compute the radial basis function kernel (see [`gaussiankernel`](@ref)) on `u`, `v`.
If `k` is not `nothing`, multiply the result by it.
"""
const multiplicativegaussiankernel = combined(*, gaussiankernel)


text:
using MinionsREPL
using Documenter

DocMeta.setdocmeta!(MinionsREPL, :DocTestSetup, :(using MinionsREPL); recursive=true)

makedocs(;
    modules=[MinionsREPL],
    authors="Roger-luo <rogerluo.rl18@gmail.com> and contributors",
    repo="https://github.com/Roger-luo/MinionsREPL.jl/blob/{commit}{path}#{line}",
    sitename="MinionsREPL.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Roger-luo.github.io/MinionsREPL.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Roger-luo/MinionsREPL.jl",
)


text:
import Pkg
println(pwd())
if !isfile("Project.toml")
    println("Project.toml が見つからない。")
    exit(1)
end
Pkg.activate(".")
using HtmlTable, DataFrames
M = [
    "Jun. 2020" 1001 5;
    "Jun. 2020" 2055 12;
    "Jul. 2020" 1201 4;
    "Aug. 2020" 1001 5;
    "Aug. 2020" 2301 9;
    "Aug. 2020" 3002 18;
    "Sep. 2020" 1001 16;
    "Oct. 2020" 1201 3;
]
df = DataFrame(M, [:MONTH, :CODE, :AMOUNT])
html = fromDataFrame(df, title="From DataFrame")
write("./test/result_df.html", html)
println("test/result_df.html が更新されました。")

max_stars_repo_name: makandat/julia

max_issues_repo_name: makandat/julia

max_forks_repo_name: makandat/julia


text:
# Collection of protected function for GP

"""Protected division"""
pdiv(x, y, undef=10e6) = ifelse(y==0 , x+undef  , div(x,y))
"""Protected exponential"""
pexp(x, undef=10e15)   = ifelse(x>=32, x+undef  , exp(x))
"""Protected natural log"""
plog(x, undef=10e6)    = ifelse(x==0 , -undef   , log(abs(x)))
"""Protected sq.root"""
psqrt(x)               = sqrt(abs(x))
"""Protected sin(x)"""
psin(x, undef=10e6)    = ifelse(isinf(x), undef, sin(x))
"""Protected cos(x)"""
pcos(x, undef=10e6)    = ifelse(isinf(x), undef, cos(x))
"""Protected exponentiation operation"""
function ppow(x, y, undef=10e6)
    if y>=10
        x+y+undef
    elseif y < 1
        pow(abs(x),y)
    else
        pow(x,y)
    end
end


text:
using ImageMagick
using FactCheck, Images, Colors, FixedPointNumbers, ZipFile

workdir = joinpath(tempdir(), "Images")
writedir = joinpath(workdir, "write")
if !isdir(workdir)
    mkdir(workdir)
end
if !isdir(writedir)
    mkdir(writedir)
end

facts("Read remote") do
    urlbase = "http://www.imagemagick.org/Usage/images/"

    function getfile(name)
        file = joinpath(workdir, name)
        if !isfile(file)
            file = download(urlbase*name, file)
        end
        file
    end

    context("Gray") do
        file = getfile("jigsaw_tmpl.png")
        img = ImageMagick.load(file)
        @fact colorspace(img) --> "Gray"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> Gray{UFixed8}
        outname = joinpath(writedir, "jigsaw_tmpl.png")
        ImageMagick.save(outname, img)
        imgc = ImageMagick.load(outname)
        @fact img.data --> imgc.data
        @fact reinterpret(UInt32, data(map(mapinfo(RGB24, img), img))) -->
            map(x->x&0x00ffffff, reinterpret(UInt32, data(map(mapinfo(ARGB32, img), img))))
        @fact mapinfo(UInt32, img) --> mapinfo(RGB24, img)
        @fact data(convert(Image{Gray{Float32}}, img)) --> float32(data(img))
        mapi = mapinfo(RGB{UFixed8}, img)
        imgrgb8 = map(mapi, img)
        @fact imgrgb8[1,1].r --> img[1].val
        open(outname, "w") do file
            @compat show(file, MIME("image/png"), img)
        end
    end

    context("Gray with alpha channel") do
        file = getfile("wmark_image.png")
        img = ImageMagick.load(file)
        @fact colorspace(img) --> "GrayA"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> Images.ColorTypes.GrayA{UFixed8}
        if is_linux()
            outname = joinpath(writedir, "wmark_image.png")
            ImageMagick.save(outname, img)
            sleep(0.2)
            imgc = ImageMagick.load(outname)
            @fact img.data --> imgc.data
            open(outname, "w") do file
                @compat show(file, MIME("image/png"), img)
            end
        end
        @fact reinterpret(UInt32, data(map(mapinfo(RGB24, img), img))) -->
            map(x->x&0x00ffffff, reinterpret(UInt32, data(map(mapinfo(ARGB32, img), img))))
        @fact mapinfo(UInt32, img) --> mapinfo(ARGB32, img)
    end

    context("RGB") do
        file = getfile("rose.png")
        img = ImageMagick.load(file)
        # Mac reader reports RGB4, imagemagick reports RGB
        @fact colorspace(img) --> "RGB"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> RGB{UFixed8}
        outname = joinpath(writedir, "rose.tiff")
        ImageMagick.save(outname, img)
        imgc = ImageMagick.load(outname)
        T = eltype(imgc)
        # Why does this one fail on OSX??
        is_apple() || @fact img.data --> imgc.data
        @fact reinterpret(UInt32, data(map(mapinfo(RGB24, img), img))) -->
            map(x->x&0x00ffffff, reinterpret(UInt32, data(map(mapinfo(ARGB32, img), img))))
        @fact mapinfo(UInt32, img) --> mapinfo(RGB24, img)
        mapi = mapinfo(RGB{UFixed8}, img)
        imgrgb8 = map(mapi, img)
        @fact data(imgrgb8) --> data(img)
        convert(Array{Gray{UFixed8}}, img)
        convert(Image{Gray{UFixed8}}, img)
        convert(Array{Gray}, img)
        convert(Image{Gray}, img)
        imgs = separate(img)
        @fact permutedims(convert(Image{Gray}, imgs), [2,1]) --> convert(Image{Gray}, img)
        # Make sure that all the operations in README will work:
        buf = Array(UInt32, size(img))
        buft = Array(UInt32, reverse(size(img)))
        uint32color(img)
        uint32color!(buf, img)
        imA = convert(Array, img)
        uint32color(imA)
        uint32color!(buft, imA)
        uint32color(imgs)
        uint32color!(buft, imgs)
        imr = reinterpret(UFixed8, img)
        uint32color(imr)
        uint32color!(buf, imr)
        is_apple() || begin
            imhsv = convert(Image{HSV}, float32(img))
            uint32color(imhsv)
            uint32color!(buf, imhsv)
            @fact pixelspacing(restrict(img)) --> [2.0,2.0]
        end
        outname = joinpath(writedir, "rose.png")
        open(outname, "w") do file
            @compat show(file, MIME("image/png"), img)
        end
    end

    context("RGBA with 16 bit depth") do
        file = getfile("autumn_leaves.png")
        img = ImageMagick.load(file)
        @fact colorspace(img) --> "BGRA"
        @fact ndims(img) --> 2
        @fact colordim(img) --> 0
        @fact eltype(img) --> Images.ColorTypes.BGRA{UFixed16}
        outname = joinpath(writedir, "autumn_leaves.png")
        is_apple() || begin
            ImageMagick.save(outname, img)
            sleep(0.2)
            imgc = ImageMagick.load(outname)
            @fact img.data --> imgc.data
            @fact reinterpret(UInt32, data(map(mapinfo(RGB24, img), img))) -->
                map(x->x&0x00ffffff, reinterpret(UInt32, data(map(mapinfo(ARGB32, img), img))))
            @fact mapinfo(UInt32, img) --> mapinfo(ARGB32, img)
        end
        open(outname, "w") do file
            @compat show(file, MIME("image/png"), img)
        end
    end

    context("Indexed") do
        file = getfile("present.gif")
        img = ImageMagick.load(file)
        @fact nimages(img) --> 1
        @fact reinterpret(UInt32, data(map(mapinfo(RGB24, img), img))) -->
            map(x->x&0x00ffffff, reinterpret(UInt32, data(map(mapinfo(ARGB32, img), img))))
        @fact mapinfo(UInt32, img) --> mapinfo(RGB24, img)
        outname = joinpath(writedir, "present.png")
        open(outname, "w") do file
            @compat show(file, MIME("image/png"), img)
        end
    end

    context("Images with a temporal dimension") do
        fname = "swirl_video.gif"
        #fname = "bunny_anim.gif"  # this one has transparency but LibMagick gets confused about its size
        file = getfile(fname)  # this also has transparency
        img = ImageMagick.load(file)
        @fact timedim(img) --> 3
        @fact nimages(img) --> 26
        outname = joinpath(writedir, fname)
        ImageMagick.save(outname, img)
        imgc = ImageMagick.load(outname)
        # Something weird happens after the 2nd image (compression?), and one starts getting subtle differences.
        # So don't compare the values.
        # Also take the opportunity to test some things with temporal images
        @fact storageorder(img) --> ["x", "y", "t"]
        @fact haskey(img, "timedim") --> true
        @fact timedim(img) --> 3
        s = getindexim(img, 1:5, 1:5, 3)
        @fact timedim(s) --> 0
        s = sliceim(img, :, :, 5)
        @fact timedim(s) --> 0
        imgt = sliceim(img, "t", 1)
        @fact reinterpret(UInt32, data(map(mapinfo(RGB24, imgt), imgt))) -->
            map(x->x&0x00ffffff, reinterpret(UInt32, data(map(mapinfo(ARGB32, imgt), imgt))))
    end

    context("Extra properties") do
        is_apple() || begin
            file = getfile("autumn_leaves.png")
            # List properties
            extraProps = ImageMagick.load(file, extrapropertynames=true)

            img = ImageMagick.load(file,extraprop=extraProps)
            props = properties(img)
            for key in extraProps
                @fact haskey(props, key) --> true
                @fact props[key] --> not(nothing)
            end
            img = ImageMagick.load(file, extraprop=extraProps[1])
            props = properties(img)
            @fact haskey(props, extraProps[1]) --> true
            @fact props[extraProps[1]] --> not(nothing)

            println("The following \"Undefined property\" warning indicates normal operation")
            img = ImageMagick.load(file, extraprop="Non existing property")
            props = properties(img)
            @fact haskey(props, "Non existing property") --> true
            @fact props["Non existing property"] --> nothing
        end
    end
end

facts("EXIF orientation") do
    function test_orientation(r, odict)
        for f in r.files
            bn = basename(f.name)
            if haskey(odict, bn)
                so = odict[bn]
                data = read(f, UInt8, f.uncompressedsize)
                img = readblob(data)
                @fact spatialorder(img) --> so
            end
        end
    end

    url = "http://www.galloway.me.uk/media/other/EXIF_Orientation_Samples.zip"
    fn = joinpath(workdir, "EXIF_Orientation_Samples.zip")
    download(url, fn)
    r = ZipFile.Reader(fn)
    test_orientation(r, Dict("up.jpg"=>["x", "y"],
                             "left-mirrored.jpg"=>["y", "x"]))
end