text:
module MonteCarloMeasurements
using LinearAlgebra, Statistics, Random, StaticArrays, RecipesBase, GenericLinearAlgebra, MacroTools
using Distributed: pmap
import StatsBase: ProbabilityWeights
using Lazy: @forward
import Base: add_sum

using Distributions, StatsBase, Requires

const DEFAUL_NUM_PARTICLES = 500
const DEFAUL_STATIC_NUM_PARTICLES = 100

const COMPARISON_FUNCTION = Ref{Function}(mean)
const USE_UNSAFE_COMPARIONS = Ref(false)

"""
    unsafe_comparisons(onoff=true; verbose=true)
Toggle the use of a comparison function without warning. By default `mean` is used to reduce particles to a floating point number for comparisons. This function can be changed, example: `set_comparison_function(median)`
"""
function unsafe_comparisons(onoff=true; verbose=true)
    USE_UNSAFE_COMPARIONS[] = onoff
    if onoff && verbose
        @info "Unsafe comparisons using the function `$(COMPARISON_FUNCTION[])` has been enabled globally. Use `@unsafe` to enable in a local expression only or `unsafe_comparisons(false)` to turn off unsafe comparisons"
    end
end
"""
    set_comparison_function(f)

Change the Function used to reduce particles to a number for comparison operators
Toggle the use of a comparison Function without warning using the Function `unsafe_comparisons`.
"""
function set_comparison_function(f)
    COMPARISON_FUNCTION[] = f
end

"""
    @unsafe expression
Activates unsafe comparisons for the provided expression only. The expression is surrounded by a try/catch block to robustly restore unsafe comparisons in case of exception.
"""
macro unsafe(ex)
    ex2 = if @capture(ex, assigned_vars__ = y_)
        if length(assigned_vars) == 1
            esc(assigned_vars[1])
        else
            esc.(assigned_vars[1].args)
        end
    else
        :(res)
    end
    quote
        previous_state = USE_UNSAFE_COMPARIONS[]
        unsafe_comparisons(true, verbose=false)
        local res
        try
            res = ($(esc(ex)))
        finally
            unsafe_comparisons(previous_state, verbose=false)
        end
        $ex2 = res
    end
end

export ±, ∓, .., AbstractParticles,Particles,StaticParticles, WeightedParticles, MvParticles, sigmapoints, transform_moments, ≲,≳, systematic_sample, outer_product, meanstd, meanvar, register_primitive, register_primitive_multi, register_primitive_single, ℝⁿ2ℝⁿ_function, ℂ2ℂ_function, ℂ2ℂ_function!, resample!, bootstrap, sqrt!, exp!, sin!, cos!, wasserstein
# Plot exports
export errorbarplot, mcplot, ribbonplot

# Statistics reexport
export mean, std, cov, var, quantile, median
# Distributions reexport
export Normal, MvNormal, Cauchy, Beta, Exponential, Gamma, Laplace, Uniform, fit, logpdf

export unsafe_comparisons, @unsafe, set_comparison_function

export @bymap, @bypmap, Workspace, with_workspace, has_particles, mean_object

include("types.jl")
include("register_primitive.jl")
include("sampling.jl")
include("particles.jl")
include("distances.jl")
include("complex.jl")
include("sigmapoints.jl")
include("resampling.jl")
include("bymap.jl")
include("deconstruct.jl")
include("diff.jl")
include("plotting.jl")
include("optimize.jl")

# This is defined here so that @bymap is loaded
Base.log(p::Matrix{<:AbstractParticles}) = @bymap log(p) # Matrix more specific than StridedMatrix used in Base.log
LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles}) = @bymap eigvals(p)

function __init__()
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
end

end

# TODO: ifelse?


text:
domain = example_domain(10,10,0.1)
globdat = example_global_data(domain)

m = 10
n = 10
EBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
for j = 1:n+1
    EBC[(j-1)*(m+1) + m+1, :] .= -1
    g[(j-1)*(m+1) + m+1, :] = rand(2)
end

domain = Domain(domain.nodes, domain.elements, 2, EBC, g, zeros(Int64, (m+1)*(n+1), 2), zeros((m+1)*(n+1), 2))

domain.state = rand(length(domain.state))
fint0, stiff0 = assembleStiffAndForce( domain, Δt)
bd_contrib = fint0 - stiff0 * domain.state[domain.eq_to_dof]

domain.state[domain.eq_to_dof] = rand(domain.neqs)
fint0, stiff0 = assembleStiffAndForce( domain, Δt)
bd_contrib0 = fint0 - stiff0 * domain.state[domain.eq_to_dof]
@info "bd contrib", norm(bd_contrib-bd_contrib0)

ext: jl


text:
@defstruct ADAMOptions <: AbstractOptimizerOptions (
  (lr           :: Real = 0.001, lr > 0),
  (grad_clip    :: Real = 0, grad_clip >= 0),
  (weight_decay :: Real = 0.00001, weight_decay >= 0),
  (beta1        :: Real = 0.9,  beta1 > 0),
  (beta2        :: Real = 0.999,  beta2 > 0),
  (epsilon      :: Real = 1e-8, epsilon > 0),
  lr_scheduler  :: Any  = nothing
)

"""
     ADAM

The solver described in Diederik Kingma, Jimmy Ba: *Adam: A Method for
Stochastic Optimization*. arXiv:1412.6980 [cs.LG].

    ADAM(; kwargs...)

* `lr::Real`: default `0.001`, learning rate.
* `lr_scheduler::AbstractLearningRateScheduler`: default `nothing`, a
       dynamic learning rate scheduler. If set, will overwrite the `lr`
       parameter.
* `beta1::Real`: default `0.9`.
* `beta2::Real`: default `0.999`.
* `epsilon::Real`: default `1e-8`.
* `grad_clip::Real`: default `0`, if positive, will clip the gradient
       into the range `[-grad_clip, grad_clip]`.
* `weight_decay::Real`: default `0.00001`, weight decay is equivalent
       to adding a global l2 regularizer for all the parameters.
"""
type ADAM <: AbstractOptimizer
  opts  :: ADAMOptions
  state :: OptimizationState

  function ADAM(; kwargs...)
    opts = ADAMOptions(;kwargs...)
    opts.lr_scheduler = get_lr_scheduler(opts.lr_scheduler, opts.lr)

    new(opts)
  end
end

type ADAMState
  current_lr :: Float64  # current learning rate
  mt         :: NDArray
  vt         :: NDArray
  beta1Power :: Float64
  beta2Power :: Float64
end

function create_state(self :: ADAM, index :: Int, weight :: NDArray)
  return ADAMState( get_learning_rate(self.opts.lr_scheduler, self.state),
                    zeros(size(weight), context(weight)),
                    zeros(size(weight), context(weight)),
                    self.opts.beta1,
                    self.opts.beta2 )
end

function update(self :: ADAM, index :: Int, weight :: NDArray, grad :: NDArray, state :: ADAMState)
  lr = state.current_lr
  grad = normalized_gradient(self.opts, self.state, weight, grad)

  state.mt = self.opts.beta1 * state.mt + (1 - self.opts.beta1) * grad
  state.vt = self.opts.beta2 * state.vt + (1 - self.opts.beta2) * (grad .* grad)

  mt = state.mt / (1 - state.beta1Power)
  vt = state.vt / (1 - state.beta2Power)

  state.beta1Power *= self.opts.beta1
  state.beta2Power *= self.opts.beta2

  @inplace weight .+= -lr * mt ./ (sqrt(vt) + self.opts.epsilon)
end


text:
"""
Given a bond (i, j) it has an environment {k} which could be
cylindrical, ellipsoidal, ...
The ball B(ri, rcut_env) must be a superset of this environment.
"""
function cutoff_env end

abstract type BondEnvelope{T} <: B1pMultiplier{T}  end

# we can't call this `evaluate` since that would return a Vector
# whereas _inner_evaluate(...) should just give us a value.
function _inner_evaluate(env::BondEnvelope, X::AbstractState)
   if X.be == :env
      return _evaluate_env(env, X)
   elseif X.be == :bond
      return _evaluate_bond(env, X)
   else
      error("invalid X.be value")
   end
end

Base.filter(env::BondEnvelope, X::AbstractState) =
         _inner_evaluate(env,X) != 0

struct CylindricalBondEnvelope{T} <: BondEnvelope{T}
   r0cut::T   # bond-length cutoff
   rcut::T    # env cutoff into r-direction
   zcut::T    #
   p0::Int
   pr::Int
   pz::Int
   floppy::Bool
   λ::T
end

CylindricalBondEnvelope(r0cut, rcut, zcut; p0 = 2, pr = 2, pz = 2, floppy = true, λ = .5) =
      CylindricalBondEnvelope(r0cut, rcut, zcut, p0, pr, pz, floppy, λ)

cutoff_env(env::CylindricalBondEnvelope) = sqrt(env.rcut^2 + (env.r0cut + env.zcut)^2)

struct EllipsoidBondEnvelope{T} <: BondEnvelope{T}
   r0cut::T
   rcut::T
   zcut::T  # must satisfy zcut >= r0cut/2
   p0::Int
   pr::Int
   floppy::Bool
   λ::T
end

EllipsoidBondEnvelope(r0cut, zcut, rcut; p0=2, pr=2, floppy=false, λ=.5) = EllipsoidBondEnvelope(r0cut, zcut, rcut, p0, pr, floppy, λ)
EllipsoidBondEnvelope(r0cut, cut; p0=2, pr=2, floppy=false, λ=.5) = EllipsoidBondEnvelope(r0cut, cut, cut, p0, pr, floppy, λ)

cutoff_env(env::EllipsoidBondEnvelope) = env.zcut + env.rcut

function _evaluate_bond(env::BondEnvelope, X::AbstractState)
   r = norm(X.rr0)
   return ((r/env.r0cut)^2 - 1)^(env.p0) * (r <= env.r0cut)
end

function _evaluate_env(env::BondEnvelope, X::AbstractState)
   # convert to cylindrical coordinates
   z, r = cat2cyl(env, X)
   zeff = _eff_zcut(env, X)
   # then return the correct cutoff
   return _eval_env_inner(env, z, r, zeff)
end

"""
`cat2cyl`: Convert cartesian coordinate to cylindrical coordinate
Given `env` that contains information of the centre and `X` the state
return Cylindrical coordinate of X as `z, r = cat2cyl(env, X)`
"""
function cat2cyl(env::BondEnvelope, X::AbstractState)
   x_centre = env.λ * X.rr0
   x̃ = X.rr - x_centre
   if norm(X.rr0) >0
      r̂b = X.rr0/norm(X.rr0)
      z̃ = dot(x̃, r̂b)
      r̃ = norm( x̃ - z̃ * r̂b )
   else
      z̃ = 0.0
      r̃ = norm(x̃)
   end
   return z̃, r̃
end

# Compute effective zcut when `floppy = true` (to enable the cutoff in z direction to grow with the bond)
_eff_zcut(env::BondEnvelope, X::AbstractState) = env.zcut + env.floppy * norm(env.λ * X.rr0)

# The envelope for environment bonds
_eval_env_inner(env::CylindricalBondEnvelope, z, r, zeff) = ( (z/zeff)^2 - 1 )^(env.pz) * (abs(z) <= zeff) *
         ( (r/env.rcut)^2 - 1 )^(env.pr) * (r <= env.rcut)

_eval_env_inner(env::EllipsoidBondEnvelope, z, r, zeff) = ( (z/zeff)^2 +(r/env.rcut)^2 - 1.0)^(env.pr) * ((z/zeff)^2 +(r/env.rcut)^2 <= 1)


text:
module GreenDeltaStar

using Elliptic

export Green

include("green.jl")

end


text:
using BAT, ValueShapes

#Model definition to generate samples from a n-dim gaussian shell
struct GaussianShellDensity<:AbstractDensity
    lambda::Vector{Float64}
    r::Float64
    sigma::Float64
    dimensions::Int64
end
ValueShapes.totalndof(model::GaussianShellDensity) = model.dimensions

#define likelihood for the Gaussian Shell
function DensityInterface.logdensityof(target::GaussianShellDensity, v::AbstractArray{Float64, 1})
    diff::Float64 = 0
    for i in eachindex(v)
        diff += (target.lambda[i] - v[i]) * (target.lambda[i] - v[i])
    end
    diff = sqrt(diff)
    expo::Float64 = exp(-(diff - target.r) * (diff - target.r) / (2 * target.sigma^2))
    return log(1.0 / sqrt(2 * pi * target.sigma^2) * expo)
end

algorithm = MCMCSampling(mcalg = MetropolisHastings())
#algorithm = MCMCSampling(mcalg = MetropolisHastings(ARPWeighting{Float64}()))

#define model and #dimensions
dim = 2
model = GaussianShellDensity(zeros(dim), 5.0, 2.0, dim)

#define boundaries
lo_bounds = [-30.0 for i = 1:dim]
hi_bounds = [ 30.0 for i = 1:dim]
bounds = BAT.HyperRectBounds(lo_bounds, hi_bounds)

#Harmonic Mean Integration
#True integral value for 2D Gaussian Shell I = 31.4411
#True integral value for 10D Gaussian Shell I = 1.1065e9

#BAT.jl samples
bat_samples = bat_sample(PosteriorDensity(model, bounds), algorithm).result
data = BAT.HMIData(bat_samples)
BAT.hm_integrate!(data)

using Plots; pyplot()
plot(data, rscale = 0.25)


text:
module SignatureV4

using AWSCore
using Base64
using Dates
using HTTP
using HTTP.Pairs
using HTTP.URIs
using HTTP: Headers, Layer, Request
using IniFile
using MbedTLS

export AWS4AuthLayer, sign_signature!

"""
    AWS4AuthLayer{Next} <: HTTP.Layer

Abstract type used by [`HTTP.request`](@ref) to add an
[AWS Signature v4](http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html)
authentication layer to the request.

Historically this layer has been placed in-between the `MessageLayer` and `RetryLayer` in the HTTP
stack. The request stack can be found
[here](https://github.com/JuliaWeb/HTTP.jl/blob/v0.8.6/src/HTTP.jl#L534).

An example of how to layer insert can be found below:
```julia
using AWSCore
using HTTP

stack = HTTP.stack()

# This will insert the AWS4AuthLayer before the RetryLayer
insert(stack, RetryLayer, AWS4AuthLayer)
```
"""
abstract type AWS4AuthLayer{Next<:Layer} <: Layer{Next} end

"""
    HTTP.request(::Type{AWS4AuthLayer}, url::HTTP.URI, req::HTTP.Request, body) -> HTTP.Response

Perform the given request, adding a layer of AWS authentication using
[AWS Signature v4](http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html).
An "Authorization" header to the request.
"""
function HTTP.request(::Type{AWS4AuthLayer{Next}}, url::URI, req::Request, body; kw...) where Next
    if !haskey(kw, :aws_access_key_id)
        creds = AWSCredentials()
        creds = (aws_access_key_id=creds.access_key_id, aws_secret_access_key=creds.secret_key)
        kw = merge(creds, kw)
    end
    sign_signature!(req.method, url, req.headers, req.body; kw...)
    return HTTP.request(Next, url, req, body; kw...)
end

# Normalize whitespace to the form required in the canonical headers.
# Note that the expected format for multiline headers seems not to be explicitly
# documented, but Amazon provides a test case for it, so we'll match that behavior.
# We replace each `\n` with a `,` and remove all whitespace around the newlines,
# then any remaining contiguous whitespace is replaced with a single space.
function _normalize_ws(s::AbstractString)
    if any(isequal('\n'), s)
        return join(map(_normalize_ws, split(s, '\n')), ',')
    end

    return replace(strip(s), r"\s+" => " ")
end

"""
    sign_signature!(method::String, url::HTTP.URI, headers::HTTP.Headers, body; kwargs...)

Add an "Authorization" header to `headers`, modifying it in place.
The header contains a computed signature based on the given credentials as well as
metadata about what was used to compute the signature.
For more information, see the AWS documentation on the
[Signature v4](http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html)
process.

# Keyword arguments

All keyword arguments to this function are optional, as they have default values.

* `body_sha256`: A precomputed SHA-256 sum of `body`
* `body_md5`: A precomputed MD5 sum of `body`
* `timestamp`: The timestamp used in request signing (defaults to now in UTC)
* `aws_service`: The AWS service for the request (determined from the URL)
* `aws_region`: The AWS region for the request (determined from the URL)
* `aws_access_key_id`: AWS access key (read from the environment)
* `aws_secret_access_key`: AWS secret access key (read from the environment)
* `aws_session_token`: AWS session token (read from the environment, or empty)
* `token_in_signature`: Use `aws_session_token` when computing the signature (`true`)
* `include_md5`: Add the "Content-MD5" header to `headers` (`true`)
* `include_sha256`: Add the "x-amz-content-sha256" header to `headers` (`true`)
"""
function sign_signature!(
    method::String,
    url::URI,
    headers::Headers,
    body::Vector{UInt8};
    body_sha256::Vector{UInt8}=digest(MD_SHA256, body),
    body_md5::Vector{UInt8}=digest(MD_MD5, body),
    t::Union{DateTime,Nothing}=nothing,
    timestamp::DateTime=now(Dates.UTC),
    aws_service::String=String(split(url.host, ".")[1]),
    aws_region::String=String(split(url.host, ".")[2]),
    aws_access_key_id::String=ENV["AWS_ACCESS_KEY_ID"],
    aws_secret_access_key::String=ENV["AWS_SECRET_ACCESS_KEY"],
    aws_session_token::String=get(ENV, "AWS_SESSION_TOKEN", ""),
    token_in_signature=true,
    include_md5=true,
    include_sha256=true,
    kw...,
)
    # ISO8601 date/time strings for time of request...
    date = Dates.format(timestamp, dateformat"yyyymmdd")
    datetime = Dates.format(timestamp, dateformat"yyyymmdd\THHMMSS\Z")

    # Authentication scope...
    scope = [date, aws_region, aws_service, "aws4_request"]

    # Signing key generated from today's scope string...
    signing_key = string("AWS4", aws_secret_access_key)
    for element in scope
        signing_key = digest(MD_SHA256, element, signing_key)
    end

    # Authentication scope string...
    scope = join(scope, "/")

    # SHA256 hash of content...
    content_hash = bytes2hex(body_sha256)

    # HTTP headers...
    rmkv(headers, "Authorization")
    setkv(headers, "host", url.host)
    setkv(headers, "x-amz-date", datetime)
    include_md5 && setkv(headers, "Content-MD5", base64encode(body_md5))
    if (aws_service == "s3" && method == "PUT") || include_sha256
        # This header is required for S3 PUT requests. See the documentation at
        # https://docs.aws.amazon.com/AmazonS3/latest/API/sig-v4-header-based-auth.html
        setkv(headers, "x-amz-content-sha256", content_hash)
    end
    if aws_session_token != ""
        setkv(headers, "x-amz-security-token", aws_session_token)
    end

    # Sort and lowercase() Headers to produce canonical form...
    unique_header_keys = Vector{String}()
    normalized_headers = Dict{String,Vector{String}}()
    for (k, v) in sort!([lowercase(k) => v for (k, v) in headers], by=first)
        # Some services want the token included as part of the signature
        if k == "x-amz-security-token" && !token_in_signature
            continue
        end
        # In Amazon's examples, they exclude Content-Length from signing. This does not
        # appear to be addressed in the documentation, so we'll just mimic the example.
        if k == "content-length"
            continue
        end
        if !haskey(normalized_headers, k)
            normalized_headers[k] = Vector{String}()
            push!(unique_header_keys, k)
        end
        push!(normalized_headers[k], _normalize_ws(v))
    end
    canonical_headers = map(unique_header_keys) do k
        string(k, ':', join(normalized_headers[k], ','))
    end
    signed_headers = join(unique_header_keys, ';')

    # Sort Query String...
    query = sort!(collect(queryparams(url.query)), by=first)

    # Paths for requests to S3 should be escaped but not normalized. See
    # http://docs.aws.amazon.com/AmazonS3/latest/API/sig-v4-header-based-auth.html#canonical-request
    # Note that escapepath escapes ~ per RFC 1738, but Amazon includes an example in their
    # signature v4 test suite where ~ remains unescaped. We follow the spec here and thus
    # deviate from Amazon's example in this case.
    path = escapepath(aws_service == "s3" ? url.path : URIs.normpath(url.path))

    # Create hash of canonical request...
    canonical_form = join([method,
                           path,
                           escapeuri(query),
                           join(canonical_headers, "\n"),
                           "",
                           signed_headers,
                           content_hash], "\n")
    canonical_hash = bytes2hex(digest(MD_SHA256, canonical_form))

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$scope\n$canonical_hash"
    signature = bytes2hex(digest(MD_SHA256, string_to_sign, signing_key))

    # Append Authorization header...
    setkv(headers, "Authorization", string(
        "AWS4-HMAC-SHA256 ",
        "Credential=$aws_access_key_id/$scope, ",
        "SignedHeaders=$signed_headers, ",
        "Signature=$signature"
    ))

    return headers
end
end # module AWS4AuthRequest


text:
############################################################################################
"""
$(TYPEDEF)

Contains information about log-likelihood, expected sample size and proposal trajectory.

# Fields
$(TYPEDFIELDS)
"""
struct PMCMCDiagnostics{T,P<:ParticleFilterDiagnostics,M<:MCMCDiagnostics} <: AbstractDiagnostics
    "Diagnostics used for all Baytes kernels"
    base::BaytesCore.BaseDiagnostics{T}
    "Particle Filter diagnostics."
    pf::P
    "MCMC diagnostics."
    mcmc::M
    function PMCMCDiagnostics(
        base::BaytesCore.BaseDiagnostics{T}, pf::P, mcmc::M
    ) where {T,P<:ParticleFilterDiagnostics,M<:MCMCDiagnostics}
        return new{T,P,M}(base, pf, mcmc)
    end
end

############################################################################################
"""
$(SIGNATURES)
Show relevant diagnostic results.

# Examples
```julia
```

"""
function generate_showvalues(diagnostics::D) where {D<:PMCMCDiagnostics}
    mcmc = generate_showvalues(diagnostics.mcmc)
    pf = generate_showvalues(diagnostics.pf)
    return function showvalues()
        return mcmc()..., pf()...
    end
end

############################################################################################
#=
function get_prediction(diagnostics::PMCMCDiagnostics)
    return get_prediction(diagnostics.pf)
end
=#
############################################################################################
#export
export PMCMCDiagnostics, generate_showvalues


text:
using Mocking

statusEndpoint = "api/$curApiVersion/status"

"""
$(SIGNATURES)
Get the status of the Synchrony service.
"""
function getStatus()::String
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(statusEndpoint)"
    response = @mock _sendRestRequest(config, HTTP.get, url)
    if(response.status != 200)
        error("Error getting server status, received $(response.status) with body '$(String(response.body))'.")
    end
    return String(response.body)
end

"""
$(SIGNATURES)
Print the current service status.
"""
function printStatus()::Nothing
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    try
        serviceStatus = getStatus()
        printstyled("Synchrony service status: ", color=:blue)
        printstyled("$(serviceStatus)\r\n", color=(serviceStatus == "UP!" ? :green : :red))
    catch ex
        println("Unable to get service status, error is:");
        showerror(stderr, ex, catch_backtrace())
    end
end


text:
function add_image_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(adjust_contrast_cb,b["contrast_min_slider"],"value-changed",Void,(),false,(handles,))
    signal_connect(adjust_contrast_cb,b["contrast_max_slider"],"value-changed",Void,(),false,(handles,))

    signal_connect(sharpen_cb,b["sharpen_window"],"value-changed",Void,(),false,(handles,))
    signal_connect(sharpen_cb,b["sharpen_reps"],"value-changed",Void,(),false,(handles,))
    signal_connect(sharpen_cb,b["sharpening_filter_type"],"changed",Void,(),false,(handles,))

    nothing
end

function adjust_contrast_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    try

        han.im_adj.contrast_min = get_gtk_property(han.b["adj_contrast_min"],:value,Int64)
        han.im_adj.contrast_max = get_gtk_property(han.b["adj_contrast_max"],:value,Int64)

        plot_image(han,han.current_frame')

    catch
        println("Error while adjusting contrast")
    end

    nothing
end

function sharpen_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    try

        han.im_adj.sharpen_filter = get_gtk_property(han.b["sharpening_filter_type"],:active,Int64) + 1
        han.im_adj.sharpen_win = get_gtk_property(han.b["adj_sharpen_window"],:value,Int64)
        han.im_adj.sharpen_reps = get_gtk_property(han.b["adj_sharpen_reps"],:value,Int64)

        plot_image(han,han.current_frame')

    catch
        println("Error while adjusting sharpening settings")
    end
    nothing
end

sharpen_mode(b::Gtk.GtkBuilder)=getproperty(b["sharpen_image_button"],:active,Bool)
anisotropic_mode(b::Gtk.GtkBuilder)=getproperty(b["aniso_button"],:active,Bool)
local_contrast_mode(b::Gtk.GtkBuilder)=getproperty(b["local_contrast_enhance"],:active,Bool)


text:
################
# NSClsdict0
################

const _NSClsdict0 = Dict{Symbol, Function}()

_NSClsdict0[:type] = nsc -> nsc.__type

_NSClsdict0[:cls] = nsc -> nsc.__cls

_NSClsdict0[:init] = nsc -> nsc.__init[1]
_NSClsdict0[:post] = nsc -> nsc.__post[1]
_NSClsdict0[:parallel] = nsc -> nsc.__code.__parallel.obj

_NSClsdict0[:cst] = nsc -> nsc.__code.cst
_NSClsdict0[:dfn] = nsc -> nsc.__code.dfn
_NSClsdict0[:req] = nsc -> nsc.__code.req
_NSClsdict0[:prp] = nsc -> nsc.__code.prp
_NSClsdict0[:fnc] = nsc -> nsc.__code.fnc
_NSClsdict0[:mth] = nsc -> nsc.__code.mth
_NSClsdict0[:var] = nsc -> nsc.__code

_NSClsdict0[:instances] = nsc -> [i for (a, k, i) ∈ nsc.__instances]
_NSClsdict0[:clr_instances!] =
    nsc -> (deleteat!(nsc.__instances, 1:length(nsc.__instances)); return)


text:
#=
    NEMO: Next Energy Modeling system for Optimization.
    https://github.com/sei-international/NemoMod.jl

    Copyright © 2019: Stockholm Environment Institute U.S.

	File description: Tests for NemoMod package. Running full suite of tests requires
        GLPK, Cbc, CPLEX, Gurobi, Mosek, and Xpress solvers. However, this file is configured
        so that tests for Cbc and proprietary solvers are skipped if solvers are not present.
=#

if !@isdefined NemoMod
    using NemoMod
end

using Test, SQLite, DataFrames, JuMP, Distributed

const TOL = 1e-3  # Default tolerance for isapprox() comparisons

"""Helper function for deleting a file after Julia has been told to release it (e.g.,
    with finalize(db); db = nothing; GC.gc())."""
function delete_file(path::String, max_del_attempts::Int)
    del_attempts::Int = 0  # Number of deletions that have been attempted

    while isfile(path) && del_attempts < max_del_attempts
        try
            del_attempts += 1
            rm(path; force = true)
            sleep(0.5)  # Without this, the next isfile() test sometimes throws an I/O error
        catch
            # Wait and continue
            sleep(0.5)
        end
    end
end  # delete_file(path::String)

@testset "Solving a scenario" begin
    include(joinpath(@__DIR__, "cbc_tests.jl"))
    include(joinpath(@__DIR__, "glpk_tests.jl"))
    include(joinpath(@__DIR__, "cplex_tests.jl"))
    include(joinpath(@__DIR__, "gurobi_tests.jl"))
    include(joinpath(@__DIR__, "mosek_tests.jl"))
    include(joinpath(@__DIR__, "xpress_tests.jl"))
end  # @testset "Solving a scenario"

@testset "Writing optimization problem for a scenario" begin
    @info "Testing function to write optimization problem for a scenario."
    optprobfile = joinpath(@__DIR__, "storage_test_prob.gz")

    # Wrap file creation operation in a function to ensure Julia releases output file for deletion
    function write_opt_prob()
        dbfile = joinpath(@__DIR__, "storage_test.sqlite")
        chmod(dbfile, 0o777)  # Make dbfile read-write. Necessary because after Julia 1.0, Pkg.add makes all package files read-only

        # Write output file for optimization problem
        NemoMod.writescenariomodel(dbfile; numprocs=1, restrictvars=true, quiet = false, writefilename = optprobfile)

        # Clean up scenario database
        db = SQLite.DB(dbfile)
        NemoMod.dropresulttables(db)
        testqry = SQLite.DBInterface.execute(db, "VACUUM")

        @test isfile(optprobfile)
    end  # write_opt_prob()

    write_opt_prob()
    GC.gc()

    # Try up to 20 times to delete file
    delete_file(optprobfile, 20)

    @test !isfile(optprobfile)
end  # @testset "Writing optimization problem for a scenario"

@testset "Other database operations" begin
    @info "Testing other database operations."

    @testset "Create a new NEMO database" begin
        # Wrap SQLite operations in a function in order to get Julia to release new DB file for deletion
        function new_nemo_db()
            db = NemoMod.createnemodb(joinpath(@__DIR__, "new_nemo.sqlite"))

            @test isfile(joinpath(@__DIR__, "new_nemo.sqlite"))
            # Test that AccumulatedAnnualDemand table exists
            @test !SQLite.done(SQLite.DBInterface.execute(db, "PRAGMA table_info('AccumulatedAnnualDemand')"))
        end  # new_nemo_db()

        new_nemo_db()
        GC.gc()

        # Try up to 20 times to delete file
        delete_file(joinpath(@__DIR__, "new_nemo.sqlite"), 20)

        @test !isfile(joinpath(@__DIR__, "new_nemo.sqlite"))
    end  # @testset "Create a new NEMO database"

    @testset "Set parameter default" begin
        # Wrap SQLite operations in a function in order to get Julia to release new DB file for deletion
        function param_default_db()
            db = NemoMod.createnemodb(joinpath(@__DIR__, "param_default.sqlite"))
            @test SQLite.done(SQLite.DBInterface.execute(db, "SELECT val FROM DefaultParams WHERE tablename = 'VariableCost'"))  # No rows in query result

            NemoMod.setparamdefault(db, "VariableCost", 1.0)

            @test DataFrame(SQLite.DBInterface.execute(db, "SELECT val FROM DefaultParams WHERE tablename = 'VariableCost'"))[1,:val] == 1.0
        end  # param_default_db()

        param_default_db()
        GC.gc()

        # Try up to 20 times to delete file
        delete_file(joinpath(@__DIR__, "param_default.sqlite"), 20)

        @test !isfile(joinpath(@__DIR__, "param_default.sqlite"))
    end  # @testset "Set parameter default"
end  # @testset "Other database operations"


text:
"""
    The model_ functions modify the JuMP model.
    They take in a forecast dictionary called f and a params dictionary called p.
"""

function model_base!(f::Dict, p::Dict)
    """
        Create base model.
        Only variable so far are variables eGrid and gGrid to contain electricity/gas imports.
    """
    T = f[:T]
    sesi = Model(solver = GurobiSolver(OutputFlag = 0))

    @variable(sesi, eGrid[1:T-1] >= 0)
    @variable(sesi, gGrid[1:T-1] >= 0)
    return sesi
end

function model_cap_thstorage!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Choose how much thermal storage to build
    """
    p[:cap_thstorage] = true
    @variable(sesi, hotTankMax >= 0)
    @variable(sesi, coldTankMax >= 0)

    @objective(sesi, Min, getobjective(sesi) + hotTankMax*p[:capCostHot]*f[:T]
        + coldTankMax*p[:capCostCold]*f[:T])
    return sesi
end

function model_cap_machines!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Choose how many machines to build
    """
    p[:cap_machines] = true
    @variables sesi begin
        nHRC, Int
        nHeater, Int
        nChiller, Int
        HRCCAP >= 0
        chillerCAP >= 0
        heaterCAP >= 0
    end

    @constraints sesi begin
        # TODO: take hard-coded numbers out of here
        HRCCAP == nHRC * p[:HRCCAP]/3
        chillerCAP == nChiller * p[:chillerCAP]/4
        heaterCAP == nHeater * p[:heaterCAP]/3
    end

    @objective(sesi, Min, getobjective(sesi)+ nHRC*p[:capCostHRC] +
        nHeater*p[:capCostHeater] + nChiller*p[:capCostChiller])
    return sesi
end

function model_water!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Add water to the model. These are hot and cold water streams.
    """
    T = f[:T];
    @variables sesi begin
        # HRC
        eHRC[1:T-1] >= 0
        hwHRC[1:T-1] >= 0
        cwHRC[1:T-1] >= 0

        # Chiller
        eChiller[1:T-1] >= 0
        cwChiller[1:T-1] >= 0

        # Heater
        eHeater[1:T-1]>=0
        gHeater[1:T-1]>=0
        hwHeater[1:T-1]>=0

        # Storage tanks
        hwTank[1:T]>=0
        cwTank[1:T]>=0

        # unmet loads
        umHot[1:T-1]>=0
        umCold[1:T-1]>=0
    end

    @constraints sesi begin

        # Definition constraints
        [t=1:T-1], eHRC[t] == p[:HRCPE] * cwHRC[t]
        [t=1:T-1], hwHRC[t] == p[:HRCHotToColdRatio] * cwHRC[t]
        [t=1:T-1], eChiller[t] == f[:chillerPE][t] * cwChiller[t]
        [t=1:T-1], eHeater[t] == p[:heaterEPE] * hwHeater[t]
        [t=1:T-1], gHeater[t] == 1/p[:heaterPE] * hwHeater[t]

        # Dynamics
        [t=1:T-1], (hwTank[t+1] == hwTank[t] + hwHRC[t] + hwHeater[t]
            - f[:loadHW][t] + umHot[t])
        [t=1:T-1], (cwTank[t+1] == cwTank[t] + cwHRC[t] + cwChiller[t]
            - f[:loadCW][t] + umCold[t])
        [t=1:T-1], umCold[t] <= f[:loadCW][t]
        [t=1:T-1], umHot[t] <= f[:loadHW][t]
        [t=1:T-1], hwTank[t+1] - hwTank[t] <= p[:S_h_dot_max]
        [t=1:T-1], hwTank[t] - hwTank[t+1] <= p[:S_h_dot_max]
        [t=1:T-1], cwTank[t+1] - cwTank[t] <= p[:S_c_dot_max]
        [t=1:T-1], cwTank[t] - cwTank[t+1] <= p[:S_c_dot_max]
    end

    # Water tank constraints
    if p[:cap_thstorage]
        @constraints sesi begin
            hotTankConstraint[t=1:T-1],( hwTank[t] <=
                sesi[:hotTankMax]*(1-p[:hotTankSafety]))
            coldTankConstraint[t=1:T-1],(cwTank[t] <=
                sesi[:coldTankMax]*(1-p[:coldTankSafety]))
            [t=1:T], hwTank[t] >= sesi[:hotTankMax]*p[:hotTankSafety]
            [t=1:T], cwTank[t] >= sesi[:coldTankMax]*p[:coldTankSafety]
        end
    else
        @constraints sesi begin
            hotTankConstraint[t=1:T],( hwTank[t] <=
                p[:hotTankMax]*(1-p[:hotTankSafety]))
            coldTankConstraint[t=1:T],(cwTank[t] <=
                p[:coldTankMax]*(1-p[:coldTankSafety]))
            [t=1:T], hwTank[t] >= p[:hotTankMax]*p[:hotTankSafety]
            [t=1:T], cwTank[t] >= p[:coldTankMax]*p[:coldTankSafety]
            # Boundary conditions - remove if we are doing capital costs,
            # in that case the horizons are long so it does not matter
            hwTank[1] == f[:hwTankInit]
            cwTank[1] == f[:cwTankInit]
            hwTank[T] >= f[:hwTankEnd]
            cwTank[T] >= f[:cwTankEnd]
        end
    end

    # Capacity constraints
    if p[:cap_machines]
        @constraints sesi begin
            [t=1:T-1], cwHRC[t] <= HRCCAP
            [t=1:T-1], cwChiller[t] <= chillerCAP
            [t=1:T-1], hwHeater[t] <= heaterCAP
        end
    else
        @constraints sesi begin
            [t=1:T-1], cwHRC[t] <= p[:HRCCAP]
            [t=1:T-1], cwChiller[t] <= p[:chillerCAP]
            [t=1:T-1], hwHeater[t] <= p[:heaterCAP]
        end
    end

    if p[:umLoadsL2]  # use L2 penalty
        @objective(sesi, Min, getobjective(sesi) + p[:umHot] * umHot'umHot
            + p[:umCold] * umCold'umCold)
    else  # use L1 penalty
        @objective(sesi, Min, getobjective(sesi) + p[:umHot] * sum(umHot)
            + p[:umCold] * sum(umCold))
    end

    return sesi
end

function model_water_switching!(sesi::JuMP.Model, f::Dict, p::Dict)
    """ Model switching constraints for heating/cooling machines.
    """
    T = f[:T]
    p[:sw] = true
    @variables sesi begin
        swHRC[1:T-2] >= 0
        swHeater[1:T-2] >= 0
        swChiller[1:T-2] >= 0
    end

     for t=1:T-2
         @constraints sesi begin
             swHRC[t] >= sesi[:eHRC][t] - sesi[:eHRC][t+1]
             swHRC[t] >= sesi[:eHRC][t+1] - sesi[:eHRC][t]
             swHeater[t] >= sesi[:eHeater][t] - sesi[:eHeater][t+1]
             swHeater[t] >= sesi[:eHeater][t+1] - sesi[:eHeater][t]
             swChiller[t] >= sesi[:eChiller][t] - sesi[:eChiller][t+1]
             swChiller[t] >= sesi[:eChiller][t+1] - sesi[:eChiller][t]
         end
     end

    # norm 2 penalty makes more sense
    @objective(sesi, Min, getobjective(sesi) + p[:swPenalty] * (swHRC'swHRC +
        swChiller'swChiller + swHeater'swHeater))

    f[:swPenalty] = p[:swPenalty]

    return sesi
end

function model_storage!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Add electrical storage to the model.
    """
    T = f[:T]

    # Electrical storage
    @variable(sesi, elecTransfer_in[1:T-1]>=0)
    @variable(sesi, elecTransfer_out[1:T-1]>=0)

    # Amount in electrical storage
    @variable(sesi, elecStored[1:T] >= 0)

    # Boundary conditions
    @constraint(sesi,elecStored[1] == f[:elecStorageInit])
    @constraint(sesi, elecStored[T] == p[:elecStorageInit])

    for t = 1:T-1
        # Operating constraints
        @constraint(sesi, elecTransfer_in[t] <= p[:elecTransferMax])
        @constraint(sesi, elecTransfer_out[t] <= p[:elecTransferMax])

        # Max storage constraint
        @constraint(sesi, elecStored[t] <= p[:elecStorageMax])

        # Dynamics
        @constraint(sesi, elecStored[t+1] == elecStored[t]
            + p[:elecTransferEta] * elecTransfer_in[t]
            - elecTransfer_out[t] / p[:elecTransferEta])
    end

    @objective(sesi, Min, getobjective(sesi) + p[:usageCostBattery]
        * sum(elecTransfer_out[t] for t=1:T-1))
    return sesi
end

function model_demand!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Add demand-related variables to the model.
        Create a monthly indexed variable to contain peak demand.
        Add constraint that peak is greater then eGrid at all hours of the month.
        Add cost of demand charge in objective function
    """
    T, M = f[:T], f[:M]
    p[:peak] = true
    @variable(sesi, peak[1:M] >= 0)
    for t = 1:T-1
        # m = getMonth(t - 1 + f["startInd"]) - getMonth(f["startInd"]) + 1
        @constraint(sesi, peak[getMonthNum(t, f)] >= sesi[:eGrid][t])
    end
    @objective(sesi, Min, getobjective(sesi)
               + sum(peak[m] * f[:demandCharge][m] for m = 1:M))

    return sesi
end

function model_carbon!(sesi::JuMP.Model, f::Dict, p::Dict)

    T = f[:T]
    p[:carbon]=true
    @objective(sesi, Min, getobjective(sesi)
        + f[:carbonPrice] * (sum(f[:carbonIntensityGrid][t]
            * sesi[:eGrid][t] for t = 1:T-1) + f[:carbonIntensityGas]
            * sum(sesi[:gGrid][t] for t = 1:T-1)))

    return sesi
end

function model_energy_price!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Add variable energy price to the model.
    """
    T = f[:T]
    @objective(sesi, Min, getobjective(sesi) + sum(sesi[:eGrid][t] * f[:elecPrice][t]
        + sesi[:gGrid][t] * f[:gasPrice] for t = 1:T-1))
    return sesi
end

function model_finalize!(sesi::JuMP.Model, f::Dict, p::Dict)
    """
        Add aggregate constraints to the model.
    """
    T = f[:T]
    for t = 1:T-1
        @expression(sesi, eGridExpr, f[:loadKWH][t])
        for var in [:eChiller, :eHeater, :eHRC, :elecTransfer_in, :elecEV_in]
            if haskey(sesi.objDict, var)
                @expression(sesi, eGridExpr, eGridExpr + sesi[var][t])
            end
        end
        if haskey(sesi.objDict, :elecTransfer_out)
            @expression(sesi, eGridExpr, eGridExpr -
			sesi[:elecTransfer_out][t])
        end
        @constraint(sesi, sesi[:eGrid][t] == eGridExpr)
    end
    for t = 1:T-1
        @expression(sesi, gGridExpr, 0)
        for var in [:gHeater]
            if haskey(sesi.objDict, var)
                @expression(sesi, gGridExpr, gGridExpr + sesi[var][t])
            end
        end
        @constraint(sesi, sesi[:gGrid][t] == gGridExpr)
    end
    return sesi
end

max_stars_repo_name: jdechalendar/sesi

max_issues_repo_name: jdechalendar/sesi

max_forks_repo_name: jdechalendar/sesi


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#
"""
`Actors` implements the classical Actor Model and is
based on the primitives defined in `ActorInterfaces.Classic`.
It provides:

- basic primitives for creating actors,
    sending messages to them and changing behavior:
    [`spawn`](@ref), [`send`](@ref), [`become`](@ref)
    with `Addr` and [`self`](@ref),
- [`onmessage`](@ref), executed by an actor on a
    received message,
- a message protocol with predefined messages,
- an API based on the protocol with primitives
    [`receive`](@ref) and [`request`](@ref) and further
    API functions [`become!`](@ref), [`call`](@ref),
    [`cast`](@ref), [`exec`](@ref), [`exit!`](@ref),
    [`init!`](@ref), [`query`](@ref), [`term!`](@ref),
    [`update!`](@ref),
- error handling with actor
    - connections: [`connect`](@ref), [`disconnect`](@ref), [`trapExit`](@ref),
    - monitors: [`monitor`](@ref), [`demonitor`](@ref),
    - supervisors: [`supervisor`](@ref), [`supervise`](@ref),
        [`unsupervise`](@ref), [`start_actor`](@ref),
        [`start_task`](@ref), [`count_children`](@ref),
        [`which_children`](@ref), [`terminate_child`](@ref),
- an actor registry: [`register`](@ref), [`unregister`](@ref),
    [`whereis`](@ref), [`registered`](@ref)

and more.

The current stable, registered version is installed with
```julia
pkg> add Actors
```

The development version is installed with:
```julia
pkg> add "https://github.com/JuliaActors/Actors.jl"
```
"""
module Actors

"Gives the package version."
const version = v"0.2.5"

using ActorInterfaces.Classic
using Distributed, .Threads, Serialization, Dates, Proquint
import ActorInterfaces.Classic: onmessage
import Base: structdiff

include("types.jl")
include("messages.jl")
include("links.jl")
include("com.jl")
include("connections.jl")
include("logging.jl")
include("errorkernel.jl")
include("supervisor.jl")
include("remote_failures.jl")
include("checkpoints.jl")
include("protocol.jl")
include("actor.jl")
include("api.jl")
include("registry.jl")
include("init.jl")
include("diag.jl")
include("utils.jl")

export
    # common types
    Msg, @msg, Request, Response, Link, Bhv,
    # -------------------------------
    # exported interface primitives
    send,
    # the following ones must be imported explicitly:
    # - newLink,
    # - spawn,
    # - _ACT,
    # - onmessage
    # -------------------------------
    # API primitives
    receive, request, send_after,
    # API
    Args, become, self, stop,
    become!, call, cast, exec, exit!, info, init!,
    query, term!, update!,
    # registry
    register, unregister, whereis, registered,
    # error handling
    connect, disconnect, monitor, demonitor, trapExit,
    # supervising
    supervisor, supervise, unsupervise,
    set_strategy, count_children, which_children,
    delete_child, start_actor, start_task, terminate_child,
    # checkpointing
    checkpointing, checkpoint, restore, get_checkpoints,
    save_checkpoints, load_checkpoints, @chkey

end


text:
@testset "vector" begin

    x = [1, 2, 3]

    @test is_transformable(x)

    @testset "apply" begin

        @testset "OneToOne" begin
            T = FakeOneToOneTransform()
            @test FeatureTransforms.apply(x, T) == ones(3)
            @test FeatureTransforms.apply(x, T; inds=[2, 3]) == ones(2)
        end

        @testset "OneToMany" begin
            T = FakeOneToManyTransform()
            @test FeatureTransforms.apply(x, T) == [ones(3) ones(3)]
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == [ones(2) ones(2)]
        end

        @testset "ManyToOne" begin
            T = FakeManyToOneTransform()
            @test FeatureTransforms.apply(x, T; dims=1) == fill(1.0)
            @test FeatureTransforms.apply(x, T; dims=1, inds=[1, 2]) == fill(1.0)
            @test_throws BoundsError FeatureTransforms.apply(x, T; dims=2)
        end

        @testset "ManyToMany" begin
            T = FakeManyToManyTransform()
            @test FeatureTransforms.apply(x, T; dims=1) == [ones(3) ones(3)]
            @test FeatureTransforms.apply(x, T; dims=1, inds=[1, 3]) == [ones(2) ones(2)]
            @test_throws BoundsError FeatureTransforms.apply(x, T; dims=2)
        end

    end

    @testset "apply!" begin
        T = FakeOneToOneTransform()

        _x = copy(x)
        FeatureTransforms.apply!(_x, T)
        @test _x == ones(3)

        # https://github.com/invenia/FeatureTransforms.jl/issues/68
        _x = copy(x)
        @test_broken FeatureTransforms.apply!(_x, T; inds=[1, 2])
        @test_broken _x == [1, 1, 3]
    end

    @testset "apply_append" begin

        @testset "OneToOne" begin
            T = FakeOneToOneTransform()
            @test FeatureTransforms.apply_append(x, T; append_dim=1) == [1, 2, 3, 1, 1, 1]
            @test FeatureTransforms.apply_append(x, T; append_dim=2) == [x ones(3)]
            @test FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=1) == [1, 2, 3, 1, 1]
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=2)
        end

        @testset "OneToMany" begin
            T = FakeOneToManyTransform()
            @test FeatureTransforms.apply_append(x, T; append_dim=2) == [x ones(3) ones(3)]
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=2)
        end

        @testset "ManyToOne" begin
            T = FakeManyToOneTransform()
            @test FeatureTransforms.apply_append(x, T; dims=1, append_dim=1) == [1, 2, 3, 1]
            @test_throws BoundsError FeatureTransforms.apply_append(x, T; dims=1, append_dim=2)
            @test_throws BoundsError FeatureTransforms.apply_append(x, T; dims=2, append_dim=1)
            @test_throws BoundsError FeatureTransforms.apply_append(x, T; dims=2, append_dim=2)
        end

        @testset "ManyToMany" begin
            T = FakeManyToManyTransform()
            @test FeatureTransforms.apply_append(x, T; append_dim=2) == [x ones(3) ones(3)]
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=2)
        end

    end

end


text:
export postwalk, prewalk, blockfilter!, blockfilter, collect_blocks, gatecount

"""
    parse_block(n, ex)

This function parse the julia object `ex` to a quantum block,
it defines the syntax of high level interfaces. `ex` can be
a function takes number of qubits `n` as input or it can be
a pair.
"""
function parse_block end

parse_block(n::Int, x::Function) = x(n)

parse_block(n::Int, ex) =
    throw(Meta.ParseError("cannot parse expression $ex, expect a pair or quantum block"))

function parse_block(n::Int, x::AbstractBlock{D}) where {D}
    n == nqudits(x) || throw(ArgumentError("number of qubits does not match: $x"))
    return x
end

# if it is a single qubit pair, parse it to put block
parse_block(n::Int, x::Pair{Int,<:AbstractBlock}) = error("got $x, do you mean put($x)?")
# infer the number of qubits if the inner function was curried
parse_block(n::Int, x::Pair{Int,<:Function}) = error("got $x, do you mean put($x)?")

# error if it is not single qubit case
function parse_block(n::Int, x::Pair)
    error(
        "please specifiy the block type of $x, consider to use concentrate for large block in local scope.",
    )
end

"""
    prewalk(f, src::AbstractBlock)

Walk the tree and call `f` once the node is visited.
"""
function prewalk(f::Base.Callable, src::AbstractBlock)
    out = f(src)
    for each in subblocks(src)
        prewalk(f, each)
    end
    return out
end

"""
    postwalk(f, src::AbstractBlock)

Walk the tree and call `f` after the children are visited.
"""
function postwalk(f::Base.Callable, src::AbstractBlock)
    for each in subblocks(src)
        postwalk(f, each)
    end
    return f(src)
end

blockfilter!(f, v::Vector, blk::AbstractBlock) = postwalk(x -> f(x) ? push!(v, x) : v, blk)

blockfilter(f, blk) = blockfilter!(f, [], blk)

"""
    collect_blocks(block_type, root)

Return a [`ChainBlock`](@ref) with all block of `block_type` in root.
"""
collect_blocks(::Type{T}, x::AbstractBlock) where {T<:AbstractBlock} =
    blockfilter!(x -> x isa T, T[], x)

#expect(op::AbstractBlock, r::AbstractRegister) = r' * apply!(copy(r), op)

#expect(op::AbstractBlock, dm::DensityMatrix) = mapslices(x->sum(mat(op).*x)[], dm.state, dims=[1,2]) |> vec
expect(op::AbstractBlock, dm::DensityMatrix{1}) =
    sum(mat(op) .* dropdims(dm.state, dims = 3))

"""
    expect(op::AbstractBlock, reg) -> Vector
    expect(op::AbstractBlock, reg => circuit) -> Vector
    expect(op::AbstractBlock, density_matrix) -> Vector

Get the expectation value of an operator, the second parameter can be a register `reg` or a pair of input register and circuit `reg => circuit`.

expect'(op::AbstractBlock, reg=>circuit) -> Pair
expect'(op::AbstractBlock, reg) -> AbstracRegister

Obtain the gradient with respect to registers and circuit parameters.
For pair input, the second return value is a pair of `gψ=>gparams`,
with `gψ` the gradient of input state and `gparams` the gradients of circuit parameters.
For register input, the return value is a register.

!!! note

    For batched register, `expect(op, reg=>circuit)` returns a vector of size number of batch as output. However, one can not differentiate over a vector loss, so `expect'(op, reg=>circuit)` accumulates the gradient over batch, rather than returning a batched gradient of parameters.
"""
function expect(op::AbstractBlock, dm::DensityMatrix)
    mop = mat(op)
    return sum(transpose(dm.state) .* mop)
end

expect(op::AbstractBlock, reg::ArrayReg) = reg' * apply!(copy(reg), op)

function expect(op::AbstractBlock, reg::BatchedArrayReg)
    B = YaoArrayRegister._asint(nbatch(reg))
    ket = apply!(copy(reg), op)
    if !(reg.state isa Transpose)
        C = conj!(reshape(ket.state, :, B))
        A = reshape(reg.state, :, B)
        dropdims(sum(A .* C, dims = 1), dims = 1) |> conj
    elseif size(reg.state, 2) == B
        Na = size(reg.state, 1)
        C = conj!(reshape(ket.state.parent, B, Na))
        A = reshape(reg.state.parent, B, Na)
        dropdims(sum(A .* C, dims = 2), dims = 2) |> conj
    else
        Na = size(reg.state, 1)
        C = conj!(reshape(ket.state.parent, :, B, Na))
        A = reshape(reg.state.parent, :, B, Na)
        dropdims(sum(A .* C, dims = (1, 3)), dims = (1, 3)) |> conj
    end
end

for REG in [:ArrayReg, :BatchedArrayReg]
    @eval function expect(op::Add, reg::$REG)
        sum(opi -> expect(opi, reg), op)
    end
    @eval function expect(op::Scale, reg::$REG)
        factor(op) * expect(content(op), reg)
    end
end

function expect(op, plan::Pair{<:AbstractRegister,<:AbstractBlock})
    expect(op, copy(plan.first) |> plan.second)
end

# obtaining Dense Matrix of a block
LinearAlgebra.Matrix(blk::AbstractBlock) = Matrix(mat(blk))

"""
    operator_fidelity(b1::AbstractBlock, b2::AbstractBlock) -> Number

Operator fidelity defined as

```math
F^2 = \\frac{1}{d^2}\\left[{\\rm Tr}(b1^\\dagger b2)\\right]
```

Here, `d` is the size of the Hilbert space. Note this quantity is independant to global phase.
See arXiv: 0803.2940v2, Equation (2) for reference.
"""
function operator_fidelity(b1::AbstractBlock, b2::AbstractBlock)
    U1 = mat(b1)
    U2 = mat(b2)
    return abs(sum(conj(U1) .* U2)) / size(U1, 1)
end

gatecount(blk::AbstractBlock) = gatecount!(blk, Dict{Type{<:AbstractBlock},Int}())
function gatecount!(
    c::Union{ChainBlock,KronBlock,PutBlock,Add,CachedBlock},
    storage::AbstractDict,
)
    (gatecount!.(c |> subblocks, Ref(storage)); storage)
end

function gatecount!(c::RepeatedBlock, storage::AbstractDict)
    k = typeof(content(c))
    n = length(c.locs)
    if haskey(storage, k)
        storage[k] += n
    else
        storage[k] = n
    end
    storage
end

function gatecount!(c::Union{PrimitiveBlock,Daggered,ControlBlock}, storage::AbstractDict)
    k = typeof(c)
    if haskey(storage, k)
        storage[k] += 1
    else
        storage[k] = 1
    end
    storage
end


text:
module Trixi2Vtk

# Include other packages
using Glob: glob
using HDF5: h5open, attributes, haskey
using ProgressMeter: @showprogress, Progress, next!
using StaticArrays: SVector
using TimerOutputs
using Trixi: Trixi, transfinite_mapping, coordinates2mapping, polynomial_interpolation_matrix,
             gauss_lobatto_nodes_weights, TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh
using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, vtk_save, paraview_collection

# Include all top-level submodule files
include("interpolate.jl")
include("io.jl")
include("pointlocators.jl")
include("vtktools.jl")

# Include top-level conversion method
include("convert.jl")

# export types/functions that define the public API of Trixi2Vtk
export trixi2vtk

end # module Trixi2Vtk


text:
function preprocessKARC(Hop, g, params::Tparams, calls, max_calls) #where T
    τ = params.τ
    ζ = params.ζ
    nshifts = params.nshifts
    shifts = params.shifts

    n = length(g)
    gNorm2 = BLAS.nrm2(n, g, 1)
    precision =  max(1e-12,min(0.5,(gNorm2^ζ)))

    #ϵ = sqrt(eps(T)) # * 100.0
    ϵ = 1e-12#sqrt(eps()) # * 100.0
    cgtol = max(ϵ, min(0.09, 0.01 * norm(g)^(1.0 + ζ)))

    (xShift, stats) = cg_lanczos_shift_seq(Hop,
                                           -g,
                                           shifts,
                                           itmax=min(max_calls-sum(calls),2*n),
                                           #τ = τ,
                                           atol = 1.0e-8,
                                           rtol = precision,
                                           verbose=false,
                                           check_curvature=true)

#    (xShift, stats) = cg_lanczos_shift_seq(Hop,
#                                           -g,
#                                           shifts,
#                                           itmax=min(max_calls-sum(calls),2*n),
#                                           #τ = τ,
#                                           atol = cgtol,
#                                           rtol = ϵ,
#                                           verbose=false,
#                                           check_curvature=true)

    positives = collect(findfirst(!, stats.flagged):length(stats.flagged))

    success = false
    good_grad = false
    if VERSION >= v"1.1.0"
        xShift = hcat(xShift...)
    end
    dirs = [ (xShift[:,i]) for i = 1 : nshifts ];
    Ndirs = map(norm, dirs);

    d = g # bidon
    return  PDataK(d, -1.0, ζ, τ, 0 , positives, xShift,   shifts, nshifts, Ndirs, true)
end

function decreaseKARC(X :: PDataK, α:: Float64, TR:: TrustRegion)
    X.indmin += 1
    p_imin = X.positives[X.indmin]
    α2 = X.norm_dirs[p_imin]^X.τ/X.shifts[p_imin]

    targetα = α*TR.decrease_factor

    # fix α to its "ideal" value to satisfy αλ=||d||^τ
    # while ensuring α decreases enough
    while α2 > targetα
        X.indmin += 1
        p_imin = X.positives[X.indmin]
        α2 = X.norm_dirs[p_imin]^X.τ/X.shifts[p_imin]
    end

    X.d = X.xShift[:,p_imin]
    X.λ = X.shifts[p_imin]

    return α2
end


text:
module MechGlueInterpolations
import Unitfu: AbstractQuantity, Quantity, ustrip, norm, unit, zero
import Unitfu: Dimensions, FreeUnits, m
import Interpolations
using Interpolations: interpolate, extrapolate, Gridded, Linear
include("matrix_interpolation.jl")
end


text:
using EffectiveWaves, ClassicalOrthogonalPolynomials
using LinearAlgebra, Statistics, Test

@testset "Discrete radial system for sphere scattering" begin

## Set parameters
    particle_medium = Acoustic(3; ρ=0.1, c=0.1);
    medium = Acoustic(3; ρ=1.0, c=1.0);

    R = 5.0; r = 1.0

    separation_ratio = 1.02
    kas = [0.01,0.2]
    ks = kas ./ r

    vol_fraction = 0.15

    basis_orders = Int.(round.(4.0 .* kas)) .+ 1
    basis_field_orders = Int.(round.(3.0 .* ks .* R)) .+ 1

    ωs = ks .* real(medium.c)

    s1 = Specie(
        particle_medium, Sphere(r);
        number_density = vol_fraction / volume(Sphere(r)),
        exclusion_distance = separation_ratio
    );

    species = [s1]
    # species = [s1,s1]

    region_shape = Sphere([0.0,0.0,0.0], R)
    material = Material(Sphere(R),species);

    sourceradial = regular_spherical_source(medium, [1.0+0.0im];
       position = [0.0,0.0,0.0], symmetry = RadialSymmetry{3}()
    );

    a12 = 2.0 * s1.exclusion_distance * outer_radius(s1)
    rs = 0.0:0.02:(R - outer_radius(s1));
    xs = [
        radial_to_cartesian_coordinates([r,0.5,0.6])
    for r in rs];

## The full discrete system does not calculate the translation matrix for |r1 - r2| <= a12 , whereas the radial method does. To get an exact match between the methods we will use a special pair-correlation which does not allow two particles on the same spherical annulus

    polynomial_order = 1
    P = Legendre()

    function gls_radial_simple(r1,r2)
        if abs(r1-r2) < a12
            zeros(polynomial_order+1)
        else
            1 ./ (1:(polynomial_order+1)) .* cos(abs(r1-r2) * pi/(2R))
        end
    end

    function pair_radial_simple(r1,r2,u)
        Pus = P[u, 1:(polynomial_order + 1)] .* (2 .* (0:polynomial_order) .+ 1) ./ (4pi)
        return sum(Pus .* gls_radial_simple(r1,r2))
    end

    function pair_corr_simple(x1,s1,x2,s2)
        if norm(x1) < 1e-12 || norm(x2) < 1e-12
            pair_radial_simple(norm(x1),norm(x2),0.0)
        else
            pair_radial_simple(norm(x1),norm(x2),dot(x1,x2) / (norm(x1)*norm(x2)))
        end
    end

    discrete_field_radials = [
        EffectiveWaves.discrete_system_radial(ωs[i], sourceradial, material, Symmetry(sourceradial,material);
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            polynomial_order = polynomial_order,
            gls_pair_radial = gls_radial_simple,
            sigma_approximation = false
        )
    for i in eachindex(ωs)];

    discrete_rad_scats = [d.coefficient_field.(xs) for d in discrete_field_radials];

    mat_dcoefs_radial = [
        material_scattering_coefficients(discrete_field_radials[i]; rtol = 1e-3,maxevals = Int(5e3))
    for i in eachindex(ωs)];

    # fully numerical method
    tmp = discrete_system(ωs[1], sourceradial, material;
        basis_order = 0, basis_field_order = 0, rtol = 1.0, maxevals = 4
    );

    # import EffectiveWaves: discrete_system
    ST = typeof(tmp);
    discrete_fields = ST[
        discrete_system(ωs[i], sourceradial, material;
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            rtol = 5e-3, maxevals = Int(1e5),
            pair_corr = pair_corr_simple
            # pair_corr = pair_corr_smooth
        )
    for i in eachindex(ωs)];

    discrete_scats = [d.coefficient_field.(xs) for d in discrete_fields];

    mat_dcoefs = [
        material_scattering_coefficients(discrete_fields[i]; rtol = 1e-3,maxevals = Int(5e3))
    for i in eachindex(ωs)];

    errors = [
        abs(mat_dcoefs_radial[i][1] - mat_dcoefs[i][1]) / abs(mat_dcoefs[i][1])
    for i in eachindex(mat_dcoefs_radial)]

    @test errors[1] < 1e-4
    @test errors[2] < 2e-3

    errors = [
        norm(discrete_rad_scats[i] - discrete_scats[i]) / norm(discrete_rad_scats[i])
    for i in eachindex(discrete_rad_scats)];

    @test errors[1] < 1e-3
    @test errors[2] < 3e-3

## For a pair correlation that depends only on inter particle distances, the two discrete methods are only approximately the same when using a high order approximation for the pair correlation.

    # use 0 basis order to speed up code. Large order was tested above
    basis_orders = basis_orders .* 0

    polynomial_order = 60
    polynomial_order = 20
    pair_corr_inf(z) = hole_correction_pair_correlation([0.0,0.0,0.0],s1, [0.0,0.0,z],s1)

    # pair_corr_inf_smooth = smooth_pair_corr_distance(
    #     pair_corr_inf, a12;
    #     smoothing = 0.0, max_distance = 2R,
    #     polynomial_order = polynomial_order
    # )

    # using Plots
    # zs = 0.0:0.01:(2R)
    # plot(pair_corr_inf_smooth,zs)
    # plot!(pair_corr_inf,zs, linestyle=:dash)

    gls_radial = gls_pair_radial_fun(pair_corr_inf, a12;
        sigma_approximation = false,
        polynomial_order = polynomial_order
    )

    # plot(abs.(gls_radial(3.0,3.0)), ylims = (0.0,2.5))

    P = Legendre{Float64}()
    function pair_radial_smooth2(r1,r2,u)
        Pus = P[u, 1:(polynomial_order + 1)] .* (2 .* (0:polynomial_order) .+ 1) ./ (4pi)

        return sum(Pus .* gls_radial(r1,r2))
    end
    #
    # plot(zs,pair_radial_smooth2.(zs,0.0,0.0))
    # plot!(zs,pair_corr_inf_smooth.(zs), linestyle=:dash)

    function pair_corr_smooth(x1,s1,x2,s2)
        if norm(x1) < 1e-12 || norm(x2) < 1e-12
            pair_radial_smooth2(norm(x1),norm(x2),0.0)
        else
            pair_radial_smooth2(norm(x1),norm(x2),dot(x1,x2) / (norm(x1)*norm(x2)))
        end
    end

## Reduced radial method
    # import EffectiveWaves: discrete_system_radial
    discrete_field_radials = [
        discrete_system_radial(ωs[i], sourceradial, material, Symmetry(sourceradial,material);
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            polynomial_order = polynomial_order,
            gls_pair_radial = gls_radial,
            sigma_approximation = false
        )
    for i in eachindex(ωs)]

    discrete_rad_scats = [
        d.coefficient_field.(xs)
    for d in discrete_field_radials];

    mat_dcoefs_radial = [
        material_scattering_coefficients(discrete_field_radials[i];
            rtol = 1e-3, maxevals = Int(5e3))
    for i in eachindex(ωs)];

## fully numerical method
    tmp = discrete_system(ωs[1], sourceradial, material;
        basis_order = 0, basis_field_order = 0, rtol = 1.0, maxevals = 4
    );

    # import EffectiveWaves: discrete_system
    ST = typeof(tmp);
    discrete_fields = ST[
        discrete_system(ωs[i], sourceradial, material;
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            rtol = 5e-3, maxevals = Int(1e5),
            pair_corr = pair_corr_smooth
            # pair_corr = hole_correction_pair_correlation
        )
    for i in eachindex(ωs)];

    discrete_scats = [
        d.coefficient_field.(xs)
    for d in discrete_fields];

    mat_dcoefs = [
        material_scattering_coefficients(discrete_fields[i]; rtol = 1e-3,maxevals = Int(5e3))
    for i in eachindex(ωs)];

    errors = [
        abs(mat_dcoefs_radial[i][1] - mat_dcoefs[i][1]) / abs(mat_dcoefs[i][1])
    for i in eachindex(mat_dcoefs_radial)]

    @test errors[1] < 1e-3
    @test errors[2] < 1e-2

    errors = [
        norm(discrete_rad_scats[i] - discrete_scats[i]) / norm(discrete_rad_scats[i])
    for i in eachindex(discrete_rad_scats)];

    # i = 2
    # fun = abs
    # plot(rs, abs.([d[1] for d in discrete_scats[i]]))
    # plot!(rs, abs.([d[1] for d in discrete_rad_scats[i]]), linestyle=:dash)

    @test errors[1] < 2e-3
    @test errors[2] < 2e-2

   #  # mat_coefs_radial = material_scattering_coefficients.(wavemodes_radial);
   # abs(mat_dcoefs_radial[i][1] - mat_dcoefs[i][1]) / abs(mat_dcoefs[i][1])
   # # abs(mat_dcoefs[i][1] - mat_coefs_radial[i][1]) / abs(mat_coefs_radial[i][1])
   # # abs(mat_dcoefs_radial[i][1] - mat_coefs_radial[i][1]) / abs(mat_coefs_radial[i][1])
   #
   # j = 1;
   # i = 2;
   # # eff_rad0s = [s[j] for s in eff_scats_radial[i]];
   #
   # df_rad0s = [s[j] for s in discrete_rad_scats[i]];
   # # df_rad1s = [s[3] for s in discrete_rad_scats[2]];
   # df0s = [d[j] for d in discrete_scats[i]];
   # # df1s = [d[3] for d in discrete_scats[2]];
   #
   # # rad_scats2 = scattered_field.(xs);
   # # df_rad2s = [s[j] for s in rad_scats2];
   #
   # zs = [x[3] for x in xs];
   # zs - norm.(xs);
   #
   # using Plots
   # gr(linewidth=2.0)
   # fun = abs;
   # plot(zs,fun.(df0s), lab = "$fun fully discrete")
   # plot!(zs,fun.(df_rad0s), linestyle = :dash, lab = "$fun discrete radial")
   # plot!(zs,fun.(eff_rad0s),
   #      linestyle = :dot,
   #      lab = "$fun eff. wave method",
   #      xlab = "radial distance",
   #      ylab = ""
   # )
end

@testset "Effective sphere scattering" begin

## Set parameters
    particle_medium = Acoustic(3; ρ=10.0, c=10.0);
    medium = Acoustic(3; ρ=1.0, c=1.0);

    R = 5.0
    r = 1.0

    separation_ratio = 1.02

    kas = [0.01,0.2]
    ks = kas ./ r

    vol_fraction = 0.12

    basis_orders = Int.(round.(4. .* kas)) .+ 1
    basis_orders = min.(basis_orders,1)

    basis_field_orders = Int.(round.(4.0 .* ks .* R)) .+ 1
    basis_field_orders = max.(basis_field_orders,2)

    ωs = ks .* real(medium.c)

    s1 = Specie(
        particle_medium, Sphere(r);
        number_density = vol_fraction / volume(Sphere(r)),
        exclusion_distance = separation_ratio
    );

    species = [s1]
    # species = [s1,s1]

    a12 = 2.0 * outer_radius(s1) * s1.exclusion_distance

## define sources and material

    psource = PlaneSource(medium, [0.0,0.0,1.0]);
    source = plane_source(medium; direction = [0.0,0.0,1.0])

    region_shape = Sphere([0.0,0.0,0.0], R)
    material = Material(Sphere(R),species);

    eff_medium = effective_medium(medium, species; numberofparticles = material.numberofparticles)
    ks_low = ωs ./ eff_medium.c

    keff_arr = [
        wavenumbers(ωs[i], medium, species;
            # num_wavenumbers = 4,
            basis_order = basis_orders[i],
            tol = 1e-7,
            numberofparticles = material.numberofparticles
        )
    for i in eachindex(ωs)]

    keffs = [ks[1] for ks in keff_arr]

## Plane wave reflection from a sphere

    ## effective waves solution
    pwavemodes_azi = [
        WaveMode(ωs[i], keffs[i], psource, material;
           basis_order = basis_orders[i],
           basis_field_order = basis_field_orders[i]
           , source_basis_field_order = basis_field_orders[i]
        )
    for i in eachindex(ωs)];

    pscat_fields = scattering_field.(pwavemodes_azi);

    ## discrete numerical solution of the average integral equations

    # increasing these parameters does lead to more accurate solutions, but convergences is slow. To increase accuracy you need to increase basis_field_order and maxevals.
    rtol = 1e-2; maxevals = Int(5e3);

    # this below is just to get the typeof tmp, to then avoid a weird unionall error which should hopefully go away when updating Julia at some point
    tmp = discrete_system(ωs[1], psource, material;
        basis_order = 0, basis_field_order = 0, rtol = 1.0, maxevals = 2
    );
    ST = typeof(tmp);
    discrete_fields = ST[
        discrete_system(ωs[i], psource, material;
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            rtol = rtol, maxevals = maxevals
        )
    for i in eachindex(ωs)];

    # avoid right next to the surface due the boundary layer
    rs = 0.0:0.01:(R - 2*outer_radius(s1));
    xs = [ radial_to_cartesian_coordinates([r,0.2,1.2]) for r in rs];

    errors = [
        norm.(discrete_fields[i].coefficient_field.(xs) - pscat_fields[i].(xs)) ./ norm.(pscat_fields[i].(xs))
    for i in eachindex(ωs)];

    @test maximum(mean.(errors)) < 0.01
    @test maximum(maximum.(errors)) < 0.02

    # The scattering coefficients from the whole sphere has a smaller error. Indicating that the errors above occur more on the higher order modes which are smaller
    mat_coefs_pwaves = material_scattering_coefficients.(pwavemodes_azi);
    mat_coefs_discretes = material_scattering_coefficients.(discrete_fields;
        rtol = rtol,
        maxevals = maxevals
    );

    errors = [
        abs(norm(mat_coefs_pwaves[i][1:length(mat_coefs_discretes[i])]) / norm(mat_coefs_discretes[i]) - 1.0)
    for i in eachindex(ωs)];

    @test errors[1] < 1e-4
    @test errors[2] < 3e-3

## Radially symmetric scattering from a sphere

   # Define two radially symmetric sources, but only dispatch for radial symmetry on one of them
    sourceradial =  regular_spherical_source(medium, [1.0+0.0im];
       position = [0.0,0.0,0.0], symmetry = RadialSymmetry{3}()
    );

    sourceazi =  regular_spherical_source(medium, [1.0+0.0im];
       position = [0.0,0.0,0.0], symmetry = AzimuthalSymmetry{3}()
    );

    wavemodes_azi = [
        WaveMode(ωs[i], keffs[i], sourceazi, material;
           basis_order = basis_orders[i],
           basis_field_order = basis_field_orders[i]
        )
    for i in eachindex(ωs)];

    # Note there is no basis_field_order for RadialSymmetry below, due to symmetry restrictions
    wavemodes_radial = [
        WaveMode(ωs[i], keffs[i], sourceradial, material;
           basis_order = basis_orders[i]
        )
    for i in eachindex(ωs)];

    scat_fields_azi = scattering_field.(wavemodes_azi);
    scat_fields_radial = scattering_field.(wavemodes_radial);

    # pscat_field = scattering_field(pwavemode)

    # res = discrete_system_residue(pscat_field, ω, source, material, AzimuthalSymmetry{3}();
    #     basis_order = basis_order, mesh_points = 5,
    #     rtol = 1e-2, maxevals = Int(1e4)
    # )

    rtol = 1e-2; maxevals = Int(1e4);
    # this below is just to get the typeof tmp, to then avoid a weird unionall error which should hopefully go away when updating Julia at some point
    tmp = discrete_system(ωs[1], sourceradial, material;
        basis_order = 0, basis_field_order = 0, rtol = 1.0, maxevals = 4
    );
    ST = typeof(tmp);
    discrete_fields = ST[
        discrete_system(ωs[i], sourceradial, material;
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            rtol = rtol, maxevals = maxevals
        )
    for i in eachindex(ωs)];

    xs = [ radial_to_cartesian_coordinates([r,0.0,0.0]) for r in rs];

    eff_scats_azi = [s.(xs) for s in scat_fields_azi];
    eff_scats_radial = [s.(xs) for s in scat_fields_radial];
    discrete_scats = [d.coefficient_field.(xs) for d in discrete_fields];

    # assuming radial or azimuthal symmetry should lead to exactly the same fields
    errors = [
        norm.(eff_scats_radial[i] - eff_scats_azi[i]) ./ norm.(eff_scats_radial[i])
    for i in eachindex(ωs)];

    @test maximum(maximum.(errors)) < 1e-12

    # the discrete method does have a difference
    errors = [
        norm.(eff_scats_radial[i] - discrete_scats[i]) ./ norm.(eff_scats_radial[i])
    for i in eachindex(ωs)];

    @test minimum(mean.(errors)) < 1e-4
    @test maximum(mean.(errors)) < 1e-3
    @test maximum(maximum.(errors)) < 2e-3

    mat_coefs_radial = material_scattering_coefficients.(wavemodes_radial);

    mat_coefs_disc_radial = material_scattering_coefficients.(discrete_fields;
        rtol = rtol,
        maxevals = maxevals
    );

    errors = norm.(mat_coefs_disc_radial - mat_coefs_radial) ./ norm.(mat_coefs_radial)

    @test errors[1] < 1e-5
    @test errors[2] < 2e-3

# Test the reduced radial discrete method

    pair_corr_inf(z) = hole_correction_pair_correlation([0.0,0.0,0.0],s1, [0.0,0.0,z],s1)

    polynomial_order = 15
    pair_corr_inf_smooth = smooth_pair_corr_distance(
        pair_corr_inf, a12;
        smoothing = 0.3, max_distance = 2R,
        polynomial_order = polynomial_order
    )

    # using Plots
    # zs = 0.0:0.01:(2R)
    # plot(pair_corr_inf_smooth,zs)
    # plot!(pair_corr_inf,zs, linestyle=:dash)

    gls_radial = gls_pair_radial_fun(pair_corr_inf_smooth, a12;
        sigma_approximation = false,
        polynomial_order = polynomial_order
    )

    discrete_field_radials = [
        discrete_system_radial(ωs[i], sourceradial, material, Symmetry(sourceradial,material);
            basis_order = basis_orders[i],
            basis_field_order = basis_field_orders[i],
            # polynomial_order = 15,
            # pair_corr_distance = pair_corr_inf
            gls_pair_radial = gls_radial
        )
    for i in eachindex(ωs)];

    discrete_rad_scats = [
        d.coefficient_field.(xs)
    for d in discrete_field_radials];

    errors = [
        norm.(discrete_rad_scats[i] - discrete_scats[i]) ./ norm.(discrete_scats[i])
    for i in eachindex(ωs)];

    @test mean(mean.(errors)) < 1e-3
    @test maximum(mean.(errors)) < 1e-3
    @test maximum(maximum.(errors)) < 4e-3

    mat_coefs_disc_radial2 = [
        material_scattering_coefficients(discrete_field_radials[i];
            rtol = rtol,
            maxevals = maxevals
        )
    for i in eachindex(ωs)];

    errors = norm.(mat_coefs_disc_radial2 - mat_coefs_disc_radial) ./ norm.(mat_coefs_disc_radial)
    @test maximum(errors) < 1e-3

# Calculate low frequency scattering
    # Linc = basis_orders[1] + basis_field_orders[1]
    # source_coefficients = regular_spherical_coefficients(sourceradial)(Linc,zeros(3),ωs[1]);
    #
    # material_low = Material(
    #     Sphere(outer_radius(material.shape) - outer_radius(s1)),
    #     species
    # );
    # #
    # effective_sphere = Particle(eff_medium, material_low.shape);
    # Tmat = MultipleScattering.t_matrix(effective_sphere, medium, ωs[1], Linc);
    # scat_coef_low = Tmat * source_coefficients;
end


text:
"""
    $(SIGNATURES)

Perform the nonlinear numerical operations to approximate the convolution with a particular user defined likelihood function (conditional), which as been prepared in the `frl` object.  This function uses root finding to enforce a non-linear function constraint.

Notes:
- remember this is a deepcopy of original sfidx, since we are generating a proposal distribution and not directly replacing the existing variable belief estimate

Future work:
- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
- improve handling of n and particleidx, especially considering future multithreading support

"""
function approxConvOnElements!(ccwl::CommonConvWrapper{T},
                               elements::Union{Vector{Int}, UnitRange{Int}}, ::Type{MultiThreaded}  ) where {T <: Union{FunctorPairwise, FunctorPairwiseMinimize}}
  Threads.@threads for n in elements
    # ccwl.thrid_ = Threads.threadid()
    ccwl.cpt[Threads.threadid()].particleidx = n
    # ccall(:jl_, Nothing, (Any,), "starting loop, thrid_=$(Threads.threadid()), partidx=$(ccwl.cpt[Threads.threadid()].particleidx)")
    numericRootGenericRandomizedFnc!( ccwl )
  end
  nothing
end

"""
    $(SIGNATURES)

Perform the nonlinear numerical operations to approximate the convolution with a particular user defined likelihood function (conditional), which as been prepared in the `frl` object.  This function uses root finding to enforce a non-linear function constraint.

Notes:
- remember this is a deepcopy of original sfidx, since we are generating a proposal distribution and not directly replacing the existing variable belief estimate

Future work:
- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
- improve handling of n and particleidx, especially considering future multithreading support

"""
function approxConvOnElements!(ccwl::CommonConvWrapper{T},
                               elements::Union{Vector{Int}, UnitRange{Int}}, ::Type{SingleThreaded}) where {T <: Union{FunctorPairwise, FunctorPairwiseMinimize}}
  #
  for n in elements
    ccwl.cpt[Threads.threadid()].particleidx = n
    numericRootGenericRandomizedFnc!( ccwl )
  end
  nothing
end

"""
    $(SIGNATURES)

Perform the nonlinear numerical operations to approximate the convolution with a particular user defined likelihood function (conditional), which as been prepared in the `frl` object.  This function uses root finding to enforce a non-linear function constraint.

Notes:
- remember this is a deepcopy of original sfidx, since we are generating a proposal distribution and not directly replacing the existing variable belief estimate

Future work:
- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
- improve handling of n and particleidx, especially considering future multithreading support

"""
function approxConvOnElements!(ccwl::CommonConvWrapper{T},
                               elements::Union{Vector{Int}, UnitRange{Int}} )  where {T <: Union{FunctorPairwise, FunctorPairwiseMinimize}}
  #
  approxConvOnElements!(ccwl, elements, ccwl.threadmodel)
end

"""
    $(SIGNATURES)

Prepare a common functor computation object `prepareCommonConvWrapper{T}` containing the user factor functor along with additional variables and information using during approximate convolution computations.
"""
function prepareCommonConvWrapper!(ccwl::CommonConvWrapper{T},
                                   Xi::Vector{DFGVariable},
                                   solvefor::Symbol,
                                   N::Int  ) where {T <: FunctorInferenceType}
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, N, solvefor)
  # should be selecting for the correct multihypothesis mode here with `gwp.params=ARR[??]`
  ccwl.params = ARR
  ccwl.measurement = getSample(ccwl.usrfnc!, maxlen) # ccwl.samplerfnc
  if ccwl.specialzDim
    ccwl.zDim = ccwl.usrfnc!.zDim[sfidx]
  else
    ccwl.zDim = size(ccwl.measurement[1],1) # TODO -- zDim aspect needs to be reviewed
  end
  ccwl.varidx = sfidx

  ccwl.xDim = size(ccwl.params[sfidx],1)
  # ccwl.xDim = size(ccwl.cpt[1].X,1)
  # info("what? sfidx=$(sfidx), ccwl.xDim = size(ccwl.params[sfidx]) = $(ccwl.xDim), size=$(size(ccwl.params[sfidx]))")
  for i in 1:Threads.nthreads()
    ccwl.cpt[i].X = ccwl.params[sfidx]
    ccwl.cpt[i].p = collect(1:size(ccwl.cpt[i].X,1))
    ccwl.cpt[i].Y = zeros(ccwl.xDim)
    ccwl.cpt[i].res = zeros(ccwl.xDim) # used in ccw functor for FunctorPairwiseMinimize
  end

  return sfidx, maxlen
end

"""
    $(SIGNATURES)

Common function to compute across a single user defined multi-hypothesis ambiguity per factor.  This function dispatches both `FunctorPairwise` and `FunctorPairwiseMinimize` factors.
"""
function computeAcrossHypothesis!(ccwl::CommonConvWrapper{T},
                                  allelements,
                                  activehypo,
                                  certainidx,
                                  sfidx) where {T <:Union{FunctorPairwise, FunctorPairwiseMinimize}}
  count = 0
  # TODO remove assert once all GenericWrapParam has been removed
  # @assert norm(ccwl.certainhypo - certainidx) < 1e-6
  for (mhidx, vars) in activehypo
    count += 1
    if sfidx in certainidx || mhidx in certainidx || mhidx == sfidx
      # hypo case mhidx, sfidx = $mhidx, $sfidx
      for i in 1:Threads.nthreads()  ccwl.cpt[i].activehypo = vars; end
      approxConvOnElements!(ccwl, allelements[count])
    # elseif mhidx == sfidx
    #   # multihypo, do conv case, mhidx == sfidx
    #   ah = sort(union([sfidx;], certainidx))
    #   @assert norm(ah - vars) < 1e-10
    #   for i in 1:Threads.nthreads()  ccwl.cpt[i].activehypo = ah; end
    #   approxConvOnElements!(ccwl, allelements[count])
    elseif mhidx != sfidx   # sfidx in uncertnidx
      # multihypo, take other value case
      # sfidx=2, mhidx=3:  2 should take a value from 3
      # sfidx=3, mhidx=2:  3 should take a value from 2
         # DEBUG sfidx=2, mhidx=1 -- bad when do something like multihypo=[0.5;0.5] -- issue 424
      ccwl.params[sfidx][:,allelements[count]] = view(ccwl.params[mhidx],:,allelements[count])
    else
      error("computeAcrossHypothesis -- not dealing with multi-hypothesis case correctly")
    end
  end
  nothing
end

"""
    $(SIGNATURES)

Prepare data required for null hypothesis cases during convolution.
"""
function assembleNullHypothesis(ccwl::CommonConvWrapper{T},
                                maxlen::Int,
                                spreadfactor::Float64 ) where {T}
  #
  nhc = rand(ccwl.usrfnc!.nullhypothesis, maxlen) .- 1
  val = ccwl.params[ccwl.varidx]
  d = size(val,1)
  var = Statistics.var(val,dims=2) .+ 1e-3
  ENT = Distributions.MvNormal(zeros(d), spreadfactor*Matrix(Diagonal(var[:])))
  allelements = 1:maxlen
  return allelements, nhc, ENT
end

"""
    $(SIGNATURES)

Do true and null hypothesis computations based on data structures prepared earlier -- specific to `FunctorPairwiseNH`.  This function will be merged into a standard case for `FunctorPairwise/Minimize` in the future.
"""
function computeAcrossNullHypothesis!(ccwl::CommonConvWrapper{T},
                                      allelements,
                                      nhc,
                                      ENT  ) where {T <: FunctorPairwiseNH}
  #
  # TODO --  Threads.@threads see area4 branch
  for n in allelements
    # ccwl.gwp(x, res)
    if nhc[n] != 0
      ccwl.cpt[Threads.threadid()].particleidx = n
      numericRootGenericRandomizedFnc!( ccwl )
    else
      ccwl.params[ccwl.varidx][:,n] += rand(ENT)
    end
  end
  nothing
end

"""
    $(SIGNATURES)

Multiple dispatch wrapper for `<:FunctorPairwise` types, to prepare and execute the general approximate convolution with user defined factor residual functions.  This method also supports multihypothesis operations as one mechanism to introduce new modality into the proposal beliefs.

Planned changes will fold null hypothesis in as a standard feature and no longer appear as a separate `InferenceType`.
"""
function evalPotentialSpecific(Xi::Vector{DFGVariable},
                               ccwl::CommonConvWrapper{T},
                               solvefor::Symbol,
                               measurement::Tuple=(zeros(0,100),);
                               N::Int=size(measurement[1],2),
                               spreadfactor::Float64=10.0,
                               dbg::Bool=false ) where {T <: FunctorPairwiseNH}
  #
  @warn "FunctorPairwiseNH will be deprecated in favor of common `nullhypo=` interface."
  # TODO -- could be constructed and maintained at addFactor! time
  sfidx, maxlen = prepareCommonConvWrapper!(ccwl, Xi, solvefor, N)
  # prepare nullhypothesis
  allelements, nhc, ENT = assembleNullHypothesis(ccwl, maxlen, spreadfactor)

  # Compute across the true or null hypothesis
  computeAcrossNullHypothesis!(ccwl, allelements, nhc, ENT )

  return ccwl.params[ccwl.varidx]
end

function evalPotentialSpecific(Xi::Vector{DFGVariable},
                               ccwl::CommonConvWrapper{T},
                               solvefor::Symbol,
                               measurement::Tuple=(zeros(0,100),);
                               N::Int=size(measurement[1],2),
                               dbg::Bool=false  ) where {T <: Union{FunctorPairwise, FunctorPairwiseMinimize}}
  #

  # Prep computation variables
  sfidx, maxlen = prepareCommonConvWrapper!(ccwl, Xi, solvefor, N)
  # check for user desired measurement values
  if 0 < size(measurement[1],1)
    ccwl.measurement = measurement
  end

  # Check which variables have been initialized
  isinit = map(x->isInitialized(x), Xi)

  # assemble how hypotheses should be computed
  _, allelements, activehypo, mhidx = assembleHypothesesElements!(ccwl.hypotheses, maxlen, sfidx, length(Xi), isinit)
  certainidx = ccwl.certainhypo

  # perform the numeric solutions on the indicated elements
  # error("ccwl.xDim=$(ccwl.xDim)")
  computeAcrossHypothesis!(ccwl, allelements, activehypo, certainidx, sfidx)

  return ccwl.params[ccwl.varidx]
end

function evalPotentialSpecific(Xi::Vector{DFGVariable},
                               ccwl::CommonConvWrapper{T},
                               solvefor::Symbol,
                               measurement::Tuple=(zeros(0,0),);
                               N::Int=size(measurement[1],2),
                               dbg::Bool=false ) where {T <: FunctorSingleton}
  #
  fnc = ccwl.usrfnc!

  nn = (N <= 0 ? size(getVal(Xi[1]),2) : N)
  ccwl.measurement = 0 < size(measurement[1],1) ? measurement : getSample(ccwl.usrfnc!, nn)
  if !ccwl.partial
    return ccwl.measurement[1]
  else
    val = deepcopy(getVal(Xi[1]))
    i = 0
    for dimnum in fnc.partial
      i += 1
      val[dimnum,:] = ccwl.measurement[1][i,:]
    end
    return val
  end
end

function evalPotentialSpecific(Xi::Vector{DFGVariable},
                               ccwl::CommonConvWrapper{T},
                               solvefor::Symbol,
                               measurement::Tuple=(zeros(0,100),);
                               N::Int=size(measurement[1],2),
                               spreadfactor::Float64=10.0,
                               dbg::Bool=false ) where {T <: FunctorSingletonNH}
  #
  @warn "FunctorSingletonNH will be deprecated in favor of common `nullhypo=` interface."
  fnc = ccwl.usrfnc!

  val = getVal(Xi[1])
  d = size(val,1)
  var = Statistics.var(val, dims=2) .+ 1e-3

  # determine amount share of null hypothesis particles
  ccwl.measurement = getSample(ccwl.usrfnc!, N)
  # values of 0 imply null hypothesis
  # ccwl.usrfnc!.nullhypothesis::Distributions.Categorical
  nhc = rand(ccwl.usrfnc!.nullhypothesis, N) .- 1

  # TODO -- not valid for manifold
  # TODO bad memory management
  ENT = Distributions.MvNormal(zeros(d), spreadfactor*Matrix(Diagonal(var[:])) )

  for i in 1:N
    if nhc[i] == 0
      ccwl.measurement[1][:,i] = val[:,i] + rand(ENT)  # TODO use view and inplace add operation
    end
  end
  # TODO -- returning to memory location inside
  return ccwl.measurement[1]
end

"""
    $(SIGNATURES)

Single entry point for evaluating factors from factor graph, using multiple dispatch to locate the correct `evalPotentialSpecific` function.
"""
function evalFactor2(dfg::AbstractDFG,
                     fct::DFGFactor,
                     solvefor::Symbol,
                     measurement::Tuple=(zeros(0,100),);
                     N::Int=size(measurement[1],2),
                     dbg::Bool=false  )
  #

  ccw = solverData(fct).fnc
  # TODO -- this build up of Xi is excessive and could happen at addFactor time
  Xi = DFGVariable[]
  count = 0
  # TODO replace fncargvID with neighbors
  variablelist = Vector{Symbol}(undef, length(solverData(fct).fncargvID))
  for id in solverData(fct).fncargvID
    count += 1
    xi = DFG.getVariable(dfg, id)
    push!(Xi, xi )

    # TODO do only once at construction time -- staring it here to be sure the code is calling factors correctly
    variablelist[count] = xi.label

    # TODO bad way to search for `solvefor`
    if xi.label == solvefor
      for i in 1:Threads.nthreads()
        ccw.cpt[i].factormetadata.solvefor = xi.label
      end
    end
  end
  for i in 1:Threads.nthreads()
    ccw.cpt[i].factormetadata.variablelist = variablelist
  end
  return evalPotentialSpecific(Xi, ccw, solvefor, measurement, N=N, dbg=dbg)
end

# import IncrementalInference: evalFactor2, approxConv
"""
    $(SIGNATURES)

Draw samples from the approximate convolution of `towards` symbol using factor `fct` relative to the other variables.  In addition the `api` can be adjusted to recover the data from elsewhere (likely to be replaced/removed in the future).
"""
function approxConv(dfg::AbstractDFG,
                    fc::DFGFactor,
                    towards::Symbol,
                    measurement::Tuple=(zeros(0,0),);
                    N::Int=size(measurement[1],2) )
  #
  v1 = getVariable(dfg, towards)
  N = N == 0 ? getNumPts(v1) : N
  return evalFactor2(dfg, fc, v1.label, measurement, N=N)
end
function approxConv(dfg::AbstractDFG,
                    fct::Symbol,
                    towards::Symbol,
                    measurement::Tuple=(zeros(0,0),);
                    N::Int=size(measurement[1],2) )
  #
  fc = getFactor(dfg, fct)
  return approxConv(dfg, fc, towards, measurement, N=N)
end

function approxConvBinary(arr::Array{Float64,2},
                          meas::FunctorInferenceType,
                          outdims::Int,
                          measurement::Tuple=(zeros(0,size(arr,2)),);
                          varidx::Int=2,
                          N::Int=size(arr,2)  )
  #
  # N = N == 0 ? size(arr,2) : N
  pts = zeros(outdims,N);
  t = Array{Array{Float64,2},1}()
  push!(t,arr)
  push!(t,pts)

  measurement = size(measurement[1],2) == 0 ? getSample(meas, N) : measurement

  zDim = size(measurement[1],1)
  ccw = CommonConvWrapper(meas, t[varidx], zDim, t, varidx=varidx, measurement=measurement)  # N=> size(measurement[1],2)

  for n in 1:N
    ccw.cpt[Threads.threadid()].particleidx = n
    numericRootGenericRandomizedFnc!( ccw )
  end
  return pts
end

"""
    $(SIGNATURES)

Compute proposal belief on `vertid` through `fct` representing some constraint in factor graph.
Always full dimension variable node -- partial constraints will only influence subset of variable dimensions.
The remaining dimensions will keep pre-existing variable values.

Notes
- fulldim is true when "rank-deficient" -- TODO swap to false (or even float)
"""
function findRelatedFromPotential(dfg::G,
                                  fct::DFGFactor,
                                  varid::Symbol,
                                  N::Int,
                                  dbg::Bool=false  )::Tuple{BallTreeDensity,Float64} where G <: AbstractDFG
  # assuming it is properly initialized TODO
  ptsbw = evalFactor2(dfg, fct, varid, N=N, dbg=dbg);
  # determine if evaluation is "dimension-deficient"

  # solvable dimension
  inferdim = getFactorSolvableDim(dfg, fct, varid)
  # zdim = getFactorDim(fct)
  # vdim = getVariableDim(DFG.getVariable(dfg, varid))

  # TODO -- better to upsample before the projection
  Ndim = size(ptsbw,1)
  Npoints = size(ptsbw,2)
  # Assume we only have large particle population sizes, thanks to addNode!
  manis = getManifolds(dfg, varid)
  # manis = getSofttype(DFG.getVariable(dfg, varid)).manifolds # older
  p = AMP.manikde!(ptsbw, manis)
  if Npoints != N # this is where we control the overall particle set size
      p = resample(p,N)
  end
  return p, inferdim
end

#


text:
function extractinfofromconfig( config_path::String,
    molecule_names::Vector{String}) where T <: Real

    N_compounds = length(molecule_names)

    # read.
    file_strings = readlines(config_path)

    cs_delta_group = Vector{Vector{Vector{Float64}}}(undef, N_compounds)
    #λ_group_labels = Vector{Vector{Vector{Int}}}(undef, N_compounds)

    for n = 1:N_compounds

        cs_delta_group[n] = extractmoleculeinfofromconfig(file_strings,
                                    molecule_names[n], config_path)
    end

    return cs_delta_group#, λ_group_labels
end

function condenseΔcsconfig(cs_delta_group::Vector{Vector{Vector{T}}}) where T
    #
    out = Vector{Vector{T}}(undef, length(cs_delta_group))

    for n = 1:length(cs_delta_group)
        out[n] = Vector{T}(undef, length(cs_delta_group[n]))

        for i = 1:length(cs_delta_group[n])
            out[n][i] = cs_delta_group[n][i][1]
        end
    end

    return out
end

function convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T where T <: Real

    return (x-a)*(d-c)/(b-a)+c
end

function displaypartitionsizes(A::NMRSpectraSimulator.FIDModelType{T,SST}) where {T,SST}

    return collect( length(A.part_inds_compound[k]) for k = 1:length(A.part_inds_compound) )
end

function forcesymmetric(A::Matrix{T})::Matrix{T} where T <: Real
    return (A+A')./2
end

# average Δc vector for each partition element.
function viewaverageΔc(A::NMRSpectraSimulator.FIDModelType{T,SST}) where {T,SST}

    N_spin_groups = length(A.part_inds_compound)
    out = Vector{Vector{Vector{T}}}(undef, N_spin_groups)

    for i = 1:N_spin_groups

        part_size = length(A.part_inds_compound[i])

        out[i] = Vector{Vector{T}}(undef, part_size)

        for (k,inds) in enumerate(A.part_inds_compound[i])

            list_of_Δc_m = A.Δc_m_compound[i][inds]

            # insert error handing here later in case list_of_Δc_m is empty.
            N_spins = length(list_of_Δc_m[1])
            out[i][k] = zeros(T, N_spins)

            N_components = length(list_of_Δc_m)
            for j = 1:N_components
                #
                out[i][k] += list_of_Δc_m[j]
            end
            out[i][k] = out[i][k] ./ N_components

            #out[i][k] = Statistics.mean(tmp)
        end
    end

    return out
end

function mergeinds(inds_set::Vector{Vector{Int}}, merge_set) where T

    # check.
    for l = 1:length(merge_set)
        @assert 0 < l <= length(inds_set)
    end

    # merge.
    inds = Vector{Int}(undef, 0)
    for i in merge_set
        push!(inds, inds_set[i]...)
    end

    unique!(inds)
    return inds
end

function runNLopt!(  opt,
    p0::Vector{T},
    obj_func,
    grad_func,
    p_lbs,
    p_ubs;
    max_iters = 10000,
    xtol_rel = 1e-12,
    ftol_rel = 1e-12,
    maxtime = Inf) where T

    @assert length(p0) == length(p_lbs) == length(p_ubs)

    opt.maxeval = max_iters
    opt.lower_bounds = p_lbs
    opt.upper_bounds = p_ubs
    opt.xtol_rel = xtol_rel
    opt.ftol_rel = ftol_rel
    opt.maxtime = maxtime

    opt.min_objective = (xx, gg)->genericNLoptcostfunc!(xx, gg, obj_func, grad_func)

    # optimize.
    (minf, minx, ret) = NLopt.optimize(opt, p0)

    N_evals = opt.numevals

    return minf, minx, ret, N_evals
end

function genericNLoptcostfunc!(x::Vector{T}, df_x, f, df)::T where T <: Real

    #
    if length(df_x) > 0
        df_x[:] = df(x)
    end

    return f(x)
end

###### tmp.
function runalignment2(shift_constants,
    U_rad_cost,
    y_cost::Vector{Complex{T}},
    LS_inds,
    Es::Vector{NMRSpectraSimulator.κCompoundFIDType{T, SST}},
    As::Vector{NMRSpectraSimulator.CompoundFIDType{T, SST}},
    fs::T,
    SW::T,
    λ_lbs, λ_ubs, κs_β_DOFs, κs_β_orderings;
    optim_algorithm::Symbol = :GN_ESCH,
    max_iters = 5000,
    xtol_rel = 1e-7,
    ftol_rel = 1e-12,
    maxtime = Inf,
    w = ones(T, length(Es)),
    κ_lb_default = 0.2,
    κ_ub_default = 50.0,
    λ_each_lb = 0.7,
    λ_each_ub = 5.0) where {T <: Real, SST}

    @assert length(U_rad_cost) == length(y_cost) == length(LS_inds)

    #
    N_d = sum( getNd(As[n]) for n = 1:length(As) )
    N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    N_λ = sum( getNλ(As[n]) for n = 1:length(As) )
    # N_vars = N_d + N_β + N_λ

    #### initial values.
    shift_lb = -ones(T, N_d)
    shift_ub = ones(T, N_d)
    shift_initial = zeros(T, N_d)

    β_lb = ones(T, N_β) .* (-π)
    β_ub = ones(T, N_β) .* (π)
    β_initial = zeros(T, N_β)

    λ_lb = λ_each_lb .* ones(T, N_λ)
    λ_ub = λ_each_ub .* ones(T, N_λ)
    λ_initial = ones(T, N_λ)

    ### initial phase estimate. Constant phase model.
    #β_est = estimateconstantphase(U_cost, y_cost)
    #βS_base = initializeFIDparameter(αS, β_est)
    #fill!(β_initial, β_est)
    #fill!(β_initial, 0.0)

    ### set up.
    p_lb = [ shift_lb; β_lb; λ_lb ]
    p_ub = [ shift_ub; β_ub; λ_ub ]
    p_initial = [shift_initial; β_initial; λ_initial]
    println("p_lb = ", p_lb)
    println("p_ub = ", p_ub)
    println("p_initial = ", p_initial)
    println("length(p_initial) = ", length(p_initial))
    println()

    q, updatedfunc, updateβfunc, updateλfunc, updateκfunc, #updatewfunc,
    E_BLS, κ_BLS, b_BLS, getshiftfunc, getβfunc, getλfunc,
    N_vars_set = setupcostcLshiftLS(Es, As, fs, SW,
        LS_inds, U_rad_cost, y_cost, shift_constants;
        w = w, κ_lb_default = κ_lb_default, κ_ub_default = κ_ub_default)
    #updatewfunc(p_initial)
    updateκfunc(p_initial)
    parseκ!(Es, κ_BLS)

    ### cost function.
    obj_func = pp->costcLshift(U_rad_cost, y_cost,
    updatedfunc, updateβfunc, updateλfunc, updateκfunc, pp, Es, E_BLS, κ_BLS, b_BLS, q)

    grad_func = xx->FiniteDiff.finite_difference_gradient(obj_func, xx)

    # force eval for debugging.
    println("Timing: obj_func")
    @time cost_initial = obj_func(p_initial)
    cost_initial = obj_func(p_initial)
    println("obj_func(p_initial) = ", cost_initial)
    println()

    opt = NLopt.Opt(optim_algorithm, length(p_initial)) # global evolutionary.

    minf, minx, ret, numevals = runNLopt!(  opt,
        p_initial,
        obj_func,
        grad_func,
        p_lb,
        p_ub;
        max_iters = max_iters,
        xtol_rel = xtol_rel,
        ftol_rel = ftol_rel,
        maxtime = maxtime)

    println("numevals = ", numevals)
    println("norm(p_initial-minx) = ", norm(p_initial-minx))

    println("Timing: obj_func")
    @time cost_final = obj_func(minx)
    println("obj_func(minx) = ", cost_final)
    println()

    return minx, q, κ_BLS, getshiftfunc, getβfunc, getλfunc,
        obj_func, N_vars_set, p_initial
end


text:
"""
    HypertextLiteral

This library provides a `@htl()` macro which implements interopolation
that is aware of hypertext escape context.
"""
module HypertextLiteral
export @htl

include("primitives.jl") # Wrap, Unwrap, EscapeProxy
include("macro.jl")      # @htl macro and `Result` object
include("notation.jl")   # @htl_str non-standard string literal
include("convert.jl")    # runtime conversion of objects
include("lexer.jl")      # interpolate string to macro expression
include("rewrite.jl")    # macro optimizations called by interpolate

end