text:
using PowerSimulationsDynamics
using Sundials

"""
Case 24:
This case study a 15-state grid following inverter against an infinite bus located at bus 1, with the inverter located at bus 2.
The perturbation increase the reference power (analogy for mechanical power) from 0.5 to 0.7.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test24.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

# time span
tspan = (0.0, 2.0);

# PSCAD benchmark data
csv_file = joinpath(TEST_FILES_DIR, "benchmarks/pscad/Test24/Test24_p.csv")
t_offset = 9.0

#Define Fault using Callbacks
Pref_change = ControlReferenceChange(1.0, case_inv, :P_ref, 0.7)

@testset "Test 24 Grid Following Inverter ResidualModel" begin
    path = (joinpath(pwd(), "test-24"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            omib_sys, # system
            path,
            tspan,
            Pref_change,
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test24_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff_val[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test24_eigvals) < 1e-3

        #Solve problem in equilibrium
        @test execute!(sim, Sundials.IDA(), dtmax = 0.001, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain frequency data
        series = get_state_series(results, ("generator-102-1", :p_oc))
        t = series[1]
        p = series[2]

        #Obtain PSCAD benchmark data
        M = get_csv_data(csv_file)
        t_pscad = M[:, 1] .- t_offset
        p_pscad = M[:, 2]

        power = PSID.get_activepower_series(results, "generator-102-1")
        rpower = PSID.get_reactivepower_series(results, "generator-102-1")
        ir = PSID.get_real_current_series(results, "generator-102-1")
        ii = PSID.get_imaginary_current_series(results, "generator-102-1")
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(ir, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(ii, Tuple{Vector{Float64}, Vector{Float64}})
        @test LinearAlgebra.norm(p - p_pscad) <= 5e-3
        @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 24 Grid Following Inverter MassMatrixModel" begin
    path = (joinpath(pwd(), "test-24"))
    !isdir(path) && mkdir(path)
    try
        #Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            omib_sys, # system
            path,
            tspan,
            Pref_change,
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test24_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff_val[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test24_eigvals) < 1e-3

        #Solve problem in equilibrium
        @test execute!(sim, Rodas4(), dtmax = 0.001, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        #Obtain frequency data
        series = get_state_series(results, ("generator-102-1", :p_oc))
        t = series[1]
        p = series[2]

        #Obtain PSCAD benchmark data
        M = get_csv_data(csv_file)
        t_pscad = M[:, 1] .- t_offset
        p_pscad = M[:, 2]

        power = PSID.get_activepower_series(results, "generator-102-1")
        rpower = PSID.get_reactivepower_series(results, "generator-102-1")
        ir = PSID.get_real_current_series(results, "generator-102-1")
        ii = PSID.get_imaginary_current_series(results, "generator-102-1")
        @test isa(power, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(rpower, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(ir, Tuple{Vector{Float64}, Vector{Float64}})
        @test isa(ii, Tuple{Vector{Float64}, Vector{Float64}})
        @test LinearAlgebra.norm(p - p_pscad) <= 5e-3
        @test LinearAlgebra.norm(t - round.(t_pscad, digits = 3)) == 0.0

    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
using MLJ
using OutlierDetection
using OutlierDetection.Examples
using Random:MersenneTwister
using Test

include("evaluator.jl")
include("detector.jl")

# Number of dimensions for training
const trainDim = 100
encoder, decoder = MLPAutoEncoder(trainDim, 5, [50,20]; bias=false);

# julia
test_detector(DNN(d=0.1))
test_detector(DNN(d=0.1, parallel=true))
test_detector(KNN())
test_detector(KNN(parallel=true, reduction=:mean, algorithm=:balltree))
test_detector(KNN(parallel=true, reduction=:median))
test_detector(LOF())
test_detector(LOF(parallel=true))
test_detector(COF())
test_detector(COF(parallel=true))
test_detector(ABOD())
test_detector(ABOD(parallel=true, enhanced=true))
test_detector(AE(encoder=encoder, decoder=decoder))
test_detector(DeepSAD(encoder=encoder, decoder=decoder))
test_detector(ESAD(encoder=encoder, decoder=decoder))

# python
test_detector(PyABOD())
test_detector(PyCBLOF(random_state=0))
test_detector(PyCOF())
test_detector(PyCOPOD())
test_detector(PyHBOS())
test_detector(PyIForest(random_state=0))
test_detector(PyKNN())
test_detector(PyLMDD(random_state=0))
test_detector(PyLODA())
test_detector(PyLOF())
test_detector(PyLOCI())
test_detector(PyMCD(random_state=0))
test_detector(PyOCSVM())
test_detector(PyPCA(random_state=0))
# The following detector is only supported in the latest version of PyOD and currently fails
# test_detector(PyROD())
test_detector(PySOD())
test_detector(PySOS())

ext: jl


text:
#==============================================================================#
# AWSLambdaWrapper.jl
#
# AWS Lambda wrapper for Julia.
#
# See http://docs.aws.amazon.com/lambda/latest/dg/API_Reference.html
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

__precompile__()

module AWSLambdaWrapper

using JSON

function invoke_jl_lambda(lambda_module, input::String, output::IOStream)

    b64in = Base64DecodePipe(IOBuffer(input))
    b64out = Base64EncodePipe(output)
    serialize(b64out, lambda_module.lambda_function(deserialize(b64in)...))
    close(b64out)
end

function invoke_lambda(lambda_module, input, output::IOStream)

    if input isa Dict{String,Any} && haskey(input, "jl_data")
        invoke_jl_lambda(lambda_module, input["jl_data"], output)
    else
        JSON.print(output, lambda_module.lambda_function_with_event(input))
    end
end

function main(lambda_module::Module)

    # Wait for python wrapper to send "\n", signalling that input is waiting.
    while !eof(STDIN)
        readavailable(STDIN)

        # Read JSON input from file.
        input = JSON.parsefile("/tmp/lambda_in", dicttype=Dict{String,Any})

        # Invoke lambda_module and save output to file.
        open("/tmp/lambda_out", "w") do output
            invoke_lambda(lambda_module, input, output)
        end

        # Signal end of output.
        write(STDOUT, "\0\n")
    end
end

context() = JSON.parsefile("/tmp/lambda_context")

precompile(main, (Module,))
precompile(invoke_jl_lambda, (Module, String, IOStream))
precompile(invoke_lambda, (Module, String, IOStream))
precompile(invoke_lambda, (Module, Dict{String,Any}, IOStream))
precompile(invoke_lambda, (Module, Vector{Any}, IOStream))

end # module AWSLambdaWrapper

#==============================================================================#
# End of file.
#==============================================================================#


text:
import .StructArrays: StructArray

# A basic chains constructor that works with the Transition struct we defined.
function AbstractMCMC.bundle_samples(
    ts::Vector{<:AbstractTransition},
    model::DensityModel,
    sampler::MHSampler,
    state,
    chain_type::Type{StructArray};
    kwargs...
)
    samples = AbstractMCMC.bundle_samples(
        ts, model, sampler, state, Vector{NamedTuple};
        kwargs...
    )
    return StructArray(samples)
end

AbstractMCMC.chainscat(c::StructArray, cs::StructArray...) = vcat(c, cs...)


text:
module RoboSimples

using PyCall
using AbstractActuators

export NRoboClient, NRoboTest
export move, moveX, moveY, moveZ, devposition, setreference
export rmove, rmoveX, rmoveY, rmoveZ
export positionX, positionY, positionZ
export setreferenceX, setreferenceY, setreferenceZ
export numaxes, axesnames, moveto

struct NRoboClient <: AbstractCartesianRobot
    devname::String
    ip::String
    port::Int32
    server::PyObject
    axes::Vector{String}
end

AbstractActuators.numaxes(dev::NRoboClient) = length(dev.axes)
AbstractActuators.axesnames(dev::NRoboClient) = dev.axes

function NRoboClient(devname, ip="192.168.0.140", port=9543; axes=["x", "y", "z"])
    xmlrpc = pyimport("xmlrpc.client")
    server = xmlrpc.ServerProxy("http://$ip:$port")
    NRoboClient(devname, ip, port, server, axes)
end

AbstractActuators.move(dev::NRoboClient, ax, mm; r=false) =
    dev.server["move"](mm, string(ax), r)

AbstractActuators.move(dev::NRoboClient, ax::Integer, mm; r=false) =
    move(dev, dev.axes[ax], mm; r=r)

function AbstractActuators.move(dev::NRoboClient, axes::AbstractVector,
                                x=AbstractVector; r=false)
    ndof = length(axes)

    for i in 1:ndof
        move(dev, axes[i], x[i]; r=r)
    end
    return
end

AbstractActuators.moveto(dev::NRoboClient, x) =
    move(dev, dev.axes, x, r=false)

AbstractActuators.moveX(dev::NRoboClient, mm) = dev.server["moveX"](mm)
AbstractActuators.moveY(dev::NRoboClient, mm) = dev.server["moveY"](mm)
AbstractActuators.moveZ(dev::NRoboClient, mm) = dev.server["moveZ"](mm)

AbstractActuators.rmoveX(dev::NRoboClient, mm) = dev.server["rmoveX"](mm)
AbstractActuators.rmoveY(dev::NRoboClient, mm) = dev.server["rmoveY"](mm)
AbstractActuators.rmoveZ(dev::NRoboClient, mm) = dev.server["rmoveZ"](mm)

import Base
function AbstractActuators.devposition(dev::NRoboClient; pulses=false)
    x = dev.server["position"]("x", pulses)
    y = dev.server["position"]("y", pulses)
    z = dev.server["position"]("z", pulses)

    return Dict{String,Float64}("x"=>x, "y"=>y, "z"=>z)
end

AbstractActuators.devposition(dev::NRoboClient, ax; pulses=false) =
    dev.server["position"](string(ax), pulses)
AbstractActuators.devposition(dev::NRoboClient, ax::Integer; pulses=false) =
    devposition(dev, dev.axes[ax]; pulses=pulses)

function AbstractActuators.devposition(dev::NRoboClient, axes::AbstractVector;
                                       pulses=false)
    ndof = length(axes)

    if ndof == 1
        return [devposition(dev, dev.axes[axes[1]]; pulses=pulse)]
    elseif ndof == 2
        return [devposition(dev, dev.axes[axes[1]]; pulses=pulse),
                devposition(dev, dev.axes[axes[2]]; pulses=pulse)]
    else
        return [devposition(dev, dev.axes[axes[1]]; pulses=pulse),
                devposition(dev, dev.axes[axes[2]]; pulses=pulse),
                devposition(dev, dev.axes[axes[3]]; pulses=pulse)]
    end
end

AbstractActuators.positionX(dev::NRoboClient; pulses=false) =
    devposition(dev, "x"; pulses=pulses)
AbstractActuators.positionY(dev::NRoboClient; pulses=false) =
    devposition(dev, "y"; pulses=pulses)
AbstractActuators.positionZ(dev::NRoboClient; pulses=false) =
    devposition(dev, "z"; pulses=pulses)

AbstractActuators.setreference(dev::NRoboClient, ax, mm=0; pulses=false) =
    dev.server["set_reference"](string(ax), mm, pulses)
AbstractActuators.setreference(dev::NRoboClient, ax::Integer, mm=0; pulses=false) =
    setreference(dev, dev.axes[ax], mm; pulses=pulses)

function AbstractActuators.setreference(dev::NRoboClient,
                                        ax::AbstractVector,
                                        mm=0; pulses=false)
    nax = length(ax)

    if length(mm) == 1
        mm = fill(mm[1], nax)
    end

    for i in 1:nax
        setreference(dev, ax[i], mm[i]; pulses=pulses)
    end
end

AbstractActuators.setreferenceX(dev::NRoboClient, mm=0 ; pulses=false) =
    setreference(dev, "x", mm; pulses=pulses)

AbstractActuators.setreferenceY(dev::NRoboClient, mm=0 ; pulses=false) =
    setreference(dev, "y", mm; pulses=pulses)

AbstractActuators.setreferenceZ(dev::NRoboClient, mm=0 ; pulses=false) =
    setreference(dev, "z", mm; pulses=pulses)

end


text:
struct AntiHermitianToeplitzStorage{T, scalar} <: AbstractStorage{T, scalar}
  data::Array{T,3}
  norb::Int
  N::Int
  M::Int

  function AntiHermitianToeplitzStorage{T,scalar}(data::AbstractArray{T,3}) where {T, scalar}
    norb = size(data,1)
    N = size(data,3)
    @assert size(data) == (norb, norb, N)
    @assert !(scalar && norb > 1)
    new(data, norb, N, N)
  end
end

function AntiHermitianToeplitzStorage(data::AbstractArray{T,3}, scalar=false) where T
  AntiHermitianToeplitzStorage{T,scalar}(data)
end

function AntiHermitianToeplitzStorage(data::AbstractArray{T,4}, scalar=false) where T
  norb = size(data,1)
  N = size(data,3)
  @assert size(data) == (norb, norb, N, N)
  col = zeros(T, norb, norb, N)
  for i in 1:N
    col[:,:,i] = data[:,:,i,1]
  end
  AntiHermitianToeplitzStorage(col, scalar)
end

function AntiHermitianToeplitzStorage(::Type{T}, N::Integer, norb=1, scalar=false) where T <: Number
  data = zeros(T, norb, norb, N)
  return AntiHermitianToeplitzStorage(data, scalar)
end

AntiHermitianToeplitzStorage(N::Integer, norb=1, scalar=false) = TimeInvariantAntiHermitian(ComplexF64, N, norb,  scalar)

function Base.getindex(X::AntiHermitianToeplitzStorage{T,scalar}, k, l, i::Int, j::Int) where {T, scalar}
  i < j && return copy(-adjoint(X[l,k,j,i]))
  n = i-j+1
  return X.data[k,l,n]
end

function Base.setindex!(X::AntiHermitianToeplitzStorage{T,scalar}, v, k, l, i::Int, j::Int) where {T, scalar}
  i < j && return X[l,k,j,i] = -adjoint(v)
  n = i-j+1
  X.data[k,l,n] = v
end


text:
struct TracedModel{S<:AbstractSampler,V<:AbstractVarInfo,M<:Model,E<:Tuple}
    model::M
    sampler::S
    varinfo::V
    evaluator::E
end

function TracedModel(
    model::Model,
    sampler::AbstractSampler,
    varinfo::AbstractVarInfo,
)
    # evaluate!!(m.model, varinfo, SamplingContext(Random.AbstractRNG, m.sampler, DefaultContext()))
    context = SamplingContext(DynamicPPL.Random.GLOBAL_RNG, sampler, DefaultContext())
    evaluator = _get_evaluator(model, varinfo, context)
    return TracedModel{AbstractSampler,AbstractVarInfo,Model,Tuple}(model, sampler, varinfo, evaluator)
end

# Smiliar to `evaluate!!` except that we return the evaluator signature without excutation.
# TODO: maybe move to DynamicPPL
@generated function _get_evaluator(
    model::Model{_F,argnames}, varinfo, context
) where {_F,argnames}
    unwrap_args = [
        :($DynamicPPL.matchingvalue(context_new, varinfo, model.args.$var)) for var in argnames
    ]
    # We want to give `context` precedence over `model.context` while also
    # preserving the leaf context of `context`. We can do this by
    # 1. Set the leaf context of `model.context` to `leafcontext(context)`.
    # 2. Set leaf context of `context` to the context resulting from (1).
    # The result is:
    # `context` -> `childcontext(context)` -> ... -> `model.context`
    #  -> `childcontext(model.context)` -> ... -> `leafcontext(context)`
    return quote
        context_new = DynamicPPL.setleafcontext(
            context, DynamicPPL.setleafcontext(model.context, DynamicPPL.leafcontext(context))
        )
        (model.f, model, DynamicPPL.resetlogp!!(varinfo), context_new, $(unwrap_args...))
    end
end

function Base.copy(trace::AdvancedPS.Trace{<:TracedModel})
    f = trace.f
    newf = TracedModel(f.model, f.sampler, deepcopy(f.varinfo))
    return AdvancedPS.Trace(newf, copy(trace.ctask))
end

function AdvancedPS.advance!(trace::AdvancedPS.Trace{<:TracedModel})
    DynamicPPL.increment_num_produce!(trace.f.varinfo)
    score = consume(trace.ctask)
    if score === nothing
        return
    else
        return score + DynamicPPL.getlogp(trace.f.varinfo)
    end
end

function AdvancedPS.delete_retained!(f::TracedModel)
    DynamicPPL.set_retained_vns_del_by_spl!(f.varinfo, f.sampler)
    return
end

function AdvancedPS.reset_model(f::TracedModel)
    newvarinfo = deepcopy(f.varinfo)
    DynamicPPL.reset_num_produce!(newvarinfo)
    return TracedModel(f.model, f.sampler, newvarinfo)
end

function AdvancedPS.reset_logprob!(f::TracedModel)
    DynamicPPL.resetlogp!!(f.varinfo)
    return
end


text:
module Rootionals
using UnicodeFun
using Primes

import Base: numerator, denominator, one, zero, show, ==,
    *, /, ^, sqrt, cbrt,
    float

mutable struct Rootional
    b::Vector{Int}
    e::Vector{Rational}
end

Rootional(z::Int) = simplify!(Rootional([z], [1]))
Rootional(q::Rational) = simplify!(Rootional([numerator(q), denominator(q)], [1, -1]))

one(::Type{Rootional}) = Rootional(1)
zero(::Type{Rootional}) = Rootional(0)

==(a::Rootional, b::Rootional) = (a.b == b.b) && (a.e == b.e)
==(a::Rootional, b::Union{Int,Rational}) = a == Rootional(b)
==(a::Union{Int,Rational}, b::Rootional) = Rootional(a) == b

*(a::Rootional, b::Rootional) =
    simplify!(Rootional(vcat(a.b,b.b), vcat(a.e, b.e)))
*(a::Rootional, b::Union{Int,Rational}) = a*Rootional(b)
*(a::Union{Int,Rational}, b::Rootional) = Rootional(a)*b

/(a::Rootional, b::Rootional) =
    simplify!(Rootional(vcat(a.b,b.b), vcat(a.e, -b.e)))
/(a::Rootional, b::Union{Int,Rational}) = a/Rootional(b)
/(a::Union{Int,Rational}, b::Rootional) = Rootional(a)/b

^(a::Rootional, e::Union{Int,Rational}) = simplify!(Rootional(a.b, e*a.e))

sqrt(a::Rootional) = a^(1//2)
cbrt(a::Rootional) = a^(1//3)

function format_factor(n,e)
    en = numerator(e)
    ed = denominator(e)
    core = if ed == 1
        "$(n)"
    else
        root = if ed == 2
            "√"
        elseif ed == 3
            "∛"
        elseif ed == 4
            "∜"
        else
            "$(to_superscript(ed))√"
        end
        "$(root)$(n)"
    end
    if en == 1
        core
    else
        es = to_superscript(en)
        if ed == 1
            "$(core)$(es)"
        else
            "($(core))$(es)"
        end
    end
end

function show(stream::IO, r::Rootional)
    ni = findall(r.e .> 0)
    di = findall(r.e .< 0)
    if length(ni) > 0
        for i in ni
            write(stream, format_factor(r.b[i],r.e[i]))
        end
    else
        write(stream, "1")
    end
    if length(di) > 0
        write(stream, "/")
        for i in di
            write(stream, format_factor(r.b[i],-r.e[i]))
        end
    end
end

function numerator(r::Rootional)
    ni = findall(r.e .> 0)
    if length(ni) > 0
        Rootional(r.b[ni], r.e[ni])
    else
        one(Rootional)
    end
end

function denominator(r::Rootional)
    di = findall(r.e .< 0)
    if length(di) > 0
        Rootional(r.b[di], -r.e[di])
    else
        one(Rootional)
    end
end

function float(r::Rootional)
    f = 1.0
    for i in eachindex(r.b)
        f *= float(r.b[i])^r.e[i]
    end
    f
end

"Simplify a rootional by cancelling and merging common bases."
function simplify!(r::Rootional)
    d = Dict{Int,Rational}()
    for i in eachindex(r.b)
        for (f,e) in factor(r.b[i])
            d[f] = get(d, f, 0) + e*r.e[i]
        end
    end
    delete!(d, 1)
    for b in keys(d)
        d[b] == 0 && delete!(d, b)
    end
    length(d) == 0 && (d[1]=1)
    r.b = sort(collect(keys(d)))
    r.e = [d[b] for b in r.b]
    r
end

export Rootional, simplify!

end # module


text:
module Binscatters

using Statistics
using DataFrames
using FixedEffectModels
using RecipesBase

include("utils.jl")
"""
    binscatter(df::Union{DataFrame, GroupedDataFrame}, f::FormulaTerm, n::Integer = 20;
                weights::Union{Symbol, Nothing} = nothing, seriestype::Symbol = :scatter,
                kwargs...)

Outputs a binned scatterplot

### Arguments
* `df`: A Table, a DataFrame or a GroupedDataFrame
* `f`: A formula created using [`@formula`](@ref). The variable(s) in the left-hand side are plotted on the y-axis. The first variable in the right-hand side is plotted on the x-axis. Add other variables for controls.
* `n`: Number of bins (default to 20).

### Keyword arguments
* `weights`: A symbol for weights
* `seriestype`:  `:scatter` (the default) plots bins, `:scatterpath` adds a line to connect the bins, `:linearfit` adds a regression line.
* `kwargs...`: Additional attributes for [`plot`](@ref).

### Examples
```julia
using DataFrames, Binscatters, RDatasets, Plots
df = dataset("plm", "Cigar")
binscatter(df, @formula(Sales ~ Price))

# Change the number of bins
binscatter(df, @formula(Sales ~ Price), 10)

# Residualize w.r.t. controls
binscatter(df, @formula(Sales ~ Price + NDI))
binscatter(df, @formula(Sales ~ Price + fe(Year)))

# Plot multiple variables on the y-axis
binscatter(df, @formula(Sales + NDI ~ Price))

# Plot binscatters within groups
df.Post70 = df.Year .>= 70
binscatter(groupby(df, :Post70), @formula(Sales ~ Price))
```
"""
binscatter

mutable struct Binscatter
    args
end
binscatter(args...; kwargs...) = RecipesBase.plot(Binscatter(args); kwargs...)
binscatter!(args...; kwargs...) = RecipesBase.plot!(Binscatter(args); kwargs...)

# User recipe
@recipe function f(bs::Binscatter; weights = nothing)
    df = bin(bs.args...; weights = weights)
    if df isa DataFrame
        cols = names(df)
        N = length(cols)
        x = df[!, end]
        Y = Matrix(df[!, 1:(end-1)])
        @series begin
            seriestype --> :scatter
            markerstrokealpha --> 0.0
            xguide --> cols[end]
            if size(Y, 2) == 1
                yguide --> cols[1]
                label --> ""
            else
                label --> reshape(cols[1:(end-1)], 1, N-1)
            end
            x, Y
        end
    elseif df isa GroupedDataFrame
        for (k, out) in pairs(df)
            cols = string.(valuecols(df))
            N = length(cols)
            x = out[!, end]
            Y = Matrix(out[!, (end-(N-1)):(end-1)])
            str = "(" * join((string(k) * " = " * string(v) for (k, v) in pairs(NamedTuple(k))), ", ") * ")"
            @series begin
                seriestype --> :scatter
                markerstrokealpha --> 0.0
                xguide --> cols[end]
                if size(Y, 2) == 1
                    yguide --> cols[1]
                    label --> str
                else
                    label --> reshape(cols[1:(end-1)], 1, N-1) .* " " .* str
                end
                x, Y
            end
        end
    end
end

function bin(df::AbstractDataFrame, @nospecialize(f::FormulaTerm), n::Integer = 20;
            weights::Union{Symbol, Nothing} = nothing)
    df = partial_out(df, _shift(f); weights = weights, align = false, add_mean = true)[1]
    cols = names(df)
    df.__cut = cut(df[!, end], n; allowempty = true)
    df = groupby(df, :__cut)
    combine(df, cols .=> mean .=> cols; keepkeys = false)
end

function bin(df::GroupedDataFrame, @nospecialize(f::FormulaTerm), n::Integer = 20;
            weights::Union{Symbol, Nothing} = nothing)
    combine(d -> bin(d, f, n; weights = weights), df; ungroup = false)
end

function bin(df, @nospecialize(f::FormulaTerm), n::Integer = 20;
            weights::Union{Symbol, Nothing} = nothing)
    bin(DataFrame(df), f, n; weights = weights)
end
export @formula, fe, binscatter, binscatter!

end


text:
module RomanNumerals

# Defines the exception the struct RomanNumeral
include("types.jl")

# Defines the functions toroman, fromroman and parse(RomanNumeral, x), to
# convert from Roman-formatted string to int and vice versa
include("parsing.jl")

# Defines new method dispatching specific to Roman numerals
include("arithmetic.jl")

export RomanNumeral, @rn_str

end # module RomanNumerals


text:
using Luna
import StaticArrays: SMatrix, SVector
import LinearAlgebra: mul!, ldiv!

a = 50e-6
gas = :Ar
pres = 5

τfwhm = 30e-15
λ0 = 400e-9
energy = 30e-6

modes = (Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0),
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=π/2))
nmodes = length(modes)

grid = Grid.EnvGrid(10e-2, λ0, (160e-9, 3000e-9), 1e-12)

energyfun = Fields.energyfuncs(grid)[1]
normfun = NonlinearRHS.norm_modal(grid.ω)

function gausspulse(t)
    It = Maths.gauss(t, fwhm=τ)
    Et = @. sqrt(It)
end

field = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)
Etlin = gausspulse(grid.t)
cenergy = energyfun(Etlin)
Etlin = sqrt(energy)/sqrt(cenergy) .* Etlin

Et = [Etlin zero(grid.t)] .+ 0im
Ew = FFTW.fft(Et, 1);
Ewd = size(Ew)
Ewp = Array{ComplexF64,2}(undef, Ewd[2], Ewd[1])
permutedims!(Ewp, Ew, [2, 1])

CP = Polarisation.rotate(Polarisation.WP(π/2), π/6 + π/8)
Ewp .= CP * Ewp
permutedims!(Ew, Ewp, [2, 1])
Et = FFTW.ifft(Ew, 1);

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),)

# dummy, we don't use these
inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, normfun, responses, inputs,
                              modes, :xy; full=false)

Eω .= Ew

statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
linop = LinearOps.make_const_linop(grid, modes, λ0)

Luna.run(Eω, grid, linop, transform, FT, output)

import FFTW
import PyPlot:pygui, plt

ω = FFTW.fftshift(grid.ω)
t = grid.t

zout = output.data["z"]
Eout = output.data["Eω"]
Eoutd = size(Eout)
Ewp = Array{ComplexF64,3}(undef, Eoutd[2], Eoutd[1], Eoutd[3])
permutedims!(Ewp, Eout, [2, 1, 3])
S = Array{Float64,3}(undef, 4, Eoutd[1], Eoutd[3])
elip = Array{Float64,2}(undef, Eoutd[1], Eoutd[3])

for i = 1:Eoutd[3]
    for j = 1:Eoutd[1]
        S[:,j,i] .= Polarisation.Stokes(Ewp[:,j,i])
        elip[j,i] = Polarisation.ellipticity(S[:,j,i])
    end
end

pygui(true)
plt.figure()
plt.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(S[3,:,:]./S[1,:,:],1))
plt.colorbar()

plt.figure()
plt.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(S[2,:,:]./S[1,:,:],1))
plt.colorbar()

plt.figure()
plt.pcolormesh(zout, ω./2π.*1e-15, FFTW.fftshift(elip,1))
plt.colorbar()

Etout = FFTW.ifft(Eout, 1)
It = abs2.(Etout)

Ilog = FFTW.fftshift(log10.(Maths.normbymax(abs2.(Eout))), 1)

for i = 1:nmodes
    plt.figure()
    plt.subplot(121)
    plt.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog[:,i,:]))
    plt.clim(-6, 0)
    plt.xlim(0,2.0)
    plt.colorbar()
    plt.subplot(122)
    plt.pcolormesh(t.*1e15, zout, transpose(It[:,i,:]))
    plt.xlim(-30.0,100.0)
    plt.colorbar()
end


text:
"""
    ResultType
"""
@enum ResultType begin
    RESULT_RAW = 0
    RESULT_NODE
    RESULT_VALUE
end
function with_result(f, type::ResultType=RESULT_VALUE; this=nothing)
    scope = type == RESULT_RAW ? nothing : open_scope()
    result = f()
    ret = if type == RESULT_RAW
        isnothing(result) ? get_undefined() : result
    elseif type == RESULT_NODE
        node_value(result)
    else
        value(result; this=this)
    end
    close_scope(scope)
    ret
end

"""
    run_script(script; context=current_context(), result=RESULT_VALUE)

Evaluate the given JavaScript script in a given `context`. If `context` is `nothing`,
the script will be run in the global scope.

If `result` is set `RESULT_RAW`, the raw `napi_value` will be returned.
Note that it is only available inside current `napi_handle_scope`.
Thus, if you did't open one, the `napi_value` would remain in the global
`napi_handle_scope`.

If `result` is `RESULT_VALUE` (which is the default), the result will be converted to
the most suitable Julia type before being returned.
"""
function run_script(
    script::AbstractString,
    result=RESULT_VALUE;
    context=current_context(),
)
    script = strip(script)
    script = if startswith(script, '{')
        "($script)"
    else
        script
    end
    if isnothing(context)
        with_result(result) do
            @napi_call napi_run_script(script::NapiValue)::NapiValue
        end
    else
        _VM.runInContext(
            script, context;
            result=result
        )
    end
end

const LOCALS_NAME = "__jlnode_locals"
const LOCALS_COUNTER = Ref(0)
@global_js_const _JS_CREATE_LOCALS = """(ctx, i) => {
    if (!("$LOCALS_NAME" in ctx)) {
        ctx.$LOCALS_NAME = {}
    }
    const locals = ctx.$LOCALS_NAME
    locals[i] = {}
    return locals[i]
}"""
# From pyeval.jl
const _INTERPOLATE_FSM = Dict(
    ('p', '\'') => '\'',
    ('p', '"') => '"',
    ('p', '`') => '`',
    ('\'', '\'') => 'p',
    ('"', '"') => 'p',
    ('`', '`') => 'p',
    ('p', '$') => '$',

    ('`', '$') => 'S',
    ('p', '/') => '/',
    ('/', '/') => 'c',
    ('/', '*') => 'C',
    ('c', '\n') => 'p',
    ('C', '*') => 'E',
    ('E', '*') => 'E',
    ('E', '/') => 'p',
)
# Shouldn't be used before `__init__()`.
function interpolate_code(code::AbstractString)
    LOCALS_COUNTER[] += 1
    locals_sym = "$LOCALS_NAME[$(LOCALS_COUNTER[])]"

    buf = IOBuffer()
    state = 'p'
    i = 1 # position in code
    locals = Dict{Union{String,Int},Any}()
    numlocals = 0
    js_interpolating = 0
    js_braces = 0
    while i <= lastindex(code)
        c = code[i]
        newstate = get(_INTERPOLATE_FSM, (state, c), '?')
        if newstate == '$'
            i += 1
            i > lastindex(code) && error("unexpected end of string after \$")
            interp_literal = false
            if code[i] == '$' # $$foo pastes the string foo into the Python code
                i += 1
                interp_literal = true
            end
            expr, i = Meta.parse(code, i, greedy=false)
            if interp_literal
                # need to save both the expression and the position
                # in the string where it should be interpolated
                locals[position(buf)+1] = expr
            else
                numlocals += 1
                localvar = "$locals_sym[$numlocals]"
                locals[string(numlocals)] = expr
                print(buf, localvar)
            end
        else
            if newstate == '?'
                if state == 'p'
                    if c == '{'
                        js_braces += 1
                        newstate = 'p'
                    elseif c == '}'
                        if js_braces > 0
                            js_braces -= 1
                            newstate = 'p'
                        elseif js_interpolating > 0
                            js_interpolating -= 1
                            newstate = '`'
                        else
                            error("Unexpected symbol }")
                        end
                    elseif c == '\\'
                        if i < lastindex(code)
                            i0 = nextind(code, i)
                            if code[i0] == '$'
                                c = '$'
                                i = i0
                            end
                        end
                        newstate = 'p'
                    else
                        newstate = 'p'
                    end
                elseif state == 'E'
                    newstate = 'C'
                elseif state == 'S'
                    if c == '{'
                        js_interpolating += 1
                        newstate = 'p'
                    else
                        newstate = '`'
                    end
                elseif state ∈ ('\'', '"', '`', '/', 'c', 'C')
                    newstate = state
                end
            end
            print(buf, c)
            state = newstate
            i = nextind(code, i)
        end
    end
    return String(take!(buf)), locals_sym, locals
end

"""
    node"...js code..."[options]

Evaluate the given JavaScript script in the current context.
Use `\$` to interpolate Julia variables.

Available `options`:

- `r`: Keep the result in the raw `NapiValue`.
- `o`: Do not try to convert the result into Julia values.
- `g`: Run the script in the global scope (instead of in a `vm.Context`).
- `f`: Do not delete locals after running (for the interpolated variables).

See docs for `run_script` for details.
"""
macro node_str(code, options...)
    result, global_, keep_locals = if length(options) == 1
        o = options[1]
        if 'r' ∈ o
            RESULT_RAW
        elseif 'o' ∈ o
            RESULT_NODE
        else
            RESULT_VALUE
        end, 'g' ∈ o, 'f' ∈ o
    else
        RESULT_VALUE, false, false
    end
    context = global_ ? :(nothing) : :(current_context())
    code, locals_sym, locals = interpolate_code(code)
    counter = LOCALS_COUNTER[]
    code_expr = Expr(:call, esc(:(Base.string)))
    i0 = firstindex(code)
    for i in sort!(collect([k for k in keys(locals) if k isa Int]))
        push!(code_expr.args, code[i0:prevind(code,i)], esc(locals[i]))
        i0 = i
    end
    push!(code_expr.args, code[i0:lastindex(code)])
    filter!(locals) do (k, _)
        k isa String
    end
    has_locals = length(locals) > 0
    assign_locals = if has_locals
        exprs = quote
            locals = NodeCall._JS_CREATE_LOCALS(
                NodeCall.get_global($context), $counter;
                context=$context, result=RESULT_RAW
            )
        end
        for (k, v) in locals
            if k isa String
                push!(exprs.args, :(
                    set!(locals, $k, $(esc(v)))
                ))
            end
        end
        exprs
    else
        nothing
    end
    delete_locals = if keep_locals || !has_locals
        nothing
    else
        delete_expr = "delete $locals_sym"
        :(run_script($delete_expr, RESULT_RAW; context=$context))
    end
    quote
        with_result($result) do
            $assign_locals
            code = $code_expr
            ret = run_script(code, RESULT_RAW; context=$context)
            $delete_locals
            ret
        end
    end
end

require(id::AbstractString; result=RESULT_VALUE) = run_script("globalThis.require('$(id)')", result)


text:
export topright, topleft, bottomright, bottomleft, topline, rightline, bottomline, leftline

Rect(p, w, h) = Rect(p, w, h, deg(0))

needed_attributes(::Type{Rect}) = needed_attributes(Circle)

function bottomleft(r::Rect)
    r.center + rotate(Point(-r.width * 0.5, -r.height * 0.5), r.angle)
end

function bottomright(r::Rect)
    r.center + rotate(Point( r.width * 0.5, -r.height * 0.5), r.angle)
end

function topleft(r::Rect)
    r.center + rotate(Point(-r.width * 0.5,  r.height * 0.5), r.angle)
end

function topright(r::Rect)
    r.center + rotate(Point( r.width * 0.5,  r.height * 0.5), r.angle)
end

topline(r::Rect) = Line(topleft(r), topright(r))
bottomline(r::Rect) = Line(bottomleft(r), bottomright(r))
leftline(r::Rect) = Line(bottomleft(r), topleft(r))
rightline(r::Rect) = Line(bottomright(r), topright(r))

function Point(r::Rect, nx::Real, ny::Real, mode::Symbol=:norm)
    if mode == :norm
        r.center + Point(r.angle) * (nx - 0.5) * r.width + Point(r.angle + deg(90)) * (ny - 0.5) * r.height
    else
        error("Mode $mode is invalid")
    end
end

function Rect(t::Txt)
    w = t.extent.width
    h = t.extent.height
    wh = P(w, h)
    ang = t.angle
    xshift = @match t.halign begin
        :c => 0
        :l => w/2
        :r => -w/2
        _ => error("Not implemented")
    end
    yshift = @match t.valign begin
        :c => 0
        :t => h/2
        :b => -h/2
        # :bl => (t.pos - rotate(wh / 2, ang)).y
        _ => error("Not implemented")
    end
    center = t.pos + rotate(P(xshift, yshift), ang)
    Rect(center, w, h, ang)
end

function Rect(t::Txt, margin::Real)
    r = Rect(t)
    Rect(r.center, r.width + 2 * margin, r.height + 2 * margin, r.angle)
end

Rect(bb::BBox) = Rect((bb.from + bb.to) / 2, bb.to.x - bb.from.x, bb.to.y - bb.from.y, deg(0))


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the equation of time.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/sun/equation_of_time.jl
# ===================================

# Function equation_of_time
# -------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Figure 3-22: Equation of Time Variation [1, p. 179].
#
# From this figure, one can see that:
#
#    |     Day     | Equation of Time |
#    |-------------|------------------|
#    | February 11 | [-15,-14]        |
#    | May 11      | [+3,+4]          |
#    | July 26     | [-7,-6]          |
#    | November 2  | [16,17]          |
#
################################################################################

@testset "Function equation_of_time" begin
    # February 11.
    eot = equation_of_time(DatetoJD(2000, 2, 11, 0, 0, 0))*12*60/pi
    @test -15 < eot < -14

    # May 11.
    eot = equation_of_time(DatetoJD(2000, 5, 11, 0, 0, 0))*12*60/pi
    @test +3 < eot < +4

    # July 26.
    eot = equation_of_time(DatetoJD(2000, 7, 26, 0, 0, 0))*12*60/pi
    @test -7 < eot < -6

    # November 2.
    eot = equation_of_time(DatetoJD(2000, 11, 2, 0, 0, 0))*12*60/pi
    @test 16 < eot < 17
end


text:
using PackageCompiler
using Libdl: Libdl

packages = [:PyCall, :DiffEqBase, :DifferentialEquations, :OrdinaryDiffEq]

sysimage_path = joinpath(@__DIR__, "sys_julia_project." * Libdl.dlext)

#create_sysimage(packages; sysimage_path=sysimage_path)

create_sysimage(packages; sysimage_path=sysimage_path,
                precompile_execution_file=joinpath(@__DIR__, "compile_exercise_script.jl"))


text:
# This file includes portions that were formerly part of Julia in modified form.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module PrimesIterator

export ModulePrimesIterator
export Primes, PrimePi, PrimeSieve

"""

* Primes, PrimePi, PrimeSieve
"""
const ModulePrimesIterator = ""

# Primes generating functions

const wheel = [4, 2, 4, 2, 4, 6, 2, 6]
const wheel_primes = [7, 11, 13, 17, 19, 23, 29, 31]
const wheel_indices = [
    0,    0,    0,    0,    0,    0,    0,    1,    1,
    1,    1,    2,    2,    3,    3,    3,    3,    4,
    4,    5,    5,    5,    5,    6,    6,    6,    6,
    6,    6,    7,    7
]

@inline function wheel_index(n)
    d, r = Base.divrem(n - 1, 30)
    return 8d + wheel_indices[r+2]
end

@inline function wheel_prime(n)
    d, r = (n - 1) >>> 3, (n - 1) & 7
    return 30d + wheel_primes[r+1]
end

#Internal function
function mask(limit::Int)
    limit < 7 && throw(ArgumentError("The condition limit ≥ 7 must be met."))
    n = wheel_index(limit)
    m = wheel_prime(n)
    sieve = ones(Bool, n)
    @inbounds for i = 1:wheel_index(isqrt(limit))
        if sieve[i]
            p = wheel_prime(i)
            q = p^2
            j = (i - 1) & 7 + 1
            while q ≤ m
                sieve[wheel_index(q)] = false
                q += wheel[j] * p
                j = j & 7 + 1
            end
        end
    end
    return sieve
end

#Internal function
function mask(lo::Int, hi::Int)
    7 ≤ lo ≤
    hi || throw(ArgumentError("The condition 7 ≤ lo ≤ hi must be met."))
    lo == 7 && return mask(hi)
    wlo, whi = wheel_index(lo - 1), wheel_index(hi)
    m = wheel_prime(whi)
    sieve = ones(Bool, whi - wlo)
    hi < 49 && return sieve
    small_sieve = mask(isqrt(hi))
    @inbounds for i = 1:length(small_sieve)  # don't use eachindex here
        if small_sieve[i]
            p = wheel_prime(i)
            j = wheel_index(2 * div(lo - p - 1, 2p) + 1)
            q = p * wheel_prime(j + 1)
            j = j & 7 + 1
            while q ≤ m
                sieve[wheel_index(q)-wlo] = false
                q += wheel[j] * p
                j = j & 7 + 1
            end
        end
    end
    return sieve
end

"""

Return the prime sieve, as a `BitArray`, of the positive integers (from `lo`, if specified) up to `hi`. Useful when working with either primes or composite numbers.
"""
function PrimeSieve(lo::Int, hi::Int)
    0 < lo ≤
    hi || throw(ArgumentError("The condition 0 < lo ≤ hi must be met."))
    sieve = falses(hi - lo + 1)
    lo ≤ 2 ≤ hi && (sieve[3-lo] = true)
    lo ≤ 3 ≤ hi && (sieve[4-lo] = true)
    lo ≤ 5 ≤ hi && (sieve[6-lo] = true)
    hi < 7 && return sieve
    wheel_sieve = mask(max(7, lo), hi)
    lsi = lo - 1
    lwi = wheel_index(lsi)
    @inbounds for i = 1:length(wheel_sieve)   # don't use eachindex here
        sieve[wheel_prime(i + lwi)-lsi] = wheel_sieve[i]
    end
    return sieve
end

function PrimeSieve(lo::T, hi::T) where {T<:Integer}
    lo ≤ hi ≤ typemax(Int) && return PrimeSieve(Int(lo), Int(hi))
    throw(ArgumentError("Both endpoints of the interval to sieve must be ≤ $(typemax(Int)), got $lo and $hi."))
end

PrimeSieve(limit::Int) = PrimeSieve(1, limit)

function PrimeSieve(n::Integer)
    n ≤ typemax(Int) && return PrimeSieve(Int(n))
    throw(ArgumentError("Requested number of primes must be ≤ $(typemax(Int)), got $n."))
end

"""

Return the collection of the prime numbers (from `lo`, if specified) up to `hi`.
"""
function Primes(lo::Int, hi::Int)
    lo ≤ hi || throw(ArgumentError("The condition lo ≤ hi must be met."))
    list = Int[]
    lo ≤ 2 ≤ hi && push!(list, 2)
    lo ≤ 3 ≤ hi && push!(list, 3)
    lo ≤ 5 ≤ hi && push!(list, 5)
    hi < 7 && return list
    lo = max(2, lo)
    sizehint!(
        list,
        5 +
        floor(Int, hi / (log(hi) - 1.12) - lo / (log(lo) - 1.12 * (lo > 7)))
    )
    sieve = mask(max(7, lo), hi)
    lwi = wheel_index(lo - 1)
    # don't use eachindex here
    @inbounds for i = 1:length(sieve)
        sieve[i] && push!(list, wheel_prime(i + lwi))
    end
    return list
end

Primes(n::Int) = Primes(1, n)

"""

Return the number of primes ``≤ n``.
"""
PrimePi(n::Int) = length(Primes(1, n))

#START-TEST-########################################################

using Test

function test()
    @testset "Primes" begin
        @test PrimePi(2^14) == 1900
        @test PrimePi(10^6) == 78498
    end
end

function demo()
    println(Primes(0, 100))
    println(Primes(90, 100))
    println(PrimePi(10000))
end

"""

primepi = PrimePi(1000000)
0.005080 seconds (5 allocations: 876.141 KiB)
"""
function perf()
    @time primepi = PrimePi(1000000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
[97]
1229

  0.005080 seconds (5 allocations: 876.141 KiB)

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# code for gauge fields

# TODO for now the wilson lines will be `Vector`s of SU(2) matrices
# TODO looks like for arbitrary dimensions we'll have lots of splatting inefficiency...
#   this should be ok because most of this stuff should just be for initialization...
#   in the future can fix splatting inefficiency with macros

const WilsonLineVar = Vector{UnitaryMatrix}
export WilsonLineVar

wilsonlineunity(N::Integer, d::Integer) = UnitaryMatrix[eye(Complex{Float64}, N) for i ∈ 1:d]
wilsonlineunity(N::Integer, lat::AbstractLattice{d}) where d = wilsonlineunity(N, d)
export wilsonlineunity

randwilsonlineSU2(d::Integer) = UnitaryMatrix[randSU2() for i ∈ 1:d]
randwilsonlineSU2(lat::AbstractLattice{d}) where d = randwilsonlineSU2(d)
export randwilsonlineSU2

"""
    GaugeFieldSU2{d} <: AbstractGaugeField

An SU(2) gauge field.  Note that we take the fundamental field variables to be Wilson lines emanating
from lattice site x in direction μ. Here `d` is the total number of dimensions `D+1`.
"""
struct GaugeFieldSU2{d} <: AbstractGaugeField{d}
    field::Array{WilsonLineVar, d}
    lattice::AbstractLattice{d}

    GaugeFieldSU2{d}(st::AbstractArray{WilsonLineVar,d}, lat::AbstractLattice) where d = new(st, lat)
end
export GaugeFieldSU2

function GaugeFieldSU2(spacetime::AbstractArray{WilsonLineVar,d}, lat::AbstractLattice{d}) where d
    GaugeFieldSU2{d}(spacetime, lat)
end

# for now this just creates a "vacuum" with A=0 everywhere
function GaugeFieldSU2(lat::AbstractLattice)
    # note inefficiency due to splatting
    spacetime = fieldsites(WilsonLineVar, lat)
    for idx ∈ eachindex(spacetime)
        spacetime[idx] = wilsonlineunity(2, lat)
    end
    GaugeFieldSU2(spacetime, lat)
end

# TODO ensure that splatting is really ok here!
# gets the gauge field at spacetime location x1, t
function (U::GaugeFieldSU2)(x::Tuple)
    idx = getindex(U.lattice, x)
    getindex(U.field, idx...)
end
(U::GaugeFieldSU2)(x::SVector) = U(Tuple(x))
(U::GaugeFieldSU2)(x...) = U(x...)

# returns spacetime function of U[μ]
getindex(U::GaugeFieldSU2, μ::Integer) = x -> U(x)[μ]

eachspacetimeindex(U::AbstractGaugeField{d}) where d = eachspacetimeindex(U.lattice)

function rand(::Type{GaugeFieldSU2}, lat::AbstractLattice)
    spacetime = fieldsites(WilsonLineVar, lat)
    for idx ∈ eachindex(spacetime)
        spacetime[idx] = randwilsonlineSU2(lat)
    end
    GaugeFieldSU2(spacetime, lat)
end

"""
    untracedplaquette(U::AbstractGaugeField, μ::Int, ν::Int, x::Tuple)

Compute the value of a plaquette before taking the trace.  In other words
\$\$ U_{\mu}(x)U_{\nu}(x+e_{\mu})U^{\dagger}_{\mu}(x+e_{\nu})U^{\dagger}_{\nu}(x) \$\$
"""
function untracedplaquette(U::AbstractGaugeField, μ::Integer, ν::Integer, x::SVector)
    U[μ](x)*U[ν](inc(x,μ))*U[μ](inc(x,ν))'*U[ν](x)'
end
export untracedplaquette

"""
    plaquette(U::AbstractGaugeField, μ::Integer, ν::Integer, x::SVector)

Compute the value of the plaquette at location `x` over the gauge field.

Note that we follow the definition in the "Lattice QCD For Novices" paper and define the paquette
to be the real part of the trace.

*TODO*: Note that ideally the normalization factor here should depend on the gauge group.
"""
plaquette(U::AbstractGaugeField, μ::Integer, ν::Integer, x::SVector) = real(tr(untracedplaquette(U,μ,ν,x)))/2.0
export plaquette

"""
    wilsonlagrangian(U::AbstractGaugeField, β::AbstractFloat, x::SVector)

Compute the lowest order Wilson plaquette Lagrangian at location `x` and inverse temperature `β`.

Note that this has a factor of a^4 in front of it relative to the usual G^2 Yang-Mills kinetic term.

Note that we neglect the constant term altogether.
"""
function wilsonlagrangian(U::AbstractGaugeField{d}, β::AbstractFloat, x::SVector) where d
    ℒ = 0.0
    for μ ∈ 1:d
        for ν ∈ 1:(μ-1)
            ℒ -= plaquette(U, μ, ν, x)
        end
    end
    β*ℒ
end
export wilsonlagrangian

"""
    wilsonaction(U::AbstractGaugeField, β::AbstractFloat)

Compute the Wilson plaquette action at inverse temperature `β`.

Again, the constant term is neglected.
"""
function wilsonaction(U::AbstractGaugeField{d}, β::AbstractFloat) where d
    S = 0.0
    for x ∈ eachspacetimeindex(U)
        S += wilsonlagrangian(U, β, x)
    end
    S/U.lattice.a^(4 - d) # this factor to fix the powers of a
end
export wilsonaction


text:
##############
### Testing CALiPPSO in the same system, with different solvers
##############
using CALiPPSO
using DelimitedFiles, Statistics

const N, d = 1024, 3 # size and dimensionality of the configuration  to be used
const L = 1.0 # size of each side of the system's volume

# Note that this value is x100 larger than the default one (that can be used with Gurobi). Some solvers need this smaller precision
const tol_Γ_convergence = 1e-12
const tol_overlap = CALiPPSO.default_tol_overlap
const tol_zero_forces = CALiPPSO.default_tol_zero_forces
const default_tol_displacements = CALiPPSO.default_tol_displacements_convergence
const tol_optimality = CALiPPSO.default_tol_optimality
const max_threads = CALiPPSO.max_threads

const solvers = [:Gurobi, :HiGHS, :Clp, :GLPK, :COSMO] # array of solver's names (or symbols)
using Gurobi, HiGHS, Clp, GLPK, COSMO #,Hypatia

println("This script will test the following solvers: ",  map(x->x*", ", string.(solvers))...)
println("But before the test is executed, `produce_jammed_configuration!` will be compiled with each of them.\n\n")

##############################
### Now let's define the arguments and attributes to be used for each solver, and create an empty model for initialization
##############################

# Gurobi:
const grb_args = Gurobi.Env()
const grb_attributes = Dict("OutputFlag" => 0, "FeasibilityTol" => tol_optimality, "OptimalityTol" => tol_optimality, "Method" => 3, "Threads" => max_threads)
precompile_main_function(Gurobi, grb_attributes, grb_args)

# HiGHS
const highs_args = nothing
const highs_attributes = Dict("small_matrix_value"=>0.1*tol_optimality, "primal_feasibility_tolerance" => tol_optimality, "dual_feasibility_tolerance" => tol_optimality, "solver" => "ipm", "ipm_optimality_tolerance"=>tol_optimality,  "threads" => max_threads, "parallel"=>"on", "output_flag"=>false)
printstyled("\n______________________________________________________________________________________\n", color=:yellow)
precompile_main_function(HiGHS, highs_attributes, highs_args)

# Clp
const clp_args = nothing
const clp_attributes = Dict("PrimalTolerance"=>tol_optimality, "DualTolerance" => tol_optimality, "LogLevel" => 0, "SolveType" => 5)
printstyled("\n______________________________________________________________________________________\n", color=:yellow)
precompile_main_function(Clp, clp_attributes, clp_args)

# GLPK
const glpk_args = (want_infeasibility_certificates=false, method=GLPK.MethodEnum(0) )
const glpk_attributes = Dict("msg_lev"=>GLPK.GLP_MSG_OFF, "tol_bnd"=>tol_optimality, "tol_dj"=>tol_optimality)
printstyled("\n______________________________________________________________________________________\n", color=:yellow)
precompile_main_function()

# COSMO
const cosmo_args = nothing
const cosmo_attributes = Dict("verbose"=>false, "max_iter"=>30000, "eps_abs"=>1e3*tol_optimality, "adaptive_rho"=>false, "eps_prim_inf"=>1e3*tol_optimality, "eps_dual_inf"=>1e3tol_optimality)
printstyled("\n______________________________________________________________________________________\n", color=:yellow)
precompile_main_function(COSMO, cosmo_attributes, cosmo_args)

# # Hypatia
# const hypa_args = (verbose = false, tol_abs_opt = tol_optimality, tol_feas = 0.1*tol_overlap, tol_infeas = 0.1*tol_overlap)
# const hypa_attributes = Dict()
# printstyled("\n______________________________________________________________________________________\n", color=:yellow)
# precompile_main_function(Hypatia, hypa_attributes, hypa_args)

#### create arrays of these variables, for easier calling
const solvers_args = [grb_args, highs_args, clp_args, glpk_args, cosmo_args]
const solvers_attributes = [grb_attributes, highs_attributes,clp_attributes ,glpk_attributes, cosmo_attributes]

## Array of packings to perform later comparison
all_packings = Vector{MonoPacking{d, Float64}}(undef, length(solvers))

# The initial conditions are configurations compressed with Lubachevsky--Stillinger up to a (reduced) pressure p>>10
#### read the initial configuration on which the test is performed
σ0 = readdlm("Centers-after-LS--N-$(N)--d-$d.dat", skipstart=2)[1, 1] # initial diameter of the particles
r0_common = 0.5*σ0 # initial radius of particles
Xs_common = PeriodicVectors(Matrix(readdlm("Centers-after-LS--N-$(N)--d-$d.dat", skipstart=3)'), L) # initial particles' positions

printstyled("\n\n O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O\n", color=:blue, bold=true)
printstyled("Applying CALiPPSO for jamming a system of N = ", N, " in d = ", d, "\t Initial φ = ", packing_fraction(d, r0_common, N, L), "\n", bold=:true, color=:blue)
printstyled("O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O --- O\n", color=:blue, bold=true)

for (ns, solver) in enumerate(solvers)
    Xs0 = 1.0.*Xs_common; r0=1*r0_common

    if  Symbol(solver)==:COSMO
        # continue
        overlap_tolerance = 1e3*tol_overlap
        zero_force = 1e-5
        Γ_conv = tol_Γ_convergence
    elseif Symbol(solver)==:Hypatia
        # continue
        overlap_tolerance = tol_overlap
        zero_force = 1e-5
        Γ_conv = 1e-10
    else
        # continue
        overlap_tolerance = tol_overlap
        zero_force = tol_zero_forces
        Γ_conv = tol_Γ_convergence
    end
    S_conv = default_tol_displacements

    printstyled("\n\n++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:magenta, bold=true)
    printstyled("\t\t\tUSING SOLVER: ", Symbol(solver), "\n", bold=:true, color=:magenta)
    println("\tConvergence criteria: tolerance √Γ-1 = ",Γ_conv, "\t tolerance |sᵢ| = ", S_conv)
    printstyled("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n", color=:magenta, bold=true)

    solver_args = solvers_args[ns]; solver_attrs = solvers_attributes[ns]

    @time jammed_packing, info_convergence, Γs_vs_t, smax_vs_t, iso_vs_t = produce_jammed_configuration!(Xs0, r0; verbose=true,
    solver=eval(solver), solver_args=solver_args, solver_attributes=solver_attrs, tol_Γ_convergence=Γ_conv, tol_S_convergence = S_conv,
    tol_overlap=overlap_tolerance, initial_monitor=30, zero_force=zero_force, max_iters=25)

    println("_______________________________________________________________________________________\n\n")
    times = info_convergence.times_LP_optim
    println("Info about LP times")
    println(times)
    println("(min, avg ± std, max) LP times:\t", minimum(times), ", \t", mean(times), " ± ", std(times), ", \t", maximum(times))
    println("_______________________________________________________________________________________\n\n\n")

    network_of_contacts(jammed_packing)

    all_packings[ns] = jammed_packing

end

Δpacks = map(x-> difference_in_packings(all_packings[1], x), all_packings)
Rref = all_packings[1].R
ΔRs = Rref .- getfield.(all_packings, :R)
printstyled("Comparing packing differences (only considering stable particles) , with respect to: ", solvers[1], "\n", color=:green)
for (ns, solver) in enumerate(solvers)
    println("Solver: ", solver, "\t mean and std Δ = ", Δpacks[ns], "\t and of difference of R = ", ΔRs[ns])
end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
function exponential_correlations(rateTimes::Vector{Float64}, longTermCorr::Float64, beta::Float64, gamma::Float64, time_::Float64)
  check_increasing_times(rateTimes)
  (longTermCorr <= 1.0 && longTermCorr >= 0.0) || error("ong term correlation outside of [0,1] interval")
  beta >= 0.0 || error("beta must be greater than or equal to zero")
  (gamma <= 1.0 && gamma >= 0.0) || error("gamma outside [0,1] interval")

  # Calcualte correlation matrix
  nbRows = length(rateTimes)-1
  correlations = zeros(nbRows, nbRows)

  for i = 1:nbRows
    # Correlation is defined only between (alive) stochastic rates
    if time_ <= rateTimes[i]
      correlations[i, i] = 1.0
      for j = 1:i
        if time_ <= rateTimes[j]
          correlations[i,j] = correlations[j, i] = longTermCorr + (1.0 - longTermCorr) *
                              exp(-beta * abs(^(rateTimes[i] - time_, gamma) - ^(rateTimes[j] - time_, gamma)))
        end
      end
    end
  end

  return correlations
end

type ExponentialForwardCorrelation <: PiecewiseConstantCorrelation
  numberOfRates::Int
  longTermCorr::Float64
  beta::Float64
  gamma::Float64
  rateTimes::Vector{Float64}
  times::Vector{Float64}
  correlations::Vector{Matrix{Float64}}
end

function ExponentialForwardCorrelation(rateTimes::Vector{Float64},
                                      longTermCorr::Float64 = 0.5,
                                      beta::Float64 = 0.2,
                                      gamma::Float64 = 1.0,
                                      times::Vector{Float64} = Vector{Float64}())
  # build
  numberOfRates = isempty(rateTimes) ? 0 : length(rateTimes) - 1
  numberOfRates > 1 || error("rate times must contain at least two values")

  check_increasing_times(rateTimes)

  # corrTimes must include all rateTimes but the last
  if isempty(times)
    times = rateTimes[1:end-1]
  else
    check_increasing_times(times)
  end

  if is_close(gamma, 1.0)
    temp = rateTimes[1:end-1]
    times == temp || error("corr times must be equal to all rate times (but the last)")

    c = exponential_correlations(rateTimes, longTermCorr, beta, 1.0, 0.0)
    correlations = evolved_matrices(TimeHomogeneousForwardCorrelation, c)
  else
    times[end] <= rateTimes[numberOfRates+1] || error("last corr time is after next-to_last rate time")

    correlations = similar(times, Matrix{Float64})
    time_ = times[1] / 2.0
    correlations[1] = exponential_correlations(rateTimes, longTermCorr, beta, gamma, time_)
    for k = 2:length(times)
      time_ = (times[k] + times[k-1]) / 2.0
      correlations[k] = exponential_correlations(rateTimes, longTermCorr, beta, gamma, time)
    end
  end

  return ExponentialForwardCorrelation(numberOfRates, longTermCorr, beta, gamma, rateTimes, times, correlations)
end

clone(efc::ExponentialForwardCorrelation) = ExponentialForwardCorrelation(efc.numberOfRates, efc.longTermCorr, efc.beta, efc.gamma, copy(efc.rateTimes), copy(efc.times), deepcopy(efc.correlations))


text:
function replicate(;
  refreshcomp::Bool = true,
  refreshtretcrsp::Bool = true,
  refreshtvolmcapcrsp::Bool = true,
  refreshdataseries::Bool = true,
  refreshportfolios::Bool = true,
  parallel::NBool = nothing,
  trialrunonly::Bool = false,
  table1=true)::Nothing

  local panel::DataFrame

  (!isnothing(parallel)) && (PARALLEL[] = parallel)
  REPLICATION_TYPE[] = trialrunonly ? :monthly : :daily

  #not strictly necessary, but helps stop doubling up on the multi-threading

  #this encompasses all pre-processing activities
  if sum([refreshcomp, refreshtretcrsp, refreshtvolmcapcrsp,
    refreshdataseries, refreshportfolios]) > 0

    panel = constructportfolios(
      refreshcomp = refreshcomp,
      refreshtvolmcapcrsp = refreshtvolmcapcrsp,
      refreshtretcrsp = refreshtretcrsp,
      refreshdataseries=refreshdataseries,
      refreshportfolios=refreshportfolios,
      trialrunonly = trialrunonly)
  end

  table1 && maketable1(trialrunonly=trialrunonly)

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function test_decimal_decimalperiod()
    @test [decimalperiod(1, i) for i in 1:10] == [0, 0, 1, 0, 0, 1, 6, 0, 1, 0]
end

function test_decimal_ispandigital()
    @test ispandigital(0) == false
    @test ispandigital(1) == true
    @test ispandigital(2) == false
    @test ispandigital(10) == false
    @test ispandigital(11) == false
    @test ispandigital(12) == true
    @test ispandigital(123) == true
    @test ispandigital(124) == false
    @test ispandigital(123456789) == true
    @test ispandigital(1234567890) == false
end

function test_decimal_all()
    print(rpad("Math.NumberTheory.Decimal...", 50, ' '))

    test_decimal_decimalperiod()
    test_decimal_ispandigital()

    println("PASS")
end


text:
#=##############################################################################################
Filename: blown_wing_template.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this is a template file. Convenience methods are provided to calculate the lift
        distribution of a blown wing over the specified alphas.
=###############################################################################################

# initialize parameters
"""
    VLM_BEM <: Parameters

Parameters for a blown wing system.

# Fields:

* `omegas::Array{TF,2}`:
* `Js::Array{TF,2}`:
* `Ts::Array{TF,2}`:
* `Qs::Array{TF,2}`:
* `CTs::Array{TF,2}`:
* `CQs::Array{TF,2}`:
* `ηs::Array{TF,2}`:
* `us::Vector{Vector{Vector{TF}}}`:
* `vs::Vector{Vector{Vector{TF}}}`:
* `wake_function::Vector{Function}`:
* `wake_shape_functions::Vector{Function}`:
* `axial_interpolations::Vector{Function}`:
* `swirl_interpolations::Vector{Function}`:
* `axial_multipliers::Vector{Function}`:
* `swirl_multipliers::Vector{Function}`:
* `CFs::Array{TF,2}`: [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{TF,2}`: [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step
* `cfs::Vector{Vector{Array{TF,2}}}`:
* `cms::Vector{Vector{Array{TF,2}}}`:
* `surfacenames::Vector{String}`:
* `rotor_names::Vector{String}`:
* `plot_directory::String`:
* `plot_base_name::String`:
* `plot_extension::String`:
* `plotstepi::TR`:

"""
struct VLM_BEM{TF,TR} <: Parameters
    omegas::Array{TF,2}
    Js::Array{TF,2}
    Ts::Array{TF,2}
    Qs::Array{TF,2}
    Ps::Array{TF,2}
    CTs::Array{TF,2}
    CQs::Array{TF,2}
    ηs::Array{TF,2}
    us::Vector{Vector{Vector{TF}}}
    vs::Vector{Vector{Vector{TF}}}
    wake_function::Vector{Function}
    wake_shape_functions::Vector{Function}
    axial_interpolations::Vector{Function}
    swirl_interpolations::Vector{Function}
    axial_multipliers::Vector{Function}
    swirl_multipliers::Vector{Function}
    CFs::Array{TF,2}
    CMs::Array{TF,2}
    cfs::Vector{Vector{Array{TF,2}}}
    cms::Vector{Vector{Array{TF,2}}}
    surfacenames::Vector{String}
    rotor_names::Vector{String}
    plot_directory::String
    plot_base_name::String
    plot_extension::String
    plotstepi::TR
end

"""
    ParamsSolveVLMBEM{TF,TFUN} <: Parameters

Bare-bones `<:Parameters` struct for the `solve_vlm_bem` action.

"""
struct ParamsSolveVLMBEM{TF} <: Parameters
    omegas::Array{TF,2}
    Js::Array{TF,2}
    Ts::Array{TF,2}
    Qs::Array{TF,2}
    Ps::Array{TF,2}
    CTs::Array{TF,2}
    CQs::Array{TF,2}
    ηs::Array{TF,2}
    us::Vector{Vector{Vector{TF}}}
    vs::Vector{Vector{Vector{TF}}}
    wake_function::Vector{Function}
    wake_shape_functions::Vector{Function}
    axial_interpolations::Vector{Function}
    swirl_interpolations::Vector{Function}
    axial_multipliers::Vector{Function}
    swirl_multipliers::Vector{Function}
    CFs::Array{TF,2}
    CMs::Array{TF,2}
    cfs::Vector{Vector{Array{TF,2}}}
    cms::Vector{Vector{Array{TF,2}}}
end

# function (rs::RotorSweepParameters{V1, V2, V3, V4, V5})(alphas, omega...)

# end

# rs  = RotorSweepParameter(args...)

# rs(alphas, omega, nblades...)

"""
    vlm_bem_template(vinfs, plotstepi, alphas,
        wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip,
        omegas, nblades, rhub, rtip, radii, chords, twists,
        airfoilcontours, airfoilnames, index,
        rotor_positions, rotor_orientation, spin_directions, Res_list, Ms_list;
            wake_developement_factor = 1.0, # fully developed by default
            swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
            surfacenames = ["default wing"],
            rotor_names = ["rotor 1"],
            plot_directory = joinpath(topdirectory, "data","plots",TODAY),
            plot_base_name = "default",
            plot_extension = ".pdf",
            step_symbol = L"\alpha ",
            kwargs...
    ) <: Template

Template function returns inputs for a VLM + BEM simulation.

# Arguments:

* `vinfs`
* `plotstepi`
* `alphas::Vector{Float64}`: angles of attack
* `wing_b::Float64`: wing span
* `wing_TR::Float64`: wing taper ratio
* `wing_AR::Float64`: wing aspect ratio
* `wing_θroot::Float64`: wing root twist
* `wing_θtip::Float64`: wing tip twist
* `omegas`
* `nblades::Int`
* `rhub`
* `rtip`
* `radii`
* `chords`
* `twists`
* `airfoilcontours`
* `airfoilnames`
* `index`
* `rotor_position`
* `rotor_orientation`
* `spin_directions`
* `Res_list`
* `Ms_list`

# Keyword Arguments:

* `wake_developement_factor` : value between 0-1 with 0 meaning not developed at all and 1 meaning fully developed
* `surfacenames` : names for each lifting surface used in plots
* `rotor_names` : names for each rotor used in plots
* `plot_directory` : directory to which plots are saved
* `plot_base_name` : first part of saved plot filenames
* `plot_extension` : saved plots file extension
* `step_symbol` : symbol in plots and terminal output describing each step

"""
function vlm_bem_template(vinfs, plotstepi, alphas, omegas,
    wing_system::VortexLatticeSystem,
    rotor_system::CCBladeSystem;
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames = ["default wing"],
        rotor_names = ["rotor 1"],
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"\alpha ",
        kwargs...
)
    # prepare subsystems
    wings = wing_system
    rotors = rotor_system
    inertia = Inertia(mass, inertia_x, inertia_y, inertia_z)
    nonliftingbodies = nothing
    structures = nothing
    motors = nothing
    batteries = nothing

    # build system struct
    aircraft = Aircraft(wings, rotors, inertia, nonliftingbodies, structures, motors, batteries)

    # compile actions
    actions = [solve_vlm_bem]

    # initialize parameters
    step_range = alphas
    params_solve_vlm_bem = solve_vlm_bem(aircraft, step_range)
    params_solve_vlm_bem[1] .*= omegas

    # set wake development parameters
    @assert wake_developement_factor <= 1 && wake_developement_factor >= 0 "`wake_developement_factor` must be between 0 and 1"
    axial_multipliers = Function[(distance2plane, Rtip) -> 2 for i in 1:length(aircraft.rotor_system.index)]
    swirl_multipliers = Function[(distance2plane, Rtip) -> 1 for i in 1:length(aircraft.rotor_system.index)]
    params_solve_vlm_bem = params_solve_vlm_bem[1:14]..., axial_multipliers, swirl_multipliers, params_solve_vlm_bem[17:end]...
    #= Contains:
        1 `omegas::Vector{Float64}` : a vector of rotational speeds in rad/s at the current step
        2 `Js::Array{Float64,2}` : each [i,j]th element is the advance ratio of the ith rotor at the jth step
        3 `Ts::Array{Float64,2}` : each [i,j]th element is the thrust of the ith rotor at the jth step
        4 `Qs::Array{Float64,2}` : each [i,j]th element is the torque of the ith rotor at the jth step
        5 `CTs::Array{Float64,2}` : each [i,j]th element is the thrust coefficient of the ith rotor at the jth step
        6 `CQs::Array{Float64,2}` : each [i,j]th element is the torque coefficient of the ith rotor at the jth step
        7 `ηs::Array{Float64,2}` : each [i,j]th element is the propulsive efficiency of the ith rotor at the jth step
        8 `us::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the axial induced velocity at ith step of the jth rotor at the kth radial section
        9 `vs::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the swirl induced velocity at ith step of the jth rotor at the kth radial section
        10 `wake_function::Vector{Function}` : function accepts position X and returns the rotor induced velocity
        11 `wake_shape_functions::Vector{Function}` : [i]th element is a function f(Rtip, x) describing the radial distance from the rotor axis to the boundary of the wake of the ith rotor
        12 `axial_interpolations::Vector{Function}` : [i]th element is a function f(rs, us, r, Rtip) that returns the axial component of rotor-induced velocity at distance r from the rotor axis based on the calculated axial induced velocities output from CCBlade of the ith rotor
        13 `swirl_interpolations::Vector{Function}` : [i]th element is a function f(rs, us, r, Rtip) that returns the swirl component of rotor-induced velocity at distance r from the rotor axis based on the calculated axial induced velocities output from CCBlade of the ith rotor
        14 `axial_multipliers::Vector{Function}` : [i]th element is a function f(distance2plane, Rtip) that is multiplied by the axial induced velocity function of the ith rotor
        15 `swirl_multipliers::Vector{Function}` : [i]th element is a function f(distance2plane, Rtip) that is multiplied by the swirl induced velocity function of the ith rotor
        16 `CFs::Array{TF,2}`: [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
        17 `CMs::Array{TF,2}`: [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth
        18 `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
        19 `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface
    =#

    params_post_plot_cf_cm_alpha_sweep = post_plot_cf_cm_alpha_sweep(aircraft, step_range)
    #= Contains:
        1 `CFs::Array{TF,2}`: [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
        2 `CMs::Array{TF,2}`: [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth
        3 `plot_directory::String` : directory where plots are saved
        4 `plot_base_name::String` : first portion of the saved figure file name
        5 `plot_extension::String` : extension of saved figure files
    =#

    params_post_plot_lift_moment_distribution = post_plot_lift_moment_distribution(aircraft, step_range)
    #= Contains:
        1 `cfs::Vector{Array{Float64,2}}` : vector with length equal to the number of lifting surfaces, each member containing an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl
        2 `cms::Vector{Array{Float64,2}}` : vector with length equal to the number of lifting surfaces, each member containing an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy,
        3 `surfacenames::Vector{String}` : names of each lifting surface to be shown in the legend
        4 `plot_directory::String` : path to the folder where plots will be saved
        5 `plot_base_name::String` : first part of saved figure file names
        6 `plot_extension::String` : extension of saved figure file names
        7 `plotstepi::Vector{Int}` : which steps at which to plot
    =#

    if !isnothing(surfacenames); params_post_plot_lift_moment_distribution[3] .= surfacenames; end

    params_post_plot_rotor_sweep = post_plot_rotor_sweep(aircraft, step_range)
    #= Contains:
        1 `Js::AbstractArray` : array of advance ratios of the ith rotor at the jth step
        2 `CTs::Array{Float64,2}` : array of thrust coefficients of the ith rotor at the jth step
        3 `CQs::Array{Float64,2}` : array of torque coefficients of the ith rotor at the jth step
        4 `ηs::Array{Float64,2}` : array of propulsive efficiencies of the ith rotor at the jth step
        5 `rotor_names::Vector{String}` : vector of rotor names for use in plot legend
        6 `plot_directory::String` : directory where plots are saved
        7 `plot_base_name::String` : first portion of the saved figure file name
        8 `plot_extension::String` : extension of saved figure files
    =#

    if !isnothing(rotor_names); params_post_plot_rotor_sweep[5] .= rotor_names; end
    # prepare plot directory
    if !isdir(plot_directory); mkpath(plot_directory); end
    println("==== MSG ====\n\tplot_directory = $plot_directory")

    # build parameters struct
    parameters = VLM_BEM(params_solve_vlm_bem..., params_post_plot_lift_moment_distribution[3], params_post_plot_rotor_sweep[5], plot_directory, plot_base_name, plot_extension, plotstepi)

    # build freestream_function
    function freestream_function(aircraft, parameters, environment, alphas, stepi)

        # calculate freestream
        vinf = vinfs[stepi] # + ti # arbitrary for lift distribution?
        alpha = alphas[stepi]
        beta = 0.0
        Omega = StaticArrays.@SVector zeros(3)
        freestream = Freestream(vinf, alpha, beta, Omega)

        # update reference
        # vwake = parameters.wake_function[stepi]
        reference = aircraft.wing_system.system.reference[1]
        S = reference.S
        c = reference.c
        b = reference.b
        r = reference.r
        newV = freestream.vinf
        newreference = VL.Reference(S,c,b,r,newV)
        aircraft.wing_system.system.reference[1] = newreference

        return freestream
    end

    # build environment_function
    function environment_function(aircraft, parameters, alphas, stepi)
        Environment() # arbitrary
    end

    # compile postactions
    postactions = [post_plot_cf_cm_alpha_sweep, post_plot_lift_moment_distribution, post_plot_rotor_sweep]

    # build objective_function
    objective_function(aircraft, parameters, freestream, environment, alphas) = 0.0

    return aircraft, parameters, actions, freestream_function, environment_function, postactions, objective_function, alphas, step_symbol
end

function vlm_bem_template(vinfs, plotstepi, alphas,
    wing_system::VortexLatticeSystem,
    omegas, nblades_list, rhub_list, rtip_list, radii_list, chords_list, twists_list,
    contour_paths_list, index,
    rotor_positions, rotor_orientations, spin_directions, Res_lists, Ms_lists;
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames = ["default wing"],
        rotor_names = ["rotor 1"],
        mach_correction = nothing,
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"\alpha ",
        kwargs...
)
    # prepare subsystems
    rotors = CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list,
        Res_lists, Ms_lists, contour_paths_list,
        index, rotor_positions, rotor_orientations, spin_directions;
            kwargs...
    )
    # rotors = CCBladeSystem(nblades, rhub, rtip, radii, chords, twists, contour_paths_list, index, rotor_positions, rotor_orientation, spin_directions, Res_list, Ms_list; kwargs...)

    return vlm_bem_template(vinfs, plotstepi, alphas, omegas,
        wing_system, rotors;
            wake_developement_factor, # fully developed by default
            swirl_recovery_factor, # as described in Veldhuis' paper
            surfacenames,
            rotor_names,
            mass,
            inertia_x,
            inertia_y,
            inertia_z,
            plot_directory,
            plot_base_name,
            plot_extension,
            step_symbol,
            kwargs...
    )
end

function vlm_bem_template(vinfs, plotstepi, alphas,
    wing_system::VortexLatticeSystem,
    omegas, nblades_list, rhub_list, rtip_list, radii_list, chords_list, twists_list,
    contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists, index,
    rotor_positions, rotor_orientations, spin_directions, Res_lists, Ms_lists;
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames = ["default wing"],
        rotor_names = ["rotor 1"],
        mach_correction = nothing,
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"\alpha ",
        kwargs...
)
    # prepare subsystems
    rotors = CCBladeSystem(nblades_list, rhub_list, rtip_list, mach_correction,
        radii_list, chords_list, twists_list,
        Res_lists, Ms_lists, contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists,
        index, rotor_positions, rotor_orientations, spin_directions;
            kwargs...
    )
    # rotors = CCBladeSystem(nblades, rhub, rtip, radii, chords, twists, contour_paths_list, index, rotor_positions, rotor_orientation, spin_directions, Res_list, Ms_list; kwargs...)

    return vlm_bem_template(vinfs, plotstepi, alphas, omegas,
        wing_system, rotors;
            wake_developement_factor, # fully developed by default
            swirl_recovery_factor, # as described in Veldhuis' paper
            surfacenames,
            rotor_names,
            mass,
            inertia_x,
            inertia_y,
            inertia_z,
            plot_directory,
            plot_base_name,
            plot_extension,
            step_symbol,
            kwargs...
    )
end

function vlm_bem_template(vinfs, plotstepi, alphas,
    wing_lexs, wing_leys, wing_lezs, wing_chords, wing_thetas, wing_phis, wing_nspanwisepanels, wing_nchordwisepanels,
    omegas, nblades_list, rhub_list, rtip_list, radii, chords_list, twists_list,
    contour_paths_list, index,
    rotor_positions, rotor_orientations, spin_directions, Res_lists, Ms_lists;
        Vref = 1.0, symmetric = true, iref = 1, static_margin = 0.10, liftingline_x_over_c = 0.25,
        cambers = [fill((xc) -> 0, length(lex)) for lex in wing_lexs],
        spanwise_spacings = fill(VL.Uniform(), length(wing_lexs)),
        chordwise_spacings = fill(VL.Uniform(), length(wing_lexs)),
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames = ["default wing"],
        rotor_names = ["rotor 1"],
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"\alpha ",
        kwargs...
)
    # prepare subsystems
    wings = VortexLatticeSystem(wing_lexs, wing_leys, wing_lezs, wing_chords, wing_thetas, wing_phis, wing_nspanwisepanels, wing_nchordwisepanels; Vref, symmetric, iref, static_margin, liftingline_x_over_c, cambers, spanwise_spacings, chordwise_spacings)

    return vlm_bem_template(vinfs, plotstepi, alphas,
        wings,
        omegas, nblades_list, rhub_list, rtip_list, radii, chords_list, twists_list,
        contour_paths_list, index,
        rotor_positions, rotor_orientations, spin_directions, Res_lists, Ms_lists;
            wake_developement_factor, # fully developed by default
            swirl_recovery_factor, # as described in Veldhuis' paper
            surfacenames,
            rotor_names,
            mass,
            inertia_x,
            inertia_y,
            inertia_z,
            plot_directory,
            plot_base_name,
            plot_extension,
            step_symbol,
            kwargs...
    )
end

function vlm_bem_template(vinfs, plotstepi, alphas,
    wing_lexs, wing_leys, wing_lezs, wing_chords, wing_thetas, wing_phis, wing_nspanwisepanels, wing_nchordwisepanels,
    omegas, nblades_list, rhub_list, rtip_list, radii, chords_list, twists_list,
    contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists, index,
    rotor_positions, rotor_orientations, spin_directions, Res_lists, Ms_lists;
        Vref = 1.0, symmetric = true, iref = 1, static_margin = 0.10, liftingline_x_over_c = 0.25,
        cambers = [fill((xc) -> 0, length(lex)) for lex in wing_lexs],
        spanwise_spacings = fill(VL.Uniform(), length(wing_lexs)),
        chordwise_spacings = fill(VL.Uniform(), length(wing_lexs)),
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames = ["default wing"],
        rotor_names = ["rotor 1"],
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"\alpha ",
        kwargs...
)
    # prepare subsystems
    wings = VortexLatticeSystem(wing_lexs, wing_leys, wing_lezs, wing_chords, wing_thetas, wing_phis, wing_nspanwisepanels, wing_nchordwisepanels; Vref, symmetric, iref, static_margin, liftingline_x_over_c, cambers, spanwise_spacings, chordwise_spacings)

    return vlm_bem_template(vinfs, plotstepi, alphas,
        wings,
        omegas, nblades_list, rhub_list, rtip_list, radii, chords_list, twists_list,
        contour_paths_list, uncorrected_polar_paths_lists, corrected_polar_paths_lists, index,
        rotor_positions, rotor_orientations, spin_directions, Res_lists, Ms_lists;
            wake_developement_factor, # fully developed by default
            swirl_recovery_factor, # as described in Veldhuis' paper
            surfacenames,
            rotor_names,
            mass,
            inertia_x,
            inertia_y,
            inertia_z,
            plot_directory,
            plot_base_name,
            plot_extension,
            step_symbol,
            kwargs...
    )
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
#Blosc1 low-level main API https://c-blosc2.readthedocs.io/en/latest/reference/blosc1.html#main-api

"""
    blosc_init()

Initialize the Blosc library environment.

You must call this previous to any other Blosc call, unless you want Blosc to be used simultaneously in a multi-threaded environment, in which case you can use the blosc2_compress_ctx blosc2_decompress_ctx pair.
"""
blosc_init() = @ccall lib.blosc_init()::Cvoid
"""
    blosc_destroy()

Destroy the Blosc library environment.

You must call this after to you are done with all the Blosc calls, unless you have not used blosc_init() before.
"""
blosc_destroy() = @ccall lib.blosc_destroy()::Cvoid

"""
    blosc_compress(level, shuffle, itemsize, srcsize, src, dest, destsize)

Compress a block of data in the src buffer and returns the size of compressed block.

C signature: `int blosc_compress(int clevel, int doshuffle, size_t typesize, size_t nbytes, const void *src, void *dest, size_t destsize)``

# Arguments

- `clevel` - The desired compression level and must be a number
            between 0 (no compression) and 9 (maximum compression).

- `doshuffle` - Specifies whether the shuffle compression preconditioner
    should be applied or not. `BLOSC_NOFILTER` means not applying filters,
    `BLOSC_SHUFFLE` means applying shuffle at a byte level and
    `BLOSC_BITSHUFFLE` at a bit level (slower but *may* achieve better
    compression).

- `typesize` - Is the number of bytes for the atomic type in binary
    src` buffer.  This is mainly useful for the shuffle preconditioner.
    For implementation reasons, only a 1 < typesize < 256 will allow the
    shuffle filter to work.  When typesize is not in this range, shuffle
    will be silently disabled.

- `nbytes` - The number of bytes to compress in the `src` buffer.

- `src` -  The buffer containing the data to compress.
- `dest`  - The buffer where the compressed data will be put,
    must have at least the size of `destsize`.

- `destsize` -  The size of the dest buffer. Blosc
    guarantees that if you set `destsize` to, at least,
    `BLOSC_MAX_OVERHEAD`, the compression will always succeed.

# Returns

The number of bytes compressed. If src buffer cannot be compressed into destsize, the return value is zero and you should discard the contents of the dest buffer.
"""
function blosc_compress(clevel, shuffle, itemsize, srcsize, src, dest, destsize)
    @ccall lib.blosc_compress(
            clevel::Cint, shuffle::Cint, itemsize::Csize_t, srcsize::Csize_t, src::Ptr{Cvoid}, dest::Ptr{Cvoid}, destsize::Csize_t
        )::Cint
end

"""
    blosc_decompress(src, dest, destsize)

Decompress a block of compressed data in src, put the result in dest and returns the size of the decompressed block.

C signature: `int blosc_decompress(const void *src, void *dest, size_t destsize)`

# Arguments

- `src`  - The buffer to be decompressed.
- `dest` - The buffer where the decompressed data will be put.
- `destsize` = The size of the @p dest buffer

# Returns
The number of bytes decompressed. If an error occurs, e.g. the compressed data is corrupted or the output buffer is not large enough, then a negative value will be returned instead.

"""
function blosc_decompress(src, dest, destsize)
    @ccall lib.blosc_decompress(src::Ptr{Cvoid}, dest::Ptr{Cvoid}, destsize::Csize_t)::Cint
end

"""
    blosc_getitem(src, start, nitems, dest)

Get `nitems` (of `typesize` size) in `src` buffer starting in `start`.
The items are returned in `dest` buffer, which has to have enough
space for storing all items.

C signature: `int blosc_getitem(const void *src, int start, int nitems, void *dest)`

#Arguments

- `src` - The compressed buffer from data will be decompressed.
- `start` - The zero-indexed position of the first item (of `typesize` size) from where data will be retrieved.
- `nitems` - The number of items (of `typesize` size) that will be retrieved.
- `dest` -  The buffer where the decompressed data retrieved will be put.

# Returns
The number of bytes copied to @p dest or a negative value if some error happens.
"""
function blosc_getitem(src, start, nitems, dest)
    @ccall lib.blosc_getitem(src::Ptr{Cvoid}, start::Cint, nitems::Cint, dest::Ptr{Cvoid})::Cint
end

"""
    blosc_get_nthreads()

Returns the current number of threads that are used for compression/decompression.
"""
function blosc_get_nthreads()
    @ccall lib.blosc_get_nthreads()::Cint
end

"""
    blosc_set_nthreads(nthreads)

Initialize a pool of threads for compression/decompression.
If nthreads is 1, then the serial version is chosen and a possible previous existing pool is ended. If this is not called, nthreads is set to 1 internally.

# Returns

The previous number of threads.
"""
function blosc_set_nthreads(nthreads)
    @ccall lib.blosc_set_nthreads(nthreads::Cint)::Cint
end

"""
    blosc_get_compressor()

Get the current compressor that is used for compression.
"""
function blosc_get_compressor()
    unsafe_string(@ccall lib.blosc_get_compressor()::Cstring)
end

"""
    blosc_set_compressor(compname::AbstractString)

Select the compressor to be used.

The supported ones are “blosclz”, “lz4”, “lz4hc”, “zlib” and “ztsd”. If this function is not called, then “blosclz” will be used.

#Returns

The code for the compressor (>=0). In case the compressor is not recognized, or there is not support for it in this build, it returns a -1.
"""
function blosc_set_compressor(compname::AbstractString)
    @ccall lib.blosc_set_compressor(compname::Cstring)::Cint
end

"""
    blosc_set_delta(dodelta)

Select the delta coding filter to be used.

This call should always succeed.

#Arguments

- `dodelta` – A value >0 will activate the delta filter. If 0, it will be de-activated

"""
function blosc_set_delta(dodelta)
    @ccall lib.blosc_set_delta(dodelta::Cint)::Cvoid
end

"""
    blosc_set_blocksize(blocksize)

Force the use of a specific blocksize.

If 0, an automatic blocksize will be used (the default).

# Warning

The blocksize is a critical parameter with important restrictions in the allowed values, so use this with care.
"""
function blosc_set_blocksize(blocksize)
    @ccall lib.blosc_set_blocksize(blocksize::Csize_t)::Cvoid
end

"""
    blosc_free_resources()

Free possible memory temporaries and thread resources.

Use this when you are not going to use Blosc for a long while.

# Returns
A 0 if succeeds, in case of problems releasing the resources, it returns a negative number.

"""
function blosc_free_resources()
    @ccall lib.blosc_free_resources()::Cint
end


text:
module Correlogram
using StatsBase

export call

function call(values)
  lags = 0:(size(values, 1)-1)
  vcat(hcat(lags...), hcat(StatsBase.autocor(values, lags)...))
end

end


text:
function area(x1,y1,z1,x2,y2,z2,x3,y3,z3)

    a = sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)
    b = sqrt((x2-x3)^2+(y2-y3)^2+(z2-z3)^2)
    c = sqrt((x3-x1)^2+(y3-y1)^2+(z3-z1)^2)
    p = (a+b+c)/2
    s = sqrt(p*(p-a)*(p-b)*(p-c))

    return s
end

max_stars_repo_path: src/single-core juSFEM/area.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/single-core juSFEM/area.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/single-core juSFEM/area.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
@testset "WordleGuess" begin
    WordleGuess = Wordle.WordleGuess

    @test all(hasfield.(Wordle.WordleGuess, [:guess, :result]))
    @test_throws ErrorException WordleGuess("", [:some, :random, :symbols])
    @test_throws ErrorException WordleGuess("", [:🟩, :🟩, :🟩])
    @test (@capture_out print(WordleGuess("", fill(Wordle.CORRECT, 5)))) == "🟩🟩🟩🟩🟩"
end


text:
function UnitSphereRandomVectorGenerator(arg0::jint)
    return UnitSphereRandomVectorGenerator((jint,), arg0)
end

function UnitSphereRandomVectorGenerator(arg0::jint, arg1::RandomGenerator)
    return UnitSphereRandomVectorGenerator((jint, RandomGenerator), arg0, arg1)
end

function next_vector(obj::UnitSphereRandomVectorGenerator)
    return jcall(obj, "nextVector", Vector{jdouble}, ())
end


text:
using Statistics
using NCDatasets

using ArgParse
using JSON

include("LinearRegression.jl")

correlation = (x1, x2) -> x1' * x2 / (sum(x1.^2)*sum(x2.^2)).^0.5

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Ocean data file. New variable will be appended."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--SST"
            help = "Variable name of SST."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["data-file"], "r") do ds

    dims = size(ds[parsed["SST"]])

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    if length(dims) == 3 # 2D + time
        rng = (:,:,beg_t:end_t)
    elseif length(dims) == 4 # 2D + time
        rng = (:,:,1,beg_t:end_t)
    else
        ErrorException("Unknown dimension") |> throw
    end

    global SST  = replace(ds[parsed["SST"]][rng...], missing=>NaN)

    global (Nx, Ny, Nt) = size(SST)

    if mod(Nt, 12) != 0

        SST = SST[:, :, 1:12*floor(Integer, Nt/12.0)]
        (Nx, Ny, Nt) = size(SST)

        #ErrorException("Time record is not multiple of 12") |> throw
    end

    global nyears = Int64(Nt / 12)
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
end

# MM   = Monthly Mean
# A    = Anomaly
# YYC  = Year-to-Year Correlation

SSTMM   = zeros(Float64, size(SST)[1:2]..., 12)
SSTA    = zeros(Float64, size(SST)...)
SSTAYYC = zeros(Float64, size(SST)[1:2]..., 12)
SSTAVAR = zeros(Float64, size(SST)[1:2]..., 12)

N = size(SST)[3]
x = collect(Float64, 1:N)
for i=1:Nx, j=1:Ny

    d = detrend(x, view(SST, i, j, :))

    if mask[i, j] != 0
        SSTMM[i, j, :]  .= NaN
        SSTA[i, j, :]   .= NaN
        SSTAYYC[i, j, :] .= NaN
        continue
    end

    SSTMM[i, j, :] = mean( reshape(d, 12, :), dims=2 )[:, 1]
    SSTA[i, j, :]  = d - repeat( SSTMM[i, j, :], outer=nyears)

    for m = 1:12
        d_yy = view(SSTA, i, j, m:12:(m+nyears*12-1))
        SSTAYYC[i, j, m] = correlation(d_yy[1:end-1], d_yy[2:end])
        SSTAVAR[i, j, m] = std(d_yy)
    end

end

Dataset(parsed["data-file"], "a") do ds
    if ! haskey(ds.dim, "months")
        defDim(ds, "months", 12)
    end
end

Dataset(parsed["data-file"], "a") do ds
    for (varname, vardata, vardim, attrib) in [
        ("SSTMM",  SSTMM, ("Nx", "Ny", "months"), Dict()),
        ("SSTA",    SSTA, ("Nx", "Ny", "time"),   Dict()),
        ("SSTAYYC",  SSTAYYC, ("Nx", "Ny", "months"), Dict()),
        ("SSTAVAR",  SSTAVAR, ("Nx", "Ny", "months"), Dict()),
    ]

        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic