text:
using Statistics
using StatsBase

"""
Bounds given metric between -1.0 and 1.0, where 1.0 is perfect fit.

Suitable for use with any metric that ranges from 1 to -∞.

# References
1. Mathevet, T., Michel, C., Andréassian, V., Perrin, C., 2006.
    A bounded version of the Nash-Sutcliffe criterion for better model
    assessment on large sets of basins.
    IAHS-AISH Publication 307, 211–219.
    https://iahs.info/uploads/dms/13614.21--211-219-41-MATHEVET.pdf

# Example
```julia
julia> import Streamfall: @bound, KGE
julia> @bound KGE([1,2], [3,2])
-0.35653767993482094
```
"""
macro bound(metric)
    tmp = :($metric)
    return :($tmp / (2.0 - $tmp))
end

"""
Normalizes given metric between 0.0 and 1.0, where 1.0 is perfect fit.

Suitable for use with any metric that ranges from 1 to -∞.

# References
1. Nossent, J., Bauwens, W., 2012.
    Application of a normalized Nash-Sutcliffe efficiency to improve the
    accuracy of the Sobol’ sensitivity analysis of a hydrological model.
    EGU General Assembly Conference Abstracts 237.

# Example
```julia
julia> import Streamfall: @normalize, KGE
julia> @normalize KGE([1,2], [3,2])
0.1111111111111111
```
"""
macro normalize(metric)
    return :(1.0 / (2.0 - $metric))
end

"""
Applies mean inverse approach to a metric.

Suitable for use with any metric that ranges from 1 to -∞.

If using with other macros such as `@normalize` or `@bound`,
these must come first.

# References
1. Garcia, F., Folton, N., Oudin, L., 2017.
    Which objective function to calibrate rainfall–runoff
        models for low-flow index simulations?
    Hydrological Sciences Journal 62, 1149–1166.
    https://doi.org/10.1080/02626667.2017.1308511

# Example
```julia
julia> import Streamfall: @normalize, @mean_inverse, KGE
julia> @normalize @mean_inverse KGE [1,2] [3,2] 1e-6
0.3193506006429825
```
"""
macro mean_inverse(metric, obs, sim, ϵ=1e-2)
    obj, o, s, ϵ = eval(metric), eval(obs), eval(sim), eval(ϵ)
    q = obj(o, s)
    q2 = obj(1.0 ./ (o .+ ϵ), 1.0 ./ (s .+ ϵ))
    return mean([q, q2])
end

"""
Applies split meta metric approach

If using with other macros such as `@normalize` or `@bound`,
these must come first.

# References
1. Fowler, K., Peel, M., Western, A., Zhang, L., 2018.
    Improved Rainfall-Runoff Calibration for Drying Climate:
    Choice of Objective Function.
    Water Resources Research 54, 3392–3408.
    https://doi.org/10.1029/2017WR022466

# Example
```julia
julia> using Statistics
julia> import Streamfall: @normalize, @split, KGE
julia> @normalize @split KGE repeat([1,2], 365) repeat([3,2], 365) 365 mean
0.3217309561946589
```
"""
macro split(metric, obs, sim, n, agg_func=mean)
    obj, o, s, t, func = eval(metric), eval(obs), eval(sim), eval(n), eval(agg_func)
    return naive_split_metric(o, s; n_members=t, metric=obj, comb_method=func)
end

"""The Nash-Sutcliffe Efficiency score"""
NSE(obs, sim) = 1.0 - sum((obs .- sim).^2) / sum((obs .- mean(obs)).^2)

"""Normalized Nash-Sutcliffe Efficiency score (bounded between 0 and 1).

# References
1. Nossent, J., Bauwens, W., 2012.
    Application of a normalized Nash-Sutcliffe efficiency to improve the accuracy of the Sobol’ sensitivity analysis of a hydrological model.
    EGU General Assembly Conference Abstracts 237.
"""
NNSE(obs, sim) = 1.0 / (2.0 - NSE(obs, sim))

"""Root Mean Square Error"""
RMSE(obs, sim) = (sum((sim .- obs).^2)/length(sim))^0.5

"""Coefficient of determination (R^2)

Aliases `NSE()`
"""
function R2(obs, sim)::Float64
    return NSE(obs, sim)
end

"""Adjusted R^2

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
- `p::Int` : number of explanatory variables
"""
function ADJ_R2(obs, sim, p::Int64)::Float64
    n = length(obs)
    adj_r2 = 1 - (1 - R2(obs, sim)) * ((n - 1) / (n - p - 1))

    return adj_r2
end

"""
Mean Absolute Error
"""
MAE(obs, sim) = mean(abs.(sim .- obs))

"""
Mean Error
"""
ME(obs, sim) = mean(sim .- obs)

"""
    PBIAS(obs::Vector, sim::Vector)::Float64

Percent bias between `sim` and `obs`

Model performance for streamflow can be determined to be
satisfactory if the Nash-Sutcliffe Efficiency (NSE)
score > 0.5, the RMSE standard deviation ratio (RSR) < 0.7
and percent bias (PBIAS) is +/- 25% (see [1]).

# References
1. Moriasi, D.N., Arnold, J.G., Liew, M.W.V., Bingner, R.L.,
    Harmel, R.D., Veith, T.L., 2007.
    Model Evaluation Guidelines for Systematic Quantification
    of Accuracy in Watershed Simulations.
    Transactions of the ASABE 50, 885–900.
    https://doi.org/10.13031/2013.23153
"""
PBIAS(obs, sim) = (sum(obs .- sim) * 100) / sum(obs)

"""
    RSR(obs::Vector, sim::Vector)::Float64

The RMSE-observations standard deviation ratio (RSR).

Varies between 0 and a large positive value, where 0
indicates an RMSE value of 0.

# References
1. Moriasi, D.N., Arnold, J.G., Liew, M.W.V., Bingner, R.L.,
    Harmel, R.D., Veith, T.L., 2007.
    Model Evaluation Guidelines for Systematic Quantification
    of Accuracy in Watershed Simulations.
    Transactions of the ASABE 50, 885–900.
    https://doi.org/10.13031/2013.23153
"""
function RSR(obs, sim)::Float64
    rmse = RMSE(obs, sim)
    σ_obs = std(obs)
    rsr = rmse / σ_obs
    return rsr
end

"""
    EV(obs, sim)

Explained Variance.

Indicates the amount of variation in the observations which the predictions
are able to explain.
"""
function EV(obs, sim)
    return 1.0 - (var(obs - sim) / obs)
end

"""
Relative Skill Score.

Provides an indication of model performance relative to a known benchmark score.

Suitable for use with least-squares approaches that provide skill scores ranging from 1 to -∞.

# Arguments
- `Sb` : Benchmark score
- `Sm` : Model score

# References
1. Knoben, W.J.M., Freer, J.E., Woods, R.A., 2019.
    Technical note: Inherent benchmark or not?
        Comparing Nash-Sutcliffe and Kling-Gupta efficiency scores (preprint).
    Catchment hydrology/Modelling approaches.
    https://doi.org/10.5194/hess-2019-327
"""
function relative_skill_score(Sb::Float64, Sm::Float64)
    return (Sm - Sb) / (1.0 - Sb)
end

"""
    NSE_logbias(obs, sim; metric::Function=NSE, bias_threshold::Float64=5.0, shape::Float64=2.5)

The NSE_logbias meta-metric provides a weighted combination of a least-squares approach and a
logarithmic function of bias. The metric penalizes predictions with an overall bias above a
threshold (defined as 5% in [1]).

It is also referred to as the Viney F score.

# Extended help
The penalty applied is non-symmetrical (or multiplicatively symmetrical) in that
predictions that are _double_ the observed are penalized identically to predictions that are
_half_ the observed volume.

# Arguments
- `obs::Vector` : Historic observations to compare against
- `sim::Vector` : Modeled time series
- `metric::Function` : least-squares method to use, defaults to NSE
- `bias_threshold::Float64` : Bias threshold after which the score given by `metric` is penalized, defaults to 5 (%)
- `shape::Float64` : Exponent value controlling shape of penalization (see Figure 2 in [1]).

# References
1. Viney, N. R., Perraud, J., Vaze, J., Chiew, F.H.S., Post, D.A., Yang, A. 2009
    The usefulness of bias constraints in model calibration for regionalisation
        to ungauged catchments
    18th World IMACS / MODSIM Congress, Cairns, Australia, 13 - 17 July 2009
    Available at:
    https://www.researchgate.net/publication/294697092_The_usefulness_of_bias_constraints_in_model_calibration_for_regionalisation_to_ungauged_catchments

2. Teng, J., Potter, N.J., Chiew, F.H.S., Zhang, L., Wang, B., Vaze, J., Evans, J.P., 2015.
    How does bias correction of regional climate model precipitation affect modelled runoff?
    Hydrology and Earth System Sciences 19, 711–728.
    https://doi.org/10.5194/hess-19-711-2015
"""
function NSE_logbias(obs, sim; metric=NSE, bias_threshold=5.0, shape=2.5)
    E_ns = metric(obs, sim)
    B = PBIAS(obs, sim) / 100.0
    return E_ns - bias_threshold * abs(log(1 + B))^shape
end

"""
    KGE(obs::Vector, sim::Vector; scaling::Tuple=nothing)::Float64

Calculate the 2009 Kling-Gupta Efficiency (KGE) metric.

Decomposes NSE into correlation (`r`), relative variability (`α`), and bias (`β`) terms.

A KGE score of 1 means perfect fit.
A score < -0.41 indicates that the mean of observations
provides better estimates (see Knoben et al., [2]).

The `scaling` argument expects a three-valued tuple
which scales `r`, `α` and `β` factors respectively.
If not specified, defaults to `1`.

Note: Although similar, NSE and KGE cannot be directly compared.

# References
1. Gupta, H.V., Kling, H., Yilmaz, K.K., Martinez, G.F., 2009.
    Decomposition of the mean squared error and NSE performance criteria:
    Implications for improving hydrological modelling.
    Journal of Hydrology 377, 80–91.
    https://doi.org/10.1016/j.jhydrol.2009.08.003

2. Knoben, W.J.M., Freer, J.E., Woods, R.A., 2019.
    Technical note: Inherent benchmark or not? Comparing Nash-Sutcliffe and Kling-Gupta efficiency scores (preprint).
    Catchment hydrology/Modelling approaches.
    https://doi.org/10.5194/hess-2019-327

3. Mizukami, N., Rakovec, O., Newman, A.J., Clark, M.P., Wood, A.W.,
    Gupta, H.V., Kumar, R., 2019.
    On the choice of calibration metrics for “high-flow”
        estimation using hydrologic models.
    Hydrology and Earth System Sciences 23, 2601–2614.
    https://doi.org/10.5194/hess-23-2601-2019
"""
function KGE(obs, sim; scaling=nothing)::Float64
    if isnothing(scaling)
        scaling = (1, 1, 1)
    end

    # Correlation
    r = Statistics.cor(obs, sim)
    if isnan(r)
        r = 0.0
    end

    # relative variance
    α = std(sim) / std(obs)

    # bias
    β = mean(sim) / mean(obs)

    rs = scaling[1]
    as = scaling[2]
    bs = scaling[3]

    kge = 1 - sqrt(rs*(r - 1)^2 + as*(α - 1)^2 + bs*(β - 1)^2)

    return kge
end

"""Bounded KGE, bounded between -1 and 1.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
"""
function BKGE(obs, sim)::Float64
    kge = KGE(obs, sim)
    return kge / (2 - kge)
end

"""Normalized KGE between 0 and 1.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
"""
function NKGE(obs, sim; scaling=nothing)::Float64
    return 1 / (2 - KGE(obs, sim; scaling=scaling))
end

"""Calculate the modified KGE metric (2012).

Also known as KGE prime (KGE').

# Extended help

It is not recommended to apply KGE' with log-transformed flows (see [2]).
Numerical instabilities arise as flow approaches values close to zero.
This is possible under extreme dry conditions or by chance when sub-sampling.

In cases where observations are constant or otherwise displays zero variance or zero
mean flow, this implementation applies a simple logistic function (ℯ⁻ˣ) to gain an
indication of simulated data's distance to zero.

This is to:
- avoid NaNs influencing subsequent calculations
- allow use with split methods which may partition streamflows into periods of 0 flows.

# Arguments
- `obs::Vector`: observations
- `sim::Vector` : modeled results
- `scaling::Tuple` : scaling factors in order of timing (r), magnitude (β), variability (γ).
                     Defaults to (1,1,1).

# References
1. Kling, H., Fuchs, M., Paulin, M., 2012.
    Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios.
    Journal of Hydrology 424–425, 264–277.
    https://doi.org/10.1016/j.jhydrol.2012.01.011

2. Santos, L., Thirel, G., Perrin, C., 2018.
    Technical note: Pitfalls in using log-transformed flows within the KGE criterion.
    Hydrology and Earth System Sciences 22, 4583–4591.
    https://doi.org/10.5194/hess-22-4583-2018
"""
function mKGE(obs, sim; scaling=nothing)::Float64
    if isnothing(scaling)
        scaling = (1,1,1)
    end

    # Timing (Pearson's correlation)
    r = Statistics.cor(obs, sim)
    if isnan(r)
        # can happen if obs or sim is a constant (std of 0)
        r = all(obs .== sim) ? 1.0 : 0.0
    end

    # Variability
    cv_s = StatsBase.variation(sim)
    if isnan(cv_s)
        cv_s = 0.0
    end

    cv_o = StatsBase.variation(obs)
    if isnan(cv_o)
        cv_o = 0.0
    end

    if cv_o == 0.0 && cv_s == 0.0
        γ = 1.0
    elseif cv_o == 0.0
        # Use logistic function to indicate distance from 0
        γ = 1.0 / exp(cv_s)
    else
        γ = cv_s / cv_o
    end

    # Magnitude
    μ_o = mean(obs)
    μ_s = mean(sim)
    if μ_o == 0.0
        # use logistic function to indicate distance from 0 (μ_o)
        β = 1.0 / exp(μ_s)
    else
        β = μ_s / μ_o
    end

    rs = scaling[1]
    βs = scaling[2]
    γs = scaling[3]

    mod_kge = 1.0 - sqrt(rs*(r - 1)^2 + βs*(β - 1)^2 + γs*(γ - 1)^2)

    return mod_kge
end

"""Bounded modified KGE between -1 and 1.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
"""
function BmKGE(obs, sim; scaling=nothing)::Float64
    mkge = mKGE(obs, sim; scaling=scaling)
    return mkge / (2 - mkge)
end

"""Normalized modified KGE between 0 and 1.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
"""
function NmKGE(obs, sim; scaling=nothing)::Float64
    return 1 / (2 - mKGE(obs, sim; scaling=scaling))
end

"""
Mean Inverse NmKGE

Said to produce better fits for low-flow indices
compared to mKGE (see [1]).

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
- `scaling::Tuple` : scaling factors for r, α, and β (defaults to 1.0)
- `ϵ::Float64` : small constant to use with inverse flow to allow consideration of periods with no flow. Defaults to 1e-2.

# References
1. Garcia, F., Folton, N., Oudin, L., 2017.
    Which objective function to calibrate rainfall–runoff
        models for low-flow index simulations?
    Hydrological Sciences Journal 62, 1149–1166.
    https://doi.org/10.1080/02626667.2017.1308511
"""
mean_NmKGE(obs, sim; scaling=nothing, ϵ=1e-2) = mean([Streamfall.NmKGE(obs, sim; scaling=scaling), Streamfall.NmKGE(1.0 ./ (obs .+ ϵ), 1.0 ./ (sim .+ ϵ); scaling=scaling)])

"""Calculate the non-parametric Kling-Gupta Efficiency (KGE) metric.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled
- `scaling::Tuple` : scaling factors for timing (s), variability (α), magnitude (β)

# References
1. Pool, S., Vis, M., Seibert, J., 2018.
    Evaluating model performance: towards a non-parametric variant of the Kling-Gupta efficiency.
    Hydrological Sciences Journal 63, 1941–1953.
    https://doi.org/10.1080/02626667.2018.1552002

"""
function npKGE(obs, sim; scaling=nothing)::Float64
    if isnothing(scaling)
        scaling = (1,1,1)
    end

    # flow duration curves
    μ_s = mean(sim)
    if μ_s == 0.0
        fdc_sim = repeat([0.0], length(sim))
    else
        x = length(sim) * μ_s
        fdc_sim = sort(sim / x)
    end

    μ_o = mean(obs)
    if μ_o == 0.0
        fdc_obs = repeat([0.0], length(obs))
    else
        x = length(obs) * μ_o
        fdc_obs = sort(obs / x)
    end

    α = 1 - 0.5 * sum(abs.(fdc_sim - fdc_obs))

    # Magnitude
    if μ_o == 0.0
        β = μ_s == 0.0 ? 1.0 : μ_s
    else
        β = μ_s / μ_o
    end

    # Timing and flow dynamics
    r = StatsBase.corspearman(fdc_obs, fdc_sim)
    if isnan(r)
        r = 1.0  # can occur if identical sequences are used (e.g., 0 flows)
    end

    rs = scaling[1]
    αs = scaling[2]
    βs = scaling[3]

    kge = 1 - sqrt(rs*(r - 1)^2 + αs*(α - 1)^2 + βs*(β - 1)^2)

    return kge
end

"""Bounded non-parametric KGE between -1 and 1.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
"""
function BnpKGE(obs, sim; scaling=nothing)::Float64
    npkge = npKGE(obs, sim; scaling=scaling)
    return npkge / (2 - npkge)
end

"""Normalized non-parametric KGE between 0 and 1.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results
"""
function NnpKGE(obs, sim; scaling=nothing)::Float64
    return 1 / (2 - npKGE(obs, sim; scaling=scaling))
end

"""Liu Mean Efficiency metric (LME).

Reformulation of the KGE metric said to be advantageous for capturing extreme
flow events.

# Arguments
- `obs::Vector` : observations
- `sim::Vector` : modeled results

# References
1. Liu, D., 2020.
    A rational performance criterion for hydrological model.
    Journal of Hydrology 590, 125488.
    https://doi.org/10.1016/j.jhydrol.2020.125488
"""
function LME(obs, sim)::Float64
    μ_o = mean(obs)
    μ_s = mean(sim)
    β = (μ_s / μ_o)

    r = Statistics.cor(obs, sim)
    σ_s = std(sim)
    σ_o = std(obs)
    α = (σ_s / σ_o)
    k_1 = r * α

    LME = 1 - sqrt((k_1 - 1)^2 + (β - 1)^2)

    return LME
end

function naive_split_metric(obs::Vector, sim::Vector, n_members::Int, metric::Function=NNSE)
    obs_chunks = collect(Iterators.partition(obs, n_members))
    sim_chunks = collect(Iterators.partition(sim, n_members))
    scores = Array{Float64, 1}(undef, length(obs_chunks))

    for (idx, h_chunk) in enumerate(obs_chunks)
        scores[idx] = metric(h_chunk, sim_chunks[idx])
    end

    return scores
end

"""Naive approach to split metrics.

Split metrics are a meta-objective optimization approach which "splits" data
into subperiods. The objective function is calculated for each subperiod and
then recombined. The approach addresses the lack of consideration of dry years
with least-squares.

In Fowler et al., [1] the subperiod is one year. The implementation offered here
is "naive" in that the data is partitioned into `N` chunks of `n_members` and
does not consider date/time.

# Arguments
- `obs::Vector` : Historic observations to compare against
- `sim::Vector` : Modeled time series
- `n_members::Int` : number of members per chunk (i.e., sub-samples), defaults to 365
- `metric::Function` : Objective function to apply, defaults to NNSE
- `comb_method::Function` : Recombination method, defaults to `mean`

# References
1. Fowler, K., Peel, M., Western, A., Zhang, L., 2018.
    Improved Rainfall-Runoff Calibration for Drying Climate: Choice of Objective Function.
    Water Resources Research 54, 3392–3408.
    https://doi.org/10.1029/2017WR022466
"""
function naive_split_metric(obs, sim; n_members::Int=365, metric::Function=NNSE, comb_method::Function=mean)
    scores = naive_split_metric(obs, sim, n_members, metric)
    return comb_method(scores)
end

"""
    inverse_metric(obs, sim; metric, comb_method::Function=mean)

A meta-objective function which combines the performance of the
given metric as applied to the discharge and the inverse of the
discharge.

By default, the combination method is to take the mean.

# Arguments
- obs : observed
- sim : modeled results
- `metric::Function` : objective function
- `comb_method::Function` : mean
- ϵ : offset value to use (enables use with zero-flow time steps), defaults to 1e-2

# References
1. Garcia, F., Folton, N., Oudin, L., 2017.
    Which objective function to calibrate rainfall–runoff models
        for low-flow index simulations?
    Hydrological Sciences Journal 62, 1149–1166.
    https://doi.org/10.1080/02626667.2017.1308511
"""
function inverse_metric(obs, sim, metric::Function; comb_method::Function=mean, ϵ=1e-2)
    return comb_method([metric(obs, sim), metric(1.0 ./ (obs .+ ϵ), 1.0 ./ (sim .+ ϵ))])
end

"""
Allows comparison of any model compared against a pre-defined benchmark, assuming both scores were obtained with
the same objective function.

Positive values indicate a model is better than the benchmark, and negative values indicate a model performs worse.

# Extended help
It is noted in Knoben et al., [1] that the skill score should always be contextualized with the original benchmark
value. Interpreting skill scores by themselves may become difficult if the benchmark score is already quite high.
A small improvement of no real practical value could be misconstrued as a large improvement.
As an example, if the benchmark has an KGE score of 0.999 and its counterpart 0.9995, then a skill score of 0.5 will
be reported.

# References
1. Knoben, W.J.M., Freer, J.E., Woods, R.A., 2019.
    Technical note: Inherent benchmark or not? Comparing Nash-Sutcliffe and Kling-Gupta efficiency scores (preprint).
    Catchment hydrology/Modelling approaches.
    https://doi.org/10.5194/hess-2019-327

2. Towner, J., Cloke, H.L., Zsoter, E., Flamig, Z., Hoch, J.M., Bazo, J., Coughlan de Perez, E., Stephens, E.M., 2019.
    Assessing the performance of global hydrological models for capturing peak river flows in the Amazon basin.
    Hydrology and Earth System Sciences 23, 3057–3080.
    https://doi.org/10.5194/hess-23-3057-2019
"""
function skill_score(model_score, benchmark_score)
    return (model_score - benchmark_score) / (1.0 - benchmark_score)
end


text:
import TOML

# Just to identify the output for a run, use the current time
df = Dates.format(Dates.now(), DateFormat("YYYY-mm-ddTHH:MM:SS"))
_default_outfile() = "deformation_" * Dates.format(Dates.now(), df) * ".h5"

# TODO: figure out how to incorporate the argparse checking
# with this TOML parsing

DEFAULT_KEYS = Dict(
    ### Files used in solving:

    "unw_stack_file" => UNW_FILENAME,
    "input_dset" => STACK_FLAT_SHIFTED_DSET,
    "outfile" => _default_outfile(),
    # This will be the dataset in the output file
    "outdset" => "velos",

    # name of .h5 file with mask data
    "mask_stack_file" => "masks.h5",
    "ignore_geo_file" => "",

    # If we want to add specific extra, list comma separated
    # TODO: implement or remove
    "ignore_geo_dates" => "",
    "max_temporal_baseline" => 500,

    ### Type of solver ###
    "constant_velocity" => true,
    # Instead of SBAS: just average all igrams by date in a stack
    "stack_average" => false,

    ### Extra solver options ###
    # Use L1 norm for SBAS cost instead of L2 least squares
    "L1" => true,
    # (SBAS) Strength of Tikhonov regularization
    "alpha" => 0.0,
)

ext: jl


text:
"""
Abstract type for models using a Gaussian/normal-inverse-gamma conjugate prior.
"""
abstract type AbstractBayesianLM end

"""
    BayesLM(d::Int; <keyword arguments>)

Construct a Bayesian linear model.

# Arguments

- `λ::Float64=1.0`: Prior scaling.
- `shape0::Float64=1e-3`: Inverse-gamma prior shape hyperparameter.
- `scale0::Float64=1e-3`: Inverse-gamma prior scale hyperparameter.
"""
mutable struct BayesLM <: AbstractBayesianLM
    shape0::Float64
    scale0::Float64
    β0::Matrix{Float64}
    Σ0::Hermitian{Float64,Matrix{Float64}}
    Λ0::Hermitian{Float64,Matrix{Float64}}
    shape::Float64
    scale::Float64
    β::Matrix{Float64}
    Σ::Hermitian{Float64,Matrix{Float64}}
    Λ::Hermitian{Float64,Matrix{Float64}}

    function BayesLM(d::Int64; λ::Float64=1.0, shape0::Float64=1e-3, scale0::Float64=1e-3)
        if d <= 0 || shape0 <= 0 || scale0 <= 0
            msg = "d, shape0 and scale0 must be strictly positive"
            throw(ArgumentError(msg))
        end
        β0, Σ0, Λ0 = zeros(d, 1),
        Hermitian(diagm(ones(d)) * λ),
        Hermitian(diagm(ones(d)) / λ)
        β, Σ, Λ = deepcopy(β0), deepcopy(Σ0), deepcopy(Λ0)
        return new(shape0, scale0, β0, Σ0, Λ0, shape0, scale0, β, Σ, Λ)
    end
end

function sherman_morrison_inv(Ainv::AbstractMatrix, u::AbstractMatrix, v::AbstractMatrix)
    if size(u) != size(v)
        throw(DimensionMismatch("u and v have mismatched dimensions"))
    elseif size(u, 2) != 1
        throw(ArgumentError("u and v should have exactly one column"))
    elseif size(Ainv) != (size(u, 1), size(u, 1))
        throw(DimensionMismatch("size of A does not match u/v"))
    end
    numer = Ainv * u * v' * Ainv
    denom = 1.0 + (v' * Ainv * u)[1, 1]
    return Ainv - numer / denom
end

function check_regression_data(X::AbstractMatrix, y::AbstractMatrix)
    if size(y, 1) != 1
        throw(ArgumentError("y should have exactly 1 row"))
    elseif size(X, 2) != size(y, 2)
        throw(DimensionMismatch("X and y don't match in dimension 2"))
    end
end

function fit!(lm::BayesLM, X::AbstractMatrix, y::AbstractMatrix)
    check_regression_data(X, y)
    if size(lm.β, 1) != size(X, 1)
        msg = "X does not have the same number of predictors as the model"
        throw(DimensionMismatch(msg))
    end
    G = Hermitian(X * X')
    Λ = G + lm.Λ
    Σ = size(X, 2) == 1 ? Hermitian(sherman_morrison_inv(lm.Σ, X, X)) : inv(Λ)
    β = Σ * (X * y' + lm.Λ * lm.β)
    shape = lm.shape + size(y, 2) / 2
    scale = lm.scale + 0.5 * (y * y' - β' * Λ * β + lm.β' * lm.Λ * lm.β)[1, 1]
    return lm.shape, lm.scale, lm.β, lm.Σ, lm.Λ = shape, scale, β, Σ, Λ
end

function (lm::BayesLM)(X::AbstractMatrix)
    if size(X, 1) != size(lm.β, 1)
        msg = "X does not have the same number of predictors as the model"
        throw(DimensionMismatch(msg))
    end
    return lm.β' * X
end

function shape_scale(lm::BayesLM)
    return lm.shape, lm.scale
end

"""
    std(model::AbstractBayesianLM)

Return the posterior mean of the inverse-gamma distribution.
"""
function std(model::AbstractBayesianLM)
    shape, scale = shape_scale(model)
    if shape <= 1
        throw(ArgumentError("normal-inverse-gamma has undefined mean for shape ≤ 1"))
    end
    return sqrt(scale / (shape - 1))
end


text:
module GreenFunc
using StaticArrays
using ..Lehmann
using ..CompositeGrids
# Write your package code here.

include("green/Green.jl")
export TimeDomain, ImTime, ReTime, ImFreq, ReFreq, DLRFreq
export Green2DLR, toTau, toMatFreq, toDLR, dynamic, instant

end


text:
## Comparison script
##
# using GitHub, JSON, PkgBenchmark
#
# baseline_cfg = BenchmarkConfig(id="master", juliacmd = `julia -O3`);
# cfg = BenchmarkConfig(juliacmd = `julia -O3`);
# results = judge("StaticPolynomials", cfg, baseline_cfg);
#
# gist_json = JSON.parse(
#     """
#     {
#     "description": "A benchmark for StaticPolynomials",
#     "public": false,
#     "files": {
#         "benchmark.md": {
#         "content": "$(escape_string(sprint(export_markdown, results)))"
#         }
#     }
#     }
#     """
# );
#
# posted_gist = create_gist(params = gist_json);
#
# url = get(posted_gist.html_url)
#
# println("Result available at:")
# println(url)
#
# if is_unix()
#     run(`open $url`)
# end

# # Benchmark script
# #
# using GitHub, JSON
# using PkgBenchmark
#
# cfg = BenchmarkConfig(juliacmd = `julia -O3`);
# results = benchmarkpkg("StaticPolynomials", cfg);
#
# gist_json = JSON.parse(
#     """
#     {
#     "description": "A benchmark for StaticPolynomials",
#     "public": false,
#     "files": {
#         "benchmark.md": {
#         "content": "$(escape_string(sprint(export_markdown, results)))"
#         }
#     }
#     }
#     """
# );
#
# posted_gist = create_gist(params = gist_json);
#
# url = get(posted_gist.html_url)
#
# println("Result available at:")
# println(url)
#
# if is_unix()
#     run(`open $url`)
# end

using BenchmarkTools
import StaticPolynomials
const SP = StaticPolynomials
using StaticArrays
using PolynomialTestSystems

function bevaluate(T, testsystem)
    F = SP.system(PolynomialTestSystems.equations(testsystem))
    n = SP.nvariables(F)
    x = SVector{n, T}(rand(T, n))
    @benchmarkable SP.evaluate($F, $x)
end

function static_bjacobian(T, testsystem)
    F = SP.system(PolynomialTestSystems.equations(testsystem))
    n = SP.nvariables(F)
    x = SVector{n, T}(rand(T, n))
    @benchmarkable SP.jacobian($F, $x)
end

function bjacobian(T, testsystem)
    F = SP.system(PolynomialTestSystems.equations(testsystem))
    n = SP.nvariables(F)
    x = rand(T, n)
    U = zeros(T, SP.npolynomials(F), n)
    @benchmarkable SP.jacobian!($U, $F, $x)
end

const SUITE = BenchmarkGroup()
SUITE["evaluate"] = BenchmarkGroup(["evaluate"])
SUITE["jacobian"] = BenchmarkGroup(["jacobian"])
SUITE["static jacobian"] = BenchmarkGroup(["static jacobian"])

systems = [
    (katsura(5), "katsura5"), (katsura(6), "katsura(6)"),
    (katsura(7), "katsura7"), (katsura(8), "katsura8"),
    (katsura(9), "katsura9"), (katsura(10), "kastura10"),
    (chandra(4), "chandra4"), (chandra(8), "chandra8"),
    (cyclic(6), "cyclic6"), (cyclic(7), "cyclic7"), (cyclic(9), "cyclic9"),
    (fourbar(), "fourbar"), (rps10(), "rps10")]

for T in [Float64, Complex{Float64}]
    SUITE["evaluate"][string(T)] = BenchmarkGroup()
    SUITE["jacobian"][string(T)] = BenchmarkGroup()
    SUITE["static jacobian"][string(T)] = BenchmarkGroup()
    for (system, name) in systems
        SUITE["evaluate"][string(T)][name] = bevaluate(T, system)
        SUITE["jacobian"][string(T)][name] = bjacobian(T, system)
        SUITE["static jacobian"][string(T)][name] = static_bjacobian(T, system)
    end
end


text:
using Persa
using Statistics
using Random

mutable struct Baseline{T<:AbstractFloat} <: Persa.Model
    γ::T
    λ::T
    μ::T
    bu::Vector{T}
    bv::Vector{T}
end

function Baseline(dataset::Dataset; γ::Float64 = 0.001, λ::Float64 = 0.02)
    u = users(dataset)
    v = items(dataset)

    μ = mean(dataset)
    bu = zeros(u)
    bv = zeros(v)

    return Baseline(γ, λ, μ, bu, bv)
end

function Persa.predict(model::Baseline, pref::Preference, user::Int, item::Int)
    model.μ + model.bu[user] + model.bv[item]
end

function Persa.objective(model::Baseline, dataset::Dataset)
    pref = preference(dataset)
    total = 0

    for (u, v, r) in dataset
        total += (r - predict(model, pref, u, v))^2
        total += 0.1 * (model.bu[u]^2 + model.bv[v]^2)
    end

    return total
end

function Persa.update!(model::Baseline, dataset::Dataset)
    γ = model.γ
    λ = model.λ
    pref = preference(dataset)
    idx = shuffle(1:length(dataset))

    for i in 1:length(dataset)
        u, v, r = dataset[idx[i]]

        e = r - predict(model, pref, u, v)

        model.bu[u] += γ * (e .- λ * model.bu[u])
        model.bv[v] += γ * (e .- λ * model.bv[v])
    end
end

data = [
    (user = 1, item = 1, rating = 4),
    (user = 1, item = 3, rating = 2),
    (user = 2, item = 1, rating = 1),
    (user = 2, item = 2, rating = 5),
    (user = 2, item = 3, rating = 3),
    (user = 3, item = 1, rating = 3),
    (user = 3, item = 2, rating = 4),
    (user = 4, item = 2, rating = 5),
    (user = 3, item = 3, rating = 4)
]

dataset = Dataset(data)

model = Baseline(dataset)

mach = machine(model, dataset)

train!(mach, dataset)

mach[1, 3]


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export AddLayer

@neuron type AddNeuron
end

@neuron forward(neuron::AddNeuron) do
    neuron.value = neuron.inputs[1][1] + neuron.inputs[2][1]
end

@neuron backward(neuron::AddNeuron) do
    neuron.∇inputs[1][1] = neuron.∇
    neuron.∇inputs[2][1] = neuron.∇
end

function AddEnsemble(name::Symbol, net::Net, shape)
    neurons = Array(AddNeuron, shape...)
    for i in 1:length(neurons)
        neurons[i] = AddNeuron()
    end
    Ensemble(net, name, neurons)
end

function AddLayer(name::Symbol, net::Net, input_1::AbstractEnsemble, input_2::AbstractEnsemble)
    @assert size(input_1) == size(input_2)
    ens = AddEnsemble(name, net, size(input_1))
    mapping = one_to_one(ndims(input_1))
    add_connections(net, input_1, ens, mapping)
    add_connections(net, input_2, ens, mapping)
    ens
end

import Base: +, *

function +(net::Net, ens1::AbstractEnsemble, ens2::AbstractEnsemble)
    AddLayer(gensym("ensemble"), net, ens1, ens2)
end

function +(net::Net, ens1::AbstractEnsemble, ens...)
    AddLayer(gensym("ensemble"), net, ens1, +(net, ens...))
end

export MulLayer

@neuron type MulNeuron
end

@neuron forward(neuron::MulNeuron) do
    neuron.value = neuron.inputs[1][1] * neuron.inputs[2][1]
end

@neuron backward(neuron::MulNeuron) do
    neuron.∇inputs[1][1] = neuron.inputs[2][1] * neuron.∇
    neuron.∇inputs[2][1] = neuron.inputs[1][1] * neuron.∇
end

function MulEnsemble(name::Symbol, net::Net, shape)
    neurons = Array(MulNeuron, shape...)
    for i in 1:length(neurons)
        neurons[i] = MulNeuron()
    end
    Ensemble(net, name, neurons)
end

function MulEnsemble(net::Net, shape)
    MulEnsemble(gensym("ensemble"), net, shape)
end

function MulLayer(name::Symbol, net::Net, input_1::AbstractEnsemble, input_2::AbstractEnsemble)
    @assert size(input_1) == size(input_2)
    ens = MulEnsemble(name, net, size(input_1))
    mapping = one_to_one(ndims(input_1))
    add_connections(net, input_1, ens, mapping)
    add_connections(net, input_2, ens, mapping)
    ens
end

function *(net::Net, ens1::AbstractEnsemble, ens2::AbstractEnsemble)
    MulLayer(gensym("ensemble"), net, ens1, ens2)
end


text:
import HarmonicBalance: drop_powers, max_power
import HarmonicBalance.Symbolics.expand

@variables a,b,c

@test max_power(a^2 + b, a) == 2
@test max_power( a*((a+b)^4 )^2 + a, a) == 9

@test isequal(drop_powers(a^2+b, a, 1) , b)
@test isequal(drop_powers( (a+b)^2, a,1), b^2)
@test isequal(drop_powers( (a+b)^2, [a,b],1), 0)

@test isequal(drop_powers((a+b)^3 + (a+b)^5 , [a,b],4), expand((a+b)^3))


text:
using SPIRV
using Base.Test

# we use CUDAnative to handle all the GPU codegen modification to Julias codegen
using CUDAnative, LLVM

function kernel_vadd(a, b, c)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    c[i] = a[i] + b[i]
    return nothing
end

const CUArray = CUDAnative.CuDeviceArray{Float32, 1}

ir = sprint() do io
    CUDAnative.code_llvm(io, kernel_vadd, Tuple{CUArray, CUArray, CUArray}; dump_module = true, cap = v"2.0")
end
mod = parse(LLVM.Module, ir)

spirv_dir = joinpath(homedir(), "gpustuff", "SPIRV-LLVM", "build", "bin")
llvm_spirv = joinpath(spirv_dir, "llvm-spirv")
llvm_dis = joinpath(spirv_dir, "llvm-dis")
@assert isfile(llvm_spirv)

open(`$llvm_dis -`, "w") do f
    buff = convert(Vector{UInt8}, mod)
    write(f, buff)
end


text:
function replace(exprn::Expr, dict::Dict{Union(Symbol, Expr), Union(Symbol, Expr)})
  queue = Array(Union(Number, Symbol, Expr), 0)
  unshift!(queue, exprn)

  while length(queue) != 0
    node = shift!(queue)

    for i in 2:length(node.args)
      if isa(node.args[i], Expr)
        if haskey(dict, node.args[i]) && get(dict, node.args[i], :()) != :()
          node.args[i] = dict[node.args[i]]
        else
          unshift!(queue, node.args[i])
        end
      elseif isa(node.args[i], Symbol) && haskey(dict, node.args[i]) && get(dict, node.args[i], :()) != :()
        node.args[i] = dict[node.args[i]]
      end
    end
  end
end


text:
using DiffEqBase
using SciMLBase: @add_kwonly, add_kwonly
using LinearAlgebra, Test

@add_kwonly function f(a, b; c=3, d=4)
  (a, b, c, d)
end
@test f(1, 2) == (1, 2, 3, 4)
@test f(a=1, b=2) == (1, 2, 3, 4)
@test_throws ErrorException f()

@add_kwonly g(a, b; c=3, d=4) = (a, b, c, d)
@test g(1, 2) == (1, 2, 3, 4)
@test g(a=1, b=2) == (1, 2, 3, 4)

@add_kwonly h(; c=3, d=4) = (c, d)
@test h() == (3, 4)

@test_throws ErrorException add_kwonly(:(i(c=3, d=4) = (c, d)))

dprob = DiscreteProblem((u,p,t)->2u,0.5,(0.0,1.0))
@test remake(dprob) == dprob
@test remake(dprob; u0 = 1.0).u0 == 1.0

oprob = ODEProblem((u,p,t)->2u,0.5,(0.0,1.0))
@test remake(oprob) == oprob
@test remake(oprob; u0 = 1.0).u0 == 1.0

sprob = SDEProblem((u,p,t)->2u,(u,p,t)->2u,0.5,(0.0,1.0))
@test remake(sprob) == sprob
@test remake(sprob; u0 = 1.0).u0 == 1.0

daeprob = DAEProblem((du,u,p,t)->du-2u,0.5,0.5,(0.0,1.0))
@test remake(daeprob) == daeprob
@test remake(daeprob; u0 = 1.0).u0 == 1.0

ddeprob = DDEProblem((du,u,h,p,t)->-2u,0.5,(p,t)->0.0,(0.0,1.0))
@test remake(ddeprob) == ddeprob
@test remake(daeprob; u0 = 1.0).u0 == 1.0

function f(du,u,p,t)
  du[1] = 0.2u[1]
  du[2] = 0.4u[2]
end
u0 = ones(2)
tspan = (0,1.0)

# Create a ODEProblem and test remake:
prob1 = SplitODEProblem(f,f,u0,tspan,Dict(),callback=nothing)
prob2 = @inferred remake(prob1; u0 = prob1.u0 .+ 1)
@test prob1.f === prob2.f
@test prob1.p === prob2.p
@test prob1.u0 .+ 1 ≈ prob2.u0
@test prob1.tspan == prob2.tspan
@test prob1.kwargs[:callback] === prob2.kwargs[:callback]
@test prob1.problem_type === prob2.problem_type

prob2 = @inferred remake(prob1; u0 = prob1.u0 .+ 1, callback = :test)
@test prob2.kwargs[:callback] == :test

# Test remake with SplitFunction:
prob1 = SplitODEProblem((u,p,t) -> u/2, (u,p,t) -> 2u, 1.0, (0.0,1.0))
prob2 = remake(prob1;  # prob1 is a ODEProblem
  f = remake(prob1.f;  # prob1.f is a SplitFunction
    f2 = (u,p,t) -> 3u))

# Test remake with NoiseProblem (a struct w/o isinplace type parameter):
struct DummyNoiseProcess <: SciMLBase.AbstractNoiseProcess{Int,1,Nothing,true}
    dummy
end
tspan1 = (0.0, 1.0)
tspan2 = (0.0, 2.0)
noise1 = NoiseProblem(DummyNoiseProcess(Dict()), tspan1);
noise2 = remake(noise1; tspan=tspan2);
@test noise1.noise === noise2.noise
@test noise1.tspan == tspan1
@test noise2.tspan == tspan2
@test noise1.tspan != noise2.tspan

# Test remake with TwoPointBVPFunction (manually defined):
f1 = SciMLBase.TwoPointBVPFunction(() -> 1)
f2 = remake(f1; bc = () -> 2)
@test f1.bc() == 1
@test f2.bc() == 2

# Testing remake for no recompile
u0 = [0;2.]
tspan = (0.0,6.3)
prob = ODEProblem{true,false}((du,u,p,t) -> 2u,u0,tspan)

# Function wrapper makes a new reference, it doesn't make it type instable
@test_broken remake(prob) == prob
@test remake(prob; u0 = [1;2]).u0 == [1;2]
@test remake(prob; p = (1,2)).p == (1,2)


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

__precompile__()

module Preprocessing

using Revise, JSON, LinearAlgebra

export initializeapp, initializeconfig, setconfig, initializemdp, preprocessing
export adddata, readconfig, readconfigfile, readMDP, readEXTXYZ, readJSON, readJSONfromFolder
export setlattice, setregion, setdomain, savebin

include("rotatelattice.jl");
include("transformcoords.jl");
include("initializeapp.jl");
include("initializeconfig.jl");
include("initializemdp.jl");
include("setunits.jl");
include("readconfig.jl");
include("readconfigfile.jl");
include("readMDP.jl");
include("readEXTXYZ.jl");
include("readJSON.jl");
include("readJSONfromFolder.jl");
include("readweight.jl");
include("datastruct.jl");
include("writeapp.jl");
include("writeconfig.jl");
include("checkconfig.jl");
include("readLAMMPS.jl");
include("cubemapping.jl");
include("boxperiodicimages.jl");
include("setlattice.jl");
include("setregion.jl");
include("setdomain.jl");
include("writelattice.jl");
include("writeregion.jl");
include("writedomain.jl");

function preprocessing(app)

if app.datafile != "" & app.dataformat != ""
    config = readconfigfile(app);
    app.nconfigs = config.nconfigs;
elseif app.configmode>0
    # read configurations from data file
    config = readconfig(app, app.configmode);

    # number of configurations
    app.nconfigs = config.nconfigs;
else
    app.nconfigs = 0;
    config = 0;
end

if app.snapcoefffile != ""
    app.snapcoeff = Main.DelimitedFiles.readdlm(app.snapcoefffile);
end

# nonbonded single potentials
if  length(app.pot1a[:]) > 0
    app.np1a = length(app.pot1a);  # number of nonbonded single potentials
    app.ncmu1a = length(app.mu1a); # length of mu1a
else
    app.np1a = 0;
    app.ncmu1a = 0;
end

# bonded single potentials
if  length(app.pot1b[:]) > 0
    app.np1b = length(app.pot1b);  # number of nonbonded single potentials
    app.ncmu1b = length(app.mu1b); # length of mu1a
else
    app.np1b = 0;
    app.ncmu1b = 0;
end

# nonbonded pair potentials
if  length(app.pot2a[:]) > 0
    app.np2a = length(app.pot2a);  # number of nonbonded single potentials
    app.ncmu2a = length(app.mu2a); # length of mu1a
else
    app.np2a = 0;
    app.ncmu2a = 0;
end

# bonded pair potentials
if  length(app.pot2b[:]) > 0
    app.np2b = length(app.pot2b);  # number of nonbonded single potentials
    app.ncmu2b = length(app.mu2b); # length of mu1a
else
    app.np2b = 0;
    app.ncmu2b = 0;
end

# two-body bond order potentials
if  length(app.pot2c[:]) > 0
    app.np2c = length(app.pot2c);  # number of nonbonded single potentials
    app.ncmu2c = length(app.mu2c); # length of mu1a
else
    app.np2c = 0;
    app.ncmu2c = 0;
end

# nonbonded triplet potentials
if  length(app.pot3a[:]) > 0
    app.np3a = length(app.pot3a);  # number of nonbonded single potentials
    app.ncmu3a = length(app.mu3a); # length of mu1a
else
    app.np3a = 0;
    app.ncmu3a = 0;
end

# bonded triplet potentials
if  length(app.pot3b[:]) > 0
    app.np3b = length(app.pot3b);  # number of nonbonded single potentials
    app.ncmu3b = length(app.mu3b); # length of mu1a
else
    app.np3b = 0;
    app.ncmu3b = 0;
end

# three-body bond order potentials
if  length(app.pot3c[:]) > 0
    app.np3c = length(app.pot3c);  # number of nonbonded single potentials
    app.ncmu3c = length(app.mu3c); # length of mu1a
else
    app.np3c = 0;
    app.ncmu3c = 0;
end

# nonbonded quadruplet potentials
if  length(app.pot4a[:]) > 0
    app.np4a = length(app.pot4a);  # number of nonbonded single potentials
    app.ncmu4a = length(app.mu4a); # length of mu1a
else
    app.np4a = 0;
    app.ncmu4a = 0;
end

# bonded quadruplet potentials
if  length(app.pot4b[:]) > 0
    app.np4b = length(app.pot4b);  # number of nonbonded single potentials
    app.ncmu4b = length(app.mu4b); # length of mu1a
else
    app.np4b = 0;
    app.ncmu4b = 0;
end

style = app.unitstyle;
if lowercase(style) == "lj"
    app.unitstylenum = 0;
elseif lowercase(style) == "real"
    app.unitstylenum = 1;
elseif lowercase(style) == "metal"
    app.unitstylenum = 2;
elseif lowercase(style) == "si"
    app.unitstylenum = 3;
elseif lowercase(style) == "cgs"
    app.unitstylenum = 4;
elseif lowercase(style) == "electron"
    app.unitstylenum = 5;
elseif lowercase(style) == "micro"
    app.unitstylenum = 6;
elseif lowercase(style) == "nano"
    app.unitstylenum = 7;
else
    error("Invalid unit style");
end

# descriptor flag: 0 -> Spherical Harmonics Bessel, 1-> snap
if (lowercase(app.descriptor) == "shp")
    app.descriptornum = 0;
elseif (lowercase(app.descriptor) == "snap")
    app.descriptornum = 1;
else
    app.descriptornum = -1;
end

ensemblemode = app.ensemblemode;
if (lowercase(ensemblemode) == "nve")
    app.ensemblemodenum = 0;
    app.runMD = 1;
elseif (lowercase(ensemblemode) == "nvelimit")
    app.ensemblemodenum = 1;
    app.runMD = 1;
elseif (lowercase(ensemblemode) == "nvt")
    app.ensemblemodenum = 2;
    app.runMD = 1;
else
    app.ensemblemodenum = -1;
    app.runMD = 0;
end

app.natomtype = maximum([length(app.atommasses), length(app.atomspecies)]);
#app.natomtype = length(app.atommasses);
if length(app.atomnumbers)>0
    app.atomnumbers = [0 app.atomnumbers];
else
    app.atomnumbers = reshape([0],1,1);
end
if length(app.atommasses)>0
    app.atommasses = [0.0 app.atommasses];
else
    app.atommasses = reshape([0.0],1,1);
end
if length(app.atomcharges)>0
    app.atomcharges = [0.0 app.atomcharges];
else
    app.atomcharges = reshape([0.0],1,1);
end

app.flag = [app.descriptornum app.spectrum app.training app.runMD app.potentialform app.neighpair app.energycal app.forcecal app.stresscal app.neighcell app.decomposition app.chemtype app.dftdata app.unitstylenum app.ensemblemodenum app.neighcheck];

if length(app.rcut2a) > 0
    rcut2a = app.rcut2a;
else
    rcut2a = reshape([0.0],1,1);
end
if length(app.rcut2b) > 0
    rcut2b = app.rcut2b;
else
    rcut2b = reshape([0.0],1,1);
end
if length(app.rcut2c) > 0
    rcut2c = app.rcut2c;
else
    rcut2c = reshape([0.0],1,1);
end
if length(app.rcut3a) > 0
    rcut3a = app.rcut3a;
else
    rcut3a = reshape([0.0],1,1);
end
if length(app.rcut3b) > 0
    rcut3b = app.rcut3b;
else
    rcut3b = reshape([0.0],1,1);
end
if length(app.rcut3c) > 0
    rcut3c = app.rcut3c;
else
    rcut3c = reshape([0.0],1,1);
end
if length(app.rcut4a) > 0
    rcut4a = app.rcut4a;
else
    rcut4a = reshape([0.0],1,1);
end
if length(app.rcut4b) > 0
    rcut4b = app.rcut4b;
else
    rcut4b = reshape([0.0],1,1);
end
rcut = [reshape([app.rcutml],1,1) rcut2a rcut2b rcut2c rcut3a rcut3b rcut3c rcut4a rcut4b];
rcutmax = maximum(rcut) + app.neighskin;
app.rcutsqmax = rcutmax^2;
app.boxoffset = [rcutmax rcutmax rcutmax];
app.simparam = [app.time app.dt];
app.solparam = [rcutmax app.neighskin];
app.snaparam = [app.snapnelem app.snapncoeff app.snaptwojmax app.snaprcutfac app.snaprfac0 app.snaprmin0 app.snapbzeroflag app.snapswitchflag app.snapquadraticflag app.snapchemflag  app.snapbnormflag app.snapwselfallflag];

app.ndims = zeros(20,1);
app.ndims[1] = app.dim;
app.ndims[2] = app.L;
app.ndims[3] = app.K;
app.ndims[4] = app.ntimesteps;
app.ndims[5] = app.nab;
app.ndims[6] = app.mpiprocs;
app.ndims[7] = app.backend;
app.ndims[8] = app.nconfigs;
app.ndims[9] = app.natomtype;
app.ndims[10] = app.nmoletype;
app.ndims[11] = app.neighmax;
app.ndims[12] = app.neighevery;
app.ndims[13] = app.neighdelay;
app.ndims[14] = app.globalfreq;
app.ndims[15] = app.peratomfreq;

if app.configmode>0
    m = 1;
    for i = 1:config.nconfigs
        B2C, C2B = cubemapping(config.a[:,i], config.b[:,i], config.c[:,i]);
        ximages = boxperiodicimages(app.pbc, config.a[:,i], config.b[:,i], config.c[:,i]);
        n = config.natom[i];
        config.x[:,m:(m+n-1)] = checkconfig(config.x[:,m:(m+n-1)], ximages, B2C, C2B);
        m = m + n;
    end
end

bindir = "exec";
if app.lattice != []
    writelattice(app.lattice, app.appname * "lattice.bin");
    mv(app.appname * "lattice.bin", app.sourcepath * bindir * "/" * app.appname * "lattice.bin", force = true);
end
if app.region != []
    writeregion(app.region, app.appname * "region.bin");
    mv(app.appname * "region.bin", app.sourcepath * bindir * "/" * app.appname * "region.bin", force = true);
end
if app.domain != []
    writedomain(app.domain, app.appname * "domain.bin");
    mv(app.appname * "domain.bin", app.sourcepath * bindir * "/" * app.appname * "domain.bin", force = true);
end

if app.nconfigs>0
    if app.training == 0
        config.we = reshape([],0,2);
        config.wf = reshape([],0,2);
    else
        config = readweight(app, config);
    end
    writeconfig(config, app.appname * "config.bin");
    mv(app.appname * "config.bin", app.sourcepath * bindir * "/" * app.appname * "config.bin", force = true);
end

writeapp(app, app.appname * "app.bin");
mv(app.appname * "app.bin", app.sourcepath * bindir * "/"  * app.appname * "app.bin", force = true);

return app, config

end

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
function test_flow_improve()

    A = connected_cliques([6,5,12],true) # get a connected ring
    R = Set{Int64}(1:5)
    push!(R,7)

    S,cval,vol,x = flow_improve(A,R)
    @assert S == Set(1:11)

    S,cval,vol,x = flow_improve(A,R;kappa=1)
    @assert S == Set(1:6)

    A = connected_cliques([5,5],false)
    R = Set{Int64}(1:5)
    A = sparse(A)
    S,cval,vol,x = flow_improve(A,R)

    @assert S == R

    R = Set{Int64}(6:10)
    S,cval,vol,x = flow_improve(A,R)

    @assert S == R

    # this graph is disconnected, but MQI won't see it, but flow improve will
    A = connected_cliques([5,6,7,2,3,6,7],false)
    A[11,12] = 0. # disconnect these cliques
    A[12,11] = 0.
    R = Set{Int64}(6:18)
    S,cval,vol,x = flow_improve(A,R)

    @assert S == Set(1:11)
    @assert cval/vol == 0.0

    # so if we set kappa large enough (e^10 seems like enough)
    S,cval,vol,x = flow_improve(A,R; kappa = 10)

    @assert S == Set(12:18)
    @assert cval/vol == 0.023255813953488372

    A[1,end] = 1.
    A[end,1] = 1. # connect the ring

    S,cval,vol,x = flow_improve(A,R)

    @assert S == Set(12:18)
    @assert cval/vol == 0.023255813953488372

end

test_flow_improve()

function test_methods()
    A = sparse([0 5. 0.; 5. 2. 3.; 0. 3. 0.])

    B,w = incidence_matrix_and_weight(A)
    @test B == [1 -1. 0; 0. 1 -1]
    @test w == [5.; 3.]

    x,val = stmincut(A,1,3)
    @test val == 3
    @test x == [1.; 1.; 0.]

    A[1,2] = 1.5;

    x,val = stmincut(A,1,3)
    @test val == 1.5
    @test x == [1.; 0.; 0.]

    # mincut test from
    # http://tcs.rwth-aachen.de/lehre/Graphentheorie/WS2013/Kiril_Mitev.pdf
    A = sparse([0 9 0 8 0.;
                9 0 9 1 2;
                0 9 0 7 4;
                8 1 7 0 3;
                0 2 4 3 0])
    x,val = stmincut(A,1,5)
    @test val == 9.
    @test x == [1.; 1.; 1.; 1.; 0.]

    # create a pair of cliques and verify that we get back the same clique
    A = [ones(5,5) zeros(5,5); zeros(5,5) ones(5,5)]
    A[5,6] = 1.
    A[6,5] = 1.
    A = A - diagm(diag(A))
    R = Set{Int64}(1:5)
    A = sparse(A)
    S,cval,vol,x = mqi(A,R)

    @test S == R

    R = Set{Int64}(6:10)
    S,cval,vol,x = mqi(A,R)

    @test S == R

    A = connected_cliques([5,6,7,2,3,6,7],false)
    A[11,12] = 0.
    A[12,11] = 0.
    R = Set{Int64}(6:18)
    S,cval,vol,x = mqi(A,R)

    @test S == Set(12:18)
    @test cval/vol == 0.023255813953488372

end

test_methods()


text:
module COBRAUtils

include("io.jl")
include("manipulations.jl")
include("fba.jl")
include("AV.jl")
include("lp-problem.jl")

include("Media.jl")

end # module


text:
function WNNM(Y::Array{Float32,3}, C::Float32, Sig::Array{Float32,1}, m::Array{Float32,3})

    F = svd(F2cplx(Y))

    lenSig = length(F.S)

    patNum = size(Y)[2]

    mean_Sig_arr = √(sum(Sig[:].^2) / length(Sig[:]))

    TempC  = Float32(C * √(patNum) * 2 * mean_Sig_arr^2)

    X_arr = ClosedQWNNM(F.S, TempC,eps())

    X = F.U * Diagonal(X_arr) * F.Vt

    return cplx2F(X) .+ m
end

function ClosedQWNNM(SigmaY::Array{Float32,1}, c::Float32, our_eps)

    temp = (SigmaY .- our_eps) .^ 2 - 4 * ((-1) * our_eps .* SigmaY .+ c)
    ind = Int32[]
    SigmaX = Float32[]
    for i in 1:length(temp)
        if temp[i] > 0
            push!(ind, i)
            push!(SigmaX, max(SigmaY[i] - our_eps + √(temp[i]) , 0) / 2)
        else
            push!(SigmaX, 0)
        end
    end

    return SigmaX
end

function F2cplx(M::Array{Float32,3})::Array{Complex{Float32},2}
    (h, w) = size(M)

    cM = adjoint(M)

    [cM[1:2*h, 1:w]; cM[1:2*h, w+1:2*w]]
end

function adjoint(A::Array{Float32,3})::Array{Complex{Float32},2}
    W = zeros(size(A)[1:2])
    X = A[:, :, 1]
    Y = A[:, :, 2]
    Z = A[:, :, 3]
    imm = Complex{Float32}(im)
    vcat(hcat(W + imm .* X, Y + imm .* Z),hcat(-Y + imm .* Z, W - imm .* X))
end

function cplx2F(M)
    (h, w) = size(M)

    cM = [M[1:h ÷ 2,1:w] M[1+h ÷ 2:h,1 : w]]

    unadjoint(cM)
end
function unadjoint(A)
    (r, c) = size(A)
    r2 = r ÷ 2
    c2 = c ÷ 2
    A1 = A[     1 : r2,      1 : c2]
    A2 = A[     1 : r2, c2 + 1 : c ]
    A3 = A[r2 + 1 : r,       1 : c2]
    A4 = A[r2 + 1 : r,  c2 + 1 : c ]
    C = zeros(Float32, r2, c2, 3)
    C[:,:,1] = Float32.(imag(A1 - A4))
    C[:,:,2] = Float32.(real(A2 - A3))
    C[:,:,3] = Float32.(imag(A2 + A3) )
    C./2
end


text:
module ModelSanitizer

import Requires

export sanitize!
export Model
export Data
export ForceSanitize

include("types.jl")

include("arrays.jl")
include("elements.jl")
include("sanitize.jl")
include("utils.jl")
include("zero.jl")

function __init__()::Nothing
    Requires.@require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" include("dataframes.jl")
    return nothing
end

end # module


text:
include("optimization.jl")

function test_preconditioner(L::Integer, N_sweeps::Integer, k_batch::Integer, equation::Function, name::String)
    h = 1/2^L
    k = sqrt(1/h)
    parameters_ = [k, ]
    BCs = ["D", "D", "D", "D"]
    symmetry = "symmetric"

    params = optimization_parameters(equation, ADAM(0.001), parameters_, 10, k_batch, N_sweeps, BCs, L, name, symmetry, "Spectral")
    A = params.equation(params.L, params.BCs, params.eq_params)
    prec, B, parameters_start = initialise_preconditioner(A, params.L, params.BCs, params.name, params.symmetry)
    history, prec, parameters = stopping_optimization(params; check_time=30, tol=1e-3)

    cond_init = measure_prec_condition_number(parameters_start, prec)
    cond_final = measure_condition_number(parameters, prec)
    rho_init = (cond_init-1)/(cond_init+1)
    rho_final = (cond_final-1)/(cond_final+1)
    println("Initial condition number ", round(cond_init, digits=3), " Final condition number ", round(cond_final, digits=3))
    println("Initial ρ number ", round(rho_init, digits=3), " Final ρ number ", round(rho_final, digits=3))
    println("Initial N ", round(-1/log10(rho_init), RoundUp), " Final N ", round(-1/log10(rho_final), RoundUp))
    println("Optimization length ", length(history))

    return history, prec, parameters
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
# Contains S-norms
# ----------------------------------

"""
    Compute the maximum value of the tmp_strengths vector

    maximum_value(tmp_strengths<:AbstractFloat)
"""
function maximum_value(tmp_strengths::Vector{<:AbstractFloat})
    return maximum(tmp_strengths)
end

"""
    Compute the algebraic sum of the tmp_strengths vector

    algebraic_sum(tmp_strengths<:AbstractFloat)
"""
function algebraic_sum(tmp_strengths::Vector{<:AbstractFloat})
    if length(tmp_strengths) > 1
        val = tmp_strengths[1]
        for i = 2:length(tmp_strengths)
            val = val + tmp_strengths[i] - (val * tmp_strengths[i])
        end
        return val
    else
        return tmp_strengths[1]
    end
end

"""
    Compute the bounded sum of the tmp_strengths vector

    bounded_sum(tmp_strengths<:AbstractFloat)
"""
function bounded_sum(tmp_strengths::Vector{<:AbstractFloat})
    if length(tmp_strengths) > 1
        val = tmp_strengths[1]
        for i = 2:length(tmp_strengths)
            val = min(1.0, val + tmp_strengths[i])
        end
        return val
    else
        return tmp_strengths[1]
    end
end

"""
    Compute the drastic sum of the tmp_strengths vector

    drastic_sum(tmp_strengths<:AbstractFloat)
"""
function drastic_sum(tmp_strengths::Vector{<:AbstractFloat})
    if length(tmp_strengths) > 1
        val = tmp_strengths[1]
        for i = 2:length(tmp_strengths)
            if min(val, tmp_strengths[i]) == 0.0
                val = max(val, tmp_strengths[i])
            else
                val = 1.0
            end
        end
        return val
    else
        return tmp_strengths[1]
    end
end

"""
    Compute the Einstein sum of the tmp_strengths vector

    einstein_sum(tmp_strengths<:AbstractFloat)
"""
function einstein_sum(tmp_strengths::Vector{<:AbstractFloat})
    if length(tmp_strengths) > 1
        val = tmp_strengths[1]
        for i = 2:length(tmp_strengths)
            val = (val * tmp_strengths[i]) / (1.0 + (val * tmp_strengths[i]))
        end
        return val
    else
        return tmp_strengths[1]
    end
end

"""
    Compute the Hamacher sum of the tmp_strengths vector

    hamacher_sum(tmp_strengths<:AbstractFloat)
"""
function hamacher_sum(tmp_strengths::Vector{<:AbstractFloat})
    if length(tmp_strengths) > 1
        val = tmp_strengths[1]
        for i = 2:length(tmp_strengths)
            val = (val + tmp_strengths[i] - (2 * val * tmp_strengths[i])) / (1 - (val * tmp_strengths[i]))
        end
        return val
    else
        return tmp_strengths[1]
    end
end


text:
"""
Add constraints of the following form to a CoreModelCoupled and return a modified one.

The arguments are same as for in-place `add_coupling_constraints!`.
"""
function add_coupling_constraints(m::CoreModelCoupled, args...)
    new_lp = deepcopy(m)
    add_coupling_constraints!(new_lp, args...)
    return new_lp
end

"""
Add constraints to a plain `CoreModel` (converts it to the coupled type)
"""
add_coupling_constraints(m::CoreModel, args...) =
    add_coupling_constraints(convert(CoreModelCoupled, m), args...)

"""
In-place add coupling constraints in form
```
    cₗ ≤ C x ≤ cᵤ
```
"""
function add_coupling_constraints!(
    m::CoreModelCoupled,
    c::V,
    cl::AbstractFloat,
    cu::AbstractFloat,
) where {V<:VecType}
    return add_coupling_constraints!(
        m,
        sparse(reshape(c, (1, length(c)))),
        sparse([cl]),
        sparse([cu]),
    )
end

function add_coupling_constraints!(
    m::CoreModelCoupled,
    C::M,
    cl::V,
    cu::V,
) where {M<:MatType,V<:VecType}

    all([length(cu), length(cl)] .== size(C, 1)) ||
        throw(DimensionMismatch("mismatched numbers of constraints"))
    size(C, 2) == n_reactions(m) ||
        throw(DimensionMismatch("mismatched number of reactions"))

    m.C = vcat(m.C, sparse(C))
    m.cl = vcat(m.cl, sparse(cl))
    m.cu = vcat(m.cu, sparse(cu))
end

"""
Remove coupling constraints from the linear model and return the modified model.

Arguments are the same as for in-place version `remove_coupling_constraints!`.
"""
function remove_coupling_constraints(m::CoreModelCoupled, args...)
    new_model = deepcopy(m)
    remove_coupling_constraints!(new_model, args...)
    return new_model
end

"""
Removes a set of coupling constraints from a CoreModelCoupled in-place.
"""
function remove_coupling_constraints!(m::CoreModelCoupled, constraint::Int)
    remove_coupling_constraints!(m, [constraint])
end

function remove_coupling_constraints!(m::CoreModelCoupled, constraints::Vector{Int})
    to_be_kept = filter(e -> e ∉ constraints, 1:n_coupling_constraints(m))
    m.C = m.C[to_be_kept, :]
    m.cl = m.cl[to_be_kept]
    m.cu = m.cu[to_be_kept]
end

"""
Change the lower and/or upper bounds ('cl' and 'cu') for given coupling constraints
"""
function change_coupling_bounds!(
    model::CoreModelCoupled,
    constraints::Vector{Int};
    cl::V = Float64[],
    cu::V = Float64[],
) where {V<:VecType}
    found = [index ∈ 1:n_coupling_constraints(model) for index in constraints]
    red_constraints = constraints[found]

    length(red_constraints) == length(unique(red_constraints)) ||
        error("`constraints` appears to contain duplicates")
    if !isempty(cl)
        length(constraints) == length(cl) ||
            throw(DimensionMismatch("`constraints` size doesn't match with `cl`"))
        model.cl[red_constraints] = cl[found]
    end

    if !isempty(cu)
        length(constraints) == length(cu) ||
            throw(DimensionMismatch("`constraints` size doesn't match with `cu`"))
        model.cu[red_constraints] = cu[found]
    end
end


text:
# Basic library loading and API calling

import Base: get

#
# API call wrapper
#

const libcuda = Ref{Ptr{Void}}()
function load_library()
    try
        return (Libdl.dlopen(@static is_windows() ? "nvcuda.dll" : "libcuda"), "NVIDIA")
    end

    error("Could not load CUDA (or any compatible) library")
end

macro apicall(f, argtypes, args...)
    # Escape the tuple of arguments, making sure it is evaluated in caller scope
    # (there doesn't seem to be inline syntax like `$(esc(argtypes))` for this)
    esc_args = [esc(arg) for arg in args]

    blk = Expr(:block)

    if !isa(f, Expr) || f.head != :quote
        error("first argument to @apicall should be a symbol")
    end

    # Print the function name & arguments
    @static if TRACE
        push!(blk.args, :(trace($(sprint(Base.show_unquoted,f.args[1])*"("); line=false)))
        i=length(args)
        for arg in args
            i-=1
            sep = (i>0 ? ", " : "")

            # TODO: we should only do this if evaluating `arg` has no side effects
            push!(blk.args, :(trace(repr_indented($(esc(arg))), $sep;
                  prefix=$(sprint(Base.show_unquoted,arg))*"=", line=false)))
        end
        push!(blk.args, :(trace(""; prefix=") =", line=false)))
    end

    # Generate the actual call
    @gensym status
    push!(blk.args, quote
        # NOTE: version symbol lookup needs to happen at runtime,
        #       because the mapping is created during __init__
        api_f = resolve($f)
        $status = ccall(Libdl.dlsym(libcuda[], api_f), Cint,
                       $(esc(argtypes)), $(esc_args...))
    end)

    # Print the results
    @static if TRACE
        push!(blk.args, :(trace(CuError($status); prefix=" ")))
    end

    # Check the return code
    push!(blk.args, quote
        if $status != SUCCESS.code
            err = CuError($status)
            throw(err)
        end
    end)

    return blk
end

#
# Initialization
#

const libcuda_vendor = Ref{String}()

function __init_library__()
    (libcuda[], libcuda_vendor[]) = load_library()
    __init_versioning__()
    init()
end

function init()
    @apicall(:cuInit, (Cint,), 0)
end

function vendor()
    return libcuda_vendor[]
end


text:
using MATLAB

type linear_solver_MATLAB <: abstract_linear_system_solver
  _SparseMatrix::SparseMatrixCSC{Float64,Int64}

	# options
	sym::Symbol

  function linear_solver_MATLAB()
      # start matlab session
      @matlab begin
      end

      return new();
  end
end

function ls_factor!(solver::linear_solver_MATLAB, n::Int64, m::Int64)
			sparse_matrix = solver. _SparseMatrix
			@mput sparse_matrix

			@matlab begin
				mat_L, mat_D, mat_P, mat_S = ldl(sparse_matrix,0.5);
			end

			@mget mat_D

			pos_eigs, neg_eigs, zero_eigs = eigsigns_of_B(mat_D, 0.0);
			return inertia_status(pos_eigs, neg_eigs, zero_eigs, n, m)
end

function ls_solve(solver::linear_solver_MATLAB, my_rhs::AbstractArray)
      @mput my_rhs

      @matlab begin
        mat_solution = mat_S * ( mat_P* (mat_L'\ (mat_D\ (mat_L\ (mat_P' * ( mat_S * my_rhs ) ) ) ) ) );
      end

      @mget mat_solution

      return mat_solution;
end

max_forks_repo_name: ohinder/Everest


text:
# """
# 	$(SIGNATURES)

# Common fields in `ModelSwitches`.
# """
# mutable struct CommonModelSwitches
#     objId :: ObjectId
# end

get_object_id(s :: ModelSwitches) = s.objId :: ObjectId;

# set_object_id!(s :: ModelSwitches)

# -----------


text:
@testset "Initial conditions" begin
    @testset "Fluid Sine Wave" begin
        h = ones(10,10)
        Swalbe.sinewave1D!(h, 2.0, 1, 1, 1)
        @test maximum(h) ≈ 2*(1 + 1) atol= 0.5
        @test minimum(h) ≈ 2*(1 - 1) atol= 0.5
    end

    @testset "Random interface" begin
        h = ones(10,10)
        Swalbe.randinterface!(h, 2.0, 0.1)
        @test maximum(h) <= 3.0
        @test minimum(h) >= 1.0
    end

    @testset "Single droplet" begin
        rad = 50
        θ = 1/3
        center = (50,50)
        height = Swalbe.singledroplet(ones(100,100), rad, θ, center)
        @test isa(height, Array{Float64,2})
        @test size(height) == (100, 100)
        @test findmax(height)[1] == rad * (1 - cospi(θ))
        @test findmax(height)[2] == CartesianIndex(center[1],center[2])
        # 1D case
        rad = 50
        θ = 1/3
        center = 100
        height = Swalbe.singledroplet(ones(200), rad, θ, center)
        @test isa(height, Vector{Float64})
        @test size(height) == (200,)
        @test findmax(height)[1] == rad * (1 - cospi(θ))
        @test findmax(height)[2] == center
    end

    @testset "Restart height" begin
        h1 = rand(10,10)
        h2 = rand(10,10)
        f = zeros(2, 100)
        f = Dict()
        f["h_1"] = vec(h1)
        f["h_2"] = vec(h2)
        bson("file.bson", f)
        save("file.jld2", f)
        h = Swalbe.restart_from_height("file.bson", kind="bson", timestep=0, size=(10,10))
        @test all(h2 .== h)
        h = Swalbe.restart_from_height("file.bson", kind="bson", timestep=1, size=(10,10))
        @test all(h1 .== h)
        h = Swalbe.restart_from_height("file.jld2", size=(10,10))
        @test all(h2 .== h)
        h = Swalbe.restart_from_height("file.jld2", kind="jld2", timestep=1, size=(10,10))
        @test all(h1 .== h)
        rm("file.bson")
        rm("file.jld2")
    end
end
@testset "Contact angle patterns" begin
    @testset "Default elliptical patch" begin
        theta, P = Swalbe.ellipsepattern(ones(100,100), 1/9) # with kw: center=(Int(lx/2), Int(ly/2)),δₐ=1/36, a=10, b=5
        @test size(theta) == (100, 100)
        @test theta[50,50] == 1/9 + 1/36
        @test theta[1,1] == 1/9
    end
    @testset "Off center elliptical patch" begin
        theta, P = Swalbe.ellipsepattern(ones(100,100), 1/9, center=(20,20)) # with kw: center=(Int(lx/2), Int(ly/2)),δₐ=1/36, a=10, b=5
        @test size(theta) == (100, 100)
        @test theta[20,20] == 1/9 + 1/36
        @test theta[100,100] == 1/9
    end
    @testset "Negative contrast elliptical patch" begin
        theta, P = Swalbe.ellipsepattern(ones(100,100), 1/9, δₐ=-1/36) # with kw: center=(Int(lx/2), Int(ly/2)),δₐ=1/36, a=10, b=5
        @test size(theta) == (100, 100)
        @test theta[50,50] == 1/9 - 1/36
        @test theta[100,100] == 1/9
    end
    @testset "Larger elliptical patch" begin
        theta, P = Swalbe.ellipsepattern(ones(100,100), 1/9, a=10, b=20) # with kw: center=(Int(lx/2), Int(ly/2)),δₐ=1/36, a=10, b=5
        @test size(theta) == (100, 100)
        @test theta[40,50] == 1/9 + 1/36
        @test theta[70,50] == 1/9 + 1/36
        @test theta[100,100] == 1/9
    end
    @testset "Polygon patch" begin
        lx, ly = 100, 100
        @testset "Equilateral Triangle" begin
            θₙ, P = Swalbe.trianglepattern(ones(100, 100), 1/9)
            @test all(θₙ .> 0.0)
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(lx/2), Int(ly/2)] == 1/9 + 1/36
            θₙ, P = Swalbe.trianglepattern(ones(100, 100), 1/9, center=(lx/3, 2*ly/3))
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(round(lx/3)), Int(round(2*ly/3))] == 1/9 # Because the side became to large
            θₙ, P = Swalbe.trianglepattern(ones(100, 100), 1/9, center=(0, 2*ly/3))
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(round(lx/3)), Int(round(2*ly/3))] == 1/9 # Because the side became to large
            @test θₙ[60, 1] == 1/9 + 1/36 # If too large stick to origin!
            θₙ, P = Swalbe.trianglepattern(ones(100, 100), 1/9, δₐ = -1/36)
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(lx/2), Int(ly/2)] == 1/9 - 1/36
            θₙ, P = Swalbe.trianglepattern(ones(100, 100), 1/9, side = 30) # Smaller triangle
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(lx/2), Int(ly/2)] == 1/9 + 1/36
            @test θₙ[50, 60] == 1/9 + 1/36
        end

        @testset "Box pattern" begin
            θₙ, P = Swalbe.boxpattern(ones(100, 100), 1/9)
            @test all(θₙ .> 0.0)
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(lx/2), Int(ly/2)] == 1/9 + 1/36
            θₙ, P = Swalbe.boxpattern(ones(100, 100), 1/9, center=(lx/3, 2*ly/3))
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(round(2*ly/3)), Int(round(lx/3))] == 1/9 + 1/36
            θₙ, P = Swalbe.boxpattern(ones(100, 100), 1/9, center=(0, 2*ly/3))
            @test θₙ[20,1] == 1/9 # the default value
            @test θₙ[60,1] == 1/9 + 1/36
            θₙ, P = Swalbe.boxpattern(ones(100, 100), 1/9, δₐ = -1/36)
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(lx/2), Int(ly/2)] == 1/9 - 1/36
            θₙ, P = Swalbe.boxpattern(ones(100, 100), 1/9, side = 30) # Smaller triangle
            @test θₙ[1,1] == 1/9 # the default value
            @test θₙ[Int(lx/2), Int(ly/2)] == 1/9 + 1/36
            @test θₙ[Int(lx/2) + Int(30/2), Int(ly/2) + Int(30/2)] == 1/9 + 1/36
        end
    end
end


text:
# Test maxproduct algorithm

@testset "MaxProduct" begin
    import SPNetworks.MAP: maxproduct!

    @testset "Discrete SPNs" begin

        @testset "Simple SPN" begin
            SPN = SumProductNetwork(
                [
                    SumNode([2,3,4],[0.2,0.5,0.3]),        # 1
                    ProductNode([5,7]),                    # 2
                    ProductNode([5,8]),                    # 3
                    ProductNode([6,8]),                    # 4
                    CategoricalDistribution(1, [0.6,0.4]), # 5
                    CategoricalDistribution(1, [0.1,0.9]), # 6
                    CategoricalDistribution(2, [0.3,0.7]), # 7
                    CategoricalDistribution(2, [0.8,0.2])  # 8
                ]
            )
            values = Vector{Float64}(undef,length(SPN))
            tree = Vector{UInt}(undef,length(SPN))
            mp = maxproduct!(values,tree,SPN,Set([1,2]),[1,1])
            expected = log.([0.24,0.42,0.48,0.72,0.6,0.9,0.7,0.8])
            for (i,v) in enumerate(values)
                @test v ≈ expected[i]
            end
            @test tree[1] == 3
            evidence = [0.0,0.0]
            mp = maxproduct!(evidence,SPN,Set([1,2]))
            @test logpdf(SPN,evidence) ≈ log(0.3)
            evidence = [0.0,NaN]
            mp = maxproduct!(evidence,SPN,Set([1]))
            @test logpdf(SPN,evidence) ≈ log(0.45)
        end

        @testset "Dichotomized version" begin
            SPN = SumProductNetwork([
                SumNode([3,2],[0.2,0.8]),
                SumNode([4,5],[0.625,0.375]),
                ProductNode([6,8]),
                ProductNode([6,9]),
                ProductNode([7,9]),
                CategoricalDistribution(1,[0.6,0.4]),
                CategoricalDistribution(1,[0.1,0.9]),
                CategoricalDistribution(2,[0.3,0.7]),
                CategoricalDistribution(2,[0.8,0.2])
            ]
            )
            evidence = [0.0,0.0]
            mp = maxproduct!(evidence,SPN,Set([1,2]))
            @test logpdf(SPN,evidence) ≈ log(0.3)
        end
        @testset "Selective SPN" begin
            SPN = SumProductNetwork(
                [
                    SumNode([2,3],[0.4,0.6]),              # 1
                    ProductNode([4,5,6]),                  # 2
                    ProductNode([7,8,9]),                  # 3
                    IndicatorFunction(1,2.0),  # 4
                    CategoricalDistribution(2,[0.3,0.7]),  # 5
                    CategoricalDistribution(3,[0.4,0.6]),  # 6
                    CategoricalDistribution(2,[0.8,0.2]),  # 7
                    CategoricalDistribution(3,[0.9,0.1]),  # 8
                    IndicatorFunction(1,1.0) # 9
                ]
            )
            lmaxv = -Inf
            for a=1:2,b=1:2,c=1:2
                v = logpdf(SPN,[a,b,c])
                lmaxv = v > lmaxv ? v : lmaxv
            end
            evidence = [0.0,0.0,0.0]
            mp = maxproduct!(evidence,SPN,Set([1,2,3]))
            @test logpdf(SPN,evidence) ≈ lmaxv
        end
        @testset "Breast-Cancer SPN" begin
            SPN = SumProductNetwork(normpath("$(@__DIR__)/../assets/breast-cancer.spn"), offset=1)
            query = Set([9,1,5])
            evidence = [0,2,NaN,NaN,0,NaN,2,NaN,0,1]
            maxproduct!(evidence,SPN,query)
            @test SPN(evidence) ≈ 0.01125 atol=1e-5
            # println("MAP: ", join(map( q -> "X$q = $(evidence[q])", collect(query)), ", ") )
            0.00194
            query = Set([2,5,3])
            evidence = [6,0,0,NaN,0,2,NaN,NaN,NaN,1]
            maxproduct!(evidence,SPN,query)
            @test SPN(evidence) ≈ 0.00194  atol=1e-5
            # println("MAP: ", join(map( q -> "X$q = $(evidence[q])", collect(query)), ", ") )
        end
    end
end # END of test set


text:
# Model.
@sugar function hmm(N::Int)
    z₀ ~ Normal(0.0, 1.0)
    z = z₀
    for i in 1 : N
        rand(:z => i, Normal(z, 1.0))
        rand(:x => i, Normal(z, 1.0))
    end
end

# Data.
data_len = 100
x = map(x) do k
    3.0 + randn()
end
obs = target(map(1 : data_len) do i
                 (:x => i, ) => x[i]
             end)

@testset "Inference library smoke test 1 - Bayesian linear regression" begin

    # Importance sampling.
    is_test = () -> begin
        n_samples = 50000
        ps, lnw = importance_sampling(obs, n_samples, bayesian_linear_regression, (x, ))
        zipped = zip(ps.calls, lnw)

        est_σ = sum(map(zipped) do (cl, w)
                        (cl[:σ]) * exp(w)
                    end)

        est_β = sum(map(zipped) do (cl, w)
                        (cl[:β]) * exp(w)
                    end)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # MH random walk kernel.
    mh_test = () -> begin
        n_iters = 50000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = mh(target([(:σ, ), (:β, )]), cl)
            i % 30 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # MH random walk kernel with proposal.
    proposal = (cl, y) -> begin
        σ = rand(:σ, Normal(6.0, 1.0))
        β = rand(:β, Normal(mean(y) / length(y), 3.0))
    end

    mh_test_with_proposal = () -> begin
        n_iters = 10000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = mh(cl, proposal, (y, ))
            i % 30 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # ADVI.
    surrogate = () -> begin
        μ₁ = learnable(:μ₁)
        μ₂ = learnable(:μ₂)
        σ₂ = learnable(:σ₂)
        σ = rand(:σ, Normal(1.0 + μ₁, 0.3))
        β = rand(:β, Normal(μ₂, σ₂ ^ 2))
    end

    advi_test = () -> begin
        ps = learnables([(:μ₁, ) => 5.0,
                         (:μ₂, ) => 1.0,
                         (:σ₂, ) => 1.0])
        @time ps, elbows, _ = advi(obs, ps,
                                   surrogate, (),
                                   bayesian_linear_regression, (x, );
                                   opt = ADAM(0.05, (0.9, 0.8)),
                                   iters = 1000)
    end

    # ADGV.
    surrogate = () -> begin
        μ₁ = learnable(:μ₁)
        μ₂ = learnable(:μ₂)
        σ₂ = learnable(:σ₂)
        σ = rand(:σ, Normal(1.0 + μ₁, 0.3))
        β = rand(:β, Normal(μ₂, σ₂ ^ 2))
    end

    adgv_test = () -> begin
        ps = learnables([(:μ₁, ) => 5.0,
                         (:μ₂, ) => 1.0,
                         (:σ₂, ) => 1.0])
        ps, elbows, _ = adgv(obs, ps, 10,
                             surrogate, (),
                             bayesian_linear_regression, (x, );
                             opt = ADAM(0.05, (0.9, 0.8)),
                             iters = 1000)
    end

    # HMC kernel.
    hmc_test = () -> begin
        n_iters = 1000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = hmc(target([(:σ, ), (:β, )]), cl)
            i % 10 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # Combination kernel.
    combo_kernel_test = () -> begin
        n_iters = 1000
        ret, cl = generate(obs, bayesian_linear_regression, x)
        calls = []
        for i in 1 : n_iters
            cl, _ = mh(cl, proposal, (y, ))
            cl, _ = hmc(target([(:σ, ), (:β, )]), cl)
            i % 10 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end) / length(calls)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end) / length(calls)
        @test isapprox(0.0, est_β - 3.0, atol=1e-1)
    end

    # Boomerang kernel.
    boomerang_test = () -> begin
        n_iters = 500
        ret, cl = generate(obs, bayesian_linear_regression, x)
        sel = get_target(cl)
        sel_values, _ = get_choice_gradients(sel, cl, 1.0)
        d = length(array(sel_values, Float64))
        flow = Boomerang(sparse(I, d, d), zeros(d), 1.0)
        θ = rand(MvNormal(d, 1.0))
        calls = []
        for i in 1 : n_iters
            cl, _ = pdmk(target([(:σ, ), (:β, )]), cl, flow, θ)
            i % 10 == 0 && begin
                push!(calls, cl)
            end
        end

        est_σ = sum(map(calls) do cl
                        (cl[:σ])
                    end)

        est_β = sum(map(calls) do cl
                        (cl[:β])
                    end)
    end

    @time is_test()
    @time mh_test()
    @time mh_test_with_proposal()
    @time hmc_test()
    @time combo_kernel_test()
    #@time boomerang_test()
    #@time advi_test()
    #@time adgv_test()

end


text:
### This is Symata level IO. The code for formatting and printing Mxpr, etc. is in output.jl

## NOTE: symprint and symprintln are defined in wrapout.jl

# function Symata_module_path()
#     joinpath(dirname(@__FILE__), "..")
# end

### Println

@sjdoc Println """
    Println(expr1, expr2, ...)

print the expressions and a newline.
"""

apprules(mx::Mxpr{:Println}) = (symprintln(margs(mx)...) ; Null)

### Print

@sjdoc Print """
    Print(expr1,expr2,...)

prints the expressions.
"""

apprules(mx::Mxpr{:Print}) = (symprint(margs(mx)...); Null)

# TODO. Read more than one expression from a line.
# TODO. Use a new exception type
# TODO. use EvaluateJuliaSyntaxSimple() I think it is simpler than the try catch stuff below

"""
    read_Symata_file(f::AbstractString, test::Symata_Test = Symata_NullTest() )

Read and evaluate Symata expressions from file `f`.
Return the last returned value. If `test` is given, interpret the contents
of the file as Symata tests.
"""
function read_Symata_file(f::AbstractString, test::Symata_Test = Symata_NullTest() )
    eline = ""
    local success_flag_return
    reading_test::Bool = false
    incomplete_flag::Bool = false
    incomplete_message = ""
    local incomplete_line_number = 0
    local first_line_number = 0  # line number at which expression being read started
    for (line_number, line) = enumerate(eachline(f))
        pline = sjpreprocess_string(line)
        if typeof(test) != Symata_NullTest && length(pline) > 1 && pline[1:2] == "T "
            if isa(test, SymataPlainPrintTest) println(pline) end
            reading_test = true
            if length(pline) > 2
                eline = eline * pline[3:end]
            end
        else
            eline = eline * pline   # is it worth using takebuf here ?
        end
        expr =
            try
                symata_syntax_deprecation_warnings(false) do
                    Meta.parse(eline)
                end
            catch e
                error("Syntax error in file $f, line $line_number: '", chomp(pline), "'")
                throw(e)
            end
        if typeof(expr) ==  Expr && expr.head == :incomplete
            incomplete_flag = true
            incomplete_message = expr.args[1]
            incomplete_line_number = line_number
            continue
        end
        success_flag_return =
            try
                incomplete_flag = false
                success_flag = Symata.symataevaluate(expr, EvaluateJuliaSyntaxSimple())
                if typeof(test) != Symata_NullTest && reading_test
                    reading_test = false
                    record_SymataTest(test, f, line_number, success_flag)
                end
                first_line_number = line_number + 1  # The first line number of the next expression
                success_flag
            catch e
                @warn("Symata: Error reading file $f, line $line_number\n")
                rethrow(e)
            finally
                eline = ""
            end
    end
    if incomplete_flag
        @isdefined(line_number) ||
            error(incomplete_message,
                  ". Incomplete expression starting near line number $first_line_number in file $f.")
        error(incomplete_message, " in file $f line $line_number.")
    end
    return success_flag_return
end

### Get

@mkapprule Get

@sjdoc Get """
    Get(filename)

read and evaluate Symata expressions from file `filename`.

!!! note
    Try putting empty lines between expressions if you get errors on reading.
"""

do_Get(mx::Mxpr{:Get}, fname::T) where {T<:AbstractString} =  read_Symata_file(fname)

### ReadString

@mkapprule ReadString  :nargs => 1

@sjdoc ReadString  """
    ReadString(filename)

reads `filename` and returns the contents as a `String`.
"""

do_ReadString(mx::Mxpr{:ReadString}, fname::T) where {T<:AbstractString} = readstring(fname)

# Write Symata expression as strings that can be used as input to define
# objects and properties.

# string to set attributes of sym
function attributes_set_string(sym::SymString)
    a = get_attributes(sym)
    a = map(string,a)
    sym = string(sym)
    if length(a) == 1
        "SetAttributes(" * sym * ", " * a[1] * ")"
    elseif length(a) > 1
        "SetAttributes(" * sym * ", [" * join(a, ", ")  * "])"
    else
        ""
    end
end

### Definition

@mkapprule Definition

@sjdoc Definition """
    Definition(sym)

prints definitions associated with `Symbol` `sym`.
"""

function maybeprint(io::IO, s)
    if length(s) > 0 symprintln(io, s)
        return true
    end
    false
end

## FIXME: The correct way is not this way. Instead  Definition(f) does not evaluate to anything.
## rather displaying the form Definition(f) prints the definitions
function do_Definition(mx::Mxpr{:Definition}, sym::SymString)
    write_definition(stdout, sym)
    Null
end

write_definition(sym) = write_definition(stdout, sym)

function write_definition(io::IO, sym)
    gotdef = write_definition_no_attributes(io,sym)
    gotdef = maybeprint(io,attributes_set_string(sym)) || gotdef
    if gotdef
        println(io)
    end
    return gotdef
end

write_definition_no_attributes(sym) = write_definition_no_attributes(stdout, sym)
function write_definition_no_attributes(io::IO, sym)
    gotdef = maybeprint(io,getdefinition(Symbol(sym)))
    dvdefs = jlistdownvaluedefs(sym)
    if length(dvdefs) > 0
        gotdef = true
        for def in dvdefs
            symprintln(io,def)
        end
    end
    uvdefs = jlistupvaluedefs(sym)
    if length(uvdefs) > 0
        gotdef = true
        for def in uvdefs
            symprintln(io,def)
        end
    end
#    gotdef = maybeprint(io,attributes_set_string(sym)) || gotdef
    if gotdef
        println(io)
    end
    return gotdef
end

### Save

@mkapprule Save

@sjdoc Save """
    Save(filename, obj)

append definitions associated with `obj` to file `filename`.

    Save(filename, list)

append definitions associated with the objects in `list` to `filename`.

`Save` writes defnitions given by `Definition`.
"""

do_Save(mx::Mxpr{:Save}, fname::T, sym::SJSym) where {T<:AbstractString} = save_definition(fname, [sym])
do_Save(mx::Mxpr{:Save}, fname::T, syms::Mxpr{:List}) where {T<:AbstractString} = save_definition(fname, margs(syms))

function save_definition(fname, syms)
    io = open(fname, "a")  # check stuff !
    for sym in syms
        write_definition(io,sym)
    end
    close(io)
end

@sjseealso_group(Get, Save, Definition)

### Open

# we should fold this behavior into mkapprule
@mkapprule Open  :nodefault => true

@sjdoc Open """
    Open(args...)

opens a file, runs a command, etc. The arguments are passed to Julia `open`.
The Julia methods that take array arguments are not implemented.
"""

do_Open(mx::Mxpr{:Open}, args...) = open(args...)

### Close

@mkapprule Close

@sjdoc Close """
    Close(str)

closes the `IO` stream `str`.
"""

do_Close(mx::Mxpr{:Close}, str::T) where {T<:IO} = (close(str); Null)

### STDOUT, STDERR.  These two symbols are set it client_symata.jl after the REPL is created.

@sjdoc STDOUT """
    STDOUT

the standard output stream.
"""

@sjdoc STDERR """
    STDERR

the standard output stream.
"""

@sjdoc DevNull """
    DevNull

is bound to Julia output stream `DevNull`.
"""

@mkapprule Out

@sjdoc Out """
    Out(n)

return the `n`th line of output. The value returned is evaluated in the current environment.

`O`

return the last output.

`OO`

returns the second to last output.

`OO..O`

`k` times gives the `k`th previous output, for `k <= 10`.

`Clear(Out)` clears all the saved output.

!!! note
    `O`, `OO`, etc. do not work in `Jupyter`.
"""

@sjseealso(Out, Clear)

function do_Out(mx::Mxpr{:Out}, n::Integer)
    if is_isymata_mode()
        doeval(Main.IJulia.Out[n])
    else
        doeval(get_output_by_line(n))
    end
end

@doap Out() = doeval(get_output_by_line(-1))

do_Out(mx::Mxpr{:Out}, x) = :Null

### In

@mkapprule In :nargs => 1

@sjdoc In """
    In(n)

returns the `n`th input cell. `In` only works in `Jupyter/IPython`.
"""
@doap function In(n::Integer)
    if is_isymata_mode()
        doeval(eval(Expr(:macrocall, Symbol("@symfull"), Meta.parse(Main.IJulia.In[n]))))
    else
        :Null
    end
end

# We should probably explicitly return Null
# from functions. But,this might catch all
# of the nothings from causing printing
# ...No. this did not help
apprules(expr::Nothing) = :Null

### TempName

@mkapprule TempName :nargs => 0

@sjdoc TempName """
    TempName()

generate a unique temporary file path.
"""
@doap TempName() = tempname()

### DeleteFile

@mkapprule DeleteFile :nargs => 1

@sjdoc DeleteFile """
    DelteFile(filename)

delete `filename`.
"""
@doap DeleteFile(file::AbstractString) = rm(file)

### Format

## Note: This is *not* Mma's format

@sjdoc Format """
    Format(n, options)

Format `n` using the Julia `format` function from `Formatting.jl`.
Optional arguments to `format` are given in Symata format.
For example `precision=3` becomes `Precision => 3`.

Options are `Precision`, `Commas`, `Stripzeros`, `Conversion`,
`Width`, `Zeropadding`,  `Alternative`, `Mixedfraction`, `Tryden`,
`Parens`,`Autoscale`, `Suffix`.
"""

@mkapprule Format :nodefault => true

Formatting.format(ws::SymataIO.AbstractWO; kws...) = Formatting.format(ws.x; kws...)

function do_Format(mx, args...)
    (nargs,kws) = separate_rules1(args...; keylowercase=true)
    Formatting.format(map(wrapout,nargs)...; kws...)
end

### FullForm

@sjdoc FullForm """
    FullForm(expr)

print the internal representation of `expr` and all sub-expressions as
`Head(arg1,arg2,...)`. Normal output may use infix notation instead. But, expressions
may always be entered in `FullForm`.
"""

@sjexamp( FullForm,
          ("Clear(a,b)",""),
          ("a+b","a+b"),
          ("FullForm(a+b)","Plus(a,b)"))

@sjexamp( FullForm,
          ("Sqrt(2)", "Power(2,1//2)", "Internally, square roots are represented as powers."),
          ("Sqrt(2)[2]", "1//2", "The second part of the expression is '1//2'"))

# FullForm is implemented in Output.jl

### Echo

@mkapprule Echo :nargs => 1:3

@doap Echo(x) = (symprintln(x); x)
@doap Echo(x,label::String) = (symprintln(label," ", x); x)

@doap function Echo(x,label::String, f)
    res = doeval(mxpr(f,x))
    symprintln(label," ", res)
    x
end

## Obsolete

# We are using read_Symata_file instead of this.
# We should probably delete this at some point.
# Read Symata expressions from a string and evaluate them, one by one.
function Symata_eval_string(s)
    s = sjpreprocess_string(s)
    i = 1
    local sjretval
    while !done(s,i)
        symata_syntax_deprecation_warnings(false) do
            expr, i = Meta.parse(s,i)
        end
        sjretval =
            try
                Symata.symataevaluate(expr)
            catch e
                symprintln("Reading file: got error ", e)
            end
    end
    sjretval
end

## Obsolete
Symata_eval_file(fname) = fname |> readstring |> Symata_eval_string


text:
# Autogenerated wrapper script for libvorbis_jll for i686-linux-musl
export libvorbis, libvorbisenc, libvorbisfile

using Ogg_jll
JLLWrappers.@generate_wrapper_header("libvorbis")
JLLWrappers.@declare_library_product(libvorbis, "libvorbis.so.0")
JLLWrappers.@declare_library_product(libvorbisenc, "libvorbisenc.so.2")
JLLWrappers.@declare_library_product(libvorbisfile, "libvorbisfile.so.3")
function __init__()
    JLLWrappers.@generate_init_header(Ogg_jll)
    JLLWrappers.@init_library_product(
        libvorbis,
        "lib/libvorbis.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libvorbisenc,
        "lib/libvorbisenc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libvorbisfile,
        "lib/libvorbisfile.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "intrinsics" begin

############################################################################################

@testset "math" begin
    buf = CuArray{Float32}(1)

    @eval function kernel_math_log10(a, i)
        a[1] = CUDAnative.log10(i)
        return nothing
    end

    @cuda (1, 1) kernel_math_log10(buf, Float32(100))
    val = Array(buf)
    @test val[1] ≈ 2.0
end

############################################################################################

@testset "I/O" begin

@testset "printing" begin
    _, out = @grab_output @on_device @cuprintf("")
    @test out == ""

    _, out = @grab_output @on_device @cuprintf("Testing...\n")
    @test out == "Testing...\n"

    _, out = @grab_output @on_device @cuprintf("Testing %ld...\n", Int64(42))
    @test out == "Testing 42...\n"

    _, out = @grab_output @on_device @cuprintf("Testing %u %u...\n", blockIdx().x, threadIdx().x)
    @test out == "Testing 1 1...\n"

    _, out = @grab_output @on_device begin
        @cuprintf("foo")
        @cuprintf("bar\n")
    end
    @test out == "foobar\n"
end

end

############################################################################################

# a composite type to test for more complex element types
@eval immutable RGB{T}
    r::T
    g::T
    b::T
end

@testset "shared memory" begin

n = 1024
types = [Int32, Int64, Float32, Float64]

@testset "constructors" begin
    # static
    @on_device @cuStaticSharedMem(Float32, 1)
    @on_device @cuStaticSharedMem(Float32, (1,2))
    @on_device @cuStaticSharedMem(Tuple{Float32, Float32}, 1)
    @on_device @cuStaticSharedMem(Tuple{Float32, Float32}, (1,2))
    @on_device @cuStaticSharedMem(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device @cuStaticSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2))

    # dynamic
    @on_device @cuDynamicSharedMem(Float32, 1)
    @on_device @cuDynamicSharedMem(Float32, (1, 2))
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, 1)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, (1,2))
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, 1)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2))

    # dynamic with offset
    @on_device @cuDynamicSharedMem(Float32, 1, 8)
    @on_device @cuDynamicSharedMem(Float32, (1,2), 8)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, 1, 8)
    @on_device @cuDynamicSharedMem(Tuple{Float32, Float32}, (1,2), 8)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, 1, 8)
    @on_device @cuDynamicSharedMem(Tuple{RGB{Float32}, UInt32}, (1,2), 8)
end

@testset "dynamic shmem" begin

@testset "statically typed" begin
    @eval function kernel_shmem_dynamic_typed(d, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(Float32, n)
        s[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return nothing
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    @cuda (1, n, n*sizeof(Float32)) kernel_shmem_dynamic_typed(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" begin
    @eval function kernel_shmem_dynamic_typevar{T}(d::CuDeviceArray{T}, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(T, n)
        s[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return nothing
    end

    for T in types
        a = rand(T, n)
        d_a = CuArray(a)

        @cuda (1, n, n*sizeof(T)) kernel_shmem_dynamic_typevar(d_a, n)
        @test reverse(a) == Array(d_a)
    end
end

end

@testset "static shmem" begin

@testset "statically typed" begin
    @eval function kernel_shmem_static_typed(d, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuStaticSharedMem(Float32, 1024)
        s2 = @cuStaticSharedMem(Float32, 1024)  # catch aliasing

        s[t] = d[t]
        s2[t] = 2*d[t]
        sync_threads()
        d[t] = s[tr]

        return nothing
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    @cuda (1, n) kernel_shmem_static_typed(d_a, n)
    @test reverse(a) == Array(d_a)
end

@testset "parametrically typed" begin
    @eval function kernel_shmem_static_typevar{T}(d::CuDeviceArray{T}, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuStaticSharedMem(T, 1024)
        s2 = @cuStaticSharedMem(T, 1024)  # catch aliasing

        s[t] = d[t]
        s2[t] = d[t]
        sync_threads()
        d[t] = s[tr]

        return nothing
    end

    for T in types
        a = rand(T, n)
        d_a = CuArray(a)

        @cuda (1, n) kernel_shmem_static_typevar(d_a, n)
        @test reverse(a) == Array(d_a)
    end
end

end

@testset "dynamic shmem consisting of multiple arrays" begin

# common use case 1: dynamic shmem consists of multiple homogeneous arrays
#                    -> split using `view`
@testset "homogeneous" begin
    @eval function kernel_shmem_dynamic_multi_homogeneous(a, b, n)
        t = threadIdx().x
        tr = n-t+1

        s = @cuDynamicSharedMem(eltype(a), 2*n)

        sa = view(s, 1:n)
        sa[t] = a[t]
        sync_threads()
        a[t] = sa[tr]

        sb = view(s, n+1:2*n)
        sb[t] = b[t]
        sync_threads()
        b[t] = sb[tr]

        return nothing
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    b = rand(Float32, n)
    d_b = CuArray(b)

    @cuda (1, n, 2*n*sizeof(Float32)) kernel_shmem_dynamic_multi_homogeneous(d_a, d_b, n)
    @test reverse(a) == Array(d_a)
    @test reverse(b) == Array(d_b)
end

# common use case 2: dynamic shmem consists of multiple heterogeneous arrays
#                    -> construct using pointer offset
@testset "heterogeneous" begin
    @eval function kernel_shmem_dynamic_multi_heterogeneous(a, b, n)
        t = threadIdx().x
        tr = n-t+1

        sa = @cuDynamicSharedMem(eltype(a), n)
        sa[t] = a[t]
        sync_threads()
        a[t] = sa[tr]

        sb = @cuDynamicSharedMem(eltype(b), n, n*sizeof(eltype(a)))
        sb[t] = b[t]
        sync_threads()
        b[t] = sb[tr]

        return nothing
    end

    a = rand(Float32, n)
    d_a = CuArray(a)

    b = rand(Int64, n)
    d_b = CuArray(b)

    @cuda (1, n, n*sizeof(Float32) + n*sizeof(Int64)) kernel_shmem_dynamic_multi_heterogeneous(d_a, d_b, n)
    @test reverse(a) == Array(d_a)
    @test reverse(b) == Array(d_b)
end

end

end

############################################################################################

@testset "data movement and conversion" begin

if capability(dev) >= v"3.0"
@testset "shuffle" begin

n = 14
types = [Int32, Int64, Float32, Float64]

@testset "down" begin
    @eval function kernel_shuffle_down{T}(d::CuDeviceArray{T}, n)
        t = threadIdx().x
        if t <= n
            d[t] += shfl_down(d[t], n÷2)
        end
        return nothing
    end

    for T in types
        a = T[i for i in 1:n]
        d_a = CuArray(a)

        threads = nearest_warpsize(dev, n)
        @cuda (1, threads) kernel_shuffle_down(d_a, n)

        a[1:n÷2] += a[n÷2+1:end]
        @test a == Array(d_a)
    end
end

end
end

end

############################################################################################

@testset "parallel synchronization and communication" begin

@testset "voting" begin

@testset "ballot" begin
    d_a = CuArray{UInt32}(1)

    @eval function kernel_vote_ballot(a, i)
        vote = vote_ballot(threadIdx().x == i)
        if threadIdx().x == 1
            a[1] = vote
        end

        return nothing
    end

    len = 4
    for i in 1:len
        @cuda (1,len) kernel_vote_ballot(d_a, i)
        @test Array(d_a) == [2^(i-1)]
    end
end

@testset "any" begin
    d_a = CuArray{UInt32}(1)

    @eval function kernel_vote_any(a, i)
        vote = vote_any(threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end

        return nothing
    end

    @cuda (1,2) kernel_vote_any(d_a, 1)
    @test Array(d_a) == [1]

    @cuda (1,2) kernel_vote_any(d_a, 2)
    @test Array(d_a) == [1]

    @cuda (1,2) kernel_vote_any(d_a, 3)
    @test Array(d_a) == [0]
end

@testset "all" begin
    d_a = CuArray{UInt32}(1)

    @eval function kernel_vote_all(a, i)
        vote = vote_all(threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end

        return nothing
    end

    @cuda (1,2) kernel_vote_all(d_a, 1)
    @test Array(d_a) == [1]

    @cuda (1,2) kernel_vote_all(d_a, 2)
    @test Array(d_a) == [0]

    @cuda (1,2) kernel_vote_all(d_a, 3)
    @test Array(d_a) == [0]
end

end

end

############################################################################################

end


text:
# Runs all unit tests
using LinearLeastSquares
using Base.Test

info("JULIA VERSION: $(VERSION)") # important as Julia versions are volatile
include("helpers.jl")
include("add_subtract.jl")
include("multiply_divide.jl")
include("stack.jl")
include("diag.jl")
include("repmat.jl")
include("index.jl")
include("reshape.jl")
include("transpose.jl")
include("sum_squares.jl")
include("compound.jl")
include("constraints.jl")
include("problems.jl")