text:
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of Intel Corporation nor the names of its contributors
#       may be used to endorse or promote products derived from this software
#       without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

using Latte

net = Net(128)
data, label  = HDF5DataLayer(net, "data/train.txt", "data/test.txt")

conv1_1 = ConvolutionLayer( :conv1_1, net, data,    64,  3, 1, 1)
relu1_1 = ReLULayer(        :relu1_1, net, conv1_1)
conv1_2 = ConvolutionLayer( :conv1_2, net, relu1_1, 64,  3, 1, 1)
relu1_2 = ReLULayer(        :relu1_2, net, conv1_2)
pool1   = MaxPoolingLayer(  :pool1,   net, relu1_2, 2,   2, 0; ceil_mode=true)

conv2_1 = ConvolutionLayer( :conv2_1, net, pool1,   128, 3, 1, 1)
relu2_1 = ReLULayer(        :relu2_1, net, conv2_1)
conv2_2 = ConvolutionLayer( :conv2_2, net, relu2_1, 128, 3, 1, 1)
relu2_2 = ReLULayer(        :relu2_2, net, conv2_2)
pool2   = MaxPoolingLayer(  :pool2,   net, relu2_2, 2,   2, 0; ceil_mode=true)

conv3_1 = ConvolutionLayer( :conv3_1, net, pool2,   256, 3, 1, 1)
relu3_1 = ReLULayer(        :relu3_1, net, conv3_1)
conv3_2 = ConvolutionLayer( :conv3_2, net, relu3_1, 256, 3, 1, 1)
relu3_2 = ReLULayer(        :relu3_2, net, conv3_2)
conv3_3 = ConvolutionLayer( :conv3_3, net, relu3_2, 256, 3, 1, 1)
relu3_3 = ReLULayer(        :relu3_3, net, conv3_3)
pool3   = MaxPoolingLayer(  :pool3,   net, relu3_3, 2,   2, 0; ceil_mode=true)

conv4_1 = ConvolutionLayer( :conv4_1, net, pool3,   512, 3, 1, 1)
relu4_1 = ReLULayer(        :relu4_1, net, conv4_1)
conv4_2 = ConvolutionLayer( :conv4_2, net, relu4_1, 512, 3, 1, 1)
relu4_2 = ReLULayer(        :relu4_2, net, conv4_2)
conv4_3 = ConvolutionLayer( :conv4_3, net, relu4_2, 512, 3, 1, 1)
relu4_3 = ReLULayer(        :relu4_3, net, conv4_3)
pool4   = MaxPoolingLayer(  :pool4,   net, relu4_3, 2,   2, 0; ceil_mode=true)

conv5_1 = ConvolutionLayer( :conv5_1, net, pool4,   512, 3, 1, 1)
relu5_1 = ReLULayer(        :relu5_1, net, conv5_1)
conv5_2 = ConvolutionLayer( :conv5_2, net, relu5_1, 512, 3, 1, 1)
relu5_2 = ReLULayer(        :relu5_2, net, conv5_2)
conv5_3 = ConvolutionLayer( :conv5_3, net, relu5_2, 512, 3, 1, 1)
relu5_3 = ReLULayer(        :relu5_3, net, conv5_3)
pool5   = MaxPoolingLayer(  :pool5,   net, relu5_3, 2,   2, 0; ceil_mode=true)

fc6     = InnerProductLayer(:fc6,     net, pool5, 512)
relu6   = ReLULayer(        :relu6,   net, fc6)
fc7     = InnerProductLayer(:fc7,     net, relu6,   10)

loss     = SoftmaxLossLayer(:loss, net, fc7, label)
accuracy = AccuracyLayer(:accuracy, net, fc7, label)

params = SolverParameters(
    lr_policy    = LRPolicy.Decay(.01f0, 5.0f-7),
    mom_policy   = MomPolicy.Fixed(0.9),
    max_epoch    = 300,
    regu_coef    = .0005,
    snapshot_dir = "no_dropout")
sgd = SGD(params)
solve(sgd, net)


text:
# # Theory

# The following theory described is meant to be minimal, with mainly the equations presented for reference when reading the code.

# ## Geometry

# ### Wing Parametrization
#
# A **half-wing** is defined in terms of a nested trapezoid. A single trapezoid is called a **section**. A section consists of two foil profiles, and their associated chord lengths and twist angles. Between them is their span length with associated _leading-edge_ dihedral and sweep angles. So a general half-wing consisting of ``n`` sections will have ``n`` entries for foils, chords, and twists, and ``n - 1`` entries for spans, dihedrals, sweeps for some ``n \in \mathbb N``. The following illustration should help visualize the concept.
#
# ![](https://godot-bloggy.xyz/post/diagrams/WingGeometry.svg)

# ## Aerodynamics

# The aerodynamic analyses in AeroMDAO mainly utilize potential flow theory and solve problems using a _boundary element method_. This essentially is the following Laplace equation problem with the following Robin (?) boundary conditions:

# ```math
# \nabla^2 \phi = 0, \quad \mathbf V \equiv \nabla \phi \cdot \hat{\mathbf n} = 0, \quad \lim_{\mathbf r \to \infty} \phi(\mathbf r) \to 0
# ```

# !!! note
#     Implementations of viscous-inviscid coupled analyses for drag prediction (á là [XFOIL](https://web.mit.edu/drela/Public/web/xfoil/)) are in progress.

# ### Doublet-Source Panel Method

# The doublet-source panel method predicts **inviscid, incompressible, irrotational, isentropic** external flow over surfaces in 2 dimensions.

# Source and doublet singularities are placed on the surface, and boundary conditions are imposed on their induced velocity to obtain a well-posed problem. The induced velocity is evaluated by the corresponding free-field Green function for each singularity.

# ```math

# ```

# The velocities are added to obtain the total induced velocity at a point $\mathbf r$.

# ```math

# ```

# ### Vortex Lattice Method

# The vortex lattice method predicts **inviscid, incompressible, irrotational, isentropic** external flow over "thin" surfaces in 3 dimensions.

# Vortex filaments are placed on the surface, and boundary conditions are imposed on their induced velocity to obtain a well-posed problem. The induced velocity is evaluated by the Biot-Savart integral for a vortex line of length $\ell$ with a constant circulation strength $\Gamma$.

# ```math
# \mathbf V(\mathbf r, \mathbf r') = \frac{\Gamma}{4\pi} \int_0^\ell \frac{d\boldsymbol\ell' \times (\mathbf r - \mathbf r')}{|\mathbf r - \mathbf r'|^3}
# ```

# The vortices can be set up in various configurations consisting of bound or semi-infinite filaments, commonly in the form of _horseshoes_ or _vortex rings_.

# 1. _Horseshoe elements:_
#    These are defined by a finite _bound leg_ and two semi-infinite _trailing legs_. AeroMDAO encodes this information in the `Horseshoe` type.
# 2. _Vortex rings:_
#    These are defined by four _bound legs_. AeroMDAO encodes this information in the `VortexRing` type.

# A quasi-steady freestream condition with velocity $\mathbf U$ and rotation $\boldsymbol\Omega$ (in the body's frame) defines an external flow. The induced velocity at a point is given by:

# ```math
# \mathbf V_{\infty}(\mathbf r) = - (\mathbf U + \boldsymbol\Omega \times \mathbf r)
# ```

# The velocities are added to obtain the total induced velocity at a point $\mathbf r$.

# ```math
# \mathbf V(\mathbf r) = \sum_i \frac{\Gamma_i}{4\pi} \int_0^{\ell_i} \frac{d\boldsymbol\ell_i' \times (\mathbf r - \mathbf r_i')}{|\mathbf r - \mathbf r_i'|^3} + \mathbf V_\infty(\mathbf r)
# ```

# Imposing the Neumann boundary condition $\mathbf V \cdot \hat{\mathbf n} = 0$ defines the problem. The construction, in essence, fundamentally results in the following linear system to be solved:

# ```math
# \mathbf A \boldsymbol\Gamma = -V_{\infty} \cdot [\hat{\mathbf n}_i ]_{i = 1, \ldots, N}
# ```

# ### Compressibility Corrections
# The following **Prandtl-Glauert** equation, represented in wind axes, is applied for a weakly compressible flow problem ($0.3 \leq M_\infty \leq 0.7$).

# ```math
# \beta_{PG}^2\frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} + \frac{\partial^2\phi}{\partial z^2} = 0, \quad \beta_{PG}^2 = \left(1 - M_\infty^2\right)
# ```

# The Prandtl-Glauert transformation $\phi(x,y,z; \beta_{PG}) \to \bar\phi(\bar x, \bar y, \bar z)$ converts this equation into an equivalent incompressible flow problem in a transformed geometric space. This "bar" map scales the coordinates $(\bar x,\bar y, \bar z) = (x,\beta_{PG} y, \beta_{PG} z)$ and the potential $\bar\phi = \beta_{PG}^2 \phi$ in sequence. Hence the transformed equation satisfies the Laplace equation with the Neumann boundary condition:
# ```math
# \begin{aligned}
#     \bar\nabla^2 \bar\phi & = 0, \\
#     \bar{\mathbf V} \cdot \hat{\bar{\mathbf n}} & = 0,
# \end{aligned}
# ```
# where $\bar\nabla$ is differentiation with respect to the transformed coordinates and $\hat{\bar{\mathbf n}} = (\beta_{PG} \hat n_x, \hat n_y, \hat n_z)$ which can be proved by computing the appropriate cross product.

# As the circulation is a scalar (hence invariant of the coordinate transformation but not the potential scaling), the inverse is also readily derived.
# ```math
# \begin{aligned}
#     \Gamma & = \int \mathbf V \cdot d\boldsymbol\ell = \int \nabla\phi \cdot d\boldsymbol\ell, \\
#     \bar{\Gamma} & = \int \bar{\mathbf V} \cdot d\bar {\boldsymbol\ell} = \int\bar\nabla\bar\phi \cdot d\bar{\boldsymbol\ell}, \\
#     \implies \Gamma & = \bar{\Gamma} / \beta_{PG}^2
# \end{aligned}
# ```
# Hence the solution of the resultant incompressible system in transformed coordinates provides the necessary quantities of interest for calculating the dynamics.

# ## Structures

# The structural analyses in AeroMDAO utilize _linear finite-element methods_.

# Particularly, a $2$-dimensional beam element model has been implemented following the standard formulation using cubic Hermite shape functions based on **Euler-Bernoulli beam theory**. These are embedded into a $3$-dimensional local coordinate system in the vortex lattice method without loss of generality.

# The linear system consisting of the stiffness matrix $\mathbf K$ and load vector $\mathbf f$ are solved to obtain the displacement vector $\boldsymbol\delta$.

# ```math
# \mathbf K \boldsymbol\delta = \mathbf f
# ```

# ## Aeroelasticity

# The vortex lattice method and beam element model are combined into a coupled system to perform static aeroelastic analyses. The analysis is made nonlinear via promotion of the angle of attack $\alpha$ to a variable by specifying the load factor $n$ with a given weight $W$ at fixed sideslip angle $\beta$.

# Define $\mathbf x = [\boldsymbol\Gamma, \boldsymbol\delta, \alpha]$ as the state vector satisfying the residual equations:

# ```math
# \begin{aligned}
#     \mathcal R_A(\mathbf x) & = \mathbf A(\boldsymbol\delta) \boldsymbol\Gamma - \mathbf V_\infty(\alpha) \cdot [\mathbf n_i(\boldsymbol\delta)]_{i = 1,\ldots, N} \\
#     \mathcal R_S(\mathbf x) & = \mathbf K \boldsymbol\delta - \mathbf f(\boldsymbol\delta, \boldsymbol\Gamma, \alpha) \\
#     \mathcal R_L(\mathbf x) & = L(\boldsymbol\delta, \boldsymbol\Gamma, \alpha) - n W \\
# \end{aligned}
# ```

# where the lift $L$ is obtained by transforming forces computed via the Kutta-Jowkowski theorem into wind axes.

# ```math
# [D, Y, L] = \mathbf R_B^W(\alpha, \beta)\left(\sum_{i = 1}^N \rho \mathbf V_i \times  \boldsymbol\Gamma_i \boldsymbol \ell_i(\boldsymbol\delta)\right)
# ```

# and the structural load vector $\mathbf f$ is obtained from conservative and consistent load averaging of the Kutta-Jowkowski forces in geometric axes.

# ```math

# ```

# ## Flight Dynamics

# ### Longitudinal Motion

# A standard three degrees-of-freedom rigid body model is used for performing flight dynamics analyses in the 2-dimensional longitudinal plane formulated as an initial-value problem. The coupled differential equations are shown in a canonical state-space representation subject to initial conditions $\mathbf x_0$, with the time-evolution of the state vector $\mathbf x$ driven by the forcing function $\mathbf f$.

# ```math
# \begin{aligned}
#     \frac{d\mathbf x}{dt} & = \mathbf f(\mathbf x, t) \\
#     \frac{d}{dt}
#     \begin{bmatrix}
#         u_b \\
#         w_b \\
#         Q \\
#         x_e \\
#         y_e \\
#         \Theta \\
#         \delta_e \\
#         m
#     \end{bmatrix} & =
#     \begin{bmatrix}
#         -Qw_b + (T - D(V_\infty)\cos\alpha + L(V_\infty, \delta_e)\sin\alpha - W\sin\Theta) / m \\
#         Qu_b + (  - D(V_\infty)\sin\alpha - L(V_\infty, \delta_e)\cos\alpha + W\cos\Theta) / m \\
#         (M_A(\alpha, \delta_e, \hat Q) - T\Delta_{zT}) / I_{yy} \\
#         u_b \cos\Theta - w_b \sin\Theta \\
#         u_b \sin\Theta + w_b \cos\Theta \\
#         Q \\
#         g(\mathbf x, t) \\
#         -c_T T
#     \end{bmatrix}
# \end{aligned}
# ```
# The fuel burn over time (viz. reduction of mass $m$) is computed using a specific fuel consumption value with a linear dependence on the thrust. A manual controller law for the elevator deflection angle $\delta_e$ can also be implemented by providing the function $g(\mathbf x, t)$.

# ### Full Space

# A standard six degrees-of-freedom rigid body model is used for performing flight dynamics analyses in 3 dimensions.

# ```math

# ```

# ## References

# 1. Mark Drela. _Flight Vehicle Aerodynamics_. MIT Press, 2014.
# 2. Joseph Katz and Allen Plotkin. _Low-Speed Aerodynamics, Second Edition_. Cambridge University Press, 2001.

ext: jl


text:
isdefined(Base, :__precompile__) && __precompile__()

module CxxWrap

using Compat

# Convert path if it contains lib prefix on windows
function lib_path(so_path::AbstractString)
  path_copy = so_path
  @static if is_windows()
    basedir, libname = splitdir(so_path)
    libdir_suffix = Sys.WORD_SIZE == 32 ? "32" : ""
    if startswith(libname, "lib") && !isfile(so_path)
      path_copy = joinpath(basedir*libdir_suffix, libname[4:end])
    end
  end
  return path_copy
end

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, CxxWrap did not build properly")
end
include(depsfile)
const cxx_wrap_path = _l_cxx_wrap

# Base type for wrapped C++ types
abstract CppAny
abstract CppDisplay <: Display
abstract CppArray{T,N} <: AbstractArray{T,N}

# C++ std::shared_ptr
type SharedPtr{T} <: CppAny
  cpp_object::Ptr{Void}
end

# C++ std::unique_ptr
type UniquePtr{T} <: CppAny
  cpp_object::Ptr{Void}
end

# Encapsulate information about a function
type CppFunctionInfo
  name::Any
  argument_types::Array{DataType,1}
  return_type::DataType
  function_pointer::Ptr{Void}
  thunk_pointer::Ptr{Void}
end

function __init__()
  @static if is_windows()
    Libdl.dlopen(cxx_wrap_path, Libdl.RTLD_GLOBAL)
  end
  ccall((:initialize, cxx_wrap_path), Void, (Any, Any, Any), CxxWrap, CppAny, CppFunctionInfo)

  Base.linearindexing(::ConstArray) = Base.LinearFast()
  Base.size(arr::ConstArray) = arr.size
end

# Load the modules in the shared library located at the given path
function load_modules(path::AbstractString)
  module_lib = Libdl.dlopen(path, Libdl.RTLD_GLOBAL)
  registry = ccall((:create_registry, cxx_wrap_path), Ptr{Void}, ())
  ccall(Libdl.dlsym(module_lib, "register_julia_modules"), Void, (Ptr{Void},), registry)
  return registry
end

function get_module_names(registry::Ptr{Void})
  ccall((:get_module_names, cxx_wrap_path), Array{AbstractString}, (Ptr{Void},), registry)
end

function get_module_functions(registry::Ptr{Void})
  ccall((:get_module_functions, cxx_wrap_path), Array{CppFunctionInfo}, (Ptr{Void},), registry)
end

function bind_types(registry::Ptr{Void}, m::Module)
  ccall((:bind_module_types, cxx_wrap_path), Void, (Ptr{Void},Any), registry, m)
end

function exported_symbols(registry::Ptr{Void}, modname::AbstractString)
  ccall((:get_exported_symbols, cxx_wrap_path), Array{AbstractString}, (Ptr{Void},AbstractString), registry, modname)
end

# Interpreted as a constructor for Julia  > 0.5
type ConstructorFname
  _type::DataType
end

# Interpreted as an operator call overload
type CallOpOverload
  _type::DataType
end

process_fname(fn::Symbol) = fn
process_fname(fn::ConstructorFname) = :(::$(Type{fn._type}))
function process_fname(fn::CallOpOverload)
  if VERSION < v"0.5-dev"
    return :call
  end
  return :(arg1::$(fn._type))
end

make_func_declaration(fn, argmap) = :($(process_fname(fn))($(argmap...)))
function make_func_declaration(fn::CallOpOverload, argmap)
  if VERSION < v"0.5-dev"
    return :($(process_fname(fn))($(argmap...)))
  end
  return :($(process_fname(fn))($((argmap[2:end])...)))
end

# Build the expression to wrap the given function
function build_function_expression(func::CppFunctionInfo)
  # Arguments and types
  argtypes = func.argument_types
  argsymbols = map((i) -> Symbol(:arg,i[1]), enumerate(argtypes))

  # Function pointer
  fpointer = func.function_pointer
  assert(fpointer != C_NULL)

  # Thunk
  thunk = func.thunk_pointer

  map_arg_type(t::DataType) = ((t <: CppAny) || (t <: CppDisplay) || (t <: Tuple)) || (t <: CppArray) ? Any : t

  # Build the types for the ccall argument list
  c_arg_types = [map_arg_type(t) for t in argtypes]
  return_type = map_arg_type(func.return_type)

  # Build the final call expression
  call_exp = nothing
  if thunk == C_NULL
    call_exp = :(ccall($fpointer, $return_type, ($(c_arg_types...),), $(argsymbols...))) # Direct pointer call
  else
    call_exp = :(ccall($fpointer, $return_type, (Ptr{Void}, $(c_arg_types...)), $thunk, $(argsymbols...))) # use thunk (= std::function)
  end
  assert(call_exp != nothing)

  # Generate overloads for some types
  overload_map = Dict([(Cint,[Int]), (Cuint,[UInt,Int]), (Float64,[Int])])
  nargs = length(argtypes)

  counters = ones(Int, nargs);
  for i in 1:nargs
    if haskey(overload_map, argtypes[i])
        counters[i] += length(overload_map[argtypes[i]])
    end
  end

  function recurse_overloads!(idx::Int, newargs, results)
    if idx > nargs
        push!(results, deepcopy(newargs))
        return
    end
    for i in 1:counters[idx]
        newargs[idx] = i == 1 ? argtypes[idx] : overload_map[argtypes[idx]][i-1]
        recurse_overloads!(idx+1, newargs, results)
    end
  end

  newargs = Array{DataType,1}(nargs);
  overload_sigs = Array{Array{DataType,1},1}();
  recurse_overloads!(1, newargs, overload_sigs);

  function_expressions = quote end
  for overloaded_signature in overload_sigs
    argmap = Expr[]
    for (t, s) in zip(overloaded_signature, argsymbols)
      push!(argmap, :($s::$t))
    end
    func_declaration = make_func_declaration(func.name, argmap)
    push!(function_expressions.args, :($func_declaration = $call_exp))
  end
  return function_expressions
end

# Wrap functions from the cpp module to the passed julia module
function wrap_functions(functions, julia_mod)
  basenames = Set([
    :getindex,
    :setindex!,
    :convert,
    :deepcopy_internal,
    :size,
    :+,
    :*,
    :(==)
  ])
  for func in functions
    if in(func.name, basenames)
      Base.eval(build_function_expression(func))
    else
      julia_mod.eval(build_function_expression(func))
    end
  end
end

# Create modules defined in the given library, wrapping all their functions and types
function wrap_modules(registry::Ptr{Void}, parent_mod=Main)
  module_names = get_module_names(registry)
  jl_modules = Module[]
  for mod_name in module_names
    modsym = Symbol(mod_name)
    jl_mod = parent_mod.eval(:(module $modsym end))
    push!(jl_modules, jl_mod)
    bind_types(registry, jl_mod)
  end

  module_functions = get_module_functions(registry)
  for (jl_mod, mod_functions) in zip(jl_modules, module_functions)
    wrap_functions(mod_functions, jl_mod)
  end

  for (jl_mod, mod_name) in zip(jl_modules, module_names)
    exps = [Symbol(s) for s in exported_symbols(registry, mod_name)]
    jl_mod.eval(:(export $(exps...)))
  end
end

# Wrap modules in the given path
function wrap_modules(so_path::AbstractString, parent_mod=Main)
  registry = CxxWrap.load_modules(lib_path(so_path))
  wrap_modules(registry, parent_mod)
end

# Place the functions and types into the current module
function wrap_module(so_path::AbstractString, parent_mod=Main)
  registry = CxxWrap.load_modules(lib_path(so_path))
  module_names = get_module_names(registry)
  mod_idx = 0
  wanted_name = string(module_name(current_module()))
  for (i,mod_name) in enumerate(module_names)
    if mod_name == wanted_name
      bind_types(registry, current_module())
      mod_idx = i
      break
    end
  end

  if mod_idx == 0
    error("Module $wanted_name not found in C++")
  end

  module_functions = get_module_functions(registry)
  wrap_functions(module_functions[mod_idx], current_module())

  exps = [Symbol(s) for s in exported_symbols(registry, wanted_name)]
  current_module().eval(:(export $(exps...)))
end

export wrap_modules, wrap_module

end # module


text:
using Gurobi
using Test

# When adding new tests, be mindful about creating new environments. Either
# re-use an existing environment in the module, or name the test function
# `test_MULTI_ENV_xxx` to trap the specific Gurobi error indicating that an
# environment could not be created.

const GRB_ENV = Gurobi.Env()

function runtests(mod)
    for name in names(mod; all = true)
        sname = "$(name)"
        if !startswith(sname, "test_")
            continue
        end
        @testset "$(name)" begin
            if startswith(sname, "test_MULTI_ENV")
                try
                    getfield(mod, name)()
                catch ex
                    if ex == ErrorException(
                        "Gurobi Error 10009: Failed to obtain a valid license"
                    )
                        @warn(
                            "Skipping a test because there was an issue " *
                            "creating multiple licenses. This is probably " *
                            "because you have a limited license."
                        )
                    else
                        rethrow(ex)
                    end
                end
            else
                getfield(mod, name)()
            end
        end
    end
end

@testset "MathOptInterface Tests" begin
    for file in readdir("MOI")
        include(joinpath("MOI", file))
    end
end

@testset "Deprecated functions" begin
    err = ErrorException(Gurobi._DEPRECATED_ERROR_MESSAGE)
    @test_throws err get_status()
    @test_throws err Gurobi.get_status_code()
end


text:
countScenarios(TS::Array{TimeSeries.TimeArray{Float64,1,DateTime,Array{Float64,1}}, 1}) = length(TS)
function countTimeSteps(TS::Array{TimeSeries.TimeArray{Float64,1,DateTime,Array{Float64,1}}, 1})
  if countScenarios(TS) == 0
    error("No scenarios are available!")
  end
  return length(TS[1])
end


text: Genie.secret_token!("75bfcfa4046be99aaf9cc4256c90590e252552ed26d6b500dc42b760f0a0f886")

max_stars_repo_name: mazarynxyz/x-haven

max_issues_repo_name: mazarynxyz/x-haven

max_forks_repo_name: mazarynxyz/x-haven


text:
"""
    is_operation(f)

Returns a single argument anonymous function predicate, that returns `true` if and only if
the argument to the predicate satisfies `istree` and `operation(x) == f`
"""
is_operation(f) = @nospecialize(x) -> istree(x) && (operation(x) == f)
export is_operation

"""
    node_count(t)
Count the nodes in a symbolic expression tree satisfying `istree` and `arguments`.
"""
node_count(t) = istree(t) ? reduce(+, node_count(x) for x in  arguments(t), init=0) + 1 : 1
export node_count


text:
# Filtering a PointCloud using a PassThrough filter
# http://www.pointclouds.org/documentation/tutorials/passthrough.php#passthrough

using PCLCommon
using PCLFilters
using Cxx

cloud = PointCloud{PointXYZ}(5, 1)
cloud_filtered = PointCloud{PointXYZ}()

for i in 0:length(cloud)-1
    cloud[i,:x] = randn()
    cloud[i,:y] = randn()
    cloud[i,:z] = randn()
end

println("Cloud before filtering")
for i in 0:length(cloud)-1
    println("    ", cloud[i,:x], " ", cloud[i,:y], " ", cloud[i,:z])
end

pass = PassThrough{PointXYZ}()
setInputCloud(pass, cloud)
setFilterFieldName(pass, "z")
setFilterLimits(pass, 0.0, 1.0)
# setFilterLimitsNegative(pass, true)
filter(pass, cloud_filtered)

println("Cloud after filtering")
for i in 0:length(cloud_filtered)-1
    println("    ", cloud_filtered[i,:x], " ", cloud_filtered[i,:y],
        " ", cloud_filtered[i,:z])
end


text:
using Findpeaks

using DelimitedFiles

using Plots
default(legend=false)

data = readdlm("../assets/example_spectrum.txt")
x = data[:, 1]
y = data[:, 2]

peaks = findpeaks(y, x, min_prom=1000.)

plot(x, y, title="Prominent peaks")
scatter!(x[peaks], y[peaks])
savefig("prom_peaks.png")

sep = 0.2

peaks = findpeaks(y, x, min_dist=sep)

plot(x, y, title="Peaks at least $sep units apart")
scatter!(x[peaks], y[peaks])
savefig("dist_peaks.png")


text:
# Test PUX by extending a function and computing gradient using Fourier method

push!(LOAD_PATH, string(pwd(),"/src"))
using ModifiedStokesSolver
using CurveDiscretization
import AnalyticDomains
using PUX
using Base.Test

@testset "PUX" begin
    srand(0)
    # Params
    alpha = 2.0
    numpanels = 10
    panelorder = 16
    Ngrid_int = 150
    pux_ep = 2.0
    pux_P = 40
    # Discretize
    curve = AnalyticDomains.starfish(n_arms=3, amplitude=0.1)
    dcurve = discretize(curve, numpanels, panelorder,
                        equal_arclength=true)
    ## PUX prep
    PUXStor, X, Y, Lgrid, Ngrid, interior = PUX.pux_precompute(curve, dcurve, Ngrid_int, pux_ep; P=pux_P)
    # Test function
    f(x, y) = x.^2 + y.^2 + x.*y
    fx(x, y) = 2*x + y
    fy(x, y) = 2*y + x
    F = f.(X,Y)
    Fxref = fx.(X,Y)
    Fyref = fy.(X,Y)
    # Extend
    FEXT = PUX.pux_eval(F, PUXStor)
    # Differentiate
    Fx, Fy = per_gradient(FEXT, Lgrid)
    Ex = abs.(Fx-Fxref)
    Ey = abs.(Fy-Fyref)
    E = max.(Ex,Ey)
    logE = log10.(E)
    Emax = norm(vec(E[interior]), Inf)
    @test Emax < 1e-6
end


text:
# Guessing game in Julia 1.2.
# Let's call the magic number λ, and a guess γ.

"Showing Julia's functional syntax & multi-dispatch."
equal(λ, x) = λ == x
equal(λ) = x -> equal(λ, x)
# @show equal(37)

"Check a guess against a magic number, and record this try."
function check(λ, γ, tries)
    push!(tries, γ)
    if equal(λ)(γ)
        println("Correct!")
        return true
    else
        if γ < λ
            println("Too small!")
        elseif γ > λ
            println("Too large!")
        end
        return false
    end
end

"Show game stats."
function stats(tries)
    println("<Game stats>")
    print("  You tried: ")
    for num in tries
        print("$num ")
    end
    println()
    println("  In total "*string(length(tries))*" tries used.")
end

"Game logic."
function play()
    println("Welcome to Guanzhou's guessing game!")

    # Generate random number.
    λ = rand(1:100)

    # Guessing loop.
    win = false
    tries = []
    while true
        print("Enter an int ∈ [1, 100], or 'q' to quit: ")
        ipt = strip(readline())
        if length(ipt) == 0     # Empty input.
            continue
        end
        if ipt[1] == 'q'        # Quit signal.
            break
        end
        if (ipt[1] < '0' || ipt[1] > '9') && ipt[1] != '-'
            println("WARN: not a valid integer. Try again...")
        else
            γ = parse(Int, ipt)
            if !(1 <= γ <= 100)
                println("WARN: valid input ∈ [1, 100]. Try again...")
            else
                print("  Your guess is... ")
                if check(λ, γ, tries)
                    win = true
                    break
                end
            end
        end
    end

    # Ending message.
    if win
        println("You win! Congrats ;)")
    else
        println("Sad to see you go :(")
    end
    stats(tries)
end

play();


text:
using Trace
using Base.Test
using Memento

const FMT_STR = "[{level}]:{name} - {msg}"
const LEVELS = Dict(
    "not_set" => 0,
    "trace" => 5,
    "debug" => 10,
    "info" => 20,
    "warn" => 30,
    "error" => 40,
)

io = IOBuffer()

handler = DefaultHandler(
    io, DefaultFormatter(FMT_STR)
)

logger = Logger(
    "Logger.example",
    Dict("Buffer" => handler),
    "info",
    LEVELS,
    DefaultRecord,
    true
)

skip_trace = median(map(1:1000) do i
    tic()
    @info(logger, "My skipped message")
    res = toq()
    @test isempty(takebuf_string(io))
    return res
end)

# enable tracing
Trace.enable()

log_trace = median(map(1:1000) do i
    tic()
    @info(logger, "My logged message")
    res = toq()
    @test contains(takebuf_string(io), "My logged message")
    return res
end)

# Compare against logging alone
skip_log = median(map(1:1000) do i
    tic()
    debug(logger, "My skipped message")
    return toq()
end)

log_log = median(map(1:1000) do i
    tic()
    info(logger, "My skipped message")
    return toq()
end)

println("Skipped trace time: $skip_trace")
println("Skipped log time: $skip_log")
println("Logged trace time: $log_trace")
println("Logged log time: $log_log")

@test skip_trace < skip_log
@test skip_trace < log_trace

set_level(logger, "trace")
@trace(logger, median(rand(1000)))


text:
     random_int(min=0, max=9999)="$(rand(min : max))"

     random_digit()= "$(rand(0:9))"

     random_digit_not_null()= "$(rand(1:9))"

     random_digit_or_empty()= rand(0:1)==1? "$(rand(0:9))" : " "

     random_digit_not_null_or_empty()= rand(0:1)==1? "$(rand(1:9))" : " "

     random_number( digits= 1 ) = "$(rand(0 : (10^digits) - 1))"

     random_letter()= rand(0:1)>=0.5? "$(rand('a':'z'))" : "$(rand('A':'Z'))"

     random_element( elements=("a", "b", "b"))= elements[rand(1: length(elements))]

     function randomize_nb_elements( number=10, le=false, ge=false)
                if le==true && ge==true
                    return "$(number)"
                end
                if ge==true
                    min = 100
                else
                    min=60
                end

                if le==true
                    max = 100
                else
                    max=140
                end

                 "$(((rand(min : max))number / 100) + 1 )"
            end

     function numerify( text = "####")
                    for i in text
                    text= replace(text,"#","$(rand(0:9))",1)
                    end
                return text
                end

     lexify( text="????")= (for i=text  text= replace(text,"?","$(random_letter())",1) end ; return text)

     bothify( text="## ??")=lexify(numerify(text))


text:
## Script to read a sample of gene trees and do some perturbation:
## - change the root position in a percentage of them
## - perform a random NNI move on a percentage of them
## This file is not meant to be used as a script, but use its parts
## into the simulations. That is, we do not necessarily want
## to do both types of perturbations at the same time.
## Claudia (June 2022)

## Loading packages
using PhyloNetworks
using PhyloPlots ## just to test the perturbations
using Random, Distributions

## Loading functions
include("gene-tree-perturbations-fns.jl")

## Global parameters
prc = 0.3 ## percentage of gene trees that will be perturbed
globalseed = 6041104 ## global seed

## Read sample of gene trees
treefile = "data/1_astral.in"
tree = readMultiTopology(treefile)
n = length(tree) ## number of gene trees

## Selecting subsample to be perturbed
Random.seed!(globalseed)
nk = Int(n*prc)
inds = sample(1:n,nk,replace=false)

## Generating specific seeds
seeds = sample(1:5555555,nk,replace=false)

# --------------------------------------
# Changing root to random place
# --------------------------------------

# Specific parameters
out = true ## see gene-tree-perturbations-fns.jl

for i in 1:nk
    changeRandomRoot!(tree[inds[i]], seeds[i], out=out)
end

# Write down the new sample of gene trees (some perturbed)
rootname = split(treefile, '.')[1]
outfile = string(rootname,"-root-changed.out")
writeMultiTopology(trees, outfile)

# --------------------------------------
# Random NNI move
# --------------------------------------

for i in 1:nk
    randomNNI!(tree[inds[i]], seeds[i])
end

# Write down the new sample of gene trees (some perturbed)
rootname = split(treefile, '.')[1]
outfile = string(rootname,"-nni.out")
writeMultiTopology(trees, outfile)

max_stars_repo_name: gtiley/bad-diamonds

max_issues_repo_name: gtiley/diamond-identifiability

max_forks_repo_name: gtiley/diamond-identifiability


text:
using EuroTokenizer
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@test 1 == 2


text:
using Oceananigans.Grids: halo_size

function run_horizontal_average_tests(arch, FT)
    topo = (Periodic, Periodic, Bounded)
    Nx = Ny = Nz = 4
    grid = RegularCartesianGrid(topology=topo, size=(Nx, Ny, Nz), extent=(100, 100, 100))
    Hx, Hy, Hz = halo_size(grid)

    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT)

    linear(x, y, z) = z
    set!(model, T=linear, w=linear)

    T̅ = Average(model.tracers.T, dims=(1, 2), with_halos=false)
    computed_T_profile = T̅(model)
    @test_skip size(computed_T_profile) == (1, 1, Nz)
    @test_skip computed_T_profile ≈ znodes(Cell, grid, reshape=true)

    T̅ = Average(model.tracers.T, dims=(1, 2), with_halos=true)
    computed_T_profile_with_halos = T̅(model)
    @test_skip size(computed_T_profile_with_halos) == (1, 1, Nz+2Hz)
    @test_skip computed_T_profile_with_halos[1+Hz:Nz+Hz] ≈ znodes(Cell, grid)

    w̅ = Average(model.velocities.w, dims=(1, 2), with_halos=false)
    computed_w_profile = w̅(model)
    @test_skip size(computed_w_profile) == (1, 1, Nz+1)
    @test_skip computed_w_profile ≈ znodes(Face, grid, reshape=true)

    w̅ = Average(model.velocities.w, dims=(1, 2), with_halos=true)
    computed_w_profile_with_halos = w̅(model)
    @test_skip size(computed_w_profile_with_halos) == (1, 1, Nz+1+2Hz)
    @test_skip computed_w_profile_with_halos[1+Hz:Nz+1+Hz] ≈ znodes(Face, grid)
end

function run_zonal_average_tests(arch, FT)
    topo = (Periodic, Bounded, Bounded)
    Nx = Ny = Nz = 4
    grid = RegularCartesianGrid(topology=topo, size=(Nx, Ny, Nz), extent=(100, 100, 100))
    Hx, Hy, Hz = halo_size(grid)

    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT)

    linear(x, y, z) = z
    set!(model, T=linear, v=linear)

    T̅ = Average(model.tracers.T, dims=1, with_halos=false)
    computed_T_slice = T̅(model)
    @test_skip size(computed_T_slice) == (1, Ny, Nz)

    computed_T_slice = dropdims(computed_T_slice, dims=1)
    zC = znodes(Cell, grid)
    @test_skip all(computed_T_slice[j, :] ≈ zC for j in 1:Ny)

    T̅ = Average(model.tracers.T, dims=1, with_halos=true)
    computed_T_slice_with_halos = T̅(model)
    @test_skip size(computed_T_slice_with_halos) == (1, Ny+2Hy, Nz+2Hz)

    computed_T_slice_with_halos = dropdims(computed_T_slice_with_halos, dims=1)
    @test_skip computed_T_slice_with_halos[1+Hy:Ny+Hy, 1+Hz:Nz+Hz] ≈ computed_T_slice

    v̅ = Average(model.velocities.v, dims=1, with_halos=false)
    computed_v_slice = v̅(model)
    @test_skip size(computed_v_slice) == (1, Ny+1, Nz)

    computed_v_slice = dropdims(computed_v_slice, dims=1)
    zC = znodes(Cell, grid)
    @test_skip all(computed_v_slice[j, :] ≈ zC for j in 1:Ny)

    v̅ = Average(model.velocities.v, dims=1, with_halos=true)
    computed_v_slice_with_halos = v̅(model)
    @test_skip size(computed_v_slice_with_halos) == (1, Ny+1+2Hy, Nz+2Hz)

    computed_v_slice_with_halos = dropdims(computed_v_slice_with_halos, dims=1)
    @test_skip computed_v_slice_with_halos[1+Hy:Ny+1+Hy, 1+Hz:Nz+Hz] ≈ computed_v_slice
end

function run_volume_average_tests(arch, FT)
    Nx = Ny = Nz = 4
    grid = RegularCartesianGrid(size=(Nx, Ny, Nz), extent=(100, 100, 100))
    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT)

    T₀(x, y, z) = z
    set!(model, T=T₀)

    T̅ = Average(model.tracers.T, dims=(1, 2, 3), time_interval=0.5second, with_halos=false)
    computed_scalar = T̅(model)
    @test_skip size(computed_scalar) == (1, 1, 1)
    @test_skip all(computed_scalar .≈ -50.0)

    T̅ = Average(model.tracers.T, dims=(1, 2, 3), time_interval=0.5second, with_halos=true)
    computed_scalar_with_halos = T̅(model)
    @test_skip size(computed_scalar_with_halos) == (1, 1, 1)
    @test_skip all(computed_scalar_with_halos .≈ -50.0)
end

function nan_checker_aborts_simulation(arch, FT)
    grid = RegularCartesianGrid(size=(4, 2, 1), extent=(1, 1, 1))
    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT)

    # It checks for NaNs in w by default.
    nc = NaNChecker(model; iteration_interval=1, fields=Dict(:w => model.velocities.w.data.parent))
    push!(model.diagnostics, nc)

    model.velocities.w[3, 2, 1] = NaN

    time_step!(model, 1, 1)
end

TestModel(::GPU, FT, ν=1.0, Δx=0.5) =
    IncompressibleModel(
          grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3Δx, 3Δx, 3Δx)),
       closure = IsotropicDiffusivity(FT, ν=ν, κ=ν),
  architecture = GPU(),
    float_type = FT
)

TestModel(::CPU, FT, ν=1.0, Δx=0.5) =
    IncompressibleModel(
          grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3Δx, 3Δx, 3Δx)),
       closure = IsotropicDiffusivity(FT, ν=ν, κ=ν),
  architecture = CPU(),
    float_type = FT
)

function max_abs_field_diagnostic_is_correct(arch, FT)
    model = TestModel(arch, FT)
    set!(model.velocities.u, rand(size(model.grid)))
    u_max = FieldMaximum(abs, model.velocities.u)
    return u_max(model) == maximum(abs, model.velocities.u.data.parent)
end

function advective_cfl_diagnostic_is_correct(arch, FT)
    model = TestModel(arch, FT)

    Δt = FT(1.3e-6)
    Δx = FT(model.grid.Δx)
    u₀ = FT(1.2)
    CFL_by_hand = Δt * u₀ / Δx

    model.velocities.u.data.parent .= u₀
    cfl = AdvectiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

function diffusive_cfl_diagnostic_is_correct(arch, FT)
    Δt = FT(1.3e-6)
    Δx = FT(0.5)
    ν = FT(1.2)
    CFL_by_hand = Δt * ν / Δx^2

    model = TestModel(arch, FT, ν, Δx)
    cfl = DiffusiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

get_iteration(model) = model.clock.iteration
get_time(model) = model.clock.time

function timeseries_diagnostic_works(arch, FT)
    model = TestModel(arch, FT)
    Δt = FT(1e-16)
    simulation = Simulation(model, Δt=Δt, stop_iteration=1)
    iter_diag = TimeSeries(get_iteration, model, iteration_interval=1)
    push!(simulation.diagnostics, iter_diag)
    run!(simulation)
    return iter_diag.time[end] == Δt && iter_diag.data[end] == 1
end

function timeseries_diagnostic_tuples(arch, FT)
    model = TestModel(arch, FT)
    Δt = FT(1e-16)
    simulation = Simulation(model, Δt=Δt, stop_iteration=2)
    timeseries = TimeSeries((iters=get_iteration, itertimes=get_time), model, iteration_interval=2)
    simulation.diagnostics[:timeseries] = timeseries
    run!(simulation)
    return timeseries.iters[end] == 2 && timeseries.itertimes[end] == 2Δt
end

function diagnostics_getindex(arch, FT)
    model = TestModel(arch, FT)
    simulation = Simulation(model, Δt=0, stop_iteration=0)
    iter_timeseries = TimeSeries(get_iteration, model)
    time_timeseries = TimeSeries(get_time, model)
    simulation.diagnostics[:iters] = iter_timeseries
    simulation.diagnostics[:times] = time_timeseries
    return simulation.diagnostics[2] == time_timeseries
end

function diagnostics_setindex(arch, FT)
    model = TestModel(arch, FT)
    simulation = Simulation(model, Δt=0, stop_iteration=0)

    iter_timeseries = TimeSeries(get_iteration, model)
    time_timeseries = TimeSeries(get_time, model)
    max_abs_u_timeseries = TimeSeries(FieldMaximum(abs, model.velocities.u), model, iteration_interval=1)

    push!(simulation.diagnostics, iter_timeseries, time_timeseries)
    simulation.diagnostics[2] = max_abs_u_timeseries

    return simulation.diagnostics[:diag2] == max_abs_u_timeseries
end

@testset "Diagnostics" begin
    @info "Testing diagnostics..."

    for arch in archs
        @testset "Average [$(typeof(arch))]" begin
            @info "  Testing averages [$(typeof(arch))]"
            for FT in float_types
                run_horizontal_average_tests(arch, FT)
                run_zonal_average_tests(arch, FT)
                run_volume_average_tests(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "NaN Checker [$(typeof(arch))]" begin
            @info "  Testing NaN Checker [$(typeof(arch))]"
            for FT in float_types
                @test_throws ErrorException nan_checker_aborts_simulation(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "Miscellaneous timeseries diagnostics [$(typeof(arch))]" begin
            @info "  Testing miscellaneous timeseries diagnostics [$(typeof(arch))]"
            for FT in float_types
                @test diffusive_cfl_diagnostic_is_correct(arch, FT)
                @test advective_cfl_diagnostic_is_correct(arch, FT)
                @test max_abs_field_diagnostic_is_correct(arch, FT)
                @test timeseries_diagnostic_works(arch, FT)
                @test timeseries_diagnostic_tuples(arch, FT)
                @test diagnostics_getindex(arch, FT)
                @test diagnostics_setindex(arch, FT)
            end
        end
    end
end


text:
function jsonify_details(details::Union{Dict{String,Any},Nothing})::String
    details_json::String = ""
    if !isnothing(details)
        details_json = JSON.json(details)
    end
    return details_json
end

function unsafe_string_and_free(cstring::Cstring, length::Union{Ref{Csize_t}, Nothing}=nothing)::String
    if isnothing(length)
        str = unsafe_string(cstring)
    else
        str = unsafe_string(pointer(cstring), length[])
    end
    ccall((:free_ptr, librdkitcffi), Cvoid, (Cstring,), cstring)
    return str
end

"""
    version()

Get RDKit version.

```julia
version = version()
```
"""
function version()::String
    val = ccall((:version, librdkitcffi), Cstring, ())
    version = unsafe_string_and_free(val)
    return version
end

"""
    enable_logging()

Enable RDKit logging.

```julia
enable_logging()
```
"""
function enable_logging()
    ccall((:enable_logging, librdkitcffi), Cvoid, ())
end

"""
    disable_logging()

Disable RDKit logging.

```julia
disable_logging()
```
"""
function disable_logging()
    ccall((:disable_logging, librdkitcffi), Cvoid, ())
end


text:
"""
    function three_eyes()

A cow with three eyes, brought to you by dpetrou@csua.berkeley.edu

# Example
```jldoctest
julia> cowsay("The better to see you with...", cow=Cowsay.three_eyes)
 _______________________________
< The better to see you with... >
 -------------------------------
        \\  ^___^
         \\ (ooo)\\_______
           (___)\\       )\\/\\
                ||----w |
                ||     ||
```

"""
function three_eyes(;eyes="oo", tongue="  ", thoughts="\\")
eye = first(eyes)
eyes = repeat(eye, 3)
the_cow =
"""
        $thoughts  ^___^
         $thoughts ($eyes)\\_______
           (___)\\       )\\/\\
            $tongue  ||----w |
                ||     ||
"""
return the_cow
end


text:
module geojl
  #include("geocommon.jl")
  include("mapquest.jl")
  include("geocoderus.jl")
  export MapQuest, GeoCoderUs
end

max_stars_repo_name: astrieanna/geojl

max_issues_repo_name: astrieanna/geojl

max_forks_repo_name: astrieanna/geojl


text:
# # dempe.mod	QQR2-MN-4-3
# # Original AMPL coding by Sven Leyffer
#
# # An MPEC from S. Dempe, "A necessary and sufficient optimality
# # condition for bilevel programming problems", Optimization 25,
# # pp. 341-354, 1992. From book Math. Progr. with Equil. Constr,
# # by Luo, Pang & Ralph, CUP, 1997, p. 354.
#
# # Number of variables:   2 + 1 multipliers
# # Number of constraints: 2
# # Nonlinear complementarity constraints
#
# var x;
# var z;
# var w >= 0;
#
# minimize f: (x - 3.5)^2 + (z + 4)^2;
#
# subject to
#     con1:  z - 3 + 2*z*w = 0;
#     con2:  0 >= z^2 - x     complements  w >= 0;
#
# data;
#
# # starting point
# let x := 1;
# let z := 1;
# let w := 1;
#
#
# let x :=   0.183193 ;
# let z := 0.428106;
# let w := 3.00379;

using JuMP, Complementarity
using Ipopt

m = Model(solver=IpoptSolver())

@variable(m, x)
@variable(m, z)
@variable(m, w>=0)

@NLobjective(m, Min, (x - 3.5)^2 + (z + 4)^2)

@NLconstraint(m, z - 3 + 2*z*w == 0)

@complements(m, 0 >= z^2 - x,  w >= 0, simple)

# setvalue(x, 50)
# setvalue(z, 50)
# setvalue(w, 1e5)

solve(m)

@show getobjectivevalue(m)
@assert isapprox(getobjectivevalue(m), 28.25, atol=1e-4)


text:
function add_λ_units(λ)
    expr = quote

        using Unitful
        using PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e

        local λ = $λ
        local ω = 2π*c_0/λ

        @unit unit_l "c/ω" UnitLength c_0/ω false
        @unit unit_t "1/ω" UnitTime 1/ω false
        @unit unit_B "mₑ ω/e" UnitMagneticField m_e * ω / e false
        @unit unit_E "mₑ c ω/e" UnitElectricField m_e * ω / e * c_0 false
        @unit unit_L "mₑ c²/ω" UnitAngularMomenta m_e * c_0^2 / ω false
        @unit pₑ "mₑ c" ElectronMomenta m_e * c_0 false

        Unitful.register(@__MODULE__)
    end
    eval(expr)
end

function add_λ_units(sim::EPOCHSimulation)
    λ = get_parameter(sim, :laser, :lambda) |> u"nm"
    add_λ_units(λ)
end


text:
using JuliaCon
using Test
using Distributed
import Dates

fakenow = Dates.DateTime("2020-07-29T16:30:00.000") # JuliaCon2020

@testset "JuliaCon.jl" begin
    @testset "juliacon2021()" begin
        @info "Local:"
        @test isnothing(juliacon2021())

        @info "Distributed:"
        withenv("JULIA_PROJECT"=>pwd()) do
            addprocs(4)
            @eval Main @everywhere using JuliaCon
            @test isnothing(@everywhere juliacon2021())
            rmprocs(workers())
        end
    end

    @testset "Schedule" begin
        @test isnothing(JuliaCon.now())
        @test isnothing(JuliaCon.now(; now=fakenow))
        @test isnothing(JuliaCon.today())
        @test isnothing(JuliaCon.today(; now=fakenow))
        @test isnothing(JuliaCon.today(; now=fakenow, track="BoF"))
    end
end


text:
@doc raw"""
    pc = PController(Kc)

Construct a Proportional (P) controller by specifying the controller gain defined under the following transfer function representation:

$$g_c(s)=K_c$$

# Arguments
* `Kc::Float64`: controller gain

# Example
```julia
pc = PController(1.0) # specify P controller gain
gc = TransferFunction(pc) # construct transfer function with this P-controller gain
```
"""
mutable struct PController
    Kc::Float64
end

TransferFunction(pc::PController) = TransferFunction([pc.Kc], [1.0])

@doc raw"""
    pic = PIController(Kc, τI)

Construct a Proportional-Integral (PI) controller by specifying the controller gain and integral time constant defined under the following transfer function representation:

$$g_c(s)=K_c \left[1+\frac{1}{\tau_I s}\right]$$

# Arguments
* `Kc::Float64`: controller gain
* `τI::Float64`: integral time constant

# Example
```julia
pic = PIController(1.0, 3.0) # specify PI controller params
gc = TransferFunction(pic) # construct transfer function with these PI-controller params
```
"""
mutable struct PIController
    Kc::Float64
    τI::Float64
end

TransferFunction(pic::PIController) = pic.Kc * (1 + 1 / (pic.τI * s))

mutable struct PIDController
    Kc::Float64
    τI::Float64
    τD::Float64
    α::Float64
end

@doc raw"""
    pidc = PIDController(Kc, τI, τD, α=0.0)

Construct a Proportional-Integral-Derivative (PID) controller by specifying the controller gain, integral time constant, derivative time constant, and derivative filter defined under the following transfer function representation:

$$g_c(s)=K_c \left[1+\frac{1}{\tau_I s}+\tau_D s \frac{1}{\alpha \tau_D s + 1}\right]$$

# Arguments
* `Kc::Float64`: controller gain
* `τI::Float64`: integral time constant
* `τD::Float64`: derivative time constant
* `α::Float64`: derivative filter

# Example
```julia
pidc = PIDController(1.0, 3.0, 0.1) # specify PID controller params
gc = TransferFunction(pidc) # construct transfer function with these PID-controller params
```
"""
PIDController(Kc::Float64, τI::Float64, τD::Float64; α::Float64=0.0) = PIDController(Kc, τI, τD, α)

TransferFunction(pidc::PIDController) = pidc.Kc * (1 + 1.0 / (pidc.τI * s) + pidc.τD * s / (pidc.α * pidc.τD * s + 1))

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
dir = @__DIR__
g1 = readg(dir*"/g150219.03200") # (+Bt, +Ip)
g2 = readg(dir*"/g133221.01150") # (-Bt, +Ip)
g3 = readg(dir*"/g159177.02700") # (+Bt, -Ip)
g4 = readg(dir*"/g153298.04400") # (-Bt, -Ip)

btip1 = (1,1)
btip2 = (-1,1)
btip3 = (1,-1)
btip4 = (-1,-1)

# COCOS for the GEQDSK files
cc1 = 1
cc2 = 5
cc3 = 3
cc4 = 7

M1 = efit(g1, clockwise_phi=false)
M2 = efit(g2, clockwise_phi=false)
M3 = efit(g3, clockwise_phi=false)
M4 = efit(g4, clockwise_phi=false)

r1 = (M1.axis[1]+0.1,M1.axis[2])
r2 = (M2.axis[1]+0.1,M2.axis[2])
r3 = (M3.axis[1]+0.1,M3.axis[2])
r4 = (M4.axis[1]+0.1,M4.axis[2])

test_data = ((cc1,g1,M1,r1,btip1), (cc2,g2,M2,r2,btip2), (cc3,g3,M3,r3,btip3), (cc4,g4,M4,r4,btip4))

@testset "EFIT Tests" begin

    @testset "Bt-Ip Signs" begin
        for (cc, g, M, r, btip) in test_data
            @test (sign(g.bcentr), sign(g.current)) == btip
            Bt, Ip = (sign(Bfield(M,r...)[2]), sign(Jfield(M,r...)[2]))
            @test (Bt, Ip) == btip
        end
    end

    @testset verbose = true "CurlB == μ₀J" begin
        for (cc, g, M, r, btip) in test_data
            @test curlB(M,r...) ≈ mu0*Jfield(M,r...) rtol=0.01
        end
    end

    @testset verbose = true "3D tests" begin
        for (cc, g, M, r, btip) in test_data
            @test Bfield(M,r...) ≈ Bfield(M,r[1],zero(r[1]),r[2]) rtol=0.01
            @test Jfield(M,r...) ≈ Jfield(M,r[1],zero(r[1]),r[2]) rtol=0.01

            @test gradB(M,r...) ≈ gradB(M,r[1],zero(r[1]),r[2]) rtol=0.01
            @test gradB(M,r[1],zero(r[1]),r[2]) ≈ Equilibrium.gradB_autodiff(M,r[1],zero(r[1]),r[2]) rtol=0.01

            @test curlB(M,r...) ≈ curlB(M,r[1],zero(r[1]),r[2]) rtol=0.01
            @test curlB(M,r[1],zero(r[1]),r[2]) ≈ Equilibrium.curlB_autodiff(M,r[1],zero(r[1]),r[2]) rtol=0.01
        end
    end

    @testset verbose = true "COCOS Identify" begin
        for (cc, g, M, r, btip) in test_data
            @test identify_cocos(g; clockwise_phi = false) == (cc,)
        end
    end

    @testset verbose = true "COCOS Consistancy" begin
        for (cc, g, M, r, btip) in test_data
            @test check_cocos(g, cc)
        end
    end

    @testset verbose = true "COCOS Transformation" begin
        for (cc_in, g, M, r, btip) in test_data
            for cc = 1:8
                _cc = cocos(cc)
                g_new = transform_cocos(g, cc_in, _cc)
                @test identify_cocos(g_new, clockwise_phi = _cc.sigma_RpZ < 0) == (cc,)
                @test check_cocos(g_new, cc)
                @test curlB(M,r...) ≈ mu0*Jfield(M,r...) rtol=0.01
            end
        end
    end

end


text:
"""
DuckDB validity mask
"""
mutable struct ValidityMask
    data::Vector{UInt64}

    function ValidityMask(data::Vector{UInt64})
        result = new(data)
        return result
    end
end

const BITS_PER_VALUE = 64;

function get_entry_index(row_idx)
    return ((row_idx - 1) ÷ BITS_PER_VALUE) + 1
end

function get_index_in_entry(row_idx)
    return (row_idx - 1) % BITS_PER_VALUE
end

function setinvalid(mask::ValidityMask, index)
    entry_idx = get_entry_index(index)
    index_in_entry = get_index_in_entry(index)
    mask.data[entry_idx] &= ~(1 << index_in_entry)
    return
end

function isvalid(mask::ValidityMask, index)::Bool
    entry_idx = get_entry_index(index)
    index_in_entry = get_index_in_entry(index)
    return (mask.data[entry_idx] & (1 << index_in_entry)) != 0
end

max_stars_repo_name: kimmolinna/duckdb

max_issues_repo_name: kimmolinna/duckdb

max_forks_repo_name: kimmolinna/duckdb


text:
## Majority system for individual dances

"""
    skating_single_dance(results)
    skating_single_dance(results, majority_from; [initial_place, initial_column, depth])

Calculates a skating tableau and placement for judgment results of a single dance contained
in the DataFrame `results`. It must contain the starter numbers of the dancers as a first
column named `:Numbers`, and the places assigned by the different judges in the following
columns, for example as in
```julia
DataFrame(Number = [11, 12, 13, 14, 15],
        JudgeA = [1, 2, 3, 4, 5],
        JudgeB = [1, 3, 2, 4, 5],
        JudgeC = [1, 3, 2, 5, 4],
        JudgeD = [5, 1, 2, 3, 4],
        JudgeE = [2, 1, 5, 3, 4])
```

Returns a DataFrame with the skating tableau represented as strings, as well as a DataFrame
with the starter numbers and places.

For calculations on cropped tables, the majority of votes, the initial placement to fill,
the initial column in the tableau where skating should start as well as the desired depth of
the calculation can be defined by keyword arguments.
"""
function skating_single_dance(results::DataFrame, majority_from::Int; initial_place::Int = 1,
    initial_column::Int = 1, depth = size(results, 1))

    calculation, sum_of_eval = prepare_sums(results, depth)

    calculation[!, :Place] .= 0.0
    calculation_text = string.(copy(calculation))

    current_place = initial_place
    max_cols = depth
    current_col = initial_column
    tmp_col = initial_column
    append_sum = false
    while current_place <= (initial_column - 1 + size(results, 1))
        if count(>=(majority_from), calculation[!, current_col+1], dims = 1)[1] == 1
            # @info "Clear Majority"
            idx = findfirst(>=(majority_from), calculation[!, current_col+1])
            write_result!(calculation, calculation_text, sum_of_eval, idx, current_place, current_col, current_col, max_cols)
            current_place += 1
            current_col += 1
        else
            # @info "Multiple Majorities"
            idx = findall(>=(majority_from), calculation[!, current_col+1])
            tmp_col = copy(current_col)
            while !isempty(idx)
                if count(==(maximum(calculation[idx, tmp_col+1])), calculation[idx, tmp_col+1]) == 1
                    # @info "Single maximal majority"
                    id = idx[findfirst(==(maximum(calculation[idx, tmp_col+1])), calculation[idx, tmp_col+1])]
                    write_result!(calculation, calculation_text, sum_of_eval, id, current_place, current_col, tmp_col, max_cols; append_sum = append_sum)
                    current_place += 1
                    append_sum = false
                    remove!(idx, id)
                else
                    # equal majority
                    # look at sum of evaluations up to current column
                    if length(findall(==(maximum(sum_of_eval[idx, tmp_col+1])), sum_of_eval[idx, tmp_col+1])) == 1
                        # @info "Single sum minority"
                        id = findfirst(==(minimum(sum_of_eval[idx, tmp_col+1])), sum_of_eval[idx, tmp_col+1])
                        write_result!(calculation, calculation_text, sum_of_eval, idx[id], current_place, current_col, tmp_col, max_cols; append_sum = true)
                        current_place += 1
                        append_sum = true
                        remove!(idx, idx[id])
                    else
                        # @info "Multiple equal sums"
                        calculation_text[idx, (tmp_col+1)] = string.(calculation[idx, (tmp_col+1)]) .* "*" .* "(" .* string.(sum_of_eval[idx, (tmp_col+1)]) .* ")"
                        tmp_col += 1
                        if tmp_col > max_cols
                            # @info "tmp_cols > max_cols"
                            places = range(current_place, length = length(idx))
                            place = mean(places)
                            calculation[idx, :Place] .= place
                            calculation_text[idx, :Place] .= string(place)
                            current_place = maximum(places) + 1
                            calculation[idx, (current_col+1):(max_cols+1)] .= 0
                            idx = []
                            break
                        end
                    end

                end
            end
            current_col += 1
        end
    end

    return calculation_text, calculation[!, [:Number, :Place]]
end

function skating_single_dance(results::DataFrame; initial_place = 1, initial_column = 1, depth = size(results, 1))
    return skating_single_dance(results::DataFrame, calc_majority(results);
        initial_place = initial_place, initial_column = initial_column, depth = depth)
end

"""
    calc_majority(results)

Returns the judge majority for a result DataFrame `results`.
"""
function calc_majority(results)
    return ceil(Int, (size(results, 2) - 2) / 2 + 1)
end

"""
    write_result!(calculation, calculation_text, sum_of_eval, id, current_place, current_col, tmp_col, max_cols; append_sum = false)

Write the result of the skating procedure to the DataFrames containing placement and summary
information. Remove the parts that are not needed anymore.
"""
function write_result!(calculation, calculation_text, sum_of_eval, index, current_place, current_col, tmp_col, max_cols;
    append_sum = false)
    calculation[index, :Place] = current_place
    calculation_text[index, :Place] = string(current_place)
    str = string(calculation[index, (tmp_col+1)]) * "*"
    if append_sum
        str = str * "(" * string(sum_of_eval[index, (tmp_col+1)]) * ")"
    end
    calculation_text[index, (tmp_col+1)] = str
    calculation_text[index, (tmp_col+2):(max_cols+1)] .= "-"
    calculation[index, (current_col+1):(max_cols+1)] .= 0
end


text:
using LLG
using Test

@testset "LLG.jl" begin
    # Write your tests here.
end

max_stars_repo_name: pgunnink/LLGTools

max_issues_repo_name: pgunnink/LLGTools

max_forks_repo_name: pgunnink/LLGTools


text:
"""
    global_move(mc::MC, m::Model, conf) -> accepted::Bool

Propose a global move for configuration `conf`.
Returns wether the global move has been accepted or not.
"""
global_move(mc::MC, m::Model, conf) = false

"""
    prepare_observables(mc::MC, m::Model) -> Dict{String, Observable}

Initializes observables and returns a `Dict{String, Observable}`.
In the latter, keys are abbreviations for the observables names and values
are the observables themselves.

See also [`measure_observables!`](@ref) and [`finish_observables!`](@ref).
"""
prepare_observables(mc::MC, m::Model) = Dict{String, Observable}()

"""
    measure_observables!(mc::MC, m::Model, obs::Dict{String,Observable}, conf)

Measures observables and updates corresponding `MonteCarloObservable.Observable`
objects in `obs`.

See also [`prepare_observables`](@ref) and [`finish_observables!`](@ref).
"""
measure_observables!(mc::MC, m::Model, obs::Dict{String,Observable}, conf) = nothing

"""
    finish_observables!(mc::MC, m::Model, obs::Dict{String,Observable})

Finish measurements of observables.

See also [`prepare_observables`](@ref) and [`measure_observables!`](@ref).
"""
finish_observables!(mc::MC, m::Model, obs::Dict{String,Observable}) = nothing

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
module ModelWrappers

############################################################################################
#Import External packages
import Base: Base, length, fill, fill!
import StatsBase: StatsBase, sample, sample!
import BaytesCore: BaytesCore, subset, update
using BaytesCore:
    AbstractModelWrapper,
    AbstractObjective,
    AbstractResult,
    BaytesCore,
    Tuple_to_Namedtuple,
    UpdateBool,
    UpdateTrue,
    UpdateFalse

using DocStringExtensions:
    DocStringExtensions, TYPEDEF, TYPEDFIELDS, FIELDS, SIGNATURES, FUNCTIONNAME
using ArgCheck: ArgCheck, @argcheck
using UnPack: UnPack, @unpack, @pack!
using Random: Random, AbstractRNG, GLOBAL_RNG

#!NOTE: These libraries are relevant for transform part
using LinearAlgebra: LinearAlgebra, Diagonal, LowerTriangular, tril!, diag, issymmetric
using Distributions: Distributions, Distribution, logpdf
using Bijectors:
    Bijectors, Bijector, logpdf_with_trans,
    TruncatedBijector, SimplexBijector, CorrBijector, PDBijector
#!NOTE: These libraries are only relevant for 'Differentiation' files - open to make separate library later on but for now easier to test changes with AD libraries included in a single library.
using ChainRulesCore, DistributionsAD, DiffResults
using ForwardDiff, ReverseDiff, Zygote

############################################################################################
# A bunch of constants used throughout the package
"Maximum value after transformation before tagged as non-finite."
const max_val = 1e+100
"Smallest decrease allowed in the log objective results before tagged as divergent."
const min_Δ = -1e+3

############################################################################################
#Import
include("Core/Core.jl")
include("Models/Models.jl")
include("Differentiation/Differentiation.jl")

############################################################################################
#export
export
    UpdateBool,
    UpdateTrue,
    UpdateFalse

end


text:
using DualNumbers
using Test

@test isempty(detect_ambiguities(DualNumbers, Base, Core))

@test checkindex(Bool, 1:3, dual(2))

# wrap in individual modules to avoid name conflicts.
module TestAutomaticDifferentiation
include("automatic_differentiation_test.jl")
end


text:
include((@__DIR__)*"/../preamble.jl")

include("absorbing_model/model_def.jl")
include("default_params.jl")

for ic_string in ["point_mass"]
    model_string = "qbdrap4"

    order = 1
    approx_data = CSV.read("hitting_times_model/hitting_times/data/"*
        ic_string*"/order_"*string(order)*"_model_"*model_string*".csv",DataFrame)
    x_vals = approx_data.t

    sim_data = CSV.read("hitting_times_model/hitting_times/data/"*ic_string*"/sims.csv",DataFrame)

    function cdf(x::AbstractVector,data::DataFrame)
        sort!(data,:t)

        x1 = data.t[data.φ.==1]
        x2 = data.t[data.φ.==2]

        c1 = 1
        c2 = 1

        F1 = zeros(length(x))
        F2 = zeros(length(x))

        for (c,xi) in enumerate(x)
            tmp1 = findfirst(y->y>xi,x1[c1:end])
            if !(tmp1===nothing)
                c1 += tmp1-1
            else
                c1 = length(x1)+1
            end
            F1[c] = c1-1

            tmp2 = findfirst(y->y>xi,x2[c2:end])
            if !(tmp2===nothing)
                c2 += tmp2-1
            else
                c2 = length(x2)+1
            end
            F2[c] = c2-1
        end
        return F1./length(data.t), F2./length(data.t)
    end

    F1,F2 = cdf(x_vals,sim_data)

    CSV.write((@__DIR__)*"/hitting_times/data/"*ic_string*"/sims_cdf_evaluated.csv",DataFrame(t=x_vals,phase_1=F1,phase_2=F2))

    function bootstrap(sim_data,x_vals,n_boot)
        rng = StableRNGs.StableRNG(1412)
        n_sims = length(sim_data.t)
        for n in 1:n_boot
            @show n
            ind = rand(rng,1:n_sims,n_sims)
            boot_data = sim_data[ind,:]
            F1,F2 = cdf(x_vals,boot_data)
            CSV.write(
                (@__DIR__)*"/hitting_times/data/"*ic_string*"/bootstrap/bootsample_"*string(n)*".csv",
                DataFrame(t=x_vals,phase_1=F1,phase_2=F2),
            )
        end
        return nothing
    end

    bootstrap(sim_data,x_vals,n_boot)

    # rng = StableRNGs.StableRNG(23122021)

    # pth = mkpath("hitting_times_model/hitting_times/data/"*ic_string*"/bootstrap")
    # for n in 1:n_boot
    #     idx = rand(rng,1:n_sims,n_sims)
    #     boot_cdf_fun = ht_cdf(sim_data[idx,:])
    #     boot_cdf_evaluated = DataFrame(t=approx_data.t,
    #                                     phase_1=boot_cdf_fun.(approx_data.t,1),
    #                                     phase_2=boot_cdf_fun.(approx_data.t,2))
    #     CSV.write(pth*"/bootsample_"*string(n)*".csv",boot_cdf_evaluated)
    #     print(string(n)*"... ")
    # end
end


text:
# Interface
Base.length(poly::Union{arb_poly_struct,ArbPoly}) =
    ccall(@libarb(arb_poly_length), Int, (Ref{arb_poly_struct},), poly)
Base.length(poly::Union{acb_poly_struct,AcbPoly}) =
    ccall(@libarb(acb_poly_length), Int, (Ref{acb_poly_struct},), poly)
Base.length(series::Union{ArbSeries,AcbSeries}) = series.degree + 1

# We have to define these separately so that it overrides the default
# correctly
degree(series::ArbSeries) = series.degree
degree(series::AcbSeries) = series.degree

Base.eltype(::Type{<:Union{ArbPoly,ArbSeries}}) = Arb
Base.eltype(::Type{<:Union{AcbPoly,AcbSeries}}) = Acb

Base.checkbounds(::Type{Bool}, poly::Union{ArbPoly,AcbPoly}, i::Integer) = i >= 0
Base.checkbounds(::Type{Bool}, series::Union{ArbSeries,AcbSeries}, i::Integer) =
    0 <= i <= degree(series)
Base.checkbounds(poly::Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries}, i::Integer) =
    checkbounds(Bool, poly, i) || throw(BoundsError(poly, i))

Base.@propagate_inbounds function Base.getindex(
    poly::T,
    i::Integer,
) where {T<:Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries}}
    @boundscheck checkbounds(poly, i)
    res = eltype(T)(prec = precision(poly))
    get_coeff!(res, poly, i)
    return res
end

Base.@propagate_inbounds function Base.setindex!(
    poly::Union{ArbPoly,ArbSeries},
    x::ArbLike,
    i::Integer,
)
    @boundscheck checkbounds(poly, i)
    set_coeff!(poly, i, x)
    return x
end

Base.@propagate_inbounds function Base.setindex!(
    poly::Union{AcbPoly,AcbSeries},
    x::AcbLike,
    i::Integer,
)
    @boundscheck checkbounds(poly, i)
    set_coeff!(poly, i, x)
    return x
end

Base.@propagate_inbounds function Base.setindex!(
    poly::Union{ArbPoly,ArbSeries,AcbPoly,AcbSeries},
    x,
    i::Integer,
)
    @boundscheck checkbounds(poly, i)
    set_coeff!(poly, i, convert(eltype(poly), x))
    return x
end

# Constructors
for TPoly in [:ArbPoly, :AcbPoly]
    @eval function $TPoly(poly::cstructtype($TPoly); prec::Integer = DEFAULT_PRECISION[])
        res = $TPoly(prec = prec)
        set!(res, poly)
        return res
    end

    @eval function $TPoly(coeff; prec::Integer = _precision(coeff))
        poly = $TPoly(prec = prec)
        poly[0] = coeff
        return poly
    end

    @eval function $TPoly(coeffs::AbstractVector; prec::Integer = _precision(first(coeffs)))
        poly = $TPoly(prec = prec)
        @inbounds for i = 1:length(coeffs)
            poly[i-1] = coeffs[i]
        end
        return poly
    end
end

for TSeries in [:ArbSeries, :AcbSeries]
    @eval function $TSeries(
        poly::cstructtype($TSeries),
        degree::Integer = degree(poly);
        prec::Integer = DEFAULT_PRECISION[],
    )
        res = $TSeries(degree, prec = prec)
        set!(res, poly)
        return res
    end

    @eval $TSeries(; prec::Integer = DEFAULT_PRECISION[]) = $TSeries(0, prec = prec)

    @eval function $TSeries(coeff, degree::Integer; prec::Integer = _precision(coeff))
        series = $TSeries(degree, prec = prec)
        series[0] = coeff
        return series
    end

    @eval function $TSeries(
        coeffs::AbstractVector,
        degree::Integer = length(coeffs) - 1;
        prec::Integer = _precision(first(coeffs)),
    )
        series = $TSeries(degree, prec = prec)
        @inbounds for i = 1:length(coeffs)
            series[i-1] = coeffs[i]
        end
        return series
    end
end

Base.zero(poly::T) where {T<:Union{ArbPoly,AcbPoly}} = T(prec = precision(poly))
Base.one(poly::T) where {T<:Union{ArbPoly,AcbPoly}} = one!(T(prec = precision(poly)))

Base.zero(series::T) where {T<:Union{ArbSeries,AcbSeries}} =
    T(degree(series), prec = precision(series))
Base.one(series::T) where {T<:Union{ArbSeries,AcbSeries}} =
    one!(T(degree(series), prec = precision(series)))

Base.zero(::Type{T}) where {T<:Union{ArbPoly,AcbPoly}} = T()
Base.one(::Type{T}) where {T<:Union{ArbPoly,AcbPoly}} = one!(T())

Base.zero(::Type{T}) where {T<:Union{ArbSeries,AcbSeries}} = T(0)
Base.one(::Type{T}) where {T<:Union{ArbSeries,AcbSeries}} = one!(T(0))


text:
#=struct TestPort
  Fret::Symbol
  Ffactors::Vector{Symbol}
end=#

const BLAS_MULTITHREAD = 1

struct FactorSpecification
    df::SubDataFrame #the data used
    Fret::Symbol
    Fdate::Symbol #date or period field
    Fport::Symbol

    #these are considered factors, generally used in the time series regressions
    Ffactors::Union{Nothing,Vector{Symbol}}

    #These are chracteristics, typically used for the cross-sectional regressions
    Fcharacteristics::Union{Nothing,Vector{Symbol}}

    Fcontrols::Vector{Symbol} #holds controls only

    #accessor variables
    eachport::GroupedDataFrame
    eachdate::GroupedDataFrame
    portidx::Dict{Union{Int,Symbol}, Int}
    dateidx::Dict{Date, Int}

    name::Symbol
end

#NOTE: Brackets index the test portfolios
Base.getindex(fs::FactorSpecification, sym::Symbol) = fs.eachport[portidx[sym]]
Base.getindex(fs::FactorSpecification, dt::Date) = fs.eachdate[dateidx[dt]]

#the below methods could probably be improved
Base.getindex(fs::FactorSpecification, syms::Symbol...) = view(
  fs.df, (s::Symbol-> s ∈ syms).(fs.df[!,fs.Fport]), :)
Base.getindex(fs::FactorSpecification, dts::Date...) = view(
  fs.df, (d::Date-> s ∈ dts).(fs.df[!,fs.Fdate]), :)

@inline function notsmallgroups(df::AbstractDataFrame, Fgroup::Symbol;
  minrowspergroup::Int = error("minrowspergroup must be assigned"))::Vector{Bool}

  dfcount::DataFrame = DataFrame(group = df[!,Fgroup])
  dfcount.count = Vector{MInt}(undef, size(dfcount,1))

  for sdfcount ∈ groupby(dfcount, :group)
    sdfcount.count .= size(sdfcount,1)
  end

  return dfcount.count .≥ minrowspergroup
end

#main constructor
function FactorSpecification(df::AbstractDataFrame, Fret::Symbol, Fdate::Symbol, Fport::Symbol;
    Ffactors::Union{Nothing, Vector{Symbol}}=nothing,
    Fcharacteristics::Union{Nothing, Vector{Symbol}}=nothing,
    bounds::Union{AbstractRange, Nothing} = nothing,
    Fcontrols = Vector{Symbol}(),
    name::NSymbol = nothing,
    minrowspergroup::Int = MIN_ROWS_PER_GROUP,
    skipcrosssectiondfcheck::Bool = false,
    skiptimeseriesdfcheck::Bool = false
    )

  local sdf::SubDataFrame
  local dfnames::Vector{Symbol} = names(df)

  #make sure the dataframe has everything we want
  isnothing(name) && error("name is a required key-word field!")
  (Fret ∈ dfnames) || error("return field $Fret not found in dataframe!")
  (Fdate ∈ dfnames) || error("datefield $Fdate not found in dataframe!")
  (Fport ∈ dfnames) || error("port field $Fport not found in dataframe!")
  (!isnothing(Ffactors)) && (
    (s->(s ∈ dfnames) || error("Factor $s not found in dataframe!")).(Ffactors))
  (!isnothing(Fcharacteristics)) && (
    (s->(s ∈ dfnames) || error("Factor $s not found in dataframe!")).(Fcharacteristics))
  (s->(s ∈ dfnames) || error("Control $s not found in dataframe!")).(Fcontrols)
  #println("num missing = $(size(df[(!).(completecases(df[!,[Fret,Fdate,Fport]])),:]))")
  (sum((!).(completecases(df[!,[Fret,Fdate,Fport]]))) == 0) || error("missing values found in identifiers")
  issorted(df, [Fport, Fdate]) || error("DataFrame must be sorted by $Fport and $Fdate")

  if isnothing(bounds) && typeof(df) == DataFrame
    sdf = view(df, :, :)
  elseif !isnothing(bounds)
    sdf = view(df, (d::MDate->
      (d ≥ minimum(bounds)) && (d≥maximum(bounds))).(df[!,Fdate]), :)
  else
    sdf = df
  end

  allfields::Vector{Symbol} = filter(s->!isnothing(s), [Fret; Fdate; Fport; Ffactors; Fcharacteristics; Fcontrols])
  sdf = view(sdf, completecases(sdf[!, allfields]), allfields)

  skiptimeseriesdfcheck || (sdf = view(sdf, notsmallgroups(sdf, Fport, minrowspergroup=minrowspergroup), :))
  skipcrosssectiondfcheck || (sdf = view(sdf, notsmallgroups(sdf, Fdate, minrowspergroup=minrowspergroup), :))

  (size(sdf,1)==0) && error("no rows left with reasonable sample size")

  #build subarrays of the most comon slices we will take for analysis
  eachport = groupby(sdf, Fport)
  eachdate = groupby(sdf, Fdate)
  portidx::Dict = Dict(eachport[i][1, Fport] => i for i ∈ 1:length(eachport))
  dateidx::Dict = Dict(eachdate[i][1, Fdate] => i for i ∈ 1:length(eachdate))

  return FactorSpecification(sdf, Fret, Fdate, Fport, Ffactors, Fcharacteristics, Fcontrols,
    eachport, eachdate, portidx, dateidx, name)
end

function rhs(fs::FactorSpecification, focalvars::Union{Vector{Symbol},Nothing};
    intercept::Bool = true)
  local rhs::Union{Symbol, Expr, Nothing}
  local rhsstring::String

  joinvec::Vector{NSymbol} = [focalvars; fs.Fcontrols;]
  joinvec = joinvec[(!isnothing).(joinvec)]

  rhsstring = join((string).(joinvec), " + ")
  (!intercept) && (rhsstring = "$rhsstring + 0")

  rhs = Meta.parse(rhsstring)

  return rhs
end

function rhsnames(fs::FactorSpecification, focalvars::Union{Vector{Symbol},Nothing};
    intercept::Bool = true)::Vector{Symbol}

  local nms::Vector{Union{Symbol,Nothing}}

  if intercept
    nms = [:intercept; focalvars; fs.Fcontrols;]
  else
    nms = [focalvars; fs.Fcontrols;]
  end

  return nms[(!isnothing).(nms)]
end

#creates the formulas
factorrhs(fs::FactorSpecification; intercept::Bool = true)::FMExpr = rhs(
  fs, fs.Ffactors, intercept=intercept)
factorrhsnames(fs::FactorSpecification; intercept::Bool = true)::Vector{Symbol} = rhsnames(
  fs, fs.Ffactors, intercept=intercept)
characteristicrhs(fs::FactorSpecification; intercept::Bool = true)::FMExpr = rhs(
  fs, fs.Fcharacteristics, intercept=intercept)
characteristicrhsnames(fs::FactorSpecification; intercept::Bool = true)::Vector{Symbol} = rhsnames(
  fs, fs.Fcharacteristics, intercept=intercept)

#run time series regressions with an arbitrary aggregation function
function timeseriesregressions(fs::FactorSpecification, ::Type{T} = FMLM;
  aggfunc::Function = (f(m::FMLM)::T=m),
  parallel::Bool = false)::Vector{T} where T

  Nreg::Int = length(fs.eachport)
  specs::FMSpecs = FMSpecs(Nreg, T, aggfunc=aggfunc)

  timeseriesrhs = factorrhs(fs)
  timeseriesrhsnames = factorrhsnames(fs)

  for i ∈ 1:length(fs.eachport)
    push!(specs, yspec = fs.Fret, xspec = timeseriesrhs, xnames = timeseriesrhsnames)
  end

  if parallel
    BLAS.set_num_threads(1)
  else
    BLAS.set_num_threads(BLAS_MULTITHREAD)
  end

  computeFMLMresults!(fs.eachport, specs, parallel=parallel, containsmissings=false)

  return specs.results
end

#run time series regressions with an arbitrary aggregation function
function crosssectionalregressions(fs::FactorSpecification, ::Type{T} = FMLM;
  aggfunc::Function = (f(m::FMLM)::T=m),
  parallel::Bool = false)::Vector{T} where T

  Nreg::Int = length(fs.eachport)
  specs::FMSpecs = FMSpecs(Nreg, T, aggfunc=aggfunc)

  crosssectionalrhs = characteristicrhs(fs)
  crosssectionalrhsnames = characteristicrhsnames(fs)

  for i ∈ 1:length(fs.eachdate)
    push!(specs, yspec = fs.Fret, xspec = crosssectionalrhs, xnames = crosssectionalrhsnames)
  end

  if parallel
    BLAS.set_num_threads(1)
  else
    BLAS.set_num_threads(BLAS_MULTITHREAD)
  end
  computeFMLMresults!(fs.eachdate, specs, parallel=parallel, containsmissings=false)

  #alldates::Vector{Date} = (df::SubDataFrame->df[1,fs.Fdate]).(fs.eachdate)

  #dateidx::Dict = Dict(zip(alldates, specs.results))

  return specs.results
end

struct FactorRegression
  λ::Union{Nothing, Matrix{Float64}}
  λᵢ::Vector{Float64}
  σᵢ::Vector{Float64}
  t::Int
  name::Symbol
  coefidx::Dict{Symbol, Int}
end

function manyFF(fs::FactorSpecification;
  parallel::Bool = false,
  errorfunction::Function = getNeweyWestFunc(2))::Vector{FactorRegression}

  local frs::Vector{FactorRegression} = Vector{FactorRegression}()

  #make the time-series regression aggregation function
  T::Type = FMLM
  @inline aggfunc(fmlm::FMLM) = fmlm

  results = timeseriesregressions(fs, T, aggfunc=aggfunc, parallel=parallel)
  coefnames::Vector{Symbol} = factorrhsnames(fs)
  coefidx = Dict{Symbol, Int}(zip(coefnames, 1:(length(coefnames))))

  #now build the factor regression objects, one for each regression
  for (i,port::SubDataFrame) ∈ enumerate(fs.eachport)
    local λᵢ::Vector{Float64} = results[i].β
    local σᵢ::Vector{Float64} = diag(errorfunction(results[i])).^0.5
    local t::Int = size(port,1)
    local name::Symbol = Symbol("$(fs.name)_col$(port[1,fs.Fport])")
    push!(frs, FactorRegression(nothing, λᵢ, σᵢ, t, name, coefidx))
  end

  return frs
end

function returnfactorFF(fs::FactorSpecification; parallel::Bool = false)::FactorRegression
  local t::Int
  local λ::Matrix{Float64}
  local λᵢ::Vector{Float64}
  local σᵢ::Vector{Float64}
  local name::Symbol
  local coefidx::Dict{Symbol, Int}

  #make the time-series regression aggregation function
  T::Type = Vector{Float64}
  @inline aggfunc(fmlm::FMLM) = fmlm.β

  results = timeseriesregressions(fs, T, aggfunc=aggfunc, parallel=parallel)
  coefnames::Vector{Symbol} = factorrhsnames(fs)
  coefidx = Dict{Symbol, Int}(zip(coefnames, 1:(length(coefnames))))

  #now build the factor regression object
  mλ = vcat((transpose).(results)...) #create rows and make result matrix
  origrows::Int = size(mλ,1)
  mλ = mλ[(i->isfinite(sum(mλ[i,:]))).(1:size(mλ,1)), :] #checks for na/missing and drops such rows
  newrows::Int = size(mλ,1)
  (origrows ≠ newrows) && (println("$(fs.name): removed " *
    "$(origrows - newrows) rows during DQ. $(newrows) remain."))
  #λ = DataFrame(mλ, coefnames)
  t = size(λ,1)
  λᵢ = (c->mean(c)).(eachcol(λ))
  σᵢ = (c->std(c)).(eachcol(λ)) ./ t^0.5
  name = fs.name

  return FactorRegression(mλ, λᵢ, σᵢ, t, name, coefidx)
end

function FMNeweyWest(fs::FactorSpecification,
  coefnames::Vector{Symbol},
  λᵢ::Vector{Float64},
  lags::Int)::Matrix{Float64}

  X::Matrix{Float64} = [ones(size(fs.df,1)) Matrix{Float64}(fs.df[!, coefnames[2:end]])]
  K::Int = length(λᵢ)
  R::Matrix{Float64} = qr(X).R
  Rinv::Matrix{Float64} = (R)\Matrix{Float64}(I,K,K) #this will save some computational time
  Sₜ::Matrix{Float64} = zeros(K,K) #holds teh central matrix

  #pre-allocate for the spectral matrix
  Σ::Matrix{Float64} = Matrix{Float64}(undef, K, K)

  temp::Matrix{Float64} = Matrix{Float64}(undef, K, K) #pre-allocate working matrix
  RRinv::Matrix{Float64} = BLAS.gemm('N', 'T', Rinv, Rinv) #this is equivelent to [X'X]^-1

  #iterate over all stocks
  for sdf ∈ fs.eachport
    #need to multiply through by the error
    T::Int = size(sdf,1)
    Xₙ::Matrix{Float64} = [ones(T) Matrix{Float64}(sdf[!, coefnames[2:end]])]
    ε::Vector{Float64} = sdf[!, fs.Fret] .- Xₙ * λᵢ
    Xₑ::Matrix{Float64} = Xₙ .* ε

    Sₜ::Matrix{Float64} += BLAS.gemm('T','N',1.0/T, Xₑ, Xₑ)

    for v::Int ∈ 1:lags
      #overwrites temp with (1/N)R'R
      BLAS.gemm!('T', 'N', 1.0/T, view(Xₑ, (v+1):T, :),view(Xₑ, 1:(T-v), :), 0.0, temp)
      Sₜ .+= (lags + 1 - v)/(lags+1.) .* (temp .+ temp')
    end

  end

  #this is [X'X]^-1S=[R'R]^-1S
  RRinvS::Matrix{Float64} = BLAS.gemm('N', 'N', RRinv, Sₜ)

  #finally we have T[X'X]^-1S[X'X]^-1
  BLAS.gemm!('N','N',Float64(size(X,1)), RRinvS, RRinv, 0.0, Σ)
  #println("$(diag(Σ .* dofCorrect))")sasa
  return Σ
end

#this version provides an intuitive check
function FMNeweyWestSlow(fs::FactorSpecification,
  coefnames::Vector{Symbol},
  λᵢ::Vector{Float64},
  lags::Int)::Matrix{Float64}

  X::Matrix{Float64} = [ones(size(fs.df,1)) Matrix{Float64}(fs.df[!, coefnames[2:end]])]

  #pre-allocate
  K::Int = size(X,2)
  Sₜ::Matrix{Float64} = zeros(K, K)
  #xxt::Matrix{Float64} = Matrix{Float64}(undef, K,K)
  temp::Matrix{Float64} = Matrix{Float64}(undef, K,K)

  #helper function using the Lars method
  @inline function Rₜ!(v::Int, X::Matrix{Float64}, ε::Vector{Float64};
      T::Int = size(X,1), out::Matrix{Float64} = Matrix{Float64}(undef, K,K))

    out .= 0.0
    for t ∈ (1+v):T
      out .+= X[t,:] * X[t-v,:]' * ε[t] * ε[t-v]
    end

    out ./= T
    return out
  end

  for sdf ∈ fs.eachport
    T::Int = size(sdf, 1)

    Xₙ::Matrix{Float64} = [ones(T) Matrix{Float64}(sdf[!, coefnames[2:end]])]
    ε::Vector{Float64} = sdf[!, fs.Fret] .- Xₙ * λᵢ
    Sₜ .+= Rₜ!(0, Xₙ, ε, out=temp)

    for v ∈ 1:lags #iterate over lags
      Rₜ!(v, Xₙ, ε, out=temp)
      Sₜ .+= (lags + 1. - v)/(lags + 1.) .* (temp .+ temp')
    end
  end

  XXInv::Matrix{Float64} = (X' * X) \ Matrix{Float64}(I,K,K)

  return size(X,1) * XXInv * Sₜ * XXInv
end

function characteristicFM(fs::FactorSpecification; parallel::Bool = false, lags=2)::FactorRegression
  local t::Int
  local λ::Matrix{Float64}
  local λᵢ::Vector{Float64}
  local σᵢ::Vector{Float64}
  local name::Symbol
  local coefidx::Dict{Symbol, Int}

  #make the time-series regression aggregation function
  T::Type = Vector{Float64}
  @inline aggfunc(fmlm::FMLM) = fmlm.β

  results = crosssectionalregressions(fs, T, aggfunc=aggfunc, parallel=parallel)
  coefnames::Vector{Symbol} = characteristicrhsnames(fs)
  coefidx = Dict{Symbol, Int}(zip(coefnames, 1:(length(coefnames))))

  #now build the factor regression object
  λ = vcat((transpose).(results)...) #create rows and make result matrix
  origrows::Int = size(λ,1)
  λ = λ[(i->isfinite(sum(λ[i,:]))).(1:size(λ,1)), :] #checks for na/missing and drops such rows
  newrows::Int = size(λ,1)
  (origrows ≠ newrows) && (println("$(fs.name): removed " *
    "$(origrows - newrows) rows during DQ. $(newrows) remain."))
  t = size(λ,1)
  λᵢ = (c->mean(c)).(eachcol(λ))

  Σ::Matrix{Float64} = FMNeweyWest(fs, coefnames, λᵢ, lags)
  σᵢ = diag(Σ) .^ 0.5

  println("FastNW: $σᵢ SlowNW: $(diag(FMNeweyWestSlow(fs, coefnames, λᵢ, lags)).^0.5) NormalSE: $((c->std(c)).(eachcol(λ)) ./ t)")

  name = fs.name

  return FactorRegression(λ, λᵢ, σᵢ, t, name, coefidx)
end

#=
86-50-31
struct FactorRegressions
  frs::Vector{FactorRegression}
  nameidx::Dict
  fssidx::Dict{}
end

function FactorRegressions=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@info "Test template pattern..."

println("Brew coffee process")
coffee = Coffee("", "")
prepare(coffee)

println()

println("Brew tea process")
tea = Tea("", "")
prepare(tea)


text:
using PyCall
@pyimport modeller

env = modeller.environ()
aln = modeller.alignment(env)
for (pdb, chain) in (("1b8p", "A"), ("1bdm", "A"), ("1civ", "A"),
                     ("5mdh", "A"), ("7mdh", "A"), ("1smk", "A"))
    m = modeller.model(env, file=pdb, model_segment=("FIRST:$chain", "LAST:$chain"))
    aln[:append_model](m, atom_files=pdb, align_codes=string(pdb, chain))
end
aln[:malign]()
aln[:malign3d]()
aln[:compare_structures]()
aln[:id_table](matrix_file="family.mat")
env[:dendrogram](matrix_file="family.mat", cluster_cut=-1.0)


text:
import Base: show
using Oceananigans.Utils: prettysummary

const DFBC = DefaultPrognosticFieldBoundaryCondition

bc_str(::FBC)     = "Flux"
bc_str(::PBC)     = "Periodic"
bc_str(::OBC)     = "Open"
bc_str(::VBC)     = "Value"
bc_str(::GBC)     = "Gradient"
bc_str(::ZFBC)    = "ZeroFlux"
bc_str(::DFBC)    = "Default"
bc_str(::Nothing) = "Nothing"

#####
##### BoundaryCondition
#####

Base.summary(bc::DFBC) = string("DefaultBoundaryCondition")
Base.summary(bc::PBC) = string("PeriodicBoundaryCondition")
Base.summary(bc::OBC) = string("OpenBoundaryCondition: ", prettysummary(bc.condition))
Base.summary(bc::FBC) = string("FluxBoundaryCondition: ", prettysummary(bc.condition))
Base.summary(bc::VBC) = string("ValueBoundaryCondition: ", prettysummary(bc.condition))
Base.summary(bc::GBC) = string("GradientBoundaryCondition: ", prettysummary(bc.condition))

show(io::IO, bc::BoundaryCondition) = print(io, summary(bc))

#####
##### FieldBoundaryConditions
#####

Base.summary(fbcs::FieldBoundaryConditions) = "FieldBoundaryConditions"

show_field_boundary_conditions(bcs::FieldBoundaryConditions, padding="") =
    string("Oceananigans.FieldBoundaryConditions, with boundary conditions", '\n',
           padding, "├── west: ",     summary(bcs.west), '\n',
           padding, "├── east: ",     summary(bcs.east), '\n',
           padding, "├── south: ",    summary(bcs.south), '\n',
           padding, "├── north: ",    summary(bcs.north), '\n',
           padding, "├── bottom: ",   summary(bcs.bottom), '\n',
           padding, "├── top: ",      summary(bcs.top), '\n',
           padding, "└── immersed: ", summary(bcs.immersed))

Base.show(io::IO, fieldbcs::FieldBoundaryConditions) = print(io, show_field_boundary_conditions(fieldbcs))


text:
_mkdir(path) = mkpath(dirname(path))

_readdir(dir; join::Bool = false, sort::Bool = true) =
    isdir(dir) ? readdir(dir; join, sort) : String[]

function _readdir(f::Function, dir; kwargs...)
    for file in _readdir(dir; kwargs...)
        f(file)
    end
end

function _rm(path)
    try; rm(path; recursive = true, force = true)
    catch err; end
end

function _foldersize(dir)
    size = 0
    for (root, _, files) in walkdir(dir)
        size += sum(filesize.(joinpath.(root, files)))
    end
    return size
end

function _cp(src::AbstractString, dst::AbstractString)
    try; cp(src, dst; force = true)
    catch err; end
end

# get the rel path from basename(startpath)
function _relbasepath(path, startpath)
	stoppath = basename(startpath)

	relpath_ = ""
	leftpath = path
	while true
		leftpath, base = splitdir(leftpath)
		(base == stoppath) && break
		relpath_ = isempty(relpath_) ? base : joinpath(base, relpath_)
		(leftpath == dirname(leftpath)) && break
	end

	return relpath_
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
include("../../CoordTrans/CoordTrans_ESMF.jl")
include("./CESMReader.jl")

using .CoordTrans_ESMF
using .CESMReader
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")

function mreplace(x)
    return replace(x, missing=>NaN)
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

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

        "--ESMF-wgt-file"
            help = "ESMF generated domain file."
            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))

wi = CoordTrans_ESMF.readWeightInfo(parsed["ESMF-wgt-file"])

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

    global mask = ds["mask"][:] |> mreplace
    global area = ds["area"][:] |> mreplace
    global lat = ds["yc"][1, :] |> mreplace

    global lat_weight = 2π * Re^2.0 * cos.(deg2rad.(lat))
    global ϕ          = deg2rad.(lat)

    global is_ocn = mask .== 0.0
end

let

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

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

    _TFLUX_DIV_implied = getData(fh, "TFLUX_DIV_implied", (parsed["beg-year"], parsed["end-year"]), (:, :))

    global Nx, Ny, Nt = size(_TFLUX_DIV_implied)

    if mod(Nt, 12) != 0
        throw(ErrorException("Time should be a multiple of 12."))
    end

    global nyears = Int64(Nt / 12)
    println(format("We got {:d} years of data.", nyears))

    # Do fraction weighting
    _ECONV = zeros(Float64, Ny, Nt)

    @time for t=1:Nt
        sum_flux = 0.0
        sum_area = 0.0
        for i=1:Nx, j=1:Ny
            b = (i-1) + Nx * (j-1) + 1
            if is_ocn[i, j]
                _ECONV[j, t] += _TFLUX_DIV_implied[i, j, t] #* wi.frac_b[b]
                sum_flux +=  _TFLUX_DIV_implied[i, j, t] * area[i, j]# * wi.frac_b[b]
                sum_area += area[i, j] #* wi.frac_b[b]
            end
        end

       println(t, " : ", sum_flux * ρc / sum_area)

    end

    _ECONV .*= ρc / Nx

    global IOET = zeros(Float64, Ny, Nt)

    for t = 1:Nt
        IOET[:, t] = integrate(ϕ, _ECONV[:, t] .* lat_weight)
    end

end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("IOET", IOET, ("Ny", "time"), Dict()),
        ("IOET_AM",    mean(IOET, dims=(2,))[:, 1], ("Ny",), Dict()),
        ("IOET_AMSTD",  std(IOET, dims=(2,))[:, 1], ("Ny",), Dict()),
    ]

        println("Doing var: ", varname)

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

        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


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule pigz_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("pigz")
JLLWrappers.@generate_main_file("pigz", UUID("1bc43ea1-30af-5bc8-a9d4-c018457e6e3e"))
end  # module pigz_jll


text:
#
# Factors Tests
#
# constructors, basic functions

@testset "Factors Main" begin

@testset "bad constructors" begin
    c = Dimension(:C, ["as", "ab", "cd", "sd"])
    o = Dimension(:O, ['2', 'J', 'Q', 'K', 'A'])
    s = Dimension(:S, 'a':2:'z')
    l1 = Dimension(:L, 3)
    l2 = Dimension(:L, [3, 4])

    @test_throws ArgumentError Factor([c, l1, o, s, l2])
    @test_throws DimensionMismatch Factor([c, o], rand(4, 6))
    @test_throws DimensionMismatch Factor([c, o], rand(4, 6, 3))
    @test_throws DimensionMismatch Factor(rand(4, 6, 3), :X=>3, :Y=>3, :Z=>4)
    @test_throws DimensionMismatch Factor(rand(4, 6, 3), :X=>3, :Y=>3)

    @test_throws ArgumentError Factor(Dimension(:X, []), Float64[])
end

@testset "basics" begin
    c = Dimension(:C, ["as", "ab", "cd", "sd"])
    o = Dimension(:O, ['2', 'J', 'Q', 'K', 'A'])
    s = Dimension(:S, 'a':2:'z')
    u = Dimension(:U, 0:15)
    l = Dimension(:L, 31)

    ft = Factor([c, o, s])
    @test ndims(ft) == 3
    @test indexin([:C], ft) == [1]
    @test indexin(:O, ft) == 2
    @test indexin(o, ft) == 2
    @test getdim(ft, :S) == s
    @test size(ft) == (4, 5, 13)
    @test all(values(ft) .== 0)
    @test eltype(ft) == (Float64)
    @test all(names(ft) .== [:C, :O, :S])
    @test length(ft) == prod(size(ft))
    @test getdim(ft, :S) == s

    ft = push(ft, l)

    @test ndims(ft) == 4
    @test size(ft) == (4, 5, 13, 31)
    @test size(ft, :S) == 13
    @test size(ft, :O) == 5
    @test size(ft, :L) == 31
    @test size(ft, :L, :S, :O) == (31, 13, 5)
    @test all(ft.potential .== 0)
    @test names(ft) == [:C, :O, :S, :L]
    @test length(ft) == prod(size(ft))
    @test getdim(ft, :L) == l
    @test indexin(:C, ft) == 1
    @test indexin([:L], ft) == [4]
    @test ft.dimensions[4] == l

    ft2 = permutedims(ft, [3, 2, 1, 4])

    @test ndims(ft2) == 4
    @test size(ft2) == (13, 5, 4, 31)
    @test all(ft2.potential .== 0)
    @test all(names(ft2) .== [:S, :O, :C, :L])
    @test length(ft2) == prod(size(ft))
    @test indexin(:S, ft2) == 1
    @test indexin(:C, ft2) == 3
    @test getdim(ft2, :S) == getdim(ft, :S)
    @test getdim(ft2, :L) == getdim(ft, :L)
    @test ft.dimensions[1] == ft2.dimensions[3]
    @test ft.dimensions[3] == ft2.dimensions[1]
    @test ft2.dimensions[2] == o

    permutedims!(ft, [3, 2, 1, 4])
    @test names(ft) == names(ft2)
end

@testset "basics 2" begin
    l1 = Dimension(:l1, 2)
    l2 = Dimension(:l2, 3)
    ft = Factor([l1, l2], 1:6)

    @test eltype(ft) == Float64
    @test ndims(ft) == 2
    @test size(ft) == (2, 3)
    @test all(ft.potential .== [1.0 3 5; 2 4 6])
    @test all(names(ft) .== [:l1, :l2])
    @test length(ft) == 6

    @test sub2ind(ft, 2, 3) == 6
    @test sub2ind(ft, 2, 1) == 2
    @test sub2ind(ft, 2, 2) == 4

    @test ind2sub(ft, 5) == (1, 3)
    @test ind2sub(ft, 3) == (1, 2)

    @test all(pattern(ft, :l1) .== [1, 2, 1, 2, 1, 2])
    @test all(pattern(ft, :l2) .== [1, 1, 2, 2, 3, 3])
end

@testset "empty" begin
    # see if any puke up errors
    Factor(2016)
    Factor(Inf)
    Factor(NaN)
    f = Factor(31)
    @test ndims(f) == 0
    @test names(f) == Symbol[]
    @test scope(f) == Dimension[]
    @test length(f) == 1
    @test size(f) == ()
end

@testset "pairs" begin
    @testset "default" begin
        ft = Factor(:S => 2, :X=>3, :Z=>4)
        @test ndims(ft) == 3
        @test size(ft) == (2, 3, 4)
        @test all(values(ft) .== 0)
        @test all(names(ft) .== [:S, :X, :Z])
        @test length(ft) == 24
    end

    @testset "empty" begin
        ft = Factor(nothing, :S => 2, :X=>3, :Z=>4)
        @test ndims(ft) == 3
        @test size(ft) == (2, 3, 4)
        @test all(names(ft) .== [:S, :X, :Z])
        @test length(ft) == 24
    end

    @testset "init" begin
        ft = Factor(2016, :S => 2, :X=>3, :Z=>4)
        @test size(ft) == (2, 3, 4)
        @test all(ft.potential .== 2016)
    end

    @testset "reshape" begin
        ft = Factor(1:24, :S => 2, :X=>3, :Z=>4)
        @test size(ft) == (2, 3, 4)
        @test all(ft.potential[:] .== collect(1:24))
    end

    @test_throws DimensionMismatch Factor(1:13, :X=>3, :Y=>4)
end

@testset "dict" begin
    @testset "default" begin
        ft = Factor(Assignment(:S => 2, :X=>3, :Z=>4))
        @test ndims(ft) == 3
        @test size(ft, :S, :X, :Z) == (2, 3, 4)
        @test all(values(ft) .== 0)
        @test length(ft) == 24
    end

    @testset "empty" begin
        ft = Factor(Assignment(:X => 3, :S=>2, :Z=>4), nothing)
        @test size(ft, :S, :X, :Z) == (2, 3, 4)
    end

    @testset "init" begin
        ft = Factor(Assignment(:X => 3, :S=>2, :Z=>4), 2016)
        @test all(ft.potential .== 2016)
    end
end

@testset "1 dimensional" begin
    @testset "default" begin
        ft = Factor(:S, 11:14)
        @test ndims(ft) == 1
        @test size(ft) == (4,)
        @test all(values(ft) .== [11, 12, 13, 14])
        @test length(ft) == 4
    end

    @testset "pairs" begin
        ft = Factor(11:14, :S => 4)
        @test ndims(ft) == 1
        @test size(ft) == (4,)
        @test all(values(ft) .== [11, 12, 13, 14])
        @test length(ft) == 4
    end

    @testset "default" begin
        ft = Factor(33, :S=>2:6)
        @test ndims(ft) == 1
        @test size(ft) == (5,)
        @test all(values(ft) .== 33)
    end
end

@testset "pattern states" begin
    c = Dimension(:C, ["bob", "bill"])
    o = Dimension(:O, ['a', 'b', 'c'])

    f = Factor([c, o])
    expected = ["bob" 'a';
                "bill" 'a';
                "bob" 'b';
                "bill" 'b';
                "bob" 'c';
                "bill" 'c']
    @test all(expected .== pattern_states(f))
end

end


text:
using StringParserPEG

calc2 = Grammar("""
  start => (number & op & number){(r,v,f,l,c) -> c[2](c[1],c[3])}

  op => plus | minus
  number => (-(space) & r([0-9]+)r) {(r,v,f,l,c) -> parse(Int,c[1].value)}
  plus => (-(space) & '+'){(a...) -> +}
  minus => (-(space) & '-'){(a...) -> -}
  space => r([ \\t\\n\\r]*)r
""")

data = "4+5"

(ast, pos, err) = parse(calc2, data)
ast