text:
# ---
# title: 1299. Replace Elements with Greatest Element on Right Side
# id: problem1299
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/replace-elements-with-greatest-element-on-right-side/description/>
# hidden: true
# ---
#
# Given an array `arr`, replace every element in that array with the greatest
# element among the elements to its right, and replace the last element with
# `-1`.
#
# After doing so, return the array.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [17,18,5,4,6,1]
#     Output: [18,6,6,6,1,-1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^4`
#   * `1 <= arr[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    NestedSamplers.Proposals

This module contains the different algorithms for proposing new points within a bounding volume in unit space.

The available implementations are
* [`Proposals.Rejection`](@ref) - samples uniformly within the bounding volume
* [`Proposals.RWalk`](@ref) - random walks to a new point given an existing one
* [`Proposals.RStagger`](@ref) - random staggering away to a new point given an existing one
* [`Proposals.Slice`](@ref) - slicing away to a new point given an existing one
* [`Proposals.RSlice`](@ref) - random slicing away to a new point given an existing one
"""
module Proposals

using ..NestedSamplers: prior_transform_and_loglikelihood
using ..Bounds

using Random
using LinearAlgebra
using Parameters

export AbstractProposal

"""
    NestedSamplers.AbstractProposal

The abstract type for live point proposal algorithms.

# Interface

Each `AbstractProposal` must have this function,
```julia
(::AbstractProposal)(::AbstractRNG, point, loglstar, bounds, loglikelihood, prior_transform)
```
which, given the input `point` with loglikelihood `loglstar` inside a `bounds`, returns a new point in unit space, prior space, the loglikelihood, and the number of function calls.
"""
abstract type AbstractProposal end

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

# Helper for checking unit-space bounds
unitcheck(us) = all(u -> 0 < u < 1, us)

"""
    Proposals.Rejection(;maxiter=100_000)

Propose a new live point by uniformly sampling within the bounding volume and rejecting samples that do not meet the likelihood constraints. This follows the original nested sampling algorithm proposed in Skilling (2004)[^1]

[^1]: John Skilling, 2004, AIP 735, 395  ["Nested Sampling"](https://aip.scitation.org/doi/abs/10.1063/1.1835238)

## Parameters
- `maxiter` is the maximum number of samples that can be rejected before giving up and throwing an error.
"""
Base.@kwdef struct Rejection <: AbstractProposal
    maxiter::Int = 100_000
end

@deprecate Uniform() Rejection()

function (prop::Rejection)(
    rng::AbstractRNG,
    point::AbstractVector,
    logl_star,
    bounds::AbstractBoundingSpace,
    model
)

    ncall = 0
    for _ in 1:prop.maxiter
        u = rand(rng, bounds)
        unitcheck(u) || continue
        v, logl = prior_transform_and_loglikelihood(model, u)
        ncall += 1
        logl ≥ logl_star && return u, v, logl, ncall
    end
    throw(ErrorException("Couldn't generate a proper point after $(prop.maxiter) attempts including $ncall likelihood calls. Bounds=$bounds, logl_star=$logl_star."))
end

Base.show(io::IO, p::Rejection) = print(io, "NestedSamplers.Proposals.Rejection")

"""
    Proposals.RWalk(;ratio=0.5, walks=25, scale=1)

Propose a new live point by random walking away from an existing live point. This follows the algorithm outlined in Skilling (2006).[^5]

[^5]: Skilling, 2006, Bayesian Anal. 1(4), ["Nested sampling for general Bayesian computation"](https://projecteuclid.org/journals/bayesian-analysis/volume-1/issue-4/Nested-sampling-for-general-Bayesian-computation/10.1214/06-BA127.full)

## Parameters
- `ratio` is the target acceptance ratio
- `walks` is the minimum number of steps to take
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct RWalk <: AbstractProposal
    ratio = 0.5
    walks = 25
    scale = 1.0

    @assert 1 / walks ≤ ratio ≤ 1 "Target acceptance ratio must be between 1/`walks` and 1"
    @assert walks > 1 "Number of steps must be greater than 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::RWalk)(
    rng::AbstractRNG,
    point::AbstractVector,
    logl_star,
    bounds::AbstractBoundingSpace,
    model;
    kwargs...
)
    # setup
    n = length(point)
    scale_init = prop.scale
    accept = reject = fail = nfail = nc = ncall = 0
    local du, u_prop, logl_prop, u, v, logl

    while nc < prop.walks || iszero(accept)
        # get proposed point
        while true
            # check scale factor to avoid over-shrinking
            prop.scale < 1e-5 * scale_init && error("Random walk sampling appears to be stuck.")
            # transform to proposal distribution
            du = randoffset(rng, bounds)
            u_prop = @. point + prop.scale * du
            # inside unit-cube
            unitcheck(u_prop) && break

            fail += 1
            nfail += 1
            # check if stuck generating bad numbers
            if fail > 100 * prop.walks
                @warn "Random number generation appears extremely inefficient. Adjusting the scale-factor accordingly"
                fail = 0
                prop.scale *= exp(-1/n)
            end
        end
        # check proposed point
        v_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop)
        if logl_prop ≥ logl_star
            u = u_prop
            v = v_prop
            logl = logl_prop
            accept += 1
        else
            reject += 1
        end
        nc += 1
        ncall += 1

        # check if stuck generating bad points
        if nc > 50 * prop.walks
            @warn "Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly"
            prop.scale *= exp(-1/n)
            nc = accept = reject = 0
        end
    end

    # update proposal scale using acceptance ratio
    update_scale!(prop, accept, reject, n)

    return u, v, logl, ncall
end

# update proposal scale using target acceptance ratio
function update_scale!(prop, accept, reject, n)
    ratio = accept / (accept + reject)
    norm = max(prop.ratio, 1 - prop.ratio) * n
    scale = prop.scale * exp((ratio - prop.ratio) / norm)
    prop.scale = min(scale, sqrt(n))
    return prop
end

"""
    Proposals.RStagger(;ratio=0.5, walks=25, scale=1)

Propose a new live point by random staggering away from an existing live point.
This differs from the random walk proposal in that the step size here is exponentially adjusted
to reach a target acceptance rate _during_ each proposal, in addition to _between_
proposals. This follows the algorithm outlined in Skilling (2006).[^5]

[^5]: Skilling, 2006, Bayesian Anal. 1(4), ["Nested sampling for general Bayesian computation"](https://projecteuclid.org/journals/bayesian-analysis/volume-1/issue-4/Nested-sampling-for-general-Bayesian-computation/10.1214/06-BA127.full)

## Parameters
- `ratio` is the target acceptance ratio
- `walks` is the minimum number of steps to take
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct RStagger <: AbstractProposal
    ratio = 0.5
    walks = 25
    scale = 1.0

    @assert 1 / walks ≤ ratio ≤ 1 "Target acceptance ratio must be between 1/`walks` and 1"
    @assert walks > 1 "Number of steps must be greater than 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::RStagger)(
    rng::AbstractRNG,
    point::AbstractVector,
    logl_star,
    bounds::AbstractBoundingSpace,
    model;
    kwargs...
)
    #setup
    n = length(point)
    scale_init = prop.scale
    accept = reject = fail = nfail = nc = ncall = 0
    stagger = 1
    local du, u_prop, logl_prop, u, v, logl

    while nc < prop.walks || iszero(accept)
        # get proposed point
        while true
            # check scale factor to avoid over-shrinking
            prop.scale < 1e-5 * scale_init && error("Random walk sampling appears to be stuck.")
            # transform to proposal distribution
            du = randoffset(rng, bounds)
            u_prop = @. point + prop.scale * stagger * du
            # inside unit-cube
            unitcheck(u_prop) && break

            fail += 1
            nfail += 1
            # check if stuck generating bad numbers
            if fail > 100 * prop.walks
                @warn "Random number generation appears extremely inefficient. Adjusting the scale-factor accordingly"
                fail = 0
                prop.scale *= exp(-1/n)
            end
        end
        # check proposed point
        v_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop)
        if logl_prop ≥ logl_star
            u = u_prop
            v = v_prop
            logl = logl_prop
            accept += 1
        else
            reject += 1
        end
        nc += 1
        ncall += 1

        # adjust _stagger_ to target an acceptance ratio of `prop.ratio`
        ratio = accept / (accept + reject)
        if ratio > prop.ratio
            stagger *= exp(1 / accept)
        elseif ratio < prop.ratio
            stagger /= exp(1 / reject)
        end

        # check if stuck generating bad points
        if nc > 50 * prop.walks
            @warn "Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly"
            prop.scale *= exp(-1 / n)
            nc = accept = reject = 0
        end
    end

    # update proposal scale using acceptance ratio
    update_scale!(prop, accept, reject, n)

    return u, v, logl, ncall
end

"""
    Proposals.Slice(;slices=5, scale=1)

Propose a new live point by a series of random slices away from an existing live point.
This is a standard _Gibbs-like_ implementation where a single multivariate slice is a combination of `slices` univariate slices through each axis. This follows the algorithm outlined in Neal (2003).[^6]

[^6]: Neal, 2003, Ann. Statist. 31(3), ["Slice Sampling"](https://projecteuclid.org/journals/annals-of-statistics/volume-31/issue-3/Slice-sampling/10.1214/aos/1056562461.full)

## Parameters
- `slices` is the minimum number of slices
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct Slice <: AbstractProposal
    slices = 5
    scale = 1.0

    @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::Slice)(
    rng::AbstractRNG,
    point::AbstractVector,
    logl_star,
    bounds::AbstractBoundingSpace,
    model;
    kwargs...
)
    # setup
    n = length(point)
    nc = nexpand = ncontract = 0
    local u, v, logl

    # modifying axes and computing lengths
    axes = Bounds.axes(bounds)
    axes = prop.scale .* axes'
    # slice sampling loop
    for it in 1:prop.slices

        # shuffle axis update order
        idxs = shuffle!(rng, collect(Base.axes(axes, 1)))

        # slice sample along a random direction
        for idx in idxs

            # select axis
            axis = axes[idx, :]

            u, v, logl, nc, nexpand, ncontract = sample_slice(
                rng, axis, point, logl_star,
                model,
                nc, nexpand, ncontract
            )
        end # end of slice sample along a random direction
    end # end of slice sampling loop

    # update slice proposal scale based on the relative size of the slices compared to the initial guess
    prop.scale = prop.scale * nexpand / (2.0 * ncontract)

    return u, v, logl, nc
end   # end of function Slice

"""
    Proposals.RSlice(;slices=5, scale=1)

Propose a new live point by a series of random slices away from an existing live point. This is a standard _random_ implementation where each slice is along a random direction based on the provided axes. This more closely matches the PolyChord implementation outlined in Handley et al. (2015a,b).[^7][^8]

[^7]: Handley, et al., 2015a, MNRAS 450(1), ["polychord: nested sampling for cosmology"](https://academic.oup.com/mnrasl/article/450/1/L61/986122)
[^8]: Handley, et al., 2015b, MNRAS 453(4), ["polychord: next-generation nested sampling"](https://academic.oup.com/mnras/article/453/4/4384/2593718)

## Parameters
- `slices` is the minimum number of slices
- `scale` is the proposal distribution scale, which will update _between_ proposals.
"""
@with_kw mutable struct RSlice <: AbstractProposal
    slices = 5
    scale = 1.0

    @assert slices ≥ 1 "Number of slices must be greater than or equal to 1"
    @assert scale ≥ 0 "Proposal scale must be non-negative"
end

function (prop::RSlice)(
    rng::AbstractRNG,
    point::AbstractVector,
    logl_star,
    bounds::AbstractBoundingSpace,
    model;
    kwargs...
)
    # setup
    n = length(point)
    nc = nexpand = ncontract = 0
    local u, v, logl
    # random slice sampling loop
    for it in 1:prop.slices

        # propose a direction on the unit n-sphere
        drhat = randn(rng, n)
        drhat /= norm(drhat)

        # transform and scale into parameter space
        axis = prop.scale .* (Bounds.axes(bounds) * drhat)
        u, v, logl, nc, nexpand, ncontract = sample_slice(
            rng, axis, point, logl_star,
            model,
            nc, nexpand, ncontract
        )
    end # end of random slice sampling loop

    # update random slice proposal scale based on the relative size of the slices compared to the initial guess
    prop.scale = prop.scale * nexpand / (2.0 * ncontract)

    return u, v, logl, nc
end # end of function RSlice

# Method for slice sampling
function sample_slice(rng, axis, u, logl_star, model, nc, nexpand, ncontract)
    # define starting window
    r = rand(rng)  # initial scale/offset
    u_l = @. u - r * axis  # left bound
    if unitcheck(u_l)
        v_l, logl_l = prior_transform_and_loglikelihood(model, u_l)
    else
        logl_l = -Inf
    end
    nc += 1
    nexpand += 1

    u_r = u_l .+ axis # right bound
    if unitcheck(u_r)
        v_r, logl_r = prior_transform_and_loglikelihood(model, u_r)
    else
        logl_r = -Inf
    end
    nc += 1
    nexpand += 1

    # stepping out left and right bounds
    while logl_l ≥ logl_star
        u_l .-= axis
        if unitcheck(u_l)
            v_l, logl_l = prior_transform_and_loglikelihood(model, u_l)
        else
            logl_l = -Inf
        end
        nc += 1
        nexpand += 1
    end

    while logl_r ≥ logl_star
        u_r .+= axis
        if unitcheck(u_r)
            v_r, logl_r = prior_transform_and_loglikelihood(model, u_r)
        else
            logl_r = -Inf
        end
        nc += 1
        nexpand += 1
    end

    # sample within limits. If the sample is not valid, shrink the limits until the `logl_star` bound is hit
    window_init = norm(u_r - u_l)  # initial window size
    while true
        # define slice and window
        u_hat = u_r - u_l
        window = norm(u_hat)

        # check if the slice has shrunk to be ridiculously small
        window < 1e-5 * window_init && error("Slice sampling appears to be stuck.")

        # propose a new position
        r = rand(rng)
        u_prop = @. u_l + r * u_hat
        if unitcheck(u_prop)
            v_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop)
        else
            logl_prop = -Inf
        end
        nc += 1
        ncontract += 1

        # if success, then move to the new position
        if logl_prop ≥ logl_star
            return u_prop, v_prop, logl_prop, nc, nexpand, ncontract
        # if fail, then check if the new point is to the left/right of the original point along the proposal axis and update the bounds accordingly
        else
            s = dot(u_prop - u, u_hat) # check sign (+/-)
            if s < 0   # left
                u_l = u_prop
            elseif s > 0  # right
                u_r = u_prop
            else # if `s = 0` something has gone wrong
                error("Slice sampler has failed to find a valid point.")
            end
        end
    end # end of sample within limits while
end

end # module Proposals

ext: jl


text:
#=
Created on Monday 6 January 2019
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Functions for unconstrained optimization.
=#

module Unconstrained

import STMO: Tracker, trace, notrack
using LinearAlgebra

"""
    backtracking_line_search(f, x, Δx, ∇f, α::Real=0.1,
                        β::Real=0.7)

Uses backtracking for finding the minimum over a line.

Inputs:
    - f: function to be searched over a line
    - x: initial point
    - Δx: direction to search
    - ∇f: gradient of f
    - α: hyperparameter
    - β: hyperparameter

Output:
    - t: suggested stepsize
"""
function backtracking_line_search(f, x, Δx, ∇f; α::Real=0.1,
                        β::Real=0.7)
    @assert 0 < α < 0.5 && 0 < β < 1 "incorrect values for α and/or β"
    t = 1.0
    while f(x + t * Δx) > f(x) + α * t * ∇f(x)' * Δx
        t *= β
    end
    return t
end

"""
    gradient_descent(f, x₀, ∇f; α::Real=0.2, β::Real=0.7,
        ν::Real=1e-3, tracker::Tracker=notrack)

General gradient descent algorithm.

Inputs:
    - f: function to be minimized
    - x₀: starting point
    - ∇f: gradient of the function to be minimized
    - α: parameter for btls
    - β: parameter for btls
    - ν: parameter to determine if the algorithm is converged
    - tracker: a structure to store the path

Outputs:
    - xstar: the found minimum
"""
function gradient_descent(f, x₀, ∇f; α::Real=0.2, β::Real=0.7,
      ν::Real=1e-5, tracker::Tracker=notrack)

    x = x₀  # initial value
    Δx = similar(x)
    while true
        Δx .= -∇f(x) # choose direction
        if norm(Δx) < ν
            break  # converged
        end
        t = backtracking_line_search(f, x, Δx, ∇f, α=α, β=β)
        x .+= t * Δx  # do a step
        trace(tracker, x)
    end
    return x
end

"""
    coordinate_descent(f, x₀, ∇f; α::Real=0.2, β::Real=0.7,
        ν::Real=1e-3, tracker::Tracker=notrack)

General coordinate descent algorithm.

Inputs:
    - f: function to be minimized
    - x₀: starting point
    - ∇f: gradient of the function to be minimized
    - α: parameter for btls
    - β: parameter for btls
    - ν: parameter to determine if the algorithm is converged
    - tracker: a structure to store the path

Outputs:
    - xstar: the found minimum
"""
function coordinate_descent(f, x₀::Vector, ∇f; α::Real=0.2, β::Real=0.7,
      ν::Real=1e-5, tracker::Tracker=notrack)
    x = x₀  # initial value
    Δx = zero(x)
    ∇fx = similar(x)
    while true
        ∇fx .= ∇f(x)
        i = argmax(abs.(∇fx))   # choose direction
        Δx[i] = -∇fx[i]
        if norm(∇fx) < ν
            break  # converged
        end
        t = backtracking_line_search(f, x, Δx, ∇f, α=α, β=β)  # BLS for optimal step size
        x .+= t * Δx  # do a step
        Δx[i] = 0.0  # reset
        trace(tracker, x)
    end
    return x
end

"""
    newtons_method(f, x₀, ∇f, ∇²f; α::Real=0.2, β::Real=0.7,
        ϵ::Real=1e-7, tracker::Tracker=notrack)

General Newton method.

Inputs:
    - f: function to be minimized
    - x₀: starting point
    - ∇f: gradient of the function to be minimized
    - ∇²f: Hessian of the function to be minimized
    - α: parameter for btls
    - β: parameter for btls
    - ϵ: parameter to determine if the algorithm is converged
    - tracker: (bool) store the path that is followed?

Outputs:
    - xstar: the found minimum

"""
function newtons_method(f, x₀, ∇f, ∇²f; α::Real=0.2, β::Real=0.7,
      ϵ::Real=1e-5, tracker::Tracker=notrack)

    x = x₀  # initial value
    Δx = similar(x)
    ∇fx = similar(x)
    while true
        ∇fx .= ∇f(x)
        Δx .= - ∇²f(x) \ ∇fx # choose direction
        λ² = - (Δx' * ∇fx)  # newton decrement
        if λ² < ϵ
            break  # converged
        end
        t = backtracking_line_search(f, x, Δx, ∇f, α=α, β=β)
        x .+= t * Δx  # do a step
        trace(tracker, x)
    end
    return x
end

export gradient_descent, coordinate_descent, newtons_method, backtracking_line_search

end


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = 1.0;
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    f = mu[1]*q;
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = 0.0;
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = [u u];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    #ui = sin(mu[2]*x[1]+mu[3]*x[2]-mu[4]*sqrt(mu[2]*mu[2]+mu[3]*mu[3])*t);
    uix = mu[2]*cos(mu[2]*x[1]+mu[3]*x[2]-mu[4]*sqrt(mu[2]*mu[2]+mu[3]*mu[3])*t);
    uiy = mu[3]*cos(mu[2]*x[1]+mu[3]*x[2]-mu[4]*sqrt(mu[2]*mu[2]+mu[3]*mu[3])*t);
    fb = [u uix*n[1]+uiy*n[2]]; # Neumann conditions
    return fb;
end
function initu(x, mu, eta)
    u0 = 0.0;
    return u0;
end
function initq(x, mu, eta)
    q0 = [0.0 0.0];
end
function initw(x, mu, eta)
    w0 = 0.0;
    return w0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
#using MathProgBase
import MathOptInterface
#using SparseArrays

########################################################
## this code is based on ModelReader in NLPModels
## and KNITRO.jl
########################################################

export OnePhaseSolver

const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities

# TODO
const SF = Union{MOI.ScalarAffineFunction{Float64},
                 MOI.ScalarQuadraticFunction{Float64}}
const VAF = MOI.VectorAffineFunction{Float64}
const VOV = MOI.VectorOfVariables

# ScalarAffineFunctions and VectorAffineFunctions
const SAF = MOI.ScalarAffineFunction{Float64}
const AF = Union{SAF, VAF}

const SS = Union{MOI.EqualTo{Float64},
                 MOI.GreaterThan{Float64},
                 MOI.LessThan{Float64},
                 MOI.Interval{Float64}}
# LinSets
const LS = Union{MOI.EqualTo{Float64},
                 MOI.GreaterThan{Float64},
                 MOI.LessThan{Float64}}
# VecLinSets
const VLS = Union{MOI.Nonnegatives,
                  MOI.Nonpositives,
                  MOI.Zeros}
##################################################
mutable struct VariableInfo
    lower_bound::Float64  # May be -Inf even if has_lower_bound == true
    has_lower_bound::Bool # Implies lower_bound == Inf
    lower_bound_dual_start::Union{Nothing, Float64}
    upper_bound::Float64  # May be Inf even if has_upper_bound == true
    has_upper_bound::Bool # Implies upper_bound == Inf
    upper_bound_dual_start::Union{Nothing, Float64}
    is_fixed::Bool        # Implies lower_bound == upper_bound and !has_lower_bound and !has_upper_bound.
    start::Union{Nothing, Float64}
end
VariableInfo() = VariableInfo(-Inf, false, nothing, Inf, false, nothing, false, nothing)

mutable struct OnePhaseProblem
    status::Symbol  # Final status
    #status::MOI.TerminationStatusCode
    x::Vector{Float64}  # Starting and final solution
    lambda::Vector{Float64}
    g::Vector{Float64}  # Final constraint values
    obj_val::Float64  # (length 1) Final objective
    solve_time::Float64

    # Custom attributes of the OnePhaseSolver
    iter::Class_iterate
    hist::Array{alg_history2,1}
    pars::Class_parameters

    function OnePhaseProblem()
        return new()
    end
end

##################################################
# EmptyNLPEvaluator for non-NLP problems.
struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end

MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x) = nothing
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    @assert length(J) == 0
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, s, mu)
    @assert length(H) == 0
    return
end
function MOI.eval_hessian_lagrangian(::Nothing, H, x, s, mu)
    @assert length(H) == 0
    return
end

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

mutable struct ConstraintInfo{F, S}
    func::F
    set::S
    dual_start::Union{Nothing, Float64}
end

ConstraintInfo(func, set) = ConstraintInfo(func, set, nothing)

mutable struct OnePhaseSolver <: MOI.AbstractOptimizer
    #inner::OnePhaseProblem
	inner::Union{OnePhaseProblem, Nothing}

    # Problem data.
	variable_info::Vector{VariableInfo}
	nlp_data::MOI.NLPBlockData
	sense :: MOI.OptimizationSense
	objective::Union{MOI.VariableIndex, MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}, Nothing}
    linear_le_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}}
    linear_ge_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}}
    linear_eq_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}}
	linear_int_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}}
    quadratic_le_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}}
    quadratic_ge_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}}
    quadratic_eq_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}}
	quadratic_int_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64}, MOI.Interval{Float64}}}
    nlp_dual_start::Union{Nothing, Vector{Float64}}

	##Constraint mappings
	constraint_mapping::Dict{MOI.ConstraintIndex, Union{Cint, Vector{Cint}}}
	# Parameters.
	silent::Bool
	options::Dict{String, Any}

	# Solution attributes.
    solve_time::Float64
end

function OnePhaseSolver(; options...)
	options_dict = Dict{String, Any}()

	for (name, value) in options
        options_dict[string(name)] = value
    end

	onePhaseSolverModel = OnePhaseSolver(
        OnePhaseProblem(),
        [],
        empty_nlp_data(),
        MOI.FEASIBILITY_SENSE,
        nothing,
        [],
        [],
        [],
        [],
        [],
        [],
		[],
		[],
        nothing,
		Dict{MOI.ConstraintIndex, Int}(),
        false,
        options_dict,
        NaN,
    )
    set_options(onePhaseSolverModel, options)

	return onePhaseSolverModel
end

function set_options(model::OnePhaseSolver, options)
    for (name, value) in options
        sname = string(name)
        MOI.set(model, MOI.RawOptimizerAttribute(sname), value)
    end
    return
end

MOI.get(::OnePhaseSolver, ::MOI.SolverName) = "OnePhaseSolver"

"""
    MOI.is_empty(model::OnePhaseSolver )
"""

function MOI.is_empty(model::OnePhaseSolver)
    return isempty(model.variable_info) &&
           model.nlp_data.evaluator isa EmptyNLPEvaluator &&
           model.sense == MOI.FEASIBILITY_SENSE &&
           isempty(model.linear_le_constraints) &&
           isempty(model.linear_ge_constraints) &&
           isempty(model.linear_eq_constraints) &&
		   isempty(model.linear_int_constraints) &&
           isempty(model.quadratic_le_constraints) &&
           isempty(model.quadratic_ge_constraints) &&
           isempty(model.quadratic_eq_constraints) &&
		   isempty(model.quadratic_int_constraints)
end

# MOI.get(model::OnePhaseSolver, ::MOI.SolveTime) = model.solve_time

function MOI.empty!(model::OnePhaseSolver)
    model.inner = nothing
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
	empty!(model.linear_int_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
	empty!(model.quadratic_int_constraints)
    model.nlp_dual_start = nothing
end

function MOI.get(model::OnePhaseSolver, ::MOI.ListOfVariableIndices)
	return [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
end

function MOI.add_variable(model::OnePhaseSolver)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end

function MOI.add_variables(model::OnePhaseSolver, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function MOI.is_valid(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return vi.value in eachindex(model.variable_info)
end

function MOI.Utilities.supports_default_copy_to(::OnePhaseSolver, copy_names::Bool)
    return !copy_names
end

function MOI.copy_to(model::OnePhaseSolver, src::MOI.ModelLike; copy_names = false)
    return MOI.Utilities.default_copy_to(model, src, copy_names)
end
# MathOptInterface constraints

##################################################
## Support constraints
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.VariableIndex}, ::Type{<:SS}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.Interval{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::OnePhaseSolver, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.Interval{Float64}}) = true

MOI.supports_constraint(::OnePhaseSolver, ::Type{<:SF}, ::Type{<:SS}) = true

function has_upper_bound(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].has_upper_bound
end

function has_lower_bound(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].has_lower_bound
end

function is_fixed(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return model.variable_info[vi.value].is_fixed
end

function MOI.add_constraint(
    model::OnePhaseSolver, v::MOI.VariableIndex, lt::MOI.LessThan{Float64},
)
	vi = v
    MOI.throw_if_not_valid(model, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(model, vi)
        throw(MOI.UpperBoundAlreadySet{typeof(lt), typeof(lt)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Float64}, typeof(lt)}(vi))
    end

	col = vi.value
    model.variable_info[col].upper_bound = lt.upper
    model.variable_info[col].has_upper_bound = true
	ci = MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}(col)
    return ci
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}},
    set::MOI.LessThan{Float64},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].upper_bound = set.upper
    return
end

function MOI.delete(
    model::OnePhaseSolver,
    ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].upper_bound = Inf
    model.variable_info[ci.value].has_upper_bound = false
    return
end

function MOI.add_constraint(
    model::OnePhaseSolver, v::MOI.VariableIndex, gt::MOI.GreaterThan{Float64},
)
	vi = v
    MOI.throw_if_not_valid(model, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        throw(MOI.LowerBoundAlreadySet{typeof(gt), typeof(gt)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{Float64}, typeof(gt)}(vi))
    end

	col = vi.value
    model.variable_info[col].lower_bound = gt.lower
    model.variable_info[col].has_lower_bound = true
	ci = MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}(col)
	return ci
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}},
    set::MOI.GreaterThan{Float64},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = set.lower
    return
end

function MOI.delete(
    model::OnePhaseSolver,
    ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = -Inf
    model.variable_info[ci.value].has_lower_bound = false
    return
end

function MOI.add_constraint(
    model::OnePhaseSolver,  v::MOI.VariableIndex, eq::MOI.EqualTo{Float64},
)
	vi = v
    MOI.throw_if_not_valid(model, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(eq.value).")
    end
    if has_lower_bound(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Float64}, typeof(eq)}(vi))
    end
    if has_upper_bound(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Float64}, typeof(eq)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.LowerBoundAlreadySet{typeof(eq), typeof(eq)}(vi))
    end

	col = vi.value
    model.variable_info[col].lower_bound = eq.value
    model.variable_info[col].upper_bound = eq.value
    model.variable_info[col].is_fixed = true
	ci = MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}(col)
	return ci
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}},
    set::MOI.EqualTo{Float64},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = set.value
    model.variable_info[ci.value].upper_bound = set.value
    return
end

function MOI.delete(
    model::OnePhaseSolver,
    ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = -Inf
    model.variable_info[ci.value].upper_bound = Inf
    model.variable_info[ci.value].is_fixed = false
    return
end

macro define_add_constraint(function_type, set_type, prefix)
    array_name = Symbol(string(prefix) * "_constraints")
    return quote
        function MOI.add_constraint(
            model::OnePhaseSolver, func::$function_type, set::$set_type,
        )
            check_inbounds(model, func)
            push!(model.$(array_name), ConstraintInfo(func, set))

			col = length(model.$(array_name))
			ci = MOI.ConstraintIndex{$function_type, $set_type}(col)
			model.constraint_mapping[ci] = convert(Cint, col)
			return ci
        end
    end
end

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, linear_le,
)

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, linear_ge,
)

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, linear_eq,
)

@define_add_constraint(
    MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}, linear_int,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, quadratic_le,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}, quadratic_ge,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}, quadratic_eq,
)

@define_add_constraint(
    MOI.ScalarQuadraticFunction{Float64}, MOI.Interval{Float64}, quadratic_int,
)

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.ObjectiveFunction,
    func::Union{
        MOI.VariableIndex,
        MOI.ScalarAffineFunction,
        MOI.ScalarQuadraticFunction,
    },
)
    check_inbounds(model, func)
    model.objective = func
    return
end

function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector)
  NLPModels.increment!(nlp, :neval_obj)
  if nlp.obj.type == "LINEAR"
    res = dot(nlp.obj.gradient, x) + nlp.obj.constant
  end
  if nlp.obj.type == "QUADRATIC"
    res =
      0.5 * coo_sym_dot(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, x, x) +
      dot(nlp.obj.gradient, x) +
      nlp.obj.constant
  end
  if nlp.obj.type == "NONLINEAR"
    res = MOI.eval_objective(nlp.eval, x)
  end
  return res
end

function append_to_hessian_sparsity!(
    ::Any,
    ::Union{MOI.VariableIndex,MOI.ScalarAffineFunction},
)
    return nothing
end

function append_to_hessian_sparsity!(
    hessian_sparsity,
    quad::MOI.ScalarQuadraticFunction,
)
    for term in quad.quadratic_terms
        push!(
            hessian_sparsity,
            # (term.variable_index_1.value, term.variable_index_2.value),
			(term.variable_1.value, term.variable_2.value),
        )
    end
end

function hessian_lagrangian_structure(model::OnePhaseSolver, nlp :: MathOptNLPModel)
    hessian_sparsity = Tuple{Int64,Int64}[]
    if model.objective !== nothing
        append_to_hessian_sparsity!(hessian_sparsity, model.obj)
    end
    for info in model.quadratic_le_constraints
        append_to_hessian_sparsity!(hessian_sparsity, info.func)
    end
    for info in model.quadratic_ge_constraints
        append_to_hessian_sparsity!(hessian_sparsity, info.func)
    end
    for info in model.quadratic_eq_constraints
        append_to_hessian_sparsity!(hessian_sparsity, info.func)
    end
    nlp_hessian_sparsity =
        MOI.hessian_lagrangian_structure(nlp.eval)
    append!(hessian_sparsity, nlp_hessian_sparsity)
    return hessian_sparsity
end

function hess_coord(nlp :: MathOptNLPModel, x :: Array{Float64};
  obj_weight :: Float64=1.0, y :: Array{Float64}=zeros(nlp.meta.ncon))
  NLPModels.increment!(nlp, :neval_hess)
  MOI.eval_hessian_lagrangian(nlp.eval, nlp.obj.hessian.vals, x, obj_weight, y)

  return (NLPModels.hess_structure(nlp)[1], NLPModels.hess_structure(nlp)[2], NLPModels.hess_coord(nlp, x, y, obj_weight=obj_weight))
end

function NLPModels.hess_structure!(
  nlp::MathOptNLPModel,
  rows::AbstractVector{<:Integer},
  cols::AbstractVector{<:Integer},
)
  if nlp.obj.type == "QUADRATIC"
    for index = 1:(nlp.obj.nnzh)
      rows[index] = nlp.obj.hessian.rows[index]
      cols[index] = nlp.obj.hessian.cols[index]
    end
  end
  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
    hesslag_struct = MOI.hessian_lagrangian_structure(nlp.eval)
    for index = (nlp.obj.nnzh + 1):(nlp.meta.nnzh)
      shift_index = index - nlp.obj.nnzh
      rows[index] = hesslag_struct[shift_index][1]
      cols[index] = hesslag_struct[shift_index][2]
    end
  end
  return rows, cols
end

############################
## END ModelReader CODE
############################

function status_One_Phase_To_JuMP(status::Symbol)
    # since our status are not equal to JuMPs we need to do a conversion
    if status == :Optimal
        return :Optimal
    elseif status == :primal_infeasible
        return :Infeasible
    elseif status == :dual_infeasible
        return :Unbounded
    elseif status == :MAX_IT || status === :MAX_TIME
        return :UserLimit
    else
        return :Error
    end
end

function create_pars_JuMP(options )
    pars = Class_parameters()
    for (param,value) in options
      what = split(String(param),"!") # we represent a parameter such as init.mu_scale as init!mu_scale because we cannot pass init.mu_scale as a parameter
      node = pars # root
      for i in 1:length(what)
          field = what[i]
          if i < length(what)
              node = getfield(node,Symbol(field))
          else # at the leaf
              setfield!(node,Symbol(field),value)
          end
      end
    end

    return pars
end

MOI.get(model::OnePhaseSolver, ::MOI.RawSolver) = model

# copy
MOI.supports_incremental_interface(solver::OnePhaseSolver) = true
function MOI.copy_to(solver::OnePhaseSolver, src::MOI.ModelLike)
    return MOI.Utilities.automatic_copy_to(model, src)
end

function MOI.optimize!(solver :: OnePhaseSolver)
    t = time()

	has_nlp_objective = false
	if !isa(solver.nlp_data.evaluator, EmptyNLPEvaluator)
		features = MOI.features_available(solver.nlp_data.evaluator)::Vector{Symbol}
		has_hessian = (:Hess in features)
		has_hessvec = (:HessVec in features)
		has_nlp_objective = solver.nlp_data.has_objective
		num_nlp_constraints = length(solver.nlp_data.constraint_bounds)
		has_nlp_constraints = (num_nlp_constraints > 0)

		init_feat = Symbol[]
		has_nlp_objective && push!(init_feat, :Grad)
		if has_hessian
			push!(init_feat, :Hess)
		elseif has_hessvec
			push!(init_feat, :HessVec)
		end

		if has_nlp_constraints
			push!(init_feat, :Jac)
		end

		MOI.initialize(solver.nlp_data.evaluator, init_feat)
	end

    pars = create_pars_JuMP(solver.options)

    iter, status, hist, t, err, timer = one_phase_solve(solver, pars)
    solver.inner = OnePhaseProblem()
    solver.inner.status = status_One_Phase_To_JuMP(status)
    solver.inner.x = get_original_x(iter)
    solver.inner.obj_val = iter.cache.fval
    solver.inner.lambda = get_y(iter)
    solver.inner.solve_time = time() - t

    # custom one phase features
    solver.inner.pars = pars
    solver.inner.iter = iter
    solver.inner.hist = hist
end

function MOI.optimize!(solver :: OnePhaseSolver, jumpModel:: Model)
    t = time()
    nlp = MathOptNLPModel(jumpModel)
    pars = create_pars_JuMP(solver.options)

    iter, status, hist, t, err, timer = one_phase_solve(nlp,pars)
    solver.inner = OnePhaseProblem()
    solver.inner.status = status_One_Phase_To_JuMP(status)
    solver.inner.x = get_original_x(iter)
    solver.inner.obj_val = iter.cache.fval
    solver.inner.lambda = get_y(iter)
    solver.inner.solve_time = time() - t

    # custom one phase features
    solver.inner.pars = pars
    solver.inner.iter = iter
    solver.inner.hist = hist
end

function MOI.get(
    model::MathOptInterface.Utilities.CachingOptimizer,
    attr::MOI.TerminationStatus,
)
    return MOI.get(model.optimizer, attr)
end

function check_inbounds(model::OnePhaseSolver, vi::MOI.VariableIndex)
    return MOI.throw_if_not_valid(model, vi)
end

function check_inbounds(model::OnePhaseSolver, aff::MOI.ScalarAffineFunction)
    for term in aff.terms
		MOI.throw_if_not_valid(model, term.variable)
    end
end

function check_inbounds(model::OnePhaseSolver, quad::MOI.ScalarQuadraticFunction)
    for term in quad.affine_terms
		MOI.throw_if_not_valid(model, term.variable)
    end
    for term in quad.quadratic_terms
		MOI.throw_if_not_valid(model, term.variable_1)
        MOI.throw_if_not_valid(model, term.variable_2)
    end
end

MOI.supports(::OnePhaseSolver, ::MOI.NLPBlock) = true

function MOI.supports(
    ::OnePhaseSolver, SF
)
    return true
end

function MOI.supports(
    ::OnePhaseSolver, ::MOI.ObjectiveFunction{MOI.VariableIndex}
)
    return true
end

function MOI.supports(
    ::OnePhaseSolver, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}
)
    return true
end

function MOI.supports(
    ::OnePhaseSolver, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}
)
    return true
end

MOI.supports(::OnePhaseSolver, ::MOI.ObjectiveSense) = true

MOI.supports(::OnePhaseSolver, ::MOI.Silent) = true

MOI.supports(::OnePhaseSolver, ::MOI.RawOptimizerAttribute) = true

function MOI.get(model::OnePhaseSolver, ::MOI.ObjectiveFunction)
    return model.objective
end

function MOI.set(model::OnePhaseSolver, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

function MOI.set(
    model::OnePhaseSolver, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense
)
    model.sense = sense
    return
end

MOI.get(model::OnePhaseSolver, ::MOI.ObjectiveSense) = model.sense

function MOI.set(model::OnePhaseSolver, ::MOI.Silent, value)
    model.silent = value
    return
end

MOI.get(model::OnePhaseSolver, ::MOI.Silent) = model.silent

function MOI.supports(
    ::OnePhaseSolver, ::MOI.VariablePrimalStart, ::Type{MOI.VariableIndex}
)
    return true
end

function MOI.set(
    model::OnePhaseSolver,
    ::MOI.VariablePrimalStart,
    vi::MOI.VariableIndex,
    value::Union{Real, Nothing},
)
    MOI.throw_if_not_valid(model, vi)
    model.variable_info[vi.value].start = value
    return
end

function MOI.get(model::MOIU.CachingOptimizer, args...)
	return MOI.get(model.model.optimizer, args)
end

function MOI.get(model::OnePhaseSolver, attr::MOI.VariablePrimal, v::VariableRef)
	return MOI.get(model, attr, v.variable)
end

function MOI.get(
    model::OnePhaseSolver, attr::MOI.VariablePrimal, vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, vi)
    return model.inner.x[vi.value]
end

function MOI.set(model::OnePhaseSolver, ::MOI.TimeLimitSec, value::Real)
    MOI.set(model, MOI.RawOptimizerAttribute(TIME_LIMIT), Float64(value))
end

function MOI.set(model::OnePhaseSolver, p::MOI.RawOptimizerAttribute, value)
    model.options[p.name] = value
    return
end

function MOI.get(model::OnePhaseSolver, p::MOI.RawOptimizerAttribute)
    if haskey(model.options, p.name)
        return model.options[p.name]
    end
    error("RawParameter with name $(p.name) is not set.")
end

function MOI.get(model::OnePhaseSolver, ::MOI.TerminationStatus)
    if model.inner === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
	status = model.inner.status
	return status
end

function MOI.get(model::OnePhaseSolver, ::MOI.RawStatusString)
	return string(model.inner.status)
end

# Ipopt always has an iterate available.
function MOI.get(model::OnePhaseSolver, ::MOI.ResultCount)
    return (model.inner !== nothing) ? 1 : 0
end

function MOI.get(model::OnePhaseSolver, attr::MOI.PrimalStatus)
    if !(1 <= attr.N <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end

	status = model.inner.status
	return status
end

function MOI.get(model::OnePhaseSolver, attr::MOI.DualStatus)
    if !(1 <= attr.N <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end
	status = model.inner.status
	return status
end

function MOI.get(model::OnePhaseSolver, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return model.inner.obj_val
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
x = 3

y = 2x

typeof(y)

f(x) = 2 + x

f

f(10)

function g(x, y)
    z = x + y
    return z^2
end

g(1, 2)

let s = 0
    for i in 1:10
        s += i    # Equivalent to s = s + i
    end
    s
end

typeof(1:10)

function mysum(n)
    s = 0
    for i in 1:n
        s += i
    end
    return s
end

mysum(100)

a = 3

a < 5

if a < 5
    "small"
else
    "big"
end

v = [1, 2, 3]

typeof(v)

v[2]

v[2] = 10

v2 = [i^2 for i in 1:10]

M = [1 2
     3 4]

typeof(M)

zeros(5, 5)

zeros(Int, 4, 5)

[i + j for i in 1:5, j in 1:6]


text:
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Flux
using Flux: Optimise.update!
using Statistics: mean

struct Critic{C}
    c1::C
    c2::C
end
(m::Critic)(s, a) = (inp = vcat(s, a); (m.c1(inp), m.c2(inp)))
Flux.@functor Critic

function act(actor, s, noisescale, noiseclamp=false)
    a = actor(s)
    noise = noisescale .* randn(Float32, size(a))
    noise = noiseclamp ? clamp.(noise, -5f-1, 5f-1) : noise
    return clamp.(a .+ noise, -1f0, 1f0)
end

function run!(actor, aopt, critic, copt, env;
        epochs=100, steps=300, maxt=100,  batchsize=128,
        noisescale=1f-1, noisescaletar=2f-1, γ=99f-2, τ=1f-3)
    rewards = zeros(Float32, epochs)
    actortar = deepcopy(actor)
    critictar = deepcopy(critic)
    memory = ReplayMemory{Float32, Float32, Float32}(
        length(env.observationspace), length(env.actionspace), 10*steps, batchsize
    )
    actorps = params(actor)
    criticps = params(critic)

    for i in 1:epochs
        episodes!(memory, env, steps, maxt) do s
            return act(actor, s, noisescale)
        end

        for z in 1:steps
            let (s1, a1, r, s2, done) = sample(memory)
                a2 = act(actortar, s2, noisescaletar, true)
                Qtar = r .+ (1 .- done) .* γ .* vec(min.(critictar(s2, a2)...))
                cgs = gradient(criticps) do
                    Q1s1a1, Q2s1a1 = critic(s1, a1)
                    return Flux.mse(vec(Q1s1a1), Qtar) + Flux.mse(vec(Q2s1a1), Qtar)
                end
                update!(copt, criticps, cgs)

                if z % 2 == 0
                    ags = gradient(actorps) do
                        a = actor(s1)
                        return -mean(critic(s1, a)[1])
                    end
                    update!(aopt, actorps, ags)

                    softupdate!(actortar, actor, τ)
                    softupdate!(critictar, critic, τ)
                end
            end
        end
        rewards[i] = test(s -> actor(s), env)
    end
    return rewards
end

env = CartPoleContinuous{Float32}()
aDim = length(env.actionspace)
oDim = length(env.observationspace)

actor = Chain(
    Dense(oDim, 32, relu),
    Dense(32, aDim, tanh)
)
critic = Critic(
    Chain(Dense(oDim+aDim, 32, relu), Dense(32, 16, relu), Dense(16, 1)),
    Chain(Dense(oDim+aDim, 32, relu), Dense(32, 16, relu), Dense(16, 1))
)

r = run!(actor, ADAM(), critic, ADAM(), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "td3")


text:
module CellLists

include("serial.jl")
include("threads.jl")
export CellList, near_neighbors, distance_condition, merge

end # module


text:
using linearstep
using Test

@time @testset "linearstep.jl" begin include("test_adam.jl") end
@time @testset "compare_adam.jl" begin include("test_compare_adam.jl") end


text:
module RobustDecisionProgramming

include("distributional_robustness.jl")
export Deviation, cross_assignment, polyhedral_uncertainty, polyhedral_uncertainty_set

end # module

max_stars_repo_name: jaantollander/RobustDecisionProgramming

max_issues_repo_name: jaantollander/RobustDecisionProgramming

max_forks_repo_name: jaantollander/RobustDecisionProgramming


text:
#############################################################################
#############################################################################
#
# This file implements factorization
#
#############################################################################
#############################################################################

"""
Factors a polynomial over the field Z_p.

Returns a vector of tuples of (irreducible polynomials (mod p), multiplicity) such that their product of the list (mod p) is f. Irreducibles are fixed points on the function factor.
"""

function factor(f::PolynomialModP, show::Bool = false)::Vector{Tuple{PolynomialModP,Int}}

    #Cantor Zassenhaus factorization

    degree(f) ≤ 1 && return [(f,1)]

    # make f primitive
    ff = prim_part(f)

     # make f square-free
    squares_poly = gcd(f, derivative(ff))

    ff = ff ÷ squares_poly
    show && println(" The No_Square Polynomial Based of f:")
    show && println(ff)
    show && println("")

    # make f monic
    old_coeff = leading(ff).coeff
    ff = ff ÷ old_coeff

    show && println(" The monic Polynomial Based of f:")
    show && println(ff)
    show && println("")

    dds = dd_factor(ff)

    ret_val = Tuple{PolynomialModP,Int}[]

    for (k,dd) in enumerate(dds)
        sp = dd_split(dd, k)
        sp = map((p)->p ÷ leading(p).coeff,sp) #makes the polynomials inside the list sp, monic
        for mp in sp
            push!(ret_val, (mp, multiplicity(f,mp)) )
        end
    end

    #Append the leading coefficient as well
    push!(ret_val, (leading(f).coeff* one(PolynomialModP,f.mod), 1) )

    degreetest = 0
    for i in ret_val
        degreetest = degreetest + degree(i[1]) * i[2]
    end

    if degreetest != degree(f) #when there is a factor with a multiplicity of p (the modulus)
                         #it does not show up in the factorisation since its derivative is 0
                         #This tries find these lost factors.
        l = pop!(ret_val)
        temp_poly = 2*l[1]
        push!(ret_val,(temp_poly,1))

        d = degree(f) - degreetest #this is the degree of this missing polynomial
        x = x_poly(f.mod)
        p = f.mod

        miss_poly =  f ÷ expand_factorization(ret_val) #Figuring out what the missing polynomial is so we can factorise it

        for i in 1:d÷p
            cand = gcd(x^(p^(i))-x,miss_poly) #looks for candidate factor on irriducable polys
            if degree(cand) == i  #if it is infact irriducable, we can use it straight away and add it to the factorisation

                push!(ret_val,(cand,p))

            elseif degree(cand) > 0 #else, it is the productt of multiple irriducable polys, ad thus we must factorise again

                a = factor(cand)
                for j in a

                    j[2] = j[2]*p
                end
                ret_val = vcat(ret_val,a)
            end

        end
    end

    return ret_val
end

"""
Expand a factorization.
"""

function expand_factorization(factorization::Vector{Tuple{PolynomialModP,Int}})::PolynomialModP
    length(factorization) == 1 && return first(factorization[1])^last(factorization[1])
    return *([first(tt)^last(tt) for tt in factorization]...)
end

"""
Compute the number of times g divides f
"""

function multiplicity(f::PolynomialModP, g::PolynomialModP)::Int
    (iszero(f) || iszero(g)) && return 0
    degree(gcd(f, g)) == 0 && return 0
    return 1 + multiplicity(f ÷ g, g)
end

"""
Distinct degree factorization.

Given a square free polynomial `f` returns a list, `g` such that `g[k]` is a product of irreducible polynomials of degree `k` for `k` in 1,...,degree(f) ÷ 2, such that the product of the list (mod `prime`) is equal to `f` (mod `prime`).
"""
function dd_factor(f::PolynomialModP)::Array{PolynomialModP}
    x = x_poly(f.mod)
    w = deepcopy(x)
    g = Array{PolynomialModP}(undef,degree(f.polynomial)) #Array of polynomials indexed by degree
    #Looping over degrees
    for k in 1:degree(f)
        w = rem(w^(f.mod), f)
        g[k] = gcd(w - x, f)

        f = f ÷ g[k]
    end

    #edge case for final factor
    f != one(PolynomialModP,f.mod) && push!(g,f)

    return g
end

"""
Distinct degree split.

Returns a list of irreducible polynomials of degree `d` so that the product of that list (mod prime) is the polynomial `f`.
"""

function dd_split(f::PolynomialModP, d::Int)::Vector{PolynomialModP}

    degree(f) == d && return [f]
    degree(f) == 0 && return []

    w = rand(PolynomialModP, f.mod, degree = d, monic = true)

    n_power = ((f.mod)^d-1) ÷ 2

    g = gcd(w^n_power - one(PolynomialModP,f.mod), f)

    ḡ = f ÷ g # g\bar + [TAB]

    return vcat(dd_split(g, d), dd_split(ḡ, d) )
end


text:
struct HH3DFarField{K}
    gamma::K
end

function HH3DFarField(;wavenumber=error("wavenumber is a required argument"))
    if iszero(real(wavenumber))
        HH3DFarField(-imag(wavenumber))
    else
        HH3DFarField(wavenumber*im)
    end
end

quaddata(op::HH3DFarField,rs,els) = quadpoints(rs,els,(3,))
quadrule(op::HH3DFarField,refspace,p,y,q,el,qdata) = qdata[1,q]

kernelvals(op::HH3DFarField,y,p) = exp(op.gamma*dot(y,cartesian(p)))
function integrand(op::HH3DFarField,krn,y,f,p)
    krn * f.value
end


text:
module DiscreteModelsTests

using Test
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Geometry: DiscreteModelMock
using Gridap.Io

domain = (0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)

model_1d = DiscreteModel(Polytope{1},model)
test_discrete_model(model_1d)

model_2d = DiscreteModel(Polytope{2},model)
test_discrete_model(model_2d)

model = DiscreteModelMock()
test_discrete_model(model)

grid = get_grid(model)
topo = get_grid_topology(model)
labeling = get_face_labeling(model)

@test num_dims(model) == 2
@test num_cell_dims(model) == 2
@test num_point_dims(model) == 2
@test num_faces(model) == 9 + 13 + 5
@test num_faces(model,0) == 9
@test num_faces(model,1) == 13
@test num_faces(model,2) == 5
@test num_cells(model) == 5
@test num_vertices(model) == 9
@test num_edges(model) == 13
@test num_facets(model) == 13
@test num_nodes(model) == 9

@test get_cell_node_ids(model) == get_cell_node_ids(grid)
@test get_node_coordinates(model) == get_node_coordinates(grid)
@test get_cell_type(model) == get_cell_type(grid)
@test get_reffes(model) == get_reffes(grid)
@test get_face_nodes(model,0) == get_faces(topo,0,0)
@test get_face_nodes(model,1) == get_faces(topo,1,0)
@test get_face_nodes(model,2) == get_faces(topo,2,0)
@test get_face_nodes(model) == get_face_vertices(topo)
@test get_face_own_nodes(model,0) == get_faces(topo,0,0)
@test get_face_own_nodes(model,1) == empty_table(num_faces(model,1))
@test get_face_own_nodes(model,2) == empty_table(num_faces(model,2))
r = vcat(get_face_own_nodes(model,0),get_face_own_nodes(model,1),get_face_own_nodes(model,2))
@test get_face_own_nodes(model) == r
@test get_vertex_node(model) == collect(1:num_vertices(model))
@test get_node_face_owner(model) == collect(1:num_nodes(model))
@test get_reffaces(ReferenceFE{0},model) == [VERTEX1]
@test get_reffaces(ReferenceFE{1},model) == [SEG2]
@test get_reffaces(ReferenceFE{2},model) == [QUAD4, TRI3]
@test get_reffaces(model) == [VERTEX1, SEG2, QUAD4, TRI3]
@test get_face_type(model,0) == get_face_type(topo,0)
@test get_face_type(model,1) == get_face_type(topo,1)
@test get_face_type(model,2) == get_face_type(topo,2)
@test get_face_type(model) == get_face_type(topo)
@test get_reffaces_offsets(model) == [0,1,2]

grid0 = Grid(ReferenceFE{0},model)
grid1 = Grid(ReferenceFE{1},model)
grid2 = Grid(ReferenceFE{2},model)
test_grid(grid0)
test_grid(grid1)
test_grid(grid2)
@test num_dims(grid0) == 0
@test num_dims(grid1) == 1
@test num_dims(grid2) == 2

grid0 = Triangulation(ReferenceFE{0},model)
grid1 = Triangulation(ReferenceFE{1},model)
grid2 = Triangulation(ReferenceFE{2},model)
test_triangulation(grid0)
test_triangulation(grid1)
test_triangulation(grid2)
@test num_dims(grid0) == 0
@test num_dims(grid1) == 1
@test num_dims(grid2) == 2

model = DiscreteModelMock()
dict = to_dict(model)
model2 = from_dict(DiscreteModel,dict)
test_discrete_model(model2)

domain = (0,1,0,1,0,1)
partition = (3,3,3)
model = CartesianDiscreteModel(domain,partition)

tmodel = simplexify(model)
test_discrete_model(tmodel)

model2 = DiscreteModel(grid,topo,labeling)
test_discrete_model(model2)

reffes = ReferenceFE(model,lagrangian,Float64,1)
@test isa(reffes,AbstractVector{<:ReferenceFE})

d = mktempdir()

filename = joinpath(d,"model.json")
to_json_file(model2,filename)
model3 = DiscreteModelFromFile(filename)
test_discrete_model(model3)

rm(d,recursive=true)

end # module


text:
## -- Common elemental masses in g/mol (AMU/atom)

    molarmass = Dict{String,Float64}("Pd" => 106.421,"Fl" => 289.0,"Nb" => 92.906372,"C" => 12.011,"P" => 30.9737619985,"Ag" => 107.86822,"Gd" => 157.253,"Si" => 28.085,"Ru" => 101.072,"At" => 210.0,"Uus" => 294.0,"Sb" => 121.7601,"Cs" => 132.905451966,"Cn" => 285.0,"Uut" => 285.0,"Be" => 9.01218315,"Ac" => 227.0,"Cf" => 251.0,"Bh" => 270.0,"Sr" => 87.621,"Ga" => 69.7231,"Ta" => 180.947882,"Te" => 127.603,"Np" => 237.0,"Lr" => 262.0,"Pu" => 244.0,"U" => 238.028913,"Kr" => 83.7982,"Y" => 88.905842,"Sg" => 271.0,"Ca" => 40.0784,"Au" => 196.9665695,"K" => 39.09831,"Rn" => 222.0,"Ra" => 226.0,"Ce" => 140.1161,"Uuo" => 294.0,"V" => 50.94151,"Fr" => 223.0,"Mo" => 95.951,"Pr" => 140.907662,"Th" => 232.03774,"Br" => 79.904,"Zn" => 65.382,"He" => 4.0026022,"Sc" => 44.9559085,"H" => 1.008,"Al" => 26.98153857,"S" => 32.06,"Uup" => 289.0,"Ar" => 39.9481,"Ge" => 72.6308,"Er" => 167.2593,"Fe" => 55.8452,"Mg" => 24.305,"F" => 18.9984031636,"La" => 138.905477,"Rf" => 267.0,"W" => 183.841,"Li" => 6.94,"Dy" => 162.5001,"O" => 15.999,"B" => 10.81,"Bi" => 208.980401,"Mn" => 54.9380443,"Re" => 186.2071,"Db" => 270.0,"Hf" => 178.492,"Cm" => 247.0,"Cl" => 35.45,"In" => 114.8181,"Ds" => 281.0,"Rb" => 85.46783,"Po" => 209.0,"Lv" => 293.0,"Am" => 243.0,"Pa" => 231.035882,"Se" => 78.9718,"Ba" => 137.3277,"Nd" => 144.2423,"Pm" => 145.0,"Rh" => 102.905502,"Ti" => 47.8671,"Tb" => 158.925352,"Hs" => 277.0,"Zr" => 91.2242,"Sm" => 150.362,"Cr" => 51.99616,"Cu" => 63.5463,"Fm" => 257.0,"Tc" => 97.0,"Tl" => 204.38,"Ne" => 20.17976,"Hg" => 200.5923,"Mt" => 276.0,"N" => 14.007,"Es" => 252.0,"Yb" => 173.0455,"Lu" => 174.96681,"Eu" => 151.9641,"Na" => 22.989769282,"No" => 259.0,"Os" => 190.233,"Ni" => 58.69344,"Ho" => 164.930332,"Co" => 58.9331944,"Md" => 258.0,"Ir" => 192.2173,"Pt" => 195.0849,"Tm" => 168.934222,"As" => 74.9215956,"Sn" => 118.7107,"Xe" => 131.2936,"I" => 126.904473,"Cd" => 112.4144,"Pb" => 207.21,"Rg" => 282.0,"Bk" => 247.0)
    export molarmass

    molarmasspercation = Dict{String,Float64}("MgO" => 40.304,"SO2" => 64.058,"TiO2" => 79.8651,"CaO" => 56.077400000000004,"P2O5" => 70.9712619985,"NiO" => 74.69244,"CoO" => 74.9321944,"SiO2" => 60.083,"Al2O3" => 50.980038570000005,"CO2" => 44.009,"Cr2O3" => 75.99466000000001,"FeO" => 71.8442,"K2O" => 47.097809999999996,"Fe2O3" => 79.8437,"H2O" => 9.0075,"MnO" => 70.9370443,"Na2O" => 30.989269282000002)
    export molarmasspercation

## ---


text:
include("../src/MechaJulia.jl")

MechaJulia.check_pr_ok()

max_stars_repo_name: MechaJulia/MechaJulia

max_issues_repo_name: MechaJulia/MechaJulia

max_forks_repo_name: MechaJulia/MechaJulia


text:
using LinearAlgebra, Jets, JetPack, Test

n1,n2 = 5,2

# no correctness test for non zero α -- JopLeakyIntegration should handle that case
@testset "JopNormalize, correctness T=$(T)" for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    F = JopNormalize(JetSpace(T,n1,n2))
    x1 = ones(domain(F))
    @test F*x1 ≈ ones(domain(F)) ./ n1
end

@testset "JopNormalize, linearity test, T=$(T)" for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    F = JopNormalize(JetSpace(T,n1,n2))
    J = jacobian(F, rand(domain(F)))
    lhs,rhs = linearity_test(J)
    @test lhs ≈ rhs
    lhs,rhs = linearity_test(J')
    @test lhs ≈ rhs
end

@testset "JopNormalize, dot product test, T=$(T)" for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    F = JopNormalize(JetSpace(T,n1,n2))
    J  = jacobian!(F, rand(domain(F)))
    lhs, rhs = dot_product_test(J, -1 .+ 2 .* rand(domain(J)), -1 .+ 2 .* rand(range(J)))
    @test lhs ≈ rhs
end

# note the key here is to increase the size of the nonlinear vector
@testset "JopNormalize, linearization test, T=$(T)" for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    F = JopNormalize(JetSpace(T,n1,n2))
    m0 = rand(domain(F))
    μ  = sqrt.([1/1,1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256,1/512,1/1024,1/2048,1/4096,1/8192])
    δm = rand(domain(F))
    observed, expected = linearization_test(F, m0, μ = μ, δm = δm)
    δ = minimum(abs, observed - expected)
    @test δ < 0.1
end


text:
if get(ENV, "CMBLENSING_TEST_CUDA", false) == "1"
    using Adapt, CUDA
    CUDA.allowscalar(false)
    maybegpu(x) = adapt(CuArray,x)
    storage = CuArray
    @info "Running tests on $(repr("text/plain", CUDA.device()))"
else
    maybegpu(x) = x
    storage = Array
    @info "Running tests on CPU"
end

##

using CMBLensing
using CMBLensing:
    @SMatrix, @SVector, AbstractCℓs, basis, Basis, BasisTuple,
    LinearInterpolation, Measurement, RK4Solver, seed!, ±

##

using Test
using SparseArrays
using LinearAlgebra
using Zygote
using AbstractFFTs

##

Nsides  = [(8,  (8,8)),     ((8,8),    (8,8)),     ((4,8),    (4,8))]
Nsides′ = [(128,(128,128)), ((128,128),(128,128)), ((196,128),(196,128))]

##

@testset "CMBLensing" begin

##

# basic printing sanity checks, which were super annoying to get right...
# see also: https://discourse.julialang.org/t/dispatching-on-the-result-of-unwrap-unionall-seems-weird/25677

@testset "Printing" begin

    # concrete types:
    for f in [maybegpu(FlatMap(rand(4,4))), maybegpu(FlatQUMap(rand(4,4),rand(4,4)))]
        @test occursin("pixels",sprint(show, MIME("text/plain"), f))
        @test occursin("pixels",sprint(show, MIME("text/plain"), [f,f]))
    end

    for m in ((), (MIME("text/plain"),))
        # unionall types: (the presence of "where" indicates printing correctly
        # forwarded to the default behavior)
        @test occursin("where",sprint(show, m..., FieldTuple{<:Any,<:NamedTuple{(:Q,:U)}}))
        @test occursin("where",sprint(show, m..., FlatMap{<:Any,<:Any,<:Matrix{Real}}))
        @test occursin("where",sprint(show, m..., FlatQUMap))
    end

end
##

@testset "FieldTuples" begin

    f = maybegpu(FlatMap(rand(4,4)))

    @testset "Constructors" begin

        # Enumerate the combinations of
        # * 1) no basis or names specified, 2) only basis specified, 3) only names specified, 4) basis and names specified
        # * 1) args/kwargs form 2) single Tuple/NamedTuple argument form

        @testset "F :: $F" for F in [
            FieldTuple,
            FieldTuple{QUMap},
            FieldTuple{QUMap,<:NamedTuple{(:Q,:U)}},
            FieldTuple{<:Basis,<:NamedTuple{(:Q,:U)}}
        ]

            @test (@inferred F(;Q=f, U=f)) isa F
            @test (@inferred F((Q=f, U=f))) isa F
            @test (@inferred F(f,f)) isa F
            @test (@inferred F((f,f))) isa F

        end

    end

    @testset "Basis conversions" begin

        # basis conversions
        for f_basistuple in [FieldTuple(f, f), FieldTuple(A=f, B=f)] # named and unnamed
            @test basis(@inferred    Fourier(f_basistuple)) <: BasisTuple{Tuple{Fourier,Fourier}}
            @test basis(@inferred        Map(f_basistuple)) <: BasisTuple{Tuple{Map,Map}}
            @test basis(@inferred DerivBasis(f_basistuple)) <: BasisTuple{Tuple{Fourier,Fourier}}
            @test basis(@inferred BasisTuple{Tuple{Fourier,Fourier}}(f_basistuple)) <: BasisTuple{Tuple{Fourier,Fourier}}
        end

        f_concretebasis = FieldTuple{QUMap, <:NamedTuple{(:Q,:U)}}(f,f)
        @test basis(@inferred    Fourier(f_concretebasis)) <: QUFourier
        @test basis(@inferred        Map(f_concretebasis)) <: QUMap
        @test basis(@inferred DerivBasis(f_concretebasis)) <: QUFourier

    end

end

##

@testset "Flat Constructors" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        θpix = 3
        kwargs = (Nside=Nside, θpix=θpix)
        P = Flat(;kwargs...)
        Ix = maybegpu(rand(Nside2D...))
        Il = rfft(Ix)

        @testset "f :: $F" for (F,args) in [
            (FlatMap,        (Ix,)),
            (FlatFourier,    (Il,)),
            (FlatQUMap,      (Ix,Ix)),
            (FlatQUFourier,  (Il,Il)),
            (FlatEBMap,      (Ix,Ix)),
            (FlatEBFourier,  (Il,Il)),
            (FlatIQUMap,     (Ix,Ix,Ix)),
            (FlatIQUFourier, (Il,Il,Il)),
            (FlatIEBMap,     (Ix,Ix,Ix)),
            (FlatIEBFourier, (Il,Il,Il))
        ]

            # the four basic construtor types with no type changing
            @test F(args...; kwargs...) isa F{P}
            @test (@inferred F{P}(args...)) isa F{P}
            @test (@inferred F{P,Float64}(args...)) isa F{P}
            @test (@inferred F{P,Float64,typeof(args[1])}(args...)) isa F{P}

            # forcing the eltype to something different
            @test real(eltype(@inferred F{P,Float32}(args...))) == Float32

            # broadcasting `real` keeps the Complex{T} eltype for Fourier objects
            @test (@inferred broadcast(real, (F{P}(args...)))) isa F{P}
            if eltype(args[1]) <: Complex
                @test (@inferred F{P}(map(real,args)...)) isa F{P}
            end

        end
    end

end

##

@testset "Algebra" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        fs = ((B0,f0),(B2,f2),(Bt,ft)) = [
            (Fourier,   maybegpu(FlatMap(rand(Nside2D...)))),
            (EBFourier, maybegpu(FlatQUMap(rand(Nside2D...),rand(Nside2D...)))), # named FieldTuple
            (Fourier,   maybegpu(FieldTuple(FlatMap(rand(Nside2D...)),FlatMap(rand(Nside2D...))))), # unnamed FieldTuple
            (BasisTuple{Tuple{Fourier,EBFourier}}, maybegpu(FlatIQUMap(rand(Nside2D...),rand(Nside2D...),rand(Nside2D...)))), # named nested FieldTuple
            (BasisTuple{Tuple{Fourier,EBFourier}}, maybegpu(FieldTuple(FlatMap(rand(Nside2D...)),FlatQUMap(rand(Nside2D...),rand(Nside2D...))))), # unnamed nested FieldTuple
            (Fourier,   maybegpu(FlatMap(rand(Nside2D...,2)))), # batched S0
            (EBFourier, maybegpu(FlatQUMap(rand(Nside2D...,2),rand(Nside2D...,2)))), # batched S2
        ]

        @testset "f :: $(typeof(f))" for (B,f) in fs

            local Ðf, Ðv, g, H

            @test similar(f) isa typeof(f)
            @test zero(f) isa typeof(f)
            @test similar(f,Float32) isa Field
            @test eltype(similar(f,Float32)) == Float32

            # broadcasting
            @test (@inferred f .+ f) isa typeof(f)
            @test (@inferred f .+ Float32.(f)) isa typeof(f)

            # promotion
            @test (@inferred f + B(f)) isa typeof(f)
            @test (@inferred f + B(Float32.(f))) isa typeof(f)

            # gradients
            @test (Ðf = @inferred ∇[1]*f) isa Field
            @test (∇[1]'*f ≈ -∇[1]*f)
            @test (-∇[1]'*f ≈ ∇[1]*f)
            @test (@inferred mul!(Ðf,∇[1],Ð(f))) isa Field
            @test (Ðv = @inferred ∇*f) isa FieldVector
            @test (@inferred mul!(Ðv,∇,Ð(f))) isa FieldVector
            @test ((g,H) = map(Ł, (@inferred gradhess(f)))) isa NamedTuple{<:Any, <:Tuple{FieldVector, FieldMatrix}}

            # Diagonal broadcasting
            @test (@inferred Diagonal(f) .* Diagonal(f) .* Diagonal(f)) isa typeof(Diagonal(f))

            # inverses
            @test (@inferred pinv(Diagonal(f))) isa Diagonal{<:Any,<:typeof(f)}
            @test_throws SingularException inv(Diagonal(0*f))

            # Field dot products
            D = Diagonal(f)
            if f isa FlatField && batchsize(f)>1 # batched fields not inferred
                @test (f' * f) isa Real
                @test (f' * B(f)) isa Real
                @test (f' * D * f) isa Real
            else
                @test (@inferred f' * f) isa Real
                @test (@inferred f' * B(f)) isa Real
                @test (@inferred f' * D * f) isa Real
            end
            @test sum(f, dims=:) ≈ sum(f[:])
            @test_throws Any sum(f, dims=1)
            @test sum(f, dims=2) == f

            if f isa FlatS0
                # FieldVector dot product
                @test (@inferred Diagonal.(g)' * g) isa typeof(g[1])
                @test (@inferred mul!(similar(g[1]), Diagonal.(g)', g)) isa typeof(g[1])

                # FieldMatrix-FieldVector product
                @test (@inferred Diagonal.(H) * g) isa FieldVector
                @test (@inferred Diagonal.(H) * Diagonal.(g)) isa FieldOrOpVector
                @test (@inferred mul!(Diagonal.(similar.(g)), Diagonal.(H), Diagonal.(g))) isa FieldOrOpVector

            end

            # Explicit vs. lazy DiagOp algebra
            @test (Diagonal(Ð(f)) + Diagonal(Ð(f))) isa DiagOp{<:Field{basis(Ð(f))}}
            @test (Diagonal(Ł(f)) + Diagonal(Ð(f))) isa LazyBinaryOp
            @test (Diagonal(Ł(f)) + Diagonal(Ł(f))) isa DiagOp{<:Field{basis(Ł(f))}}

            # tuple adjoints
            f0b = identity.(batch(f0,batchsize(f)))
            v = similar.(@SVector[f0b,f0b])
            @test (@inferred mul!(f0b, tuple_adjoint(f), f)) isa Field{<:Any,S0}
            @test (@inferred mul!(v, tuple_adjoint(f), @SVector[f,f])) isa FieldVector{<:Field{<:Any,S0}}

        end

        # mixed-spin
        @test (@inferred f0 .* f2) isa typeof(f2)
        @test (@inferred f0 .* ft) isa typeof(ft)

        # matrix type promotion
        @test_broken (@inferred FlatMap(rand(Float32,2,2)) + FlatMap(spzeros(Float64,2,2))) isa FlatMap{<:Any,Float64,<:Matrix}

    end

end

##

@testset "FlatS2" begin
    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides
        C = maybegpu(Diagonal(EBFourier(FlatEBMap(rand(Nside2D...), rand(Nside2D...)))))
        f = maybegpu(FlatQUMap(rand(Nside2D...), rand(Nside2D...)))
        @test C*f ≈ FlatQUFourier(C[:QQ]*f[:Q]+C[:QU]*f[:U], C[:UU]*f[:U]+C[:UQ]*f[:Q])
    end
end

##

@testset "FlatS02" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        ΣTT, ΣTE, ΣEE, ΣBB = [Diagonal(Fourier(maybegpu(FlatMap(rand(Nside2D...))))) for i=1:4]
        L = FlatIEBCov(@SMatrix([ΣTT ΣTE; ΣTE ΣEE]), ΣBB)
        f = maybegpu(IEBFourier(FlatIEBMap(rand(Nside2D...),rand(Nside2D...),rand(Nside2D...))))

        @test (sqrt(L) * @inferred(@inferred(sqrt(L)) * f)) ≈ (L * f)
        @test (L * @inferred(@inferred(pinv(L)) * f)) ≈ f
        @test @inferred(L * L) isa FlatIEBCov
        @test @inferred(L + L) isa FlatIEBCov
        @test L * Diagonal(f) isa FlatIEBCov
        @test Diagonal(f) * L isa FlatIEBCov
        @test_broken @inferred L * Diagonal(f)
        @test @inferred(diag(L)) isa FlatIEBFourier
        @test @inferred(L + I) isa FlatIEBCov
        @test @inferred(2 * L) isa FlatIEBCov
        @test @inferred(similar(L)) isa FlatIEBCov
        @test (L .= 2L) isa FlatIEBCov

    end

end

##

@testset "BatchedReal" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        Nside2D = Nside .* (1,1)

        r  = 1.
        rb = batch([1.,2])

        @testset "f :: $(typeof(f))" for (f,fb) in [
            (maybegpu(FlatMap(rand(Nside2D...))),                    maybegpu(FlatMap(rand(Nside2D...,2)))),
            (maybegpu(FlatQUMap(rand(Nside2D...),rand(Nside2D...))), maybegpu(FlatQUMap(rand(Nside2D...,2),rand(Nside2D...,2))))
        ]

            @test @inferred(r * f)  == f
            @test @inferred(r * fb) == fb
            @test unbatch(@inferred(rb * f)) == [f, 2f]
            @test unbatch(@inferred(rb * fb)) == [batchindex(fb,1), 2batchindex(fb,2)]

        end

    end

end

##

@testset "Gradients" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        @test (@inferred ∇[1] * maybegpu(FlatMap(rand(Nside2D...), ∂mode=fourier∂))) isa FlatFourier
        @test (@inferred ∇[1] * maybegpu(FlatQUMap(rand(Nside2D...), rand(Nside2D...), ∂mode=fourier∂))) isa FlatQUFourier
        @test (@inferred ∇[1] * maybegpu(FlatIQUMap(rand(Nside2D...), rand(Nside2D...), rand(Nside2D...), ∂mode=fourier∂))) isa FlatIQUFourier

        @test (@inferred ∇[1] * Fourier(FlatMap(rand(Nside2D...), ∂mode=map∂))) isa FlatMap
        @test (@inferred ∇[1] * QUFourier(FlatQUMap(rand(Nside2D...), rand(Nside2D...), ∂mode=map∂))) isa FlatQUMap
        @test (@inferred ∇[1] * BasisTuple{Tuple{Fourier,QUFourier}}(FlatIQUMap(rand(Nside2D...), rand(Nside2D...), rand(Nside2D...), ∂mode=map∂))) isa FlatIQUMap

    end

end

##

@testset "Misc" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        f = maybegpu(FlatMap(rand(Nside2D...)))

        @test                  @inferred(MidPass(100,200) .* Diagonal(Fourier(f))) isa Diagonal
        @test_throws Exception           MidPass(100,200) .* Diagonal(        f)

    end

end

##

@testset "Cℓs" begin

    @test InterpolatedCℓs(1:100, rand(100))      isa AbstractCℓs{Float64}
    @test InterpolatedCℓs(1:100, rand(100) .± 1) isa AbstractCℓs{Measurement{Float64}}
    @test (InterpolatedCℓs(1:100, 1:100) * ℓ²)[50] == 50^3

end

##

@testset "ParamDependentOp" begin

    D = Diagonal(maybegpu(FlatMap(rand(4,4))))

    @test ParamDependentOp((;x=1, y=1)->x*y*D)() ≈ D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)(z=2) ≈ D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)(x=2) ≈ 2D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)((x=2,y=2)) ≈ 4D # tuple calling form
    @test_throws MethodError ParamDependentOp((;x=1, y=1)->x*y*D)(2) # only Tuple unnamed arg is OK

end

##

@testset "Chains" begin

    chains = CMBLensing.wrap_chains([
        [Dict(:i=>1, :b=>2), Dict(:i=>2       ), Dict(:i=>3, :b=>2)],
        [Dict(:i=>1, :b=>3), Dict(:i=>2, :b=>3), Dict(:i=>3, :b=>3)],
    ])

    # basic
    @test chains[1, 1, :i] == 1
    @test chains[:, 1, :i] == [1,1]
    @test chains[:, :, :i] == [[1, 2, 3], [1, 2, 3]]

    # slices
    @test chains[1, 1:2, :i] == [1, 2]
    @test chains[:, 1:2, :i] == [[1,2], [1,2]]

    # implied : in first dims
    @test chains[:i] == [[1,2,3],[1,2,3]]
    @test chains[1,:i] == [1,2,3]

    # missing
    @test all(chains[1,:b] .=== [2, missing, 2])

end;

##

@testset "Zygote" begin

    @testset "Nside = $Nside" for (Nside,Nside2D) in Nsides

        @testset "$(typeof(f))" for (f,g,h) in [
            @repeated(maybegpu(FlatMap(rand(Nside2D...))),3),
            @repeated(maybegpu(FlatQUMap(rand(Nside2D...),rand(Nside2D...))),3)
        ]

            v = @SVector[f,f]
            D = Diagonal(f)

            @testset "Fields" begin

                @testset "sum" begin
                    @test ((δ = gradient(f -> sum(Map(f)),     Map(f))[1]); basis(δ)==basis(Map(f))     && δ ≈ one(Map(f)))
                    @test ((δ = gradient(f -> sum(Map(f)), Fourier(f))[1]); basis(δ)==basis(Fourier(f)) && δ ≈ Fourier(one(Map(f))))
                end

                @testset "B1=$B1, B2=$B2, B3=$B3" for B1=[Map,Fourier], B2=[Map,Fourier], B3=[Map,Fourier]

                    @test gradient(f -> dot(B1(f),B2(f)), B3(f))[1] ≈ 2f
                    @test gradient(f -> norm(B1(f)), B3(f))[1] ≈ f/norm(f)
                    @test gradient(f -> B1(f)' * B2(g), B3(f))[1] ≈ g
                    @test gradient(f -> sum(Diagonal(Map(f)) * B2(g)), B3(f))[1] ≈ g
                    @test gradient(f -> sum(Diagonal(Map(∇[1]*f)) * B2(g)), B3(f))[1] ≈ ∇[1]'*g
                    @test gradient(f -> B1(f)'*(D\B2(f)), B3(f))[1] ≈ D\f + D'\f
                    @test gradient(f -> (B1(f)'/D)*B2(f), B3(f))[1] ≈ D\f + D'\f
                    @test gradient(f -> B1(f)'*(D*B2(f)), B3(f))[1] ≈ D*f + D'*f
                    @test gradient(f -> (B1(f)'*D)*B2(f), B3(f))[1] ≈ D*f + D'*f
                    if B2==Map
                        @test gradient(f -> B1(f)'*Diagonal(B2(f))*f, B3(f))[1] ≈ @. 3*$B2(f)^2
                    else
                        @test_broken gradient(f -> B1(f)'*Diagonal(B2(f))*f, B3(f))[1] ≈ @. 3*$B2(f)^2
                    end

                end

                @testset "Broadcasting" begin
                    @test gradient(f -> sum(@. f*f + 2*f + 1), f)[1] ≈ 2*f+2
                    @test gradient(f -> sum(@. f^2 + 2*f + 1), f)[1] ≈ 2*f+2
                end

            end

            @testset "FieldVectors" begin

                @test gradient(f -> Map(∇[1]*f)' *     Map(v[1]) + Map(∇[2]*f)' *     Map(v[2]), f)[1] ≈ ∇' * v
                @test gradient(f -> Map(∇[1]*f)' * Fourier(v[1]) + Map(∇[2]*f)' * Fourier(v[2]), f)[1] ≈ ∇' * v
                @test gradient(f -> sum(Diagonal(Map(∇[1]*f)) * v[1] + Diagonal(Map(∇[2]*f)) * v[2]), f)[1] ≈ ∇' * v

                @test gradient(f -> @SVector[f,f]' * Map.(@SVector[g,g]), f)[1] ≈ 2g
                @test gradient(f -> @SVector[f,f]' * Fourier.(@SVector[g,g]), f)[1] ≈ 2g

                @test gradient(f -> sum(Diagonal.(Map.(∇*f))' * Fourier.(v)), f)[1] ≈ ∇' * v
                @test gradient(f -> sum(Diagonal.(Map.(∇*f))' * Map.(v)), f)[1] ≈ ∇' * v

                @test gradient(f -> sum(sum(@SVector[f,f])),                            f)[1] ≈ 2*one(f)
                @test gradient(f -> sum(sum(@SVector[f,f]      .+ @SVector[f,f])),      f)[1] ≈ 4*one(f)
                @test gradient(f -> sum(sum(@SMatrix[f f; f f] .+ @SMatrix[f f; f f])), f)[1] ≈ 8*one(f)

                @test gradient(f -> sum(sum(Diagonal.(@SMatrix[f f; f f]) * @SVector[f,f])), f)[1] ≈ 8*f
                @test_broken gradient(f -> sum(sum(@SMatrix[f f] * @SMatrix[f f; f f])), f)[1] ≈ 8*f

            end

            @testset "OuterProdOp" begin

                @test OuterProdOp(f,g) * h ≈ f*(g'*h)
                @test OuterProdOp(f,g)' * h ≈ g*(f'*h)
                @test diag(OuterProdOp(f,g)) ≈ f .* conj.(g)
                @test diag(OuterProdOp(f,g)') ≈ conj.(f) .* g
                @test diag(OuterProdOp(f,g) + OuterProdOp(f,g)) ≈ 2 .* f .* conj.(g)

            end

            if f isa FlatS0

                @testset "logdet" begin
                    @test gradient(x->logdet(x*Diagonal(Map(f))),     1)[1] ≈ size(Map(f))[1]
                    @test gradient(x->logdet(x*Diagonal(Fourier(f))), 1)[1] ≈ size(Map(f))[1]
                    L = ParamDependentOp((;x=1)->x*Diagonal(Fourier(f)))
                    @test gradient(x->logdet(L(x=x)), 1)[1] ≈ size(Map(f))[1]
                end

                @test gradient(x -> norm(x*Fourier(f)), 1)[1] ≈ norm(f)
                @test gradient(x -> norm(x*Map(f)), 1)[1] ≈ norm(f)

                L₀ = Diagonal(Map(f))
                @test gradient(x -> norm((x*L₀)*f), 1)[1] ≈ norm(L₀*f)
                L₀ = Diagonal(Fourier(f))
                @test gradient(x -> norm((x*L₀)*f), 1)[1] ≈ norm(L₀*f)

                @test gradient(x -> norm((x*Diagonal(Map(f)))*f), 1)[1] ≈ norm(Diagonal(Map(f))*f)
                @test gradient(x -> norm((x*Diagonal(Fourier(f)))*f), 1)[1] ≈ norm(Diagonal(Fourier(f))*f)

            end

        end

    end

    @testset "LinearInterpolation" begin
        @test gradient(x->LinearInterpolation([1,2,3],[1,2,3])(x), 2)[1] == 1
        @test gradient(x->LinearInterpolation([1,2,3],[1,x,3])(2), 2)[1] == 1
        @test gradient(x->LinearInterpolation([1,x,3],[1,2,3])(2), 2)[1] == -1
    end

end

##

@testset "Lensing" begin

    local f,ϕ

    Cℓ = camb().unlensed_total
    seed!(0)

    @testset "Nside = $Nside" for (Nside,) in Nsides′

        @testset "T :: $T" for T in (Float32, Float64)

            ε = sqrt(eps(T))
            Cϕ = maybegpu(Cℓ_to_Cov(Flat(Nside=Nside), T, S0, Cℓ.ϕϕ))
            @test (ϕ = @inferred simulate(Cϕ)) isa FlatS0
            Lϕ = LenseFlow(ϕ)

            ## S0
            Cf = maybegpu(Cℓ_to_Cov(Flat(Nside=Nside), T, S0, Cℓ.TT))
            @test (f = @inferred simulate(Cf)) isa FlatS0
            @test (@inferred Lϕ*f) isa FlatS0
            # adjoints
            f,g = simulate(Cf),simulate(Cf)
            @test f' * (Lϕ * g) ≈ (f' * Lϕ) * g
            # gradients
            δf, δϕ = simulate(Cf), simulate(Cϕ)
            @test FieldTuple(gradient((f′,ϕ) -> f'*(LenseFlow(ϕ)*f′), f, ϕ))' * FieldTuple(δf,δϕ) ≈
                (f'*((LenseFlow(ϕ+ε*δϕ)*(f+ε*δf))-(LenseFlow(ϕ-ε*δϕ)*(f-ε*δf)))/(2ε)) rtol=1e-2

            # S2 lensing
            Cf = maybegpu(Cℓ_to_Cov(Flat(Nside=Nside), T, S2, Cℓ.EE, Cℓ.BB))
            @test (f = @inferred simulate(Cf)) isa FlatS2
            @test (@inferred Lϕ*f) isa FlatS2
            # adjoints
            f,g = simulate(Cf),simulate(Cf)
            @test f' * (Lϕ * g) ≈ (f' * Lϕ) * g
            # gradients
            δf, δϕ = simulate(Cf), simulate(Cϕ)
            @test FieldTuple(gradient((f′,ϕ) -> f'*(LenseFlow(ϕ)*f′), f, ϕ))' * FieldTuple(δf,δϕ) ≈
                (f'*((LenseFlow(ϕ+ε*δϕ)*(f+ε*δf))-(LenseFlow(ϕ-ε*δϕ)*(f-ε*δf)))/(2ε)) rtol=1e-2

        end

    end

end

##

@testset "Posterior" begin

    Cℓ = camb()
    L = LenseFlow{RK4Solver{7}}
    T = Float64

    @testset "Nside = $Nside" for (Nside,) in Nsides′

        @testset "pol = $pol" for pol in (:I,:P)

            @unpack f,f̃,ϕ,ds,ds₀ = load_sim(
                seed     = 0,
                Cℓ       = Cℓ,
                θpix     = 3,
                Nside    = Nside,
                T        = T,
                beamFWHM = 3,
                pol      = pol,
                L        = L,
                storage  = storage,
                pixel_mask_kwargs = (edge_padding_deg=1,)
            )
            @unpack Cf,Cϕ,D = ds₀
            f° = L(ϕ)*D*f

            @test lnP(0,f,ϕ,ds) ≈ lnP(1,    f̃,  ϕ ,ds) rtol=1e-4
            @test lnP(0,f,ϕ,ds) ≈ lnP(:mix, f°, ϕ, ds) rtol=1e-4

            ε = sqrt(eps(T))
            seed!(1)
            δf,δϕ = simulate(Cf),simulate(Cϕ)

            @test FieldTuple(gradient((f,ϕ)->lnP(0,f,ϕ,ds),f,ϕ))'*FieldTuple(δf,δϕ) ≈
                (lnP(0,f+ε*δf,ϕ+ε*δϕ,ds)-lnP(0,f-ε*δf,ϕ-ε*δϕ,ds))/(2ε)  rtol=3e-2
            @test FieldTuple(gradient((f̃,ϕ)->lnP(1,f̃,ϕ,ds),f̃,ϕ))'*FieldTuple(δf,δϕ) ≈
                (lnP(1,f̃+ε*δf,ϕ+ε*δϕ,ds)-lnP(1,f̃-ε*δf,ϕ-ε*δϕ,ds))/(2ε)  rtol=3e-2
            @test FieldTuple(gradient((f°,ϕ)->lnP(:mix,f°,ϕ,ds),f°,ϕ))'*FieldTuple(δf,δϕ) ≈
                (lnP(:mix,f°+ε*δf,ϕ+ε*δϕ,ds)-lnP(:mix,f°-ε*δf,ϕ-ε*δϕ,ds))/(2ε)  rtol=3e-2

        end

    end

end

##

end


text:
using GLMakie, GeometryBasics, RPRMakie, RadeonProRender
using Colors, FileIO
using Colors: N0f8

f = (u, v) -> cos(v) * (6 - (5 / 4 + sin(3 * u)) * sin(u - 3 * v))
g = (u, v) -> sin(v) * (6 - (5 / 4 + sin(3 * u)) * sin(u - 3 * v))
h = (u, v) -> -cos(u - 3 * v) * (5 / 4 + sin(3 * u));
u = range(0; stop=2π, length=150)
v = range(0; stop=2π, length=150)
radiance = 500
lights = [EnvironmentLight(1.0, load(RPR.assetpath("studio026.exr"))),
          PointLight(Vec3f(10), RGBf(radiance, radiance, radiance * 1.1))]

fig = Figure(; resolution=(1500, 1000))
ax = LScene(fig[1, 1]; show_axis=false, scenekw=(lights=lights,))
screen = RPRMakie.RPRScreen(size(ax.scene); plugin=RPR.Tahoe)
material = RPR.UberMaterial(screen.matsys)

surface!(ax, f.(u, v'), g.(u, v'), h.(u, v'); ambient=Vec3f(0.5), diffuse=Vec3f(1), specular=0.5,
         colormap=:balance, material=material)

function Input(fig, val::RGB)
    hue = Slider(fig; range=1:380, width=200)
    lightness = Slider(fig; range=LinRange(0, 1, 100), width=200)
    labels = [Label(fig, "hue"; halign=:left), Label(fig, "light"; halign=:left)]
    layout = grid!(hcat(labels, [hue, lightness]))
    hsl = HSL(val)
    set_close_to!(hue, hsl.h)
    set_close_to!(lightness, hsl.l)
    color = map((h, l) -> RGB(HSL(h, 0.9, l)), hue.value, lightness.value)
    return color, layout
end

function Input(fig, val::Vec4)
    s = Slider(fig; range=LinRange(0, 1, 100), width=200)
    set_close_to!(s, first(val))
    return map(x -> Vec4f(x), s.value), s
end

function Input(fig, val::Bool)
    toggle = Toggle(fig; active=val)
    return toggle.active, toggle
end

sliders = (reflection_color=Input(fig, RGB(0, 0, 0)), reflection_weight=Input(fig, Vec4(0)),
           reflection_roughness=Input(fig, Vec4(0)), reflection_anisotropy=Input(fig, Vec4(0)),
           reflection_anisotropy_rotation=Input(fig, Vec4(0)), reflection_mode=Input(fig, Vec4(0)),
           reflection_ior=Input(fig, Vec4(0)), reflection_metalness=Input(fig, Vec4(0)),
           refraction_color=Input(fig, RGB(0, 0, 0)), refraction_weight=Input(fig, Vec4(0)),
           refraction_roughness=Input(fig, Vec4(0)), refraction_ior=Input(fig, Vec4(0)),
           refraction_absorption_color=Input(fig, RGB(0, 0, 0)),
           refraction_absorption_distance=Input(fig, Vec4(0)), refraction_caustics=Input(fig, true),
           sss_scatter_color=Input(fig, RGB(0, 0, 0)), sss_scatter_distance=Input(fig, Vec4(0)),
           sss_scatter_direction=Input(fig, Vec4(0)), sss_weight=Input(fig, Vec4(0)),
           sss_multiscatter=Input(fig, false), backscatter_weight=Input(fig, Vec4(0)),
           backscatter_color=Input(fig, RGB(0, 0, 0)))

labels = []
inputs = []
refresh = Observable(nothing)
for (key, (obs, input)) in pairs(sliders)
    push!(labels, Label(fig, string(key); align=:left))
    push!(inputs, input)
    on(obs) do value
        @show key value
        setproperty!(material, key, value)
        return notify(refresh)
    end
end

fig[1, 2] = grid!(hcat(labels, inputs); width=500)
GLMakie.activate!()

cam = cameracontrols(ax.scene)
cam.eyeposition[] = Vec3f(22, 0, 17)
cam.lookat[] = Vec3f(0, 0, -1)
cam.upvector[] = Vec3f(0, 0, 1)
cam.fov[] = 30

display(fig)

context, task = RPRMakie.replace_scene_rpr!(ax.scene, screen; refresh=refresh)

# Change light parameters interactively
begin
    lights[1].intensity[] = 1.5
    lights[2].radiance[] = RGBf(1000, 1000, 1000)
    lights[2].position[] = Vec3f(3, 10, 10)
    notify(refresh)
end


text:
# vim: set ft=julia sw=4 ts=4:
#
#
# A set of (micro-) benchmarks for the Julia programming language
#

module Benchmark

include("Fibonacci.jl")
#import Fibonacci

println("Fibonacci numbers")
println("-----------------")
@time println(lpad("fib_naive(35)", 16) * " = $(Fibonacci.fib_naive(35))")
@time println(lpad("fib(35)", 16) * " = $(Fibonacci.fib(35))")
@time println(lpad("fib(1000)", 16) * " = $(Fibonacci.fib(1000))")
println()

println("Perfect numbers")
println("---------------")
println(" todo...")
println()

println("Mandelbrot set")
println("--------------")
println(" todo...")
println()

println(" done.")

end # module

max_stars_repo_name: kkirstein/proglang-playground

max_issues_repo_name: kkirstein/proglang-playground

max_forks_repo_name: kkirstein/proglang-playground


text:
using CompScienceMeshes
using Test

fn = joinpath(dirname(@__FILE__),"mesa.msh")
m = CompScienceMeshes.read_gmsh_mesh(fn)

@test numcells(m) == 1132
@test numvertices(m) == 626
@test norm(vertices(m,200) - [0.454487,-0.0982151,-0.1]) < 1e-3


text:
@rule Transition(:out, Marginalisation) (q_in::Categorical, q_a::MatrixDirichlet) = begin
    a = clamp.(exp.(logmean(q_a) * probvec(q_in)), tiny, Inf)
    return Categorical(a ./ sum(a))
end

@rule Transition(:out, Marginalisation) (m_in::Categorical, q_a::MatrixDirichlet) = begin
    a = clamp.(exp.(logmean(q_a)) * probvec(m_in), tiny, Inf)
    return Categorical(a ./ sum(a))
end


text:
function parsecomment(s)

p = split(s, '=', keepempty=false);
n = length(p);

keys = Array{Any}(undef,n-1);
values = Array{Any}(undef,n-1);
for i = 1:n-1
    ii = findlast("\"", p[i]);
    if ii === nothing
        q = split(p[i], ' ', keepempty=false);
        if length(q) == 1
            keys[i] = q[1];
        elseif length(q)==2
            keys[i] = q[2];
            values[i-1] = q[1];
        end
    else
        t = split(p[i], '"', keepempty=false);
        q = split(t[1], ' ', keepempty=false);
        values[i-1] = map(x -> parse(Float64, x), q);
        keys[i] = t[2][2:end];
    end
end
t = split(p[n], '"', keepempty=false);
values[n-1] = t[1];

# handle energy
for i = 1:n-1
    if lowercase(string(keys[i])) == "energy"
        values[i] = parse(Float64, values[i]);
    end
end

# handle pbc
for i = 1:n-1
    if lowercase(string(keys[i])) == "pbc"
        values[i] = replace(values[i], "T" => "1");
        values[i] = replace(values[i], "F" => "0");
        q = split(values[i], ' ', keepempty=false);
        values[i] = map(x -> parse(Int64, x), q);
    end
end

# handle properties
indp = 0;
for i = 1:n-1
    if lowercase(string(keys[i])) == "properties"
        values[i] = split(values[i], ':', keepempty=false);
        indp = i;
        break;
    end
end

return keys, values, indp

end

function readEXTXYZ(filename, species)

config = initializeconfig(0);

f = readlines(filename);
n = length(f);

natomall = 0;
nframe = 0;
m = 0;
while (m<n)
    m = m + 1;
    nframe = nframe + 1
    natom = parse(Int64, f[m])
    config.natom = hcat(config.natom, reshape([natom],1,1))

    m = m + 1
    keys, values, indp = parsecomment(f[m]);
    for i = 1:length(keys)
        mystr = lowercase(string(keys[i]));
        if mystr == "lattice"
            d = reshape(values[i], length(values[i]), 1);
            if nframe == 1
                config.lattice = d;
            else
                config.lattice = hcat(config.lattice, d);
            end
            config.ncl = length(values[i]);
        elseif mystr == "energy"
            d = reshape([values[i]], 1, 1);
            if nframe == 1
                config.e = d;
            else
                config.e = hcat(config.e, d);
            end
            config.nce = 1;
        elseif mystr == "stress"
            d = reshape(values[i], length(values[i]), 1);
            if nframe == 1
                config.stress = d;
            else
                config.stress = hcat(config.stress, d);
            end
            config.ncs = length(values[i]);
        elseif mystr == "pbc"
            d = reshape(values[i], length(values[i]), 1);
            if nframe == 1
                config.pbc = d;
            else
                config.pbc = hcat(config.pbc, d);
            end
            config.ncp = length(values[i]);
        end
    end

    s = String.(reduce(vcat, split.(f[(m+1):(m+natom)])))
    s = reshape(s, (Int64(length(s)/natom), natom))
    m = m + natom

    # display(keys)
    # display(values)
    # display(indp)
    # display(values[indp])
    # handle properties
    prop = reshape(values[indp],(3, Int64(length(values[indp])/3)));
    nfield = size(prop,2);
    fieldtype = String.(prop[2,:]);
    fieldlength = parse.(Int64, prop[3,:]);
    fieldlength = Int64.(cumsum(fieldlength[:]));

    frame = Array{Any}(undef,nfield);
    for j = 1:nfield
        if j == 1
            k1 = 1;
            k2 = fieldlength[1];
        else
            k1 = fieldlength[j-1]+1;
            k2 = fieldlength[j];
        end
        if lowercase(fieldtype[j]) == "s" || lowercase(fieldtype[j]) == "l"
            frame[j] = s[k1,:];
            b = ones(Int64,1,natom);
            if lowercase(fieldtype[j]) == "l"
                ind = findall(frame[j][:] .== "F");
                b[ind] .= 0;
                frame[j] = b;
            else
                for k = 1:length(species)
                    ind = findall(frame[j][:] .== species[k]);
                    b[ind] .= k;
                end
                frame[j] = b;
            end
        else
           frame[j] = parse.(Float64, s[k1:k2,:]);
        end
    end

    for j = 1:nfield
        mystr = lowercase(prop[1,j]);
        if mystr == "species"
            if nframe == 1
                config.t = frame[j];
            else
                config.t = hcat(config.t, frame[j]);
            end
            config.nct = size(config.t,1);
        elseif mystr == "pos"
            if nframe == 1
                config.x = frame[j];
            else
                config.x = hcat(config.x, frame[j]);
            end
            config.ncx = size(config.x,1);
            config.dim = size(config.x,1);
        elseif mystr == "vel"
            if nframe == 1
                config.v = frame[j];
            else
                config.v = hcat(config.v, frame[j]);
            end
            config.ncv = size(config.v,1);
        elseif mystr == "forces"
            if nframe == 1
                config.f = frame[j];
            else
                config.f = hcat(config.f, frame[j]);
            end
            config.ncf = size(config.f,1);
        elseif mystr == "charge"
            if nframe == 1
                config.q = frame[j];
            else
                config.q = hcat(config.q, frame[j]);
            end
            config.ncq = size(config.q,1);
        elseif mystr == "mass"
            if nframe == 1
                config.mass = frame[j];
            else
                config.mass = hcat(config.mass, frame[j]);
            end
            config.ncm = size(config.mass,1);
        elseif mystr == "tags"
            if nframe == 1
                config.tags = frame[j];
            else
                config.tags = hcat(config.tags, frame[j]);
            end
            config.nci = size(config.tags,1);
        elseif mystr == "groups"
            if nframe == 1
                config.group = frame[j];
            else
                config.group = hcat(config.group, frame[j]);
            end
            config.ncg = size(config.group,1);
        elseif mystr == "move_mask"
            if nframe == 1
                config.move = frame[j];
            else
                config.move = hcat(config.move, frame[j]);
            end
            config.nco = size(config.move,1);
        elseif mystr == "z"
            if nframe == 1
                config.Z = frame[j];
            else
                config.Z = hcat(config.Z, frame[j]);
            end
            config.ncz = size(config.Z,1);
        end
    end
end

config.nconfigs = length(config.natom)

return config;

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using URLParse
using Base.Test

#############################
# generate test code
# To make it easy to write regression tests:
# Add a URL into urls
# Run test.jl (this file) with gen in the command line and redirect output to "do_gen_tests.jl"
# Verify that the generated test code is ok.
# Run test.jl anytime without any parameters to run regression tests
function gen_test(url::String)
    up = urlparse(url)
    println("    println(\"testing url: $(url)\")")
    println("    up = urlparse(\"$(url)\")")
    println("    @test up.scheme == \"$(up.scheme)\"")
    println("    @test up.netloc == \"$(up.netloc)\"")
    println("    @test up.url == \"$(up.url)\"")
    println("    @test up.params == \"$(up.params)\"")
    println("    @test up.query == \"$(up.query)\"")
    println("    @test up.fragment == \"$(up.fragment)\"")
    println("    @test string(hostname(up)) == \"$(hostname(up))\"")
    println("    @test string(port(up)) == \"$(port(up))\"")
    println("    @test string(username(up)) == \"$(username(up))\"")
    println("    @test string(password(up)) == \"$(password(up))\"")
    println("    @test urlunparse(up) == \"$(url)\"")
    println("")
end

const urls = ["hdfs://user:password@hdfshost:9000/root/folder/file.csv",
    "https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag",
    "ftp://ftp.is.co.za/rfc/rfc1808.txt", "http://www.ietf.org/rfc/rfc2396.txt",
    "ldap://[2001:db8::7]/c=GB?objectClass?one", "mailto:John.Doe@example.com",
    "news:comp.infosystems.www.servers.unix", "tel:+1-816-555-1212", "telnet://192.0.2.16:80/",
    "urn:oasis:names:specification:docbook:dtd:xml:4.1.2"]

function gen_tests()
    println("function do_gen_tests()")
    for url in urls
        gen_test(url)
    end
    println("end");
end

include("do_gen_tests.jl")

const do_gen = isempty(ARGS) ? false : ("gen" == ARGS[1])

if(do_gen)
    gen_tests()
else
    do_gen_tests()
end