text:
using DataInterpolations, Test
using FiniteDifferences
using DataInterpolations: derivative

function test_derivatives(func, tspan, name::String)
  trange = range(minimum(tspan), maximum(tspan), length=32)[2:end-1]
  @testset "$name" begin
    for t in trange
      # Linearly spaced points might lead to evaluations outside
      # trange
      cdiff = central_fdm(5, 1; geom=true)(_t -> func(_t), t)
      adiff = derivative(func, t)
      @test isapprox(cdiff, adiff, atol=1e-8)
    end
  end
end

# Linear Interpolation
u = 2.0collect(1:10)
t = 1.0collect(1:10)
A = LinearInterpolation(u,t)

test_derivatives(A, t, "Linear Interpolation (Vector)")

u = vcat(2.0collect(1:10)', 3.0collect(1:10)')
A = LinearInterpolation(u,t)

test_derivatives(A, t, "Linear Interpolation (Matrix)")

# Quadratic Interpolation
u = [1.0, 4.0, 9.0, 16.0]
t = [1.0, 2.0, 3.0, 4.0]
A = QuadraticInterpolation(u,t)

test_derivatives(A, t, "Quadratic Interpolation (Vector)")

u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
A = QuadraticInterpolation(u,t)

test_derivatives(A, t, "Quadratic Interpolation (Matrix)")

# Lagrange Interpolation
u = [1.0, 4.0, 9.0]
t = [1.0, 2.0, 3.0]
A = LagrangeInterpolation(u,t)

test_derivatives(A, t, "Lagrange Interpolation (Vector)")

# Lagrange Interpolation
u = [1.0 4.0 9.0; 1.0 2.0 3.0]
t = [1.0, 2.0, 3.0]
A = LagrangeInterpolation(u,t)

test_derivatives(A, t, "Lagrange Interpolation (Matrix)")

# Lagrange Interpolation
u = [[1.0, 4.0, 9.0], [3.0, 7.0, 4.0], [5.0, 4.0, 1.0]]
t = [1.0, 2.0, 3.0]
A = LagrangeInterpolation(u,t)

test_derivatives(A, t, "Lagrange Interpolation (Vector of Vectors)")

# Lagrange Interpolation
u = [[3.0 1.0 4.0; 1.0 5.0 9.0], [2.0 6.0 5.0; 3.0 5.0 8.0], [9.0 7.0 9.0; 3.0 2.0 3.0]]
t = [1.0, 2.0, 3.0]
A = LagrangeInterpolation(u,t)

test_derivatives(A, t, "Lagrange Interpolation (Vector of Matrices)")

# Akima Interpolation
u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
t = collect(0.0:10.0)
A = AkimaInterpolation(u, t)

test_derivatives(A, t, "Akima Interpolation")

# QuadraticSpline Interpolation
u = [0.0, 1.0, 3.0]
t = [-1.0, 0.0, 1.0]

A = QuadraticSpline(u,t)

test_derivatives(A, t, "Quadratic Interpolation (Vector)")

# QuadraticSpline Interpolation
u = [[1.0, 2.0, 9.0], [3.0, 7.0, 5.0], [5.0, 4.0, 1.0]]
t = [-1.0, 0.0, 1.0]

A = QuadraticSpline(u,t)

test_derivatives(A, t, "Quadratic Interpolation (Vector of Vectors)")

# QuadraticSpline Interpolation
u = [[1.0 4.0 9.0; 5.0 9.0 2.0], [3.0 7.0 4.0; 6.0 5.0 3.0], [5.0 4.0 1.0; 2.0 3.0 8.0]]
t = [-1.0, 0.0, 1.0]

A = QuadraticSpline(u,t)

test_derivatives(A, t, "Quadratic Interpolation (Vector of Matrices)")

# CubicSpline Interpolation
u = [0.0, 1.0, 3.0]
t = [-1.0, 0.0, 1.0]

A = CubicSpline(u,t)

test_derivatives(A, t, "Cubic Spline Interpolation (Vector)")

# CubicSpline Interpolation
u = [[1.0, 2.0, 9.0], [3.0, 7.0, 5.0], [5.0, 4.0, 1.0]]
t = [-1.0, 0.0, 1.0]

A = CubicSpline(u,t)

test_derivatives(A, t, "Cubic Spline Interpolation (Vector of Vectors)")

# CubicSpline Interpolation
u = [[1.0 4.0 9.0; 5.0 9.0 2.0], [3.0 7.0 4.0; 6.0 5.0 3.0], [5.0 4.0 1.0; 2.0 3.0 8.0]]
t = [-1.0, 0.0, 1.0]

A = CubicSpline(u,t)

test_derivatives(A, t, "Cubic Spline Interpolation (Vector of Matrices)")

# BSpline Interpolation and Approximation
t = [0,62.25,109.66,162.66,205.8,252.3]
u = [14.7,11.51,10.41,14.95,12.24,11.22]

A = BSplineInterpolation(u,t,2,:Uniform,:Uniform)

test_derivatives(A, t, "BSpline Interpolation (Uniform, Uniform)")

A = BSplineInterpolation(u,t,2,:ArcLen,:Average)

test_derivatives(A, t, "BSpline Interpolation (Arclen, Average)")

A = BSplineApprox(u,t,3,4,:Uniform,:Uniform)

test_derivatives(A, t, "BSpline Approx (Uniform, Uniform)")


text:
"""
    ConcreteStructs

ConcreteStructs exports the macro `@concrete`, which can be used to concretely parameterize
structs.
"""
module ConcreteStructs

export @concrete

"""
    @concrete struct
    @concrete mutable struct
    @concrete terse struct
    @concrete terse mutable struct

The `@concrete` macro makes non-concrete structs concrete, saving the boilerplate of making
type parameters. The `terse` keyword causes the types to show without their parameters while
in `:compact => true` mode of an `IOContext`.

## Examples
```julia
julia> using ConcreteStructs

julia> @concrete struct AB
           a
           b
       end

julia> ab = AB("hi", 1+im)
AB{String,Complex{Int64}}("hi", 1 + 1im)

julia> @concrete terse mutable struct CDE{D} <: Number
            d::D
            c
            e::Symbol
        end

julia> cde = CDE(1f0, (1,2.0), :yo)
CDE(1.0f0, (1, 2.0), :yo)

julia> typeof(cde)
CDE{Float32,Tuple{Int64,Float64}}

julia> @concrete terse struct FGH{T,N,G<:AbstractArray{T,N}} <: AbstractArray{T,N}
           f
           g::G
           h::T
       end

julia> Base.size(x::FGH) = size(x.g); Base.getindex(x::FGH, i...) = getindex(x.g[i...])

julia> fgh = FGH(nothing, [1,2,3], 4)
3-element FGH:
 1
 2
 3

julia> typeof(fgh)
FGH{Int64,1,Array{Int64,1},Nothing}
```
"""
macro concrete(expr)
    expr, _, _ = _concretize(expr)
    return esc(expr)
end

macro concrete(terse, expr)
    terse isa Symbol && terse == :terse || error("Invalid usage of @concrete")

    expr, struct_name, type_params = _concretize(expr)
    struct_name = string(struct_name)
    num_params = length(type_params)

    terse_string = if num_params == 0
        struct_name
    else
        :($(struct_name) * "{" * join(T.parameters[1:$num_params], ",") * "}")
    end

    return quote
        $(Base).@__doc__ $expr
        function Base.show(io::IO, T::Type{<:$(Symbol(struct_name))})
            if T isa UnionAll
                print(io, $struct_name)
            else
                print(io, $(struct_name) * "{" * join(T.parameters[1:$num_params], ",") * "}")
            end
        end
        function Base.show(io::IO, ::MIME"text/plain", T::Type{<:$(Symbol(struct_name))})
            if T isa UnionAll
                return print(io, $struct_name)
            else
                return print(io, $(struct_name) * "{" * join(T.parameters, ",") * "}")
            end
        end
    end |> esc
end

# Parse whole struct definition for the @concrete macro
function _concretize(expr)
    expr isa Expr && expr.head == :struct || error("Invalid usage of @concrete")

    is_mutable = expr.args[1]
    struct_name, type_params, super = _parse_head(expr.args[2])
    line_tuples = _parse_line.(expr.args[3].args)
    lines = first.(line_tuples)
    type_params_full = (type_params..., filter(x -> x!==nothing, last.(line_tuples))...)

    struct_type = if length(type_params_full) == 0
        struct_name
    else
        Expr(:curly, struct_name, type_params_full...)
    end

    head = Expr(:(<:), struct_type, super)
    constructor_expr = _make_constructor(struct_name, type_params, type_params_full, lines)
    body = Expr(:block, lines..., constructor_expr)
    struct_expr = Expr(:struct, is_mutable, head, body)

    return struct_expr, struct_name, type_params
end

# Make the inner constructor function
function _make_constructor(struct_name, type_params, type_params_full, lines)
    lines = map(line->line isa Expr && line.head==:(=) ? line.args[1] : line, lines)
    field_lines = filter(line -> ((line isa Expr) && (line.head === :(::))), lines)
    args = map(x->x.args, field_lines)
    vars = first.(args)
    var_types = last.(args)
    constructor_params = _get_constructor_params(type_params, var_types)
    new_params = _strip_super(type_params_full)

    if length(type_params) == length(type_params_full) && all(type_params .== type_params_full)
        return Expr(:block)
    elseif length(constructor_params)==0
        return :(
            function $struct_name($(field_lines...)) where {$(type_params_full...)}
                return new{$(new_params...)}($(vars...))
            end
        )
    else
        return :(
            function $struct_name{$(constructor_params...)}($(field_lines...)) where {$(type_params_full...)}
                return new{$(new_params...)}($(vars...))
            end
        )
    end
end

# Get the parameters that are unmatched to variables and need to be annoted in the constructor
function _get_constructor_params(type_params, var_types)
    subparams = _get_subparams(type_params)
    type_params = _strip_super(type_params)
    var_types = [subparams; _strip_super(var_types)]
    return setdiff(type_params, var_types)
end

# Strip supertype annotations
_strip_super(x) = x
_strip_super(x::Union{Tuple, AbstractVector}) = vcat(_strip_super.(x)...)
_strip_super(x::Expr) = x.head == :(<:) ? x.args[1] : x

# Get the subparameters of supertypes of subtype parameters (sorry)
_get_subparams(x) = []
_get_subparams(x::Union{Tuple, AbstractVector}) = vcat(_get_subparams.(x)...)
function _get_subparams(x::Expr)
    if x.head === :curly
        return x.args[2:end]
    elseif x.head === :(<:)
        return _get_subparams(x.args[2:end])
    end
end

# Parse the top line of the struct definition
_parse_head(head::Symbol) = (_parse_struct_def(head)..., :(Any))
function _parse_head(head::Expr)
    if head.head === :curly
        super = :(Any)
        struct_name, type_params = _parse_struct_def(head)
    elseif head.head === :(<:)
        super = head.args[2]
        struct_name, type_params = _parse_head(head.args[1])
    end

    return (struct_name, type_params, super)
end

# Parse the struct name and parameters
_parse_struct_def(struct_def::Symbol) = (struct_def, [])
_parse_struct_def(struct_def::Expr) = (struct_def.args[1], struct_def.args[2:end])

# Parse a line of the body of the struct def. Returns the line and the type parameter to be
# included in the struct header
_parse_line(line) = (line, nothing)
function _parse_line(line::Expr)
    assignment = line.head === :(=)
    annotation = nothing
    if assignment
        val = line.args[2]
        line, annotation = _parse_line(line.args[1])
    end

    out = if line isa Expr && line.head === :(<:)
        field = line.args[1]
        T = line.args[2]
        sym = Symbol(:__T_, field)
        (:($field::$sym), :($sym<:$T))
    else
        (line, annotation)
    end

    return if assignment
        (:($(out[1])=$val), out[2])
    else
        out
    end

    # if (T isa Symbol && isdefined(Base.Main, T) && !isconcretetype(Base.eval(Base.Main, T))) || T isa Expr
    #     sym = Symbol("__T_" * string(field))
    #     return (:($field::$sym), :($sym<:$T))
    # else
    #     return (line, nothing)
    # end
end
function _parse_line(line::Symbol)
    T = Symbol("__T_" * string(line))
    return (:($line::$T), T)
end

end

ext: jl


text:
module PostgreSQL
    export  Postgres,
            executemany,
            copy_from,
            escapeliteral

    using BinDeps
    @BinDeps.load_dependencies

    include("libpq_interface.jl")
    using .libpq_interface
    using DBI
    using DataFrames
    using DataArrays

    include("types.jl")
    include("dbi_impl.jl")
end

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
__precompile__()

module NLsolve

using Distances
using NLSolversBase
using LineSearches
using LinearAlgebra
using Printf

import Base.show,
       Base.push!,
       Base.getindex,
       Base.setindex!

import NLSolversBase: OnceDifferentiable, InplaceObjective, NotInplaceObjective,
       only_fj, only_fj!

using Reexport
@reexport using LineSearches
using LinearAlgebra

export OnceDifferentiable,
       n_ary,
       nlsolve,
       mcpsolve,
       converged,
       only_fj,
       only_fj!,
       fixedpoint

abstract type AbstractSolverCache end

struct IsFiniteException <: Exception
  indices::Vector{Int}
end
show(io::IO, e::IsFiniteException) = print(io,
  "During the resolution of the non-linear system, the evaluation" *
  " of the following equation(s) resulted in a non-finite number: $(e.indices)")

include("objectives/helpers.jl")

include("solvers/newton.jl")
include("solvers/broyden.jl")
include("solvers/trust_region.jl")
include("solvers/anderson.jl")
include("solvers/mcp_func_defs.jl")
include("solvers/mcp.jl")

include("nlsolve/solver_state_results.jl")
include("nlsolve/nlsolve.jl")
include("nlsolve/utils.jl")
include("nlsolve/fixedpoint.jl")

end # module


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

__precompile__()

module Optimization

using Revise

export setoptim, printerrors, linearfit, linearfit2, optimize, optimize2, validate, validate2
export mpolyeval, mpolyinterp, gradientdescent, num2string, replacestring

function deletenothing(a)
    if a !== nothing
        n = length(a)
        if n > 0
            m = 0
            for i = n:-1:1
                if a[i] === nothing
                    m = m + 1
                    deleteat!(a, i)
                end
            end
            if m==n
                a = nothing
            end
        end
    end
    return a
end

include("printerrors.jl");
include("OptimStruct.jl");
include("tensorpolyfit.jl");
include("tensorpolyopt.jl");
include("linearfit.jl");
include("maespace.jl");

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
"""
A wrapper of `TransientFEOperator` that transforms it to `ODEOperator`, i.e.,
takes A(t,uh,∂tuh,vh) and returns A(t,uF,∂tuF) where uF and ∂tuF represent the
free values of the `EvaluationFunction` uh and ∂tuh.
"""
struct ODEOpFromFEOp{C} <: ODEOperator{C}
  feop::TransientFEOperator{C}
end

function allocate_cache(op::ODEOpFromFEOp)
  Ut = get_trial(op.feop)
  Vt = ∂t(Ut)
  U = allocate_trial_space(Ut)
  V = allocate_trial_space(Vt)
  fecache = allocate_cache(op.feop)
  ode_cache = (U,V,Ut,Vt,fecache)
  ode_cache
end

function update_cache!(ode_cache,op::ODEOpFromFEOp,t::Real)
  U,V,Ut,Vt,fecache = ode_cache
  U = evaluate!(U,Ut,t)
  V = evaluate!(V,Vt,t)
  fecache = update_cache!(fecache,op.feop,t)
  (U,V,Ut,Vt,fecache)
end

function allocate_residual(op::ODEOpFromFEOp,uhF::AbstractVector,ode_cache)
  U,V,Ut,Vt,fecache = ode_cache
  uh = EvaluationFunction(U,uhF)
  allocate_residual(op.feop,uh,fecache)
end

function allocate_jacobian(op::ODEOpFromFEOp,uhF::AbstractVector,ode_cache)
  U,V,Ut,Vt,fecache = ode_cache
  uh = EvaluationFunction(U,uhF)
  allocate_jacobian(op.feop,uh,fecache)
end

function residual!(b::AbstractVector,op::ODEOpFromFEOp,t::Real,uhF::AbstractVector,uhtF::AbstractVector,ode_cache)
  Uh,Uht, = ode_cache
  uh = EvaluationFunction(Uh,uhF)
  uht = EvaluationFunction(Uht,uhtF)
  residual!(b,op.feop,t,uh,uht,ode_cache)
end

function jacobian!(A::AbstractMatrix,op::ODEOpFromFEOp,t::Real,uhF::AbstractVector,uhtF::AbstractVector,ode_cache)
  Uh,Uht, = ode_cache
  uh = EvaluationFunction(Uh,uhF)
  uht = EvaluationFunction(Uht,uhtF)
  jacobian!(A,op.feop,t,uh,uht,ode_cache)
end

function jacobian_t!(J::AbstractMatrix,op::ODEOpFromFEOp,t::Real,uhF::AbstractVector,uhtF::AbstractVector,dut_u::Real,ode_cache)
  Uh,Uht, = ode_cache
  uh = EvaluationFunction(Uh,uhF)
  uht = EvaluationFunction(Uht,uhtF)
  jacobian_t!(J,op.feop,t,uh,uht,dut_u,ode_cache)
end

function jacobian_and_jacobian_t!(J::AbstractMatrix,op::ODEOpFromFEOp,t::Real,uhF::AbstractVector,uhtF::AbstractVector,dut_u::Real,ode_cache)
  Uh,Uht, = ode_cache
  uh = EvaluationFunction(Uh,uhF)
  uht = EvaluationFunction(Uht,uhtF)
  jacobian_and_jacobian_t!(J,op.feop,t,uh,uht,dut_u,ode_cache)
end


text:
module NumberSuffix

    export SuffixIt

    function SuffixIt(num)::String
        if (typeof(num)!=Int64)
            throw(TypeError(num,"cannot suffix this type of argument",Int64))
        end
        if (num<0)
            throw(DomainError(x, "argument must be nonnegative"))
        end
        suffixedNumber=string(num)
        lastDigit=num%10
        if (lastDigit==0 || lastDigit>=4)
            suffixedNumber=suffixedNumber*"th"
        elseif (lastDigit==1)
            suffixedNumber=suffixedNumber*"st"
        elseif (lastDigit==2)
            suffixedNumber=suffixedNumber*"nd"
        else
            suffixedNumber=suffixedNumber*"rd"
        end
        return suffixedNumber
    end

end


text: # @test  res[3] |> success


text:
# Code for accessing/indexing KeyedDataset components
"""
    getproperty(ds::KeyedDataset, sym::Symbol)

Extract [`KeyedDataset`](@ref) fields, dimension keys or components in that order.
Shared axis keys are wrapped in a `ReadOnlyArray` for safety.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(zeros(3, 2); time=1:3, obj=[:a, :b]),
           :val2 => KeyedArray(ones(3, 2) .+ 1.0; time=1:3, loc=[:x, :y]),
       );

julia> collect(keys(ds.data))
2-element Vector{Tuple}:
 (:val1,)
 (:val2,)

julia> ds.time
3-element ReadOnlyArrays.ReadOnlyArray{Int64, 1, UnitRange{Int64}}:
 1
 2
 3

julia> dimnames(ds.val1)
(:time, :obj)
```
"""
function Base.getproperty(ds::KeyedDataset, sym::Symbol)
    symkey = (sym,)
    if sym in fieldnames(KeyedDataset)
        return getfield(ds, sym)
    elseif sym in dimnames(ds)
        _keys = axiskeys(ds, sym)
        if length(_keys) == 1
            return ReadOnlyArray(first(_keys))
        else
            throw(ArgumentError(
                "$sym is an ambiguous dimension in the dataset. " *
                "More than one set of values exist for this dimension: $_keys"
            ))
        end
    elseif haskey(ds.data, symkey)
        return ds[symkey]
    else
        throw(ArgumentError("type KeyedDataset has no field $sym"))
    end
end

"""
    getindex(ds::KeyedDataset, key)

Lookup [`KeyedDataset`](@ref) component by its `Tuple` key, or `Symbol` for keys of depth 1.
Shared axis keys for the returned `KeyedArray` are wrapped in a `ReadOnlyArray` for safety.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset;

julia> ds = KeyedDataset(
           :val1 => KeyedArray(zeros(3, 2); time=1:3, obj=[:a, :b]),
           :val2 => KeyedArray(ones(3, 2) .+ 1.0; time=1:3, loc=[:x, :y]),
       );

julia> ds[:val1]
2-dimensional KeyedArray(NamedDimsArray(...)) with keys:
 ↓   time ∈ 3-element ReadOnlyArrays.ReadOnlyArray{Int64,...}
  →   obj ∈ 2-element ReadOnlyArrays.ReadOnlyArray{Symbol,...}
  And data, 3×2 Array{Float64,2}:
       (:a)  (:b)
  (1)   0.0   0.0
  (2)   0.0   0.0
  (3)   0.0   0.0
"""
Base.getindex(ds::KeyedDataset, key::Symbol) = getindex(ds, (key,))
function Base.getindex(ds::KeyedDataset, key::Tuple)
    A = ds.data[key]
    names = dimnames(A)
    keys = map(zip(names, axiskeys(A))) do (n, k)
        p = (key..., n)
        any(c -> p in c, ds.constraints) ? ReadOnlyArray(k) : k
    end

    kw = NamedTuple{dimnames(A)}(keys)
    return KeyedArray(parent(A); kw...)
end

"""
    setindex!(ds::KeyedDataset{T}, val, key) -> T

Store the new `val` in the [`KeyedDataset`](@ref). If any new dimension names don't any
existing constraints then `Pattern(:__, <dimname>)` is used by default.
If the axis values of the new `val` doesn't meet the existing constraints in the dataset
then an error will be throw.

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset, constraintmap;

julia> ds = KeyedDataset(:a => KeyedArray(zeros(3); time=1:3));

julia> ds[:b] = KeyedArray(ones(3, 2); time=1:3, lag=[-1, -2]);

julia> collect(constraintmap(ds))
2-element Vector{Pair{AxisSets.Pattern, Set{Tuple}}}:
 Pattern((:__, :time)) => Set([(:b, :time), (:a, :time)])
  Pattern((:__, :lag)) => Set([(:b, :lag)])

julia> ds[:c] = KeyedArray(ones(3, 2); time=2:4, lag=[-1, -2])
ERROR: KeyAlignmentError: Misaligned dimension keys on constraint Pattern((:__, :time))
  Tuple[(:b, :time), (:a, :time)] ∈ 3-element UnitRange{Int64}
  Tuple[(:c, :time)] ∈ 3-element UnitRange{Int64}
```
"""
Base.setindex!(ds::KeyedDataset, val, key::Symbol) = setindex!(ds, val, (key,))
function Base.setindex!(ds::KeyedDataset, val, key::Tuple)
    ds.data[key] = val

    for d in dimnames(val)
        dimpath = (key..., d)

        # Similar to construction, if our dimpath isn't present in any existing constraints
        # then we introduce a new one that's `Pattern(:__, dim)`
        if all(c -> !in(dimpath, c), ds.constraints)
            push!(ds.constraints, Pattern(:__, d))
        end
    end

    validate(ds)
end

"""
    (ds::KeyedDataset)(key) -> KeyedDataset

A collable syntax for selecting of filtering a subset of a [`KeyedDataset`](@ref).

# Example
```jldoctest
julia> using AxisKeys; using AxisSets: KeyedDataset, flatten;

julia> ds = KeyedDataset(
           flatten([
               :g1 => [
                   :a => KeyedArray(zeros(3); time=1:3),
                   :b => KeyedArray(ones(3, 2); time=1:3, loc=[:x, :y]),
                ],
                :g2 => [
                    :a => KeyedArray(ones(3); time=1:3),
                    :b => KeyedArray(zeros(3, 2); time=1:3, loc=[:x, :y]),
                ]
            ])...
       );

julia> collect(keys(ds(:__, :a).data))
2-element Vector{Tuple}:
 (:g1, :a)
 (:g2, :a)

julia> collect(keys(ds(:g1, :__).data))
2-element Vector{Tuple}:
 (:g1, :a)
 (:g1, :b)
```
"""
(ds::KeyedDataset)(args...) = _filterset(ds, args...)
_filterset(ds::KeyedDataset, key...) = _filterset(ds, key)
_filterset(ds::KeyedDataset, key::Symbol) = ds[key]
_filterset(ds::KeyedDataset, key::Tuple) = _filterset(ds, Pattern(key))
_filterset(ds::KeyedDataset, key::Pattern) = _filterset(ds, in(key))
function _filterset(ds::KeyedDataset, f::Function)
    data = filter(p -> f(first(p)), pairs(ds.data))
    paths = collect(Iterators.flatten(((k..., d) for d in dimnames(v)) for (k, v) in data))
    constraints = filter(c -> any(in(c), paths), ds.constraints)
    return KeyedDataset(constraints, data)
end


text:
using NetCDF, LinearAlgebra
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
#pygui(:qt)
using PyPlot
using MultivariateStats
# eofs = pyimport("eofs")

path = "/Users/milan/Downloads/cams/ieee"
filelist = filter(x->endswith(x,"no2.grib"),readdir(path))
gribfile = xr.open_dataset(joinpath(path,filelist[1]),engine="cfgrib")
no2 = gribfile.no2.values
level = 1:size(no2)[1]

no2m = mean(no2,dims=2)[:,1]
no2mi = median(no2,dims=2)[:,1]
no2min = minimum(no2,dims=2)[:,1]
no2max = maximum(no2,dims=2)[:,1]
no2p10 = [percentile(vec(no2[i,:]),10) for i in level]
no2p90 = [percentile(vec(no2[i,:]),90) for i in level]

## FIND FIRST SMALLER
function ispower2(x::Integer)
    while true
        if x == 2
            return true
        elseif x % 2 == 1
            return false
        else
            x = x ÷ 2
        end
    end
    return false
end

## FIND FIRST SMALLER
function findFirstSmaller(x::Float64,v::Array{Float64,1})

    l = length(v)-1
    @boundscheck ispower2(l) || throw(BoundsError())
    n = Int(log2(l))-1
    idx = l ÷ 2

    if x > v[end-1]     # round to max
        return l
    end

    # binary tree search
    for i in 1:n
        if x < v[idx]
            idx -= 2^(n-i)
        else
            idx += 2^(n-i)
        end
    end

    # split off the i = n+1 case
    if x >= v[idx]
        idx += 1
    end

    return idx-1
end

## LIN PACKING
Δlin24 = (Float64(maximum(no2max)) - Float64(minimum(no2min)))/(2^24-1)
Δlin16 = (Float64(maximum(no2max)) - Float64(minimum(no2min)))/(2^16-1)

lin24 = Array(minimum(no2min):Δlin24:maximum(no2max))
lin16 = Array(minimum(no2min):Δlin16:maximum(no2max))

lin24bounds = [(lin24[i+1]+lin24[i])/2 for (i,_) in enumerate(lin24[1:end-1])]
lin24bounds = vcat(minimum(no2min),lin24bounds,maximum(no2max))

lin16bounds = [(lin16[i+1]+lin16[i])/2 for (i,_) in enumerate(lin16[1:end-1])]
lin16bounds = vcat(minimum(no2min),lin16bounds,maximum(no2max))

no2lin24 = similar(no2)
no2lin16 = similar(no2)

for i in eachindex(no2)
    no2lin24[i] = Float32(lin24[findFirstSmaller(Float64(no2[i]),lin24bounds)])
    no2lin16[i] = Float32(lin16[findFirstSmaller(Float64(no2[i]),lin16bounds)])
end

## stats
no24m = mean(no2lin24,dims=2)[:,1]
no24mi = median(no2lin24,dims=2)[:,1]
no24min = minimum(no2lin24,dims=2)[:,1]
no24max = maximum(no2lin24,dims=2)[:,1]
no24p10 = [percentile(vec(no2lin24[i,:]),10) for i in level]
no24p90 = [percentile(vec(no2lin24[i,:]),90) for i in level]

no16m = mean(no2lin16,dims=2)[:,1]
no16mi = median(no2lin16,dims=2)[:,1]
no16min = minimum(no2lin16,dims=2)[:,1]
no16max = maximum(no2lin16,dims=2)[:,1]
no16p10 = [percentile(vec(no2lin16[i,:]),10) for i in level]
no16p90 = [percentile(vec(no2lin16[i,:]),90) for i in level]

## error 24bit
de24 = abs.(log10.(no2 ./ no2lin24))

de24m = mean(de24,dims=2)[:,1]
de24mi = median(de24,dims=2)[:,1]
de24min = minimum(de24,dims=2)[:,1]
de24max = maximum(de24,dims=2)[:,1]
de24p10 = [percentile(vec(de24[i,:]),10) for i in level]
de24p90 = [percentile(vec(de24[i,:]),90) for i in level]

de16 = abs.(log10.(no2 ./ no2lin16))

de16m = mean(de16,dims=2)[:,1]
de16mi = median(de16,dims=2)[:,1]
de16min = minimum(de16,dims=2)[:,1]
de16max = maximum(de16,dims=2)[:,1]
de16p10 = [percentile(vec(de16[i,:]),10) for i in level]
de16p90 = [percentile(vec(de16[i,:]),90) for i in level]

## plot
fig,(ax,ax2,ax3) = subplots(1,3,figsize=(10,6),sharey=true)
ax.invert_yaxis()

ax.semilogx(no24m,level,"C0",lw=3,label="mean")
ax.plot(no24mi,level,"C0",lw=1.5,label="median")
ax.fill_betweenx(level,no24min,no24max,color="C0",alpha=0.2,label="Min-max range")
ax.fill_betweenx(level,no24p10,no24p90,color="C0",alpha=0.5,label="10-90% range")

ax2.semilogx(no16m,level,"C0",lw=3,label="mean")
ax2.plot(no16mi,level,"C0",lw=1.5,label="median")
ax2.fill_betweenx(level,no16min,no16max,color="C0",alpha=0.2,label="Min-max range")
ax2.fill_betweenx(level,no16p10,no16p90,color="C0",alpha=0.5,label="10-90% range")

for i in level
    # 24-bit
    xvals = no2[i,(de24[i,:] .>= 1f-2) .* (de24[i,:] .< 1f-1)]
    ax.scatter(xvals,i*ones(length(xvals)),10,"yellow",alpha=0.05)

    xvals = no2[i,(de24[i,:] .>= 1f-1) .* (de24[i,:] .< 1f0)]
    ax.scatter(xvals,i*ones(length(xvals)),10,"C1",alpha=0.1)

    # 16-bit
    xvals = no2[i,(de16[i,:] .>= 2f-1) .* (de16[i,:] .< 0.5)]
    ax2.scatter(xvals,i*ones(length(xvals)),10,"yellow",alpha=0.05)

    xvals = no2[i,(de16[i,:] .>= 0.5) .* (de16[i,:] .< 1f0)]
    ax2.scatter(xvals,i*ones(length(xvals)),10,"C1",alpha=0.1)
end

ax.scatter([1,2,3],[0,0,0],10,"yellow",alpha=0.5,label="error > 1%")
ax.scatter([1,2,3],[0,0,0],10,"C1",alpha=0.5,label="error > 10%")

ax2.scatter([1,2,3],[0,0,0],10,"yellow",alpha=0.5,label="error > 20%")
ax2.scatter([1,2,3],[0,0,0],10,"C1",alpha=0.5,label="error > 50%")

ax3.semilogx(de24m,level,"C2",lw=3,label="24bit")
ax3.plot(de24mi,level,"C2",lw=1.5)
ax3.fill_betweenx(level,de24min,de24max,color="C2",alpha=0.2)
ax3.fill_betweenx(level,de24p10,de24p90,color="C2",alpha=0.5)

ax3.semilogx(de16m,level,"C1",lw=3,label="16bit")
ax3.plot(de16mi,level,"C1",lw=1.5)
ax3.fill_betweenx(level,de16min,de16max,color="C1",alpha=0.2)
ax3.fill_betweenx(level,de16p10,de16p90,color="C1",alpha=0.5)

ax.set_title(L"NO$_2$: 24-bit linear packing")
ax2.set_title(L"NO$_2$: 16-bit linear packing")
ax3.set_title("L1 log-error")
ax.set_ylabel("model level")
ax.set_xlabel("mixing ratio kg/kg")
ax2.set_xlabel("mixing ratio kg/kg")
ax.set_ylim(137,1)
ax.legend(loc=2,scatterpoints=3)
ax2.legend(loc=2,scatterpoints=3)
ax3.legend(loc=1)

ax.set_xlim(1e-14,1e-6)
ax2.set_xlim(1e-14,1e-6)

#ax2.set_xticks(lin16[1:20])

tight_layout()


text:
"""
    calculateSWI(data::Data, options::Options=Options())

Returns the calculated SWI using 'data' and 'options'.

# Examples
```julia-repl
julia> TEs = [4,8,12]
julia> data = Data(mag, phase, header(mag), TEs);
julia> swi = calculateSWI(data);
```
# With Options
$(@doc Options)

# Examples
```julia-repl
julia> TEs = [4,8,12]
julia> data = Data(mag, phase, header(mag), TEs);
julia> options = Options(phase_hp_σ=[10,10,5], mag_softplus=false)
julia> swi = calculateSWI(data, options);
```
"""
function calculateSWI(data::Data, options::Options=Options())
    if !isnothing(options.writesteps) mkpath(options.writesteps) end
    getswimag(data, options) .* getswiphase(data, options)
end

"""
    createMIP(S::AbstractArray{<:Number,3}, d=7)

Creates minimum intensity projection of `S` over `d` slices.

# Examples
```julia-repl
julia> TEs = [4,8,12]
julia> data = Data(mag, phase, header(mag), TEs);
julia> swi = calculateSWI(data);
julia> mip = createMIP(swi);
```

See also [`createIntensityProjection`](@ref)
"""
createMIP(S::AbstractArray{<:Number,3}, d=7) = createIntensityProjection(S, minimum, d)

"""
    createIntensityProjection(S::AbstractArray{T,3}, func, d=7) where T

Calculates an intensity projection over the 3rd dimension of `S` with the function `func` over `d` slices.
Good function are `minimum`, `maximum`, `mean`, `std`, ...

# Examples
```julia-repl
julia> using Statistics
julia> a = rand(Float64, (64,64,20))
julia> createIntensityProjection(a, std)
```
"""
function createIntensityProjection(S::AbstractArray{T,3}, func, d=7) where T
    [func(S[x,y,z:z+d-1]) for x in 1:size(S,1), y in 1:size(S,2), z in 1:size(S,3)-d+1]
end


text:
# using SymPy

# x, y = @vars x y
# u = x*(1-x)*y*(1-y)
# v = x*(1-x)*y*(1-y)^2

# p = x*(1-x)*y*(1-y)

# ux = diff(u,x)
# uy = diff(u,y)
# vx = diff(v,x)
# vy = diff(v,y)
# px = diff(p,x)
# py = diff(p,y)
# f = u*ux + v*uy - (diff(ux,x)+diff(uy,y)) + px
# g = u*vx + v*vy - (diff(vx,x)+diff(vy,y)) + py
# h = diff(u, x) + diff(v, y)
# println(replace(replace(sympy.julia_code(simplify(f)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(g)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(h)), ".^"=>"^"), ".*"=>"*"))

function u_exact(x,y)
    x*(1-x)*y*(1-y)
end

function v_exact(x,y)
    x*(1-x)*y*(1-y)
end

function p_exact(x,y)
    x*(1-x)*y*(1-y)
end

function t_exact(x,y)
    x*(1-x)*y*(1-y)
end

function k_exact(x,y)
    1 + x^2 + x / (1 + y^2)
end

function k_nn(x, y)
    out = fc([x y], [20,20,20,1])^2 + 0.5 # N x 1
    squeeze(out)
    # 1 + x^2 + x / (1 + y^2)
end

function ffunc_(x, y)
    x*y*(1 - x)*(1 - y)*(-x*y*(1 - x) + x*(1 - x)*(1 - y)) + x*y*(1 - x)*(1 - y)*(-x*y*(1 - y) + y*(1 - x)*(1 - y)) - x*y*(1 - y) + 0.02*x*(1 - x) + y*(1 - x)*(1 - y) + 0.02*y*(1 - y)
end

function gfunc_(x, y)
    x*y*(1 - x)*(1 - y)*(-x*y*(1 - x) + x*(1 - x)*(1 - y)) + x*y*(1 - x)*(1 - y)*(-x*y*(1 - y) + y*(1 - x)*(1 - y)) - x*y*(1 - x) + x*(1 - x)*(1 - y) + 0.02*x*(1 - x) + 0.02*y*(1 - y)
end

function hfunc_(x,y)
    -x*y*(1 - x) - x*y*(1 - y) + x*(1 - x)*(1 - y) + y*(1 - x)*(1 - y)
end

function heat_source_func(x, y)
    # x^2*y*(x - 1)^2*(y - 1)*(2*y - 1) + x*y^2*(x - 1)*(2*x - 1)*(y - 1)^2 - 2.0*x*(x - 1) - 2.0*y*(y - 1) # for k=1.0
    (2*x^2*y*(x - 1)*(2*y - 1) + x*y*(x - 1)*(y - 1)*(y^2 + 1)^2*(x*(x - 1)*(2*y - 1) + y*(2*x - 1)*(y - 1)) - (y^2 + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))/(y^2 + 1)^2
end

using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

m = 20
n = 20
h = 1/n
nu = 0.01

F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc_, m, n, h), m, n, h)
F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc_, m, n, h), m, n, h)
H = h^2*eval_f_on_fvm_pts(hfunc_, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
Laplace = nu * constant(compute_fem_laplace_matrix1(m, n, h))
heat_source = eval_f_on_gauss_pts(heat_source_func, m, n, h)
heat_source = constant(compute_fem_source_term1(heat_source, m, n, h))

xy = fem_nodes(m, n, h)
x, y = xy[:, 1], xy[:, 2]
k = @. k_nn(x, y); k=stack(k)
kgauss = fem_to_gauss_points(k, m, n, h)
LaplaceK = compute_fem_laplace_matrix1(kgauss, m, n, h)
# kgauss = eval_f_on_gauss_pts(k_exact, m, n, h)
# LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))

# The temperature term is
# T_xx + T_yy + q = u T_x + v T_y
#
function compute_residual(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    f5 = -F1
    F = f1 + f2 + f3 + f4 + f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    g5 = -F2
    G = g1 + g2 + g3 + g4 + g5

    H0 = -B * [u;v] + H

    # T_xx + T_yy = u T_x + v T_y - heat_source
    T0 = LaplaceK * T + compute_fem_advection_matrix1(constant(ugauss),constant(vgauss), m, n, h) * T - heat_source

    R = [F;G;H0;T0]
    return R
end

function compute_jacobian(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]

    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    # M = constant(compute_fem_laplace_matrix1(m, n, h)) +
        #  compute_fem_advection_matrix1(ugauss,vgauss, m, n, h)

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, m, n, h))

    gradT = eval_grad_on_gauss_pts1(T, m, n, h)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)

    J0 = [Fu Fv                         # Jacobian of [u; v]
          Gu Gv]

    J1 = [J0 -B'                        # Jacobian of [u; v; p]
          -B spdiag(zeros(size(B,1)))]

    N = 2*(m+1)*(n+1) + m*n                                     # length of [u; v; p]
    J = [J1 SparseTensor(spzeros(N,(m+1)*(n+1)))                # Jacobian of [u; v; p; T]
        DU_TX DV_TY SparseTensor(spzeros((m+1)*(n+1), m*n)) M]
end

NT = 5    # number of iterations for Newton's method

bd = bcnode("all", m, n, h)
bd = [bd; bd .+ (m+1)*(n+1);
    2*(m+1)*(n+1)+(n-1)*m+1:2*(m+1)*(n+1)+(n-1)*m+2;
    (2*(m+1)*(n+1)+m*n) .+ bd]
# only apply Dirichlet to velocity; set left bottom two points to zero to fix rank deficient problem for pressure

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    residual = scatter_update(residual, bd, zeros(length(bd)))    # residual[bd] .= 0.0 in Tensorflow syntax

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = @. u_exact(x,y)
v0 = @. v_exact(x,y)
t0 = @. t_exact(x,y)

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @. p_exact(x,y)

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(m*n+3*(m+1)*(n+1)))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n+(m+1)*(n+1)))

V_computed = S[end, :]

V_data = matread("CHTcoupled_data.mat")["V"]

sample_size = 40
idx = rand(1:(m+1)*(n+1), sample_size)
idx = [idx; (m+1)*(n+1) .+ idx; 2*(m+1)*(n+1)+m*n .+ idx]
observed_data = V_data[idx]

noise = true
noise_level = 0.05
if noise
    noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
    observed_data = observed_data .* noise_ratio
end

loss = mean((V_computed[idx] .- observed_data)^2)
loss = loss * 1e10
# ---------------------------------------------------
# create a session and run
max_iter = 100
sess = Session(); init(sess)
loss_ = BFGS!(sess, loss, max_iter)
figure(); semilogy(loss_); savefig("chtcoupled_loss.png")

xy = fem_nodes(m, n, h)
x, y = xy[:, 1], xy[:, 2]
k0 = @. k_exact(x,y)

figure(figsize=(14,4));
subplot(131)
visualize_scalar_on_fem_points(k0, m, n, h); title("conductivity exact")
subplot(132)
visualize_scalar_on_fem_points(run(sess, k), m, n, h); title("conductivity prediction")
subplot(133)
visualize_scalar_on_fem_points(k0.-run(sess, k), m, n, h); title("conductivity difference")
savefig("chtcoupled_k.png")


text:
function Annotate(x::Real, y::Real, pix_offset::ImVec2, fmt::String...)
    ccall((:ImPlot_Annotate_Str, libcimgui), Cvoid, (Cdouble,Cdouble,ImVec2,Cstring...),
          x, y, pix_offset, fmt...)
end

function AnnotateClamped(x::Real, y::Real, pix_offset::ImVec2, fmt::String...)
    ccall((:ImPlot_AnnotateClamped_Str, libcimgui), Cvoid, (Cdouble,Cdouble,ImVec2,Cstring...),
          x, y, pix_offset, fmt...)
end

function Annotate(x::Real, y::Real, pix_offset::ImVec2, color::ImVec4, fmt::String...)
    ccall((:ImPlot_Annotate_Vec4, libcimgui), Cvoid, (Cdouble,Cdouble,ImVec2,ImVec4,Cstring...),
          x, y, pix_offset, color, fmt...)
end

function AnnotateClamped(x::Real, y::Real, pix_offset::ImVec2, color::ImVec4, fmt::String...)
    ccall((:ImPlot_AnnotateClamped_Vec4, libcimgui), Cvoid, (Cdouble,Cdouble,ImVec2,ImVec4,Cstring...),
          x, y, pix_offset, color, fmt...)
end


text:
using Documenter
using Phylo

makedocs(modules = [Phylo, Phylo.API],
         sitename = "Phylo.jl")

deploydocs(repo = "github.com/EcoJulia/Phylo.jl.git",
           devbranch = "dev",
           deps = Deps.pip("pygments",
                           "mkdocs",
                           "mkdocs-material",
                           "python-markdown-math"))


text:
## Counstruction
struct Operator{T,D}
    type::Val{T}
    data::Dict{Symbol,D}
end
Operator(t::Symbol) = Operator(Val(t),Dict{Symbol,Any}())
function Operator(t::Symbol,d::Pair{Symbol,D}...) where D<:Any
    return Operator(Val(t),Dict(d))
end

## General Functions
push!(op::Operator,d::Pair{Symbol,D}...) where D<:Any = push!(op.data,d...)
@inline getproperty(op::Operator,f::Symbol) = hasfield(Operator,f) ? getfield(op,f) : getfield(op,:data)[f]
@inline function setproperty!(op::Operator,f::Symbol,x)
    getfield(op,:data)[f] = x
end

@inline function (op::Operator)(aps::Vector{T},gps::Vector{S},k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where {T<:AbstractElement,S<:AbstractElement}
    for i in 1:length(aps)
        @inbounds op(aps[i],gps[i],k,f)
    end
end

@inline function (op::Operator)(aps::Vector{T},k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    for ap in aps
        op(ap,k,f)
    end
end
@inline function (op::Operator)(aps::Vector{T},k::AbstractMatrix{Float64}) where T<:AbstractElement
    for ap in aps
        op(ap,k)
    end
end
@inline function (op::Operator)(aps::Vector{T},f::AbstractVector{Float64}) where T<:AbstractElement
    for ap in aps
        op(ap,f)
    end
end

@inline function (op::Operator)(aps::Vector{T},s::Symbol) where T<:AbstractElement
    for ap in aps
        op(ap,s)
    end
end
@inline function (op::Operator)(aps::Vector{T}) where T<:AbstractElement
    for ap in aps
        op(ap)
    end
end

function prescribe!(ap::T,s::Symbol,f::Function) where T<:AbstractElement
    𝓖 = ap.𝓖
    for ξ in 𝓖
        𝒙 = get𝒙(ap,ξ)
        if applicable(f,𝒙...)
            v = f(𝒙...)
        else
            𝒏 = get𝒏(ap)
            v = f(𝒙...,𝒏...)
        end
        setproperty!(ξ,s,v)
    end
end

function prescribe!(aps::Vector{T},s::Symbol,f::Function) where T<:AbstractElement
    for ap in aps
        prescribe!(ap,s,f)
    end
end

##
function (op::Operator{:𝑓𝑣})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            f[I] += N[i]*ξ.u*𝑤
        end
    end
end

## Potential Problem
function (op::Operator{:∫∇v∇uvbdΩ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N,B₁,B₂,B₃ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[I,J] += op.k*(B₁[i]*B₁[j] + B₂[i]*B₂[j] + B₃[i]*B₃[j])*𝑤
            end
            f[I] += N[i]*ξ.b*𝑤
        end
    end
end

function (op::Operator{:∫∇v∇udΩ})(ap::T,k::AbstractMatrix{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        _,B₁,B₂,B₃ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[I,J] += op.k*(B₁[i]*B₁[j] + B₂[i]*B₂[j] + B₃[i]*B₃[j])*𝑤
            end
        end
    end
end

function (op::Operator{:∫vbdΩ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            f[I] += N[i]*ξ.b*𝑤
        end
    end
end

function (op::Operator{:∫vtdΓ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            f[I] += N[i]*ξ.t*𝑤
        end
    end
end

function (op::Operator{:∫vgdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[I,J] += op.α*N[i]*N[j]*𝑤
            end
            f[I] += op.α*N[i]*ξ.g*𝑤
        end
    end
end

function (op::Operator{:∫λgdΓ})(ap1::T,ap2::S,g::AbstractMatrix{Float64},q::AbstractVector{Float64}) where {T<:AbstractElement,S<:AbstractElement}
    for ξ in ap1.𝓖
        𝑤 = get𝑤(ap1,ξ)
        N = get𝝭(ap1,ξ)
        N̄ = get𝝭(ap2,ξ)
        for k in 1:length(ap2.𝓒)
            K = ap2.𝓒[k].id
            for i in 1:length(ap1.𝓒)
                I = ap1.𝓒[i].id
                g[I,K] -= N[i]*N̄[k]*𝑤
            end
            q[K] -= N̄[k]*ξ.g*𝑤
        end
    end
end

function (op::Operator{:∫∇𝑛vgdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    for ξ in ap.𝓖
        N,B₁,B₂,B₃ = get∇𝝭(ap,ξ)
        n₁ = ξ.n₁
        n₂ = ξ.n₂
        n₃ = ξ.n₃
        for i in 1:length(ap.𝓒)
            I = ap.𝓒[i].id
            for j in 1:length(ap.𝓒)
                J = ap.𝓒[j].id
                k[I,J] -= op.k*((B₁[i]*n₁+B₂[i]*n₂+B₃[i]*n₃)*N[j] + N[i]*(B₁[j]*n₁+B₂[j]*n₂+B₃[j]*n₃))*ξ.w
            end
            f[I] -= op.k*(B₁[i]*n₁+B₂[i]*n₂+B₃[i]*n₃)*ξ.g*ξ.w
        end
    end
end

function (op::Operator{:∫∇̄𝑛vgdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N,B₁,B₂,B₃ = get∇̄𝝭(ap,ξ)
        n₁ = ξ.n₁
        n₂ = ξ.n₂
        n₃ = ξ.n₃
        g = ξ.g
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                k[I,J] += (B₁[i]*n₁+B₂[i]*n₂+B₃[i]*n₃)*N[j]*𝑤
            end
            f[I] += (B₁[i]*n₁+B₂[i]*n₂+B₃[i]*n₃)*ξ.g*𝑤
        end
    end
end

function (op::Operator{:g})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64};dof::Symbol=:d) where T<:AbstractElement{:Poi1}
    x = ap.𝓒[1]
    j = x.id
    g = getproperty(x,dof)
    for i in 1:length(f)
        f[i] -= k[i,j]*g
    end
    k[j,:] .= 0.
    k[:,j] .= 0.
    k[j,j] = 1.
    f[j] = g
end

## Plane Strain
function (op::Operator{:∫∫εᵢⱼσᵢⱼvᵢbᵢdxdy})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N,B₁,B₂ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵢᵢᵢᵢ = E/(1-ν^2)
        Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
        Cᵢⱼᵢⱼ = E/2/(1+ν)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] += (Cᵢᵢᵢᵢ*B₁[i]*B₁[j] + Cᵢⱼᵢⱼ*B₂[i]*B₂[j])*𝑤
                k[2*I-1,2*J]   += (Cᵢᵢⱼⱼ*B₁[i]*B₂[j] + Cᵢⱼᵢⱼ*B₂[i]*B₁[j])*𝑤
                k[2*I,2*J-1]   += (Cᵢᵢⱼⱼ*B₂[i]*B₁[j] + Cᵢⱼᵢⱼ*B₁[i]*B₂[j])*𝑤
                k[2*I,2*J]     += (Cᵢᵢᵢᵢ*B₂[i]*B₂[j] + Cᵢⱼᵢⱼ*B₁[i]*B₁[j])*𝑤
            end
            f[2*I-1] += N[i]*ξ.b₁*𝑤
            f[2*I]   += N[i]*ξ.b₂*𝑤
        end
    end
end

function (op::Operator{:∫∫εᵢⱼσᵢⱼdxdy})(ap::T,k::AbstractMatrix{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        _,B₁,B₂ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵢᵢᵢᵢ = E/(1-ν^2)
        Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
        Cᵢⱼᵢⱼ = E/2/(1+ν)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] += (Cᵢᵢᵢᵢ*B₁[i]*B₁[j] + Cᵢⱼᵢⱼ*B₂[i]*B₂[j])*𝑤
                k[2*I-1,2*J]   += (Cᵢᵢⱼⱼ*B₁[i]*B₂[j] + Cᵢⱼᵢⱼ*B₂[i]*B₁[j])*𝑤
                k[2*I,2*J-1]   += (Cᵢᵢⱼⱼ*B₂[i]*B₁[j] + Cᵢⱼᵢⱼ*B₁[i]*B₂[j])*𝑤
                k[2*I,2*J]     += (Cᵢᵢᵢᵢ*B₂[i]*B₂[j] + Cᵢⱼᵢⱼ*B₁[i]*B₁[j])*𝑤
            end
        end
    end
end

function (op::Operator{:∫∫εᵛᵢⱼσᵛᵢⱼdxdy})(ap::T,k::AbstractMatrix{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        _,B₁,B₂ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵛ = E/(1-2*ν)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] += Cᵛ/3*B₁[i]*B₁[j]*𝑤
                k[2*I-1,2*J]   += Cᵛ/3*B₁[i]*B₂[j]*𝑤
                k[2*I,2*J-1]   += Cᵛ/3*B₂[i]*B₁[j]*𝑤
                k[2*I,2*J]     += Cᵛ/3*B₂[i]*B₂[j]*𝑤
            end
        end
    end
end

function (op::Operator{:∫∫εᵈᵢⱼσᵈᵢⱼdxdy})(ap::T,k::AbstractMatrix{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        _,B₁,B₂ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵈ = E/(1+ν)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] += Cᵈ*( 2/3*B₁[i]*B₁[j]+1/2*B₂[i]*B₂[j])*𝑤
                k[2*I-1,2*J]   += Cᵈ*(-1/3*B₁[i]*B₂[j]+1/2*B₂[i]*B₁[j])*𝑤
                k[2*I,2*J-1]   += Cᵈ*(-1/3*B₂[i]*B₁[j]+1/2*B₁[i]*B₂[j])*𝑤
                k[2*I,2*J]     += Cᵈ*( 2/3*B₂[i]*B₂[j]+1/2*B₁[i]*B₁[j])*𝑤
            end
        end
    end
end

function (op::Operator{:∫∫vᵢbᵢdxdy})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N = get𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            f[2*I-1] += N[i]*ξ.b₁*𝑤
            f[2*I]   += N[i]*ξ.b₂*𝑤
        end
    end
end

function (op::Operator{:∫vᵢtᵢds})(ap::T,f::Vector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N = get𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            f[2*I-1] += N[i]*ξ.t₁*𝑤
            f[2*I]   += N[i]*ξ.t₂*𝑤
        end
    end
end

function (op::Operator{:∫λᵢgᵢds})(ap1::T,ap2::S,g::AbstractMatrix{Float64},q::AbstractVector{Float64}) where {T<:AbstractElement,S<:AbstractElement}
    for ξ in ap1.𝓖
        𝑤 = get𝑤(ap1,ξ)
        N = get𝝭(ap1,ξ)
        N̄ = get𝝭(ap2,ξ)
        g₁ = ξ.g₁
        g₂ = ξ.g₂
        n₁₁ = ξ.n₁₁
        n₂₂ = ξ.n₂₂
        n₁₂ = ξ.n₁₂
        for k in 1:length(ap2.𝓒)
            K = ap2.𝓒[k].id
            for i in 1:length(ap1.𝓒)
                I = ap1.𝓒[i].id
                N̄ₖNᵢ = N̄[k]*N[i]
                g[2*K-1,2*I-1] -= n₁₁*N̄ₖNᵢ*𝑤
                g[2*K-1,2*I]   -= n₁₂*N̄ₖNᵢ*𝑤
                g[2*K,2*I-1]   -= n₁₂*N̄ₖNᵢ*𝑤
                g[2*K,2*I]     -= n₂₂*N̄ₖNᵢ*𝑤
            end
            q[2*K-1] -= N̄[k]*(n₁₁*g₁+n₁₂*g₂)*𝑤
            q[2*K]   -= N̄[k]*(n₁₂*g₁+n₂₂*g₂)*𝑤
        end
    end
end

function (op::Operator{:∫σᵢⱼnⱼgᵢds})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N,B₁,B₂ = get∇𝝭(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵢᵢᵢᵢ = E/(1-ν^2)
        Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
        Cᵢⱼᵢⱼ = E/2/(1+ν)
        n₁₁ = ξ.n₁₁
        n₁₂ = ξ.n₁₂
        n₂₂ = ξ.n₂₂
        g₁ = ξ.g₁
        g₂ = ξ.g₂
        n₁ = ξ.n₁
        n₂ = ξ.n₂
        C₁₁₁ = Cᵢᵢᵢᵢ*n₁*n₁₁+Cᵢᵢⱼⱼ*n₂*n₁₂
        C₁₁₂ = Cᵢᵢᵢᵢ*n₁*n₁₂+Cᵢᵢⱼⱼ*n₂*n₂₂
        C₂₂₁ = Cᵢᵢⱼⱼ*n₁*n₁₁+Cᵢᵢᵢᵢ*n₂*n₁₂
        C₂₂₂ = Cᵢᵢⱼⱼ*n₁*n₁₂+Cᵢᵢᵢᵢ*n₂*n₂₂
        C₁₂₁ = Cᵢⱼᵢⱼ*(n₁*n₁₂+n₂*n₁₁)
        C₁₂₂ = Cᵢⱼᵢⱼ*(n₂*n₁₂+n₁*n₂₂)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] -= (C₁₁₁*(N[i]*B₁[j]+B₁[i]*N[j]) + C₁₂₁*(N[i]*B₂[j]+B₂[i]*N[j]))*ξ.w
                k[2*I-1,2*J]   -= (C₁₂₁*N[i]*B₁[j] + C₁₁₂*B₁[i]*N[j] + C₂₂₁*N[i]*B₂[j] + C₁₂₂*B₂[i]*N[j])*ξ.w
                k[2*I,2*J-1]   -= (C₁₁₂*N[i]*B₁[j] + C₁₂₁*B₁[i]*N[j] + C₁₂₂*N[i]*B₂[j] + C₂₂₁*B₂[i]*N[j])*ξ.w
                k[2*I,2*J]     -= (C₁₂₂*(N[i]*B₁[j]+B₁[i]*N[j]) + C₂₂₂*(N[i]*B₂[j]+B₂[i]*N[j]))*ξ.w
            end
            f[2*I-1] -= ((C₁₁₁*B₁[i]+C₁₂₁*B₂[i])*g₁ + (C₁₁₂*B₁[i]+C₁₂₂*B₂[i])*g₂)*ξ.w
            f[2*I]   -= ((C₁₂₁*B₁[i]+C₂₂₁*B₂[i])*g₁ + (C₁₂₂*B₁[i]+C₂₂₂*B₂[i])*g₂)*ξ.w
        end
    end
end

function (op::Operator{:∫σᵢⱼnⱼgᵢvᵢgᵢds})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N,B₁,B₂ = get∇𝝭(ap,ξ)
        E = op.E
        ν = op.ν
        α = op.α
        Cᵢᵢᵢᵢ = E/(1-ν^2)
        Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
        Cᵢⱼᵢⱼ = E/2/(1+ν)
        n₁₁ = ξ.n₁₁
        n₁₂ = ξ.n₁₂
        n₂₂ = ξ.n₂₂
        g₁ = ξ.g₁
        g₂ = ξ.g₂
        n₁ = ξ.n₁
        n₂ = ξ.n₂
        C₁₁₁ = Cᵢᵢᵢᵢ*n₁*n₁₁+Cᵢᵢⱼⱼ*n₂*n₁₂
        C₁₁₂ = Cᵢᵢᵢᵢ*n₁*n₁₂+Cᵢᵢⱼⱼ*n₂*n₂₂
        C₂₂₁ = Cᵢᵢⱼⱼ*n₁*n₁₁+Cᵢᵢᵢᵢ*n₂*n₁₂
        C₂₂₂ = Cᵢᵢⱼⱼ*n₁*n₁₂+Cᵢᵢᵢᵢ*n₂*n₂₂
        C₁₂₁ = Cᵢⱼᵢⱼ*(n₁*n₁₂+n₂*n₁₁)
        C₁₂₂ = Cᵢⱼᵢⱼ*(n₂*n₁₂+n₁*n₂₂)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] -= (C₁₁₁*(N[i]*B₁[j]+B₁[i]*N[j]) + C₁₂₁*(N[i]*B₂[j]+B₂[i]*N[j]) - α*N[i]*n₁₁*N[j])*ξ.w
                k[2*I-1,2*J]   -= (C₁₂₁*N[i]*B₁[j] + C₁₁₂*B₁[i]*N[j] + C₂₂₁*N[i]*B₂[j] + C₁₂₂*B₂[i]*N[j] -α*N[i]*n₁₂*N[j])*ξ.w
                k[2*I,2*J-1]   -= (C₁₁₂*N[i]*B₁[j] + C₁₂₁*B₁[i]*N[j] + C₁₂₂*N[i]*B₂[j] + C₂₂₁*B₂[i]*N[j] -α*N[i]*n₁₂*N[j])*ξ.w
                k[2*I,2*J]     -= (C₁₂₂*(N[i]*B₁[j]+B₁[i]*N[j]) + C₂₂₂*(N[i]*B₂[j]+B₂[i]*N[j]) -α*N[i]*n₂₂*N[j])*ξ.w
            end
            f[2*I-1] -= ((C₁₁₁*B₁[i]+C₁₂₁*B₂[i]-α*N[i]*n₁₁)*g₁ + (C₁₁₂*B₁[i]+C₁₂₂*B₂[i]-α*N[i]*n₁₂)*g₂)*ξ.w
            f[2*I]   -= ((C₁₂₁*B₁[i]+C₂₂₁*B₂[i]-α*N[i]*n₁₂)*g₁ + (C₁₂₂*B₁[i]+C₂₂₂*B₂[i]-α*N[i]*n₂₂)*g₂)*ξ.w
        end
    end
end

function (op::Operator{:∫σ̄ᵢⱼnⱼgᵢds})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N,B₁,B₂ = get∇̄𝝭(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵢᵢᵢᵢ = E/(1-ν^2)
        Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
        Cᵢⱼᵢⱼ = E/2/(1+ν)
        n₁₁ = ξ.n₁₁
        n₁₂ = ξ.n₁₂
        n₂₂ = ξ.n₂₂
        g₁ = ξ.g₁
        g₂ = ξ.g₂
        n₁ = ξ.n₁
        n₂ = ξ.n₂
        C₁₁₁ = Cᵢᵢᵢᵢ*n₁*n₁₁+Cᵢᵢⱼⱼ*n₂*n₁₂
        C₁₁₂ = Cᵢᵢᵢᵢ*n₁*n₁₂+Cᵢᵢⱼⱼ*n₂*n₂₂
        C₂₂₁ = Cᵢᵢⱼⱼ*n₁*n₁₁+Cᵢᵢᵢᵢ*n₂*n₁₂
        C₂₂₂ = Cᵢᵢⱼⱼ*n₁*n₁₂+Cᵢᵢᵢᵢ*n₂*n₂₂
        C₁₂₁ = Cᵢⱼᵢⱼ*(n₁*n₁₂+n₂*n₁₁)
        C₁₂₂ = Cᵢⱼᵢⱼ*(n₂*n₁₂+n₁*n₂₂)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] += (C₁₁₁*B₁[i]*N[j] + C₁₂₁*B₂[i]*N[j])*ξ.w
                k[2*I-1,2*J]   += (C₁₁₂*B₁[i]*N[j] + C₁₂₂*B₂[i]*N[j])*ξ.w
                k[2*I,2*J-1]   += (C₁₂₁*B₁[i]*N[j] + C₂₂₁*B₂[i]*N[j])*ξ.w
                k[2*I,2*J]     += (C₁₂₂*B₁[i]*N[j] + C₂₂₂*B₂[i]*N[j])*ξ.w
            end
            f[2*I-1] += ((C₁₁₁*B₁[i] + C₁₂₁*B₂[i])*g₁ + (C₁₁₂*B₁[i] + C₁₂₂*B₂[i])*g₂)*ξ.w
            f[2*I]   += ((C₁₂₁*B₁[i] + C₂₂₁*B₂[i])*g₁ + (C₁₂₂*B₁[i] + C₂₂₂*B₂[i])*g₂)*ξ.w
        end
    end
end

function (op::Operator{:∫σ̃̄ᵢⱼnⱼgᵢds})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        N,B₁,B₂ = get∇𝝭(ap,ξ)
        _,B̄₁,B̄₂ = get∇̄𝝭(ap,ξ)
        E = op.E
        ν = op.ν
        Cᵢᵢᵢᵢ = E/(1-ν^2)
        Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
        Cᵢⱼᵢⱼ = E/2/(1+ν)
        n₁₁ = ξ.n₁₁
        n₁₂ = ξ.n₁₂
        n₂₂ = ξ.n₂₂
        g₁ = ξ.g₁
        g₂ = ξ.g₂
        n₁ = ξ.n₁
        n₂ = ξ.n₂
        C₁₁₁ = Cᵢᵢᵢᵢ*n₁*n₁₁+Cᵢᵢⱼⱼ*n₂*n₁₂
        C₁₁₂ = Cᵢᵢᵢᵢ*n₁*n₁₂+Cᵢᵢⱼⱼ*n₂*n₂₂
        C₂₂₁ = Cᵢᵢⱼⱼ*n₁*n₁₁+Cᵢᵢᵢᵢ*n₂*n₁₂
        C₂₂₂ = Cᵢᵢⱼⱼ*n₁*n₁₂+Cᵢᵢᵢᵢ*n₂*n₂₂
        C₁₂₁ = Cᵢⱼᵢⱼ*(n₁*n₁₂+n₂*n₁₁)
        C₁₂₂ = Cᵢⱼᵢⱼ*(n₂*n₁₂+n₁*n₂₂)
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            ΔB₁ = B₁[i]-B̄₁[i]
            ΔB₂ = B₂[i]-B̄₂[i]
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] -= (C₁₁₁*(N[i]*B₁[j]+ΔB₁*N[j]) + C₁₂₁*(N[i]*B₂[j]+ΔB₂*N[j]))*ξ.w
                k[2*I-1,2*J]   -= (C₁₂₁*N[i]*B₁[j] + C₁₁₂*ΔB₁*N[j] + C₂₂₁*N[i]*B₂[j] + C₁₂₂*ΔB₂*N[j])*ξ.w
                k[2*I,2*J-1]   -= (C₁₁₂*N[i]*B₁[j] + C₁₂₁*ΔB₁*N[j] + C₁₂₂*N[i]*B₂[j] + C₂₂₁*ΔB₂*N[j])*ξ.w
                k[2*I,2*J]     -= (C₁₂₂*(N[i]*B₁[j]+ΔB₁*N[j]) + C₂₂₂*(N[i]*B₂[j]+ΔB₂*N[j]))*ξ.w
            end
            f[2*I-1] -= ((C₁₁₁*ΔB₁+C₁₂₁*ΔB₂)*g₁ + (C₁₁₂*ΔB₁+C₁₂₂*ΔB₂)*g₂)*ξ.w
            f[2*I]   -= ((C₁₂₁*ΔB₁+C₂₂₁*ΔB₂)*g₁ + (C₁₂₂*ΔB₁+C₂₂₂*ΔB₂)*g₂)*ξ.w
        end
    end
end

function (op::Operator{:∫vᵢgᵢds})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        n₁₁ = ξ.n₁₁
        n₂₂ = ξ.n₂₂
        n₁₂ = ξ.n₁₂
        g₁ = ξ.g₁
        g₂ = ξ.g₂
        for i in 1:length(𝓒)
            I = 𝓒[i].id
            for j in 1:length(𝓒)
                J = 𝓒[j].id
                k[2*I-1,2*J-1] += op.α*N[i]*n₁₁*N[j]*𝑤
                k[2*I,2*J-1]   += op.α*N[i]*n₁₂*N[j]*𝑤
                k[2*I-1,2*J]   += op.α*N[i]*n₁₂*N[j]*𝑤
                k[2*I,2*J]     += op.α*N[i]*n₂₂*N[j]*𝑤
            end
            f[2*I-1] += op.α*N[i]*(n₁₁*ξ.g₁+n₁₂*ξ.g₂)*𝑤
            f[2*I]   += op.α*N[i]*(n₁₂*ξ.g₁+n₂₂*ξ.g₂)*𝑤
        end
    end
end

## error estimates
function (op::Operator{:L₂})(ap::T) where T<:AbstractElement
    Δu²= 0
    ū² = 0
    for ξ in ap.𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        ūᵢ = ξ.u
        uᵢ = 0
        for i in 1:length(ap.𝓒)
            xᵢ = ap.𝓒[i]
            I = xᵢ.id
            uᵢ += N[i]*xᵢ.d
        end
        Δu² += (uᵢ - ūᵢ)^2*𝑤
        ū²  += ūᵢ^2*𝑤
    end
    return Δu², ū²
end

function (op::Operator{:L₂})(aps::Vector{T}) where T<:AbstractElement
    L₂Norm_Δu²= 0
    L₂Norm_ū² = 0
    for ap in aps
        Δu², ū² = op(ap)
        L₂Norm_Δu² += Δu²
        L₂Norm_ū²  += ū²
    end
    return (L₂Norm_Δu²/L₂Norm_ū²)^0.5
end

function (op::Operator{:H₁})(ap::T) where T<:AbstractElement
    Δ∇u²= 0
    ∇ū² = 0
    Δu²= 0
    ū² = 0
    for ξ in ap.𝓖
        𝑤 = get𝑤(ap,ξ)
        N,B₁,B₂,B₃ = get∇𝝭(ap,ξ)
        ūᵢ = ξ.u
        ∂ūᵢ∂x = ξ.∂u∂x
        ∂ūᵢ∂y = ξ.∂u∂y
        ∂ūᵢ∂z = ξ.∂u∂z
        uᵢ = 0.
        ∂uᵢ∂x = 0.
        ∂uᵢ∂y = 0.
        ∂uᵢ∂z = 0.
        for i in 1:length(ap.𝓒)
            xᵢ = ap.𝓒[i]
            I = xᵢ.id
            uᵢ += N[i]*xᵢ.d
            ∂uᵢ∂x += B₁[i]*xᵢ.d
            ∂uᵢ∂y += B₂[i]*xᵢ.d
            ∂uᵢ∂z += B₃[i]*xᵢ.d
        end
        Δ∇u² += ((∂uᵢ∂x - ∂ūᵢ∂x)^2 + (∂uᵢ∂y - ∂ūᵢ∂y)^2 + (∂uᵢ∂z - ∂ūᵢ∂z)^2)*𝑤
        ∇ū² += (∂ūᵢ∂x^2 + ∂ūᵢ∂y^2 + ∂ūᵢ∂z^2)*𝑤
        Δu² += (uᵢ - ūᵢ)^2*𝑤
        ū² += ūᵢ^2*𝑤
    end
    return Δ∇u², ∇ū², Δu², ū²
end

function (op::Operator{:H₁})(aps::Vector{T}) where T<:AbstractElement
    H₁Norm_Δu²= 0
    H₁Norm_ū² = 0
    L₂Norm_Δu²= 0
    L₂Norm_ū² = 0
    for ap in aps
        Δ∇u², ∇ū², Δu², ū² = op(ap)
        H₁Norm_Δu² += Δu² + Δ∇u²
        H₁Norm_ū²  += ū² + ∇ū²
        L₂Norm_Δu² += Δu²
        L₂Norm_ū²  += ū²
    end
    return (H₁Norm_Δu²/H₁Norm_ū²)^0.5, (L₂Norm_Δu²/L₂Norm_ū²)^0.5
end

function (op::Operator{:Hₑ_PlaneStress})(ap::T) where T<:AbstractElement
    ΔW²= 0
    W̄² = 0
    Δu²= 0
    ū² = 0
    E = op.E
    ν = op.ν
    Cᵢᵢᵢᵢ = E/(1-ν^2)
    Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
    Cᵢⱼᵢⱼ = E/2/(1+ν)
    for ξ in ap.𝓖
        N,B₁,B₂ = get∇𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        ū₁ = ξ.u
        ū₂ = ξ.v
        ∂ū₁∂x = ξ.∂u∂x
        ∂ū₁∂y = ξ.∂u∂y
        ∂ū₂∂x = ξ.∂v∂x
        ∂ū₂∂y = ξ.∂v∂y
        ε̄₁₁ = ∂ū₁∂x
        ε̄₂₂ = ∂ū₂∂y
        ε̄₁₂ = ∂ū₁∂y + ∂ū₂∂x
        σ̄₁₁ = Cᵢᵢᵢᵢ*ε̄₁₁ + Cᵢᵢⱼⱼ*ε̄₂₂
        σ̄₂₂ = Cᵢᵢᵢᵢ*ε̄₂₂ + Cᵢᵢⱼⱼ*ε̄₁₁
        σ̄₁₂ = Cᵢⱼᵢⱼ*ε̄₁₂
        u₁ = 0.
        u₂ = 0.
        ε₁₁ = 0.
        ε₂₂ = 0.
        ε₁₂ = 0.
        for i in 1:length(ap.𝓒)
            xᵢ = ap.𝓒[i]
            I = xᵢ.id
            u₁ += N[i]*xᵢ.d₁
            u₂ += N[i]*xᵢ.d₂
            ε₁₁ += B₁[i]*xᵢ.d₁
            ε₂₂ += B₂[i]*xᵢ.d₂
            ε₁₂ += B₂[i]*xᵢ.d₁ + B₁[i]*xᵢ.d₂
        end
        σ₁₁ = Cᵢᵢᵢᵢ*ε₁₁ + Cᵢᵢⱼⱼ*ε₂₂
        σ₂₂ = Cᵢᵢᵢᵢ*ε₂₂ + Cᵢᵢⱼⱼ*ε₁₁
        σ₁₂ = Cᵢⱼᵢⱼ*ε₁₂
        ΔW² += 0.5*((σ₁₁-σ̄₁₁)*(ε₁₁-ε̄₁₁) + (σ₂₂-σ̄₂₂)*(ε₂₂-ε̄₂₂) + (σ₁₂-σ̄₁₂)*(ε₁₂-ε̄₁₂))*𝑤
        W̄² += 0.5*(σ₁₁*ε₁₁ + σ₂₂*ε₂₂ + σ₁₂*ε₁₂)*𝑤
        Δu² += ((u₁ - ū₁)^2 + (u₂ - ū₂)^2)*𝑤
        ū² += (ū₁^2 + ū₂^2)*𝑤
    end
    return ΔW², W̄², Δu², ū²
end

function (op::Operator{:Hₑ_PlaneStress})(aps::Vector{T}) where T<:AbstractElement
    HₑNorm_ΔW²= 0.0
    HₑNorm_W̄² = 0.0
    L₂Norm_Δu²= 0.0
    L₂Norm_ū² = 0.0
    for ap in aps
        ΔW², W̄², Δu², ū² = op(ap)
        HₑNorm_ΔW² += ΔW²
        HₑNorm_W̄²  += W̄²
        L₂Norm_Δu² += Δu²
        L₂Norm_ū²  += ū²
    end
    return (HₑNorm_ΔW²/HₑNorm_W̄²)^0.5, (L₂Norm_Δu²/L₂Norm_ū²)^0.5
end

function set∇𝑢!(ap::T) where T<:AbstractElement
    for ξ in ap.𝓖
        N,B₁,B₂,B₃ = get∇𝝭(ap,ξ)
        𝒙 = get𝒙(ap,ξ)
        u = 0.
        ∂u∂x = 0.
        ∂u∂y = 0.
        ∂u∂z = 0.
        for i in 1:length(ap.𝓒)
            x = ap.𝓒[i]
            I = x.id
            u += N[i]*x.d
            ∂u∂x += B₁[i]*x.d
            ∂u∂y += B₂[i]*x.d
            ∂u∂z += B₃[i]*x.d
        end
        ξ.x = 𝒙[1]
        ξ.y = 𝒙[2]
        ξ.z = 𝒙[3]
        ξ.u = u
        ξ.∂u∂x = ∂u∂x
        ξ.∂u∂y = ∂u∂y
        ξ.∂u∂z = ∂u∂z
    end
end

function (op::Operator{:H₃})(aps::Vector{T}) where T<:AbstractElement
    H₃Norm_Δu²= 0
    H₃Norm_ū² = 0
    H₂Norm_Δu²= 0
    H₂Norm_ū² = 0
    H₁Norm_Δu²= 0
    H₁Norm_ū² = 0
    L₂Norm_Δu²= 0
    L₂Norm_ū² = 0
    for ap in aps
        Δ∇³u², ∇³ū²,Δ∇²u², ∇²ū²,Δ∇u², ∇ū², Δu², ū² = op(ap)
        H₃Norm_Δu² += Δu² + Δ∇u² + Δ∇²u² + Δ∇³u²
        H₃Norm_ū²  += ū² + ∇ū² + ∇²ū² + ∇³ū²
        H₂Norm_Δu² += Δu² + Δ∇u² + Δ∇²u²
        H₂Norm_ū²  += ū² + ∇ū² + ∇²ū²
        H₁Norm_Δu² += Δu² + Δ∇u²
        H₁Norm_ū²  += ū² + ∇ū²
        L₂Norm_Δu² += Δu²
        L₂Norm_ū²  += ū²
    end
    return (H₃Norm_Δu²/H₃Norm_ū²)^0.5, (H₂Norm_Δu²/H₂Norm_ū²)^0.5, (H₁Norm_Δu²/H₁Norm_ū²)^0.5, (L₂Norm_Δu²/L₂Norm_ū²)^0.5
end
function (op::Operator{:H₃})(ap::T) where T<:AbstractElement
    Δ∇³u²= 0
    ∇³ū² = 0
    Δ∇²u²= 0
    ∇²ū² = 0
    Δ∇u²= 0
    ∇ū² = 0
    Δu²= 0
    ū² = 0
    for ξ in ap.𝓖
        𝑤 = get𝑤(ap,ξ)
        N,B₁,B₂,B₁₁,B₁₂,B₂₂,B₁₁₁,B₁₁₂,B₁₂₂,B₂₂₂ = get∇³𝝭(ap,ξ)
        ūᵢ = ξ.u
        ∂ūᵢ∂x = ξ.∂u∂x
        ∂ūᵢ∂y = ξ.∂u∂y
        ∂²ūᵢ∂x² = ξ.∂²u∂x²
        ∂²ūᵢ∂x∂y = ξ.∂²u∂x∂y
        ∂²ūᵢ∂y² = ξ.∂²u∂y²
        ∂³ūᵢ∂x³ = ξ.∂³u∂x³
        ∂³ūᵢ∂x²∂y = ξ.∂³u∂x²∂y
        ∂³ūᵢ∂x∂y² = ξ.∂³u∂x∂y²
        ∂³ūᵢ∂y³ = ξ.∂³u∂y³
        uᵢ = 0.
        ∂uᵢ∂x = 0.
        ∂uᵢ∂y = 0.
        ∂²uᵢ∂x² = 0.
        ∂²uᵢ∂x∂y = 0.
        ∂²uᵢ∂y² = 0.
        ∂³uᵢ∂x³ = 0.
        ∂³uᵢ∂x²∂y = 0.
        ∂³uᵢ∂x∂y² = 0.
        ∂³uᵢ∂y³ = 0.
        for i in 1:length(ap.𝓒)
            xᵢ = ap.𝓒[i]
            I = xᵢ.id
            uᵢ += N[i]*xᵢ.d
            ∂uᵢ∂x += B₁[i]*xᵢ.d
            ∂uᵢ∂y += B₂[i]*xᵢ.d
            ∂²uᵢ∂x² += B₁₁[i]*xᵢ.d
            ∂²uᵢ∂x∂y += B₁₂[i]*xᵢ.d
            ∂²uᵢ∂y² += B₂₂[i]*xᵢ.d
            ∂³uᵢ∂x³ += B₁₁₁[i]*xᵢ.d
            ∂³uᵢ∂x²∂y += B₁₁₂[i]*xᵢ.d
            ∂³uᵢ∂x∂y² += B₁₂₂[i]*xᵢ.d
            ∂³uᵢ∂y³ += B₂₂₂[i]*xᵢ.d
        end
        Δ∇³u² += ((∂³uᵢ∂x³ - ∂³ūᵢ∂x³)^2 + (∂³uᵢ∂x²∂y - ∂³ūᵢ∂x²∂y)^2 + (∂³uᵢ∂x∂y² - ∂³ūᵢ∂x∂y²)^2 + (∂³uᵢ∂y³ - ∂³ūᵢ∂y³)^2)*𝑤
        ∇³ū² += (∂³ūᵢ∂x³^2 + ∂³ūᵢ∂x²∂y^2  + ∂³ūᵢ∂x∂y²^2+ ∂³ūᵢ∂y³^2)*𝑤
        Δ∇²u² += ((∂²uᵢ∂x² - ∂²ūᵢ∂x²)^2 + (∂²uᵢ∂x∂y - ∂²ūᵢ∂x∂y)^2 + (∂²uᵢ∂y² - ∂²ūᵢ∂y²)^2)*𝑤
        ∇²ū² += (∂²ūᵢ∂x²^2 + ∂²ūᵢ∂x∂y^2 + ∂²ūᵢ∂y²^2)*𝑤
        Δ∇u² += ((∂uᵢ∂x - ∂ūᵢ∂x)^2 + (∂uᵢ∂y - ∂ūᵢ∂y)^2)*𝑤
        ∇ū² += (∂ūᵢ∂x^2 + ∂ūᵢ∂y^2)*𝑤
        Δu² += (uᵢ - ūᵢ)^2*𝑤
        ū² += ūᵢ^2*𝑤
    end
    return Δ∇³u², ∇³ū², Δ∇²u², ∇²ū², Δ∇u², ∇ū², Δu², ū²
end

function set∇𝑢!(aps::Vector{T}) where T<:AbstractElement
    for ap in aps
        set∇𝑢!(ap)
    end
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
using BinaryBuilder, Pkg

name = "HPTT"
version = v"1.0.5"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/springer13/hptt.git", "a55c2a927d5462e81abeb12081fd345024caf5f6"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
atomic_patch -p1 $WORKSPACE/srcdir/patches/clang_compatibility.patch
mkdir -p ${libdir}
mkdir -p ${includedir}
export hpttdir=${WORKSPACE}/srcdir/hptt
cp ${hpttdir}/include/* ${includedir}
export CXXFLAGS="-O3 -std=c++11 -DNDEBUG -fopenmp -fPIC"
if [[ ${proc_family} == intel ]]; then
    export CXXFLAGS="$CXXFLAGS -mavx -DHPTT_ARCH_AVX"
elif [[ ${proc_family} == power ]]; then
    export CXXFLAGS="$CXXFLAGS -DHPTT_ARCH_IBM -maltivec -mabi=altivec";
## specific arm optimizations seem to be broken in library
# elif [[ ${target} == arm* ]]; then
#     export CXXFLAGS="$CXXFLAGS -mfpu=neon -DHPTT_ARCH_ARM"
fi
for f in ${hpttdir}/src/*.cpp; do
    $CXX $CXXFLAGS -I ${includedir} -c $f -o ${f%.cpp}.o
done
$CXX ${hpttdir}/src/*.o $CXXFLAGS -o ${libdir}/libhptt.$dlext -shared
install_license ${hpttdir}/LICENSE.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libhptt", :libhptt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"); platforms=filter(!Sys.isbsd, platforms))
    Dependency(PackageSpec(name="LLVMOpenMP_jll", uuid="1d63c593-3942-5779-bab2-d838dc0a180e"); platforms=filter(Sys.isbsd, platforms))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0", julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
# module TestSubsetting

using Test, Dates, Temporal

@testset "Subsetting" begin
    @testset "Identity" begin
        @test X[] == X
        @test X[:] == X
        @test X[:,:] == X
    end
    @testset "Single Row" begin
        @test X[1].values == data[1,:]'
        @test Y[1].values == data[1,:]'
        @test X[1,1].values[1,1] == data[1,1]
        @test Y[1,1].values[1,1] == data[1,1]
        @test X[1,1:2].values == data[1,1:2]'
        @test Y[1,1:2].values == data[1,1:2]'
        @test X[1,[1,2]].values == data[1,[1,2]]'
        @test Y[1,[1,2]].values == data[1,[1,2]]'
    end
    @testset "String Row Indexing" begin
        @testset "Basic" begin
            @test X["/"] == X[""] == X
            @test Y["/"] == Y[""] == Y
        end
        @testset "Single Observation" begin
            @testset "Date" begin
                for i in Temporal.DATE_STRING_LENGTHS
                    s = string(X.index[1])[1:i]
                    @test X[s] != 0
                end
            end
            @testset "Time" begin
                for i in Temporal.DATETIME_STRING_LENGTHS
                    s = string(Y.index[1])[1:i]
                    @test Y[s] != 0
                end
            end
        end
        @testset "Range" begin
            from_string = string(Y.index[1])
            thru_string = string(Y.index[end])
            @testset "Date" begin
                for i in Temporal.DATE_STRING_LENGTHS
                    s1 = from_string[1:i]
                    s2 = thru_string[1:i]
                    @test X["$s1/$s2"] != 0
                end
            end
            @testset "Time" begin
                for i in Temporal.DATETIME_STRING_LENGTHS
                    s1 = from_string[1:i]
                    s2 = thru_string[1:i]
                    @test Y["$s1/$s2"] != 0
                end
            end
        end
    end
    @testset "Single Date/DateTime Row Indexing" begin
        @test X[dates[1]].values == data[1,:]'
        @test Y[times[1]].values == data[1,:]'
        @test X[dates[1],1].values[:] == [data[1,1]]
        @test Y[times[1],1].values[:] == [data[1,1]]
        @test X[dates[1],1:2].values == data[1,1:2]'
        @test Y[times[1],1:2].values == data[1,1:2]'
        @test X[dates[1],[1,2]].values == data[1,1:2]'
        @test Y[times[1],[1,2]].values == data[1,1:2]'
    end
    @testset "Numeric Range Row Indexing" begin
        @test X[1:2].values == data[1:2,:]
        @test Y[1:2].values == data[1:2,:]
        @test X[1:2,1].values[:] == data[1:2,1]
        @test Y[1:2,1].values[:] == data[1:2,1]
        @test X[1:2,1:2].values == data[1:2,1:2]
        @test Y[1:2,1:2].values == data[1:2,1:2]
        @test X[1:2,[1,2]].values == data[1:2,1:2]
        @test Y[1:2,[1,2]].values == data[1:2,1:2]

    end
    @testset "Date/DateTime Range Row Indexing" begin
        @test X[dates[1:2]].values == data[1:2,:]
        @test Y[times[1:2]].values == data[1:2,:]
        @test X[dates[1:2],1].values[:] == data[1:2,1]
        @test Y[times[1:2],1].values[:] == data[1:2,1]
        @test X[dates[1:2],1:2].values == data[1:2,1:2]
        @test Y[times[1:2],1:2].values == data[1:2,1:2]
        @test X[dates[1:2],[1,2]].values == data[1:2,1:2]
        @test Y[times[1:2],[1,2]].values == data[1:2,1:2]
    end
    @testset "Numeric Array Row Indexing" begin
        @test X[[1,2]].values == data[1:2,:]
        @test Y[[1,2]].values == data[1:2,:]
        @test X[[1,2],1].values[:] == data[1:2,1]
        @test Y[[1,2],1].values[:] == data[1:2,1]
        @test X[[1,2],1:2].values == data[1:2,1:2]
        @test Y[[1,2],1:2].values == data[1:2,1:2]
        @test X[[1,2],[1,2]].values == data[1:2,1:2]
        @test Y[[1,2],[1,2]].values == data[1:2,1:2]
    end
    @testset "Date/DateTime Array Row Indexing" begin
        @test X[collect(dates[1:2])].values == data[1:2,:]
        @test Y[collect(times[1:2])].values == data[1:2,:]
        @test X[collect(dates[1:2]),1].values[:] == data[1:2,1]
        @test Y[collect(times[1:2]),1].values[:] == data[1:2,1]
        @test X[collect(dates[1:2]),1:2].values == data[1:2,1:2]
        @test Y[collect(times[1:2]),1:2].values == data[1:2,1:2]
        @test X[collect(dates[1:2]),[1,2]].values == data[1:2,1:2]
        @test Y[collect(times[1:2]),[1,2]].values == data[1:2,1:2]
    end
    @testset "Date/DateTime Step Range Row Indexing" begin
        @test X[dates[1:2]].values == data[1:2,:]
        @test Y[times[1:2]].values == data[1:2,:]
        @test X[dates[1:2],1].values[:] == data[1:2,1]
        @test Y[times[1:2],1].values[:] == data[1:2,1]
        @test X[dates[1:2],1:2].values == data[1:2,1:2]
        @test Y[times[1:2],1:2].values == data[1:2,1:2]
    end
    @testset "Symbol Column Indexing" begin
        @test X[:A] == X[:,1] == X[:,:A]
        @test X[1,:A] == X[1,1]
        @test X[[:A,:B]] == X[:,1:2] == X[:,[:A,:B]]
    end
end

# end


text:
#!/usr/bin/env julia

#Start Test Script
using ConjugateGradients
using Base.Test

function test_cg()
    tA = sprandn(100,100,.1) + 10.0*speye(100)
    A = tA'*tA
    b = rand(100)
    true_x = A\b
    x, exit_code, num_iters = cg((x,y) -> A_mul_B!(x, A, y) , b)
    if norm(true_x - x) < 1e-6
        return true
    else
        return false
    end
end

function test_bicgstab()
    A = sprandn(100,100,.1) + 10.0*speye(100)
    b = rand(100)
    true_x = A\b
    x, exit_code, num_iters = bicgstab((x,y) -> A_mul_B!(x, A, y) , b)
    if norm(true_x - x) < 1e-6
        return true
    else
        return false
    end
end

@test test_cg()
@test test_bicgstab()


text:
module Soss

import Base.rand
using Random
using Reexport: @reexport

@reexport using StatsFuns
using MeasureTheory

export For

using NamedTupleTools
using SampleChains
using SymbolicCodegen

using SymbolicUtils: Symbolic
const MaybeSym{T} = Union{T, Symbolic{T}}

# MeasureTheory.For(f, dims::MaybeSym{<: Integer}...) = ProductMeasure(mappedarray(i -> f(Tuple(i)...), CartesianIndices(dims)))
# MeasureTheory.For(f, dims::MaybeSym{<: AbstractUnitRange}...) = ProductMeasure(mappedarray(i -> f(Tuple(i)...), CartesianIndices(dims)))

import MacroTools: prewalk, postwalk, @q, striplines, replace, @capture
import MacroTools
import MLStyle
# import MonteCarloMeasurements
# using MonteCarloMeasurements: Particles, StaticParticles, AbstractParticles

using Requires
using ArrayInterface: StaticInt

using TransformVariables: as, asℝ, as𝕀, asℝ₊
import TransformVariables
const TV = TransformVariables

using SimplePosets: SimplePoset

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

const ∞ = MeasureTheory.∞

"""
we use this to avoid introduce static type parameters
for generated functions
"""
_unwrap_type(a::Type{<:Type}) = a.parameters[1]

include("noted.jl")
include("core/models/abstractmodel.jl")
include("core/statement.jl")
include("core/models/model.jl")
# include("core/models/jointdistribution.jl")
include("core/canonical.jl")
include("core/dependencies.jl")
include("core/toposort.jl")
include("core/weighted.jl")
include("core/utils.jl")
include("core/models/conditional.jl")

# include("distributions/dist.jl")
# include("distributions/for.jl")
include("distributions/iid.jl")
# include("distributions/mix.jl")
# include("distributions/flat.jl")
# include("distributions/markovchain.jl")

include("primitives/rand.jl")
include("primitives/simulate.jl")
include("primitives/logdensity.jl")
include("primitives/xform.jl")
include("primitives/likelihood-weighting.jl")
# include("primitives/gg.jl")
# @init @require Bijectors="76274a88-744f-5084-9051-94815aaf08c4" begin
#     include("primitives/bijectors.jl")
# end

include("primitives/basemeasure.jl")
include("primitives/testvalue.jl")
# include("primitives/entropy.jl")

include("transforms/predict.jl")
include("transforms/markovblanket.jl")
include("transforms/utils.jl")
include("transforms/basictransforms.jl")
include("transforms/withmeasures.jl")

include("symbolic/symcall.jl")
include("symbolic/symify.jl")
include("symbolic/rules.jl")
include("symbolic/symbolic.jl")
include("symbolic/codegen.jl")

# include("particles.jl")
include("plots.jl")

include("inference/rejection.jl")
# include("inference/dynamicHMC.jl")
# include("inference/advancedhmc.jl")
include("inference/power-posterior.jl")
# include("inference/Δlogdensity.jl")

#
# # include("graph.jl")
# # # include("optim.jl")
include("importance.jl")
#

function __init__()
    @require SampleChainsDynamicHMC = "6d9fd711-e8b2-4778-9c70-c1dfb499d4c4" begin
        include("samplechains/dynamichmc.jl")
    end
end

# # # include("sobols.jl")
# # # include("fromcube.jl")
# # # include("tocube.jl")
end # module


text:
const ENTRY_SYMBOLS = "HDRATJCIM*"

"""
    foreachentry(f, path)

Iterates over each entry in AAindex file at `path` and applies function `f` to it.
"""
function foreachentry(f, path)
    open(path, "r") do io
        while !eof(io)
            # read and parse indicies until end of file
            f(readuntil(io, "//\n"))
        end
    end
end

"""
    parse(record)

Parses given AAindex `entry`.
"""
function parse(record::String)::AbstractAAIndex
    pairs = Dict()
    lines = map(String, split(record, '\n', keepempty=false))

    while !isempty(lines) && first(lines) != "//"
        line = popfirst!(lines)
        tag, value = line[1], line[3:end]

        while !isempty(lines) && isspace(first(first(lines)))
            value *= " " * popfirst!(lines)
        end

        push!(pairs, tag => value)
    end

    if 'R' in keys(pairs)
        pairs['R'] = map(String, split(pairs['R'], " "))
    end

    metadata = Metadata(
        pairs['H'],
        pairs['D'],
        'R' in keys(pairs) ? pairs['R'] : String[],
        pairs['J'],
        pairs['T'],
        pairs['A'],
        '*' in keys(pairs) ? pairs['*'] : ""
    )

    if 'M' in keys(pairs)
        return AMatrix(_parse_matrix(pairs['M'])..., metadata)
    elseif 'I' in keys(pairs)
        return Index(
            _parse_index(pairs['I']),
            _parse_correlations(pairs['C']),
            metadata
        )
    end
end

function parse_id(record::String)::String
    only(match(r"(\w{4}\d{6})", record).captures)
end

function _parse_index(data::AbstractString)
    data = replace(data, "NA" => "NaN")
    data = split(data, r"\s+", keepempty=false)

    values = filter(x -> !isnothing(x), map(x -> tryparse(Float64, x), data))

    SVector{length(values)}(values)
end

function _parse_matrix(data::AbstractString)
    header_idx = findfirst(r"^[A-Za-z\s\-=,]+\s", data)
    header, data = data[header_idx], data[header_idx.stop+1:end]
    rowids, columnids = [header[idx] for idx in findall(r"[A-Z\-]+", header)]

    data = replace(data, "NA" => "NaN")
    data = replace(data, r"\s-\s" => "NaN")
    data = split(data, r"\s+", keepempty=false)

    values = filter(x -> !isnothing(x), map(x -> tryparse(Float64, x), data))

    m, n = length(rowids), length(columnids)
    data = zeros(m, n)

    # check if matrix is lower triangular
    if length(values) < m * n
        indices = [(x, y) for x in 1:m for y in 1:x]

        for (k, (i, j)) in enumerate(indices)
            data[i, j] = values[k]
        end

        data = SHermitianCompact{m}(data)
    else
        indices = [(x, y) for x in 1:m for y in 1:n]

        for (k, (i, j)) in enumerate(indices)
            data[i, j] = values[k]
        end

        data = SMatrix{m,n}(data)
    end

    return rowids, columnids, data
end

function _parse_correlations(data::AbstractString)
    tokens = split(data, r"\s", keepempty=false)
    values = []

    for i = 1:2:length(tokens)
        push!(values, String(tokens[i]) => tryparse(Float16, tokens[i+1]))
    end

    Dict(values)
end


text:
@testset "Test QFT module" begin
    @testset "Test forward QFT" begin
        num_qubits = 5
        num_gates = num_qubits*(num_qubits+1)/2
        cct1 = QXZoo.Circuit.Circ(num_qubits)
        cct2 = QXZoo.Circuit.Circ(num_qubits)
        QXZoo.QFT.apply_qft!(cct1, collect(1:num_qubits))

        cct2 << DefaultGates.h(1)
        cct2 << DefaultGates.c_r_phase(1, 2, π/2)
        cct2 << DefaultGates.c_r_phase(1, 3, π/4)
        cct2 << DefaultGates.c_r_phase(1, 4, π/8)
        cct2 << DefaultGates.c_r_phase(1, 5, π/16)
        cct2 << DefaultGates.h(2)
        cct2 << DefaultGates.c_r_phase(2, 3, π/2)
        cct2 << DefaultGates.c_r_phase(2, 4, π/4)
        cct2 << DefaultGates.c_r_phase(2, 5, π/8)
        cct2 << DefaultGates.h(3)
        cct2 << DefaultGates.c_r_phase(3, 4, π/2)
        cct2 << DefaultGates.c_r_phase(3, 5, π/4)
        cct2 << DefaultGates.h(4)
        cct2 << DefaultGates.c_r_phase(4, 5, π/2)
        cct2 << DefaultGates.h(5)
        cct2 << CompositeGates.swap(1,5)
        cct2 << CompositeGates.swap(2,4)

        num_swaps = Int(floor(num_qubits/2))
        @test cct1.circ_ops.len == num_gates + num_swaps

        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head
        h2 = cct2.circ_ops.head
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end
    @testset "Test inverse QFT" begin
        num_qubits = 5
        num_gates = num_qubits*(num_qubits+1)/2

        cct1 = QXZoo.Circuit.Circ(num_qubits)
        cct2 = QXZoo.Circuit.Circ(num_qubits)

        QXZoo.QFT.apply_iqft!(cct1, collect(1:num_qubits))

        cct2 << CompositeGates.swap(1,5)
        cct2 << CompositeGates.swap(2,4)

        cct2 << DefaultGates.h(5)
        cct2 << DefaultGates.c_r_phase(4, 5, -π/2)
        cct2 << DefaultGates.h(4)
        cct2 << DefaultGates.c_r_phase(3, 5, -π/4)
        cct2 << DefaultGates.c_r_phase(3, 4, -π/2)
        cct2 << DefaultGates.h(3)
        cct2 << DefaultGates.c_r_phase(2, 5, -π/8)
        cct2 << DefaultGates.c_r_phase(2, 4, -π/4)
        cct2 << DefaultGates.c_r_phase(2, 3, -π/2)
        cct2 << DefaultGates.h(2)
        cct2 << DefaultGates.c_r_phase(1, 5, -π/16)
        cct2 << DefaultGates.c_r_phase(1, 4, -π/8)
        cct2 << DefaultGates.c_r_phase(1, 3, -π/4)
        cct2 << DefaultGates.c_r_phase(1, 2, -π/2)
        cct2 << DefaultGates.h(1)

        num_swaps = Int(floor(num_qubits/2))
        @test cct1.circ_ops.len == num_gates + num_swaps

        @test cct1.gate_set == cct2.gate_set

        h1 = cct1.circ_ops.head
        h2 = cct2.circ_ops.head
        for i in 1:cct1.circ_ops.len
            @test h1.data === h2.data
            h1 = h1.next
            h2 = h2.next
        end
    end

end


text:
@generated function Base.convert(::Type{T}, x::Infinite) where {T<:Real}
  Infinite <: T && return :(x)
  pinf = posinf(T)
  ninf = neginf(T)
  mkpinf = pinf===nothing ? :(throw(InexactError(:convert,T,x))) : pinf
  mkninf = ninf===nothing ? :(throw(InexactError(:convert,T,x))) : ninf
  :(x.signbit ? $mkninf : $mkpinf)
end
Base.convert(::Type{Infinite}, x::Real) = isinf(x) ? Infinite(signbit(x)) : throw(InexactError(:convert,Infinite,x))
Base.convert(::Type{Infinite}, x::Infinite) = x

(::Type{T})(x::Infinite) where {T<:AbstractFloat} = convert(T, x)

"""
    UnknownReal()

Represents an unknown real number.
"""
struct UnknownReal <: Real end

Base.one(::Type{Infinite}) = UnknownReal
Base.zero(::Type{Infinite}) = UnknownReal
Base.float(::Type{UnknownReal}) = float(Int)
Base.float(::Type{Infinite}) = float(Int)


text:
#################### PDMatDistribution ####################

typealias PDMatDistribution Union{InverseWishart, Wishart}

function unlist(d::PDMatDistribution, X::AbstractArray)
  n = dim(d)
  y = similar(X, Int(n * (n + 1) / 2))
  k = 0
  for i in 1:n, j in i:n
    k += 1
    y[k] = X[i, j]
  end
  y
end

function relistlength(d::PDMatDistribution, X::AbstractArray)
  n = dim(d)
  Y = similar(X, n, n)
  k = 0
  for i in 1:n, j in i:n
    k += 1
    Y[i, j] = Y[j, i] = X[k]
  end
  (Y, k)
end

function link(d::PDMatDistribution, X::DenseMatrix)
  n = dim(d)
  Y = zeros(n, n)
  U = chol(X)
  for i in 1:n
    Y[i, i] = log(U[i, i])
  end
  for i in 1:n, j in (i + 1):n
    Y[i, j] = U[i, j]
  end
  Y
end

function invlink(d::PDMatDistribution, X::DenseMatrix)
  n = dim(d)
  U = zeros(n, n)
  for i in 1:n
    U[i, i] = exp(X[i, i])
  end
  for i in 1:n, j in (i + 1):n
    U[i, j] = X[i, j]
  end
  At_mul_B(U, U)
end

function logpdf(d::PDMatDistribution, X::DenseMatrix, transform::Bool)
  lp = logpdf(d, X)
  if transform && isfinite(lp)
    U = chol(X)
    n = dim(d)
    for i in 1:n
      lp += (n - i + 2) * log(U[i, i])
    end
    lp += n * log(2)
  end
  lp
end


text:
using Test
using ExprParsers
using Documenter

@test isempty(detect_ambiguities(ExprParsers))

# apparently the output of different kinds of Base.Expr change quite a lot between julia versions
if v"1.4" <= VERSION < v"1.5"
  @testset "doctest" begin
    doctest(ExprParsers)
  end
end

"""
  testing that `parser(parser(expr)) == parser(expr)`

which is the property of a closure operator, hence the name
"""
function test_closure(parser, expr)
  parsed_expr′ = parse_expr(parser, expr)
  expr′ = to_expr(parsed_expr′)
  parsed_expr′′ = parse_expr(parser, expr′)
  expr′′ = to_expr(parsed_expr′′)
  @test parsed_expr′ == parsed_expr′′
  @test expr′ == expr′′
end

@testset "syntax" begin
  include("syntax.jl")
end

@testset "expr_parsers_core" begin
  include("expr_parsers_core.jl")
end

@testset "expr_parsers_with_parsed" begin
  include("expr_parsers_with_parsed.jl")
end

@testset "expr_parsers_meta" begin
  include("expr_parsers_meta.jl")
end


text:
@testset "WindowedAssociativeOp" begin
    @testset "empty" begin
        for T in [Float64, Bool, Array{Float64,1}]
            for op in [+, -, *]
                state = WindowedAssociativeOp{T,op}()
                @test state.previous_cumsum == T[]
                @test state.ri_previous_cumsum == 0
                @test state.values == T[]
            end
        end
    end
end


text:
module SetOps

import Base: (+), (-), (*), (|), (&)
export powerset, ($)

"""
For sets `A` and `B`, `A-B` is the set of all elements of `A`
that are not in `B`.
"""
(-)(A::Set, B::Set)::Set = setdiff(A,B)

"""
For sets `A` and `B`, `A+B` is the union of `A` and `B`.
"""
(+)(A::Set, B::Set)::Set= union(A,B)

"""
For sets `A` and `B`, `A|B` is the union of `A` and `B`.
"""
(|)(A::Set, B::Set)::Set = union(A,B)

"""
For sets `A` and `B`, `A&B` is the intersection of `A` and `B`.
"""
(&)(A::Set, B::Set)::Set = intersect(A,B)

"""
For sets `A` and `B`, `A\$B` is the symmetric difference of `A` and `B`.
"""
($)(A::Set, B::Set)::Set = symdiff(A,B)

"""
For sets `A` and `B`, `A*B` [also `A×B`] is the Cartesian product, i.e., the
set of all ordered pairs `(a,b)` with `a` chosen from `A` and
`b` chosen from `B`.
"""
(*)(A::Set{S}, B::Set{T}) where {S,T}=Set{Tuple{S,T}}((a,b) for a in A for b in B)
(×)(A::Set,B::Set) = A*B

"""
`powerset(A)` is the set of all subsets of the set `A`.
"""
function powerset(A::Set{T})::Set{Set{T}} where T
    if length(A) == 0
        return Set([Set{T}()])
    end

    a = first(A)
    AA = deepcopy(A)
    delete!(AA,a)

    P1 = powerset(AA)
    P2 = Set( union(X,Set([a])) for X in P1 )

    return union(P1,P2)
end

end


text:
#=
# Float32 vs Float64

Compare training performance using `Float32` vs. `Float64`.
=#

import MKL
import MLDatasets
import Makie
import CairoMakie
import RestrictedBoltzmannMachines as RBMs
using RestrictedBoltzmannMachines: cd!
using ValueHistories: MVHistory

# Using Float32

train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float32}(train_x[:, :, train_y .== 0] .≥ 0.5)
rbm = RBMs.BinaryRBM(Float32, (28,28), 128)
RBMs.initialize!(rbm, train_x)
history32 = MVHistory()
time_0 = time()
@time RBMs.cd!(
    rbm, train_x; epochs=100, batchsize=128, steps=1,
    callback = function(@nospecialize(args...); @nospecialize(kw...))
        push!(history32, :t, time() - time_0)
    end
)
nothing #hide

# Using Float64

train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float64}(train_x[:, :, train_y .== 0] .≥ 0.5)
rbm = RBMs.BinaryRBM(Float64, (28,28), 128)
RBMs.initialize!(rbm, train_x)
history64 = MVHistory()
time_0 = time()
@time RBMs.cd!(
    rbm, train_x; epochs=100, batchsize=128, steps=1,
    callback = function(@nospecialize(args...); @nospecialize(kw...))
        push!(history64, :t, time() - time_0)
    end
)
nothing #hide

# Compare

fig = Makie.Figure(resolution=(600, 400))
ax = Makie.Axis(fig[1,1], xlabel="epoch", ylabel="seconds")
Makie.lines!(ax, get(history32, :t)..., label="32")
Makie.lines!(ax, get(history64, :t)..., label="64")
Makie.ylims!(ax, low=0)
Makie.axislegend(ax, position=:rt)
fig


text:
## this script contain functions related to the reference set
#   exchange strategy. It is a multi-point exchange strategy, that
#   encourages alternation in error signs.
# given:
#   - an array, a, of pairs of real numbers (x,y)
#   - a positive number h
#   - an integer N such that N+2 <= length of a.
# find the N+2 entries in that encourages:
#   - the y values alternate in sign, when ordered by x.
#   - the y values are among the largest.

function getcandidatesfromextremas( 𝑒::Function,
                                    L::Int,
                                    𝓧::Vector{T},
                                    sample_position_template::Vector{T},
                                    tol_𝓧_spacing, # was defaulted to convert(T,1e-6/length(𝓟_in)
                                    tol_derivative_zero,
                                    h) where T

    ### get roots.
    extremas, d_set, X_set = getextremafromrootsparallel( 𝑒, L, 𝓧,
                                            sample_position_template,
                                            tol_𝓧_spacing, tol_derivative_zero)
    #
    # if not enough extrema, there might be a numerical issue,
    #   or that the interpolant degree is not high enough, or that the extrema is
    #   at one of the boundary points in 𝓟 (i.e. in the reference set 𝓧).
    # add the reference set.

    # enforce the extremas positions must have an absolute error that is larger than h.
    N_samples = 5
    𝓟, 𝑒_𝓟 = refinecandidates!( extremas,
                                𝑒,
                                h,
                                N_samples)

    return 𝓟, 𝑒_𝓟, d_set, X_set
end

# Adds uniformly spaced positions between the closest candidate to -1. Similarly for 1.
# Removes candidates that have an error that is less than the level error h.
# These are incorrect results, since the local extremas
#   should be greater or equal to the level error, which is non-zero
#   unless the target function is in the space of Chebyshev approximations.
#
# ϵ is used for comparing two floating point numbers.
# Keep only the entries in a that are larger than the level error h.
# add new entries that are between the kept entires.
# modifies candidate_array.
# f here is the error function.
# The output is sorted in terms of the level error, not the candidate position.
function refinecandidates!( candidate_array::Vector{T},
                            f::Function,
                            h,
                            N_samples::Int,
                            ϵ::T = eps(T)*2) where T
    #
    unique!(candidate_array)

    # add endpoints and equally spaced points between each end point to the
    #   closest candidate, to the candidate list.
    if !isnumericallyclose(-one(T), candidate_array[1])
        tmp = collect( LinRange(-one(T), candidate_array[1], N_samples) )
        push!(candidate_array, tmp[1:end-1]...)
    end

    if !isnumericallyclose(one(T), candidate_array[end])

        tmp = collect( LinRange(candidate_array[end], one(T), N_samples) )
        push!(candidate_array, tmp[2:end]...)
    end

    # Evaluate the error for each candidate.
    error_array = SharedArray{Float64}(length(candidate_array))
    @sync @distributed for i = 1:length(candidate_array)
        error_array[i] = convert(Float64, f(candidate_array[i]))
    end

    ## discard the candidates that have an absolute error less than the level error.
    sort_ind = sortperm(abs.(error_array), rev=true) # descending order.
    candidate_𝓧_array = candidate_array[sort_ind]
    candidate_𝑒_array = error_array[sort_ind]

    k = 1
    while k <= length(candidate_𝑒_array) && abs(candidate_𝑒_array[k]) > h
        k += 1
    end
    k -= 1

    if k > 0
        # among the candidates, at least one is as large as the level error.
        # Keep the ones that are larger or equal to the level error.
        candidate_𝑒_array = convert(Vector{T}, candidate_𝑒_array[1:k])
        candidate_𝓧_array = candidate_𝓧_array[1:k]

        return candidate_𝓧_array, candidate_𝑒_array
    end

    # case: no candidates that have absolute error that is > the level error.
    return Vector{T}(undef,0), Vector{T}(undef,0)
end

# Use the candidates and the current reference set to update the reference set.
function processreferenceset!(  𝓧::Vector{T},
                                h,
                                𝑒::Function,
                                𝓟::Vector{T},
                                𝑒_𝓟::Vector{T},
                                L∞_value::Float64,
                                L∞_position::T,
                                sign_𝑒_at_L∞_position::Bool,
                                tol_𝓧_spacing,
                                #tol_h,
                                iter,
                                verbose_flag::Bool = true,
                                L∞_allowed_porpotion::Float64 = 0.95)::Bool where T
#
    @assert length(𝓟) == length(𝑒_𝓟)

    ## sanity-check: 𝑒(𝓧) must be alternating in sign.
    𝑒_𝓧_array = SharedArray{Float64}(length(𝓧))
    @sync @distributed for i = 1:length(𝓧)
        𝑒_𝓧_array[i] = convert(Float64, 𝑒(𝓧[i]) )
    end

    for i = 2:length(𝑒_𝓧_array)
        @assert !xor( sign(𝑒_𝓧_array[i-1]) > 1.0, sign(𝑒_𝓧_array[i]) > 1.0 )
    end

    ### Decide if doing multi-point exchange or one point exchange.
    # prune one entry so that the number of elements to be removed becomes even.

    if findfirst(abs.(𝑒_𝓟).>=L∞_allowed_porpotion*L∞_value) == nothing
        println("Using single exchange")
        # single exchange: replace an entry of the reference set with the
        #   L∞_position.
        # The entry should be chosen such that it is closest to L∞_position, but
        #   has the same error sign as the L∞_value

        min_dist_unused, j = findmin(abs.( 𝓧 .- L∞_position ))

        # find second minimum. Use the fact that 𝓧 is sorted in ascending order.
        j2 = j+1
        if abs( 𝓧[j-1] - L∞_position ) < abs( 𝓧[j+1] - L∞_position )
            j2 = j-1
        end

        sign_j = sign(𝑒_𝓧_array[j]) > 0
        sign_j2 = sign(𝑒_𝓧_array[j2]) > 0

        if sign_𝑒_at_L∞_position == sign_j
            𝓧[j] = L∞_position

            println("sign match as is")

        else

            # sanity check: the errors of 𝓧[j] and 𝓧[j2] should be alternating,
            #   since they are adjacent positions.
            @assert sign_𝑒_at_L∞_position == sign_j2

            𝓧[j2] = L∞_position

            println("sign match, use closest neighbour with same sign.")
        end

        sort!(𝓧)
        return true
    end

    ### Do multi-point exchange.
    println("Using multiple exchange")
    push!(𝓟, 𝓧...)

    ## discard entries that are too close in proximity.
    makenotclose!(𝓟, convert(T,1e-14))
    #unique!(𝓟)
    @assert length(𝓟) >= length(𝓧)

    if length(𝓟) > length(𝓧)
        𝓐 = removepointsAlg7!(𝓟, 𝑒, length(𝓧))
    end
    @assert length(𝓐) == length(𝓧)

    # check exchange condition.
    #@assert all(abs.(𝑒.(𝓐)).>=h)

    ### sanity-check: 𝑒(𝓧) must be alternating in sign.
    𝑒_𝓐_array = SharedArray{Float64}(length(𝓐))
    @sync @distributed for i = 1:length(𝓐)
        𝑒_𝓐_array[i] = convert(Float64, 𝑒(𝓐[i]) )
    end

    for i = 2:length(𝑒_𝓐_array)
        if xor( sign(𝑒_𝓐_array[i-1]) > 1.0, sign(𝑒_𝓐_array[i]) > 1.0 )

            # The errors of 𝓐 should be alternating.
            if verbose_flag
                Printf.@printf("Error (iter %d): 𝑒(𝓐) is not alternating! Exit.", iter )
            end

            return false
        end
    end

    ### sanity-check:
    if findfirst(abs.(𝑒_𝓐_array).>=L∞_allowed_porpotion*L∞_value) == nothing
        if verbose_flag
            Printf.@printf("Error (iter %d): none of the entries in 𝑒(𝓧_new) is near max(𝑒([-1,1])). Exit.", iter )
        end

        return false
    end

    ### safe to update reference set.
    𝓧[:] = 𝓐
    sort!(𝓧)

    return true
end

function processreferenceset(   𝓧::Vector{T},
                                h,
                                𝑒,
                                𝓟::Vector{T},
                                𝑒_𝓟::Vector{T},
                                L∞_value::Float64,
                                L∞_position::T,
                                sign_𝑒_at_L∞_position::Bool,
                                tol_𝓧_spacing,
                                #tol_h,
                                iter,
                                verbose_flag::Bool = true) where T
    𝓧_new = copy(𝓧)
    status_flag = processreferenceset!(𝓧_new, h, 𝑒, 𝓟, 𝑒_𝓟,
                                        L∞_value,
                                        L∞_position,
                                        sign_𝑒_at_L∞_position,
                                        tol_𝓧_spacing,
                                        #tol_h,
                                        iter, verbose_flag)

    return 𝓧_new, status_flag
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataFrames
using CSV

function plant_emissions_report(solution)::DataFrame
    df = DataFrame()
    df."plant type" = String[]
    df."location name" = String[]
    df."year" = Int[]
    df."emission type" = String[]
    df."emission amount (tonne)" = Float64[]
    T = length(solution["Energy"]["Plants (GJ)"])
    for (plant_name, plant_dict) in solution["Plants"]
        for (location_name, location_dict) in plant_dict
            for (emission_name, emission_amount) in location_dict["Emissions (tonne)"]
                for year = 1:T
                    push!(
                        df,
                        [
                            plant_name,
                            location_name,
                            year,
                            emission_name,
                            round(emission_amount[year], digits = 2),
                        ],
                    )
                end
            end
        end
    end
    return df
end

write_plant_emissions_report(solution, filename) =
    CSV.write(filename, plant_emissions_report(solution))

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG