text:
include("../Random-Pertubation-on-Constrained-Manifolds/RandPertOnConstraindManifolds.jl")
include("../../example_cases/BasinStabilityForPD/pd_basin_stability.jl")
op = find_steady_state(powergrid)
rpg = rhs(powergrid)

P, Q = SimplePowerFlow(powergrid, op)

ω_idx = findall(:ω .∈ symbolsof.(powergrid.nodes))
θ_idx = findall(:θ .∈ symbolsof.(powergrid.nodes))

PerturbedState, dz = RandPertWithConstrains(powergrid, op, 5,[-π, π],[-100,100] ,Q[3])

result = sim(powergrid, PerturbedState, (0.0, 20.0), true, rpg)

plot_res(result, powergrid, 3)

plot(PowerGridSolution(result.dqsol, powergrid), θ_idx, :θ, legend = false)

#ND.jl gives us an ODEFunction so we have to construct one manually
#rpg_ode = rhs(powergrid)

##
throws = 1000
node = 4
v_vec = Array{Float64}(undef, throws)
θ_vec = Array{Float64}(undef, throws)
ω_vec = Array{Float64}(undef, throws)
i_vec = Array{Float64}(undef, throws)

for i in 1:throws
        PerturbedState, dz  =  RandPertWithConstrains(powergrid, op, node,[-π, π],[-100,100] ,Q[node])
        v_vec[i] = abs(dz[3])
        θ_vec[i] = abs(dz[1])
        ω_vec[i] = PerturbedState[node,:ω]
        i_vec[i] = PerturbedState[node,:iabs]
end

scatter(θ_vec, v_vec, marker = (:circle, 8), label = "Optim Projection")
xlabel!("θ [rad]")
ylabel!("v [p.u.]")

png("Pertubations_Optim_neu_4")

op[3,:φ]
a = powergrid.nodes[2]

println(a.H)
#=
rr = RootRhs(rhs(powergrid))
rr = (dx,x) -> rpg(dx,x, nothing, 0)
xtol = 1E-9
nl_res = nlsolve(rr, op.vec)
dx = similar(nl_res.zero)

# Define a differential equation system
@parameters t
@variables θ(t) i(t) u(t) g(t)
@derivatives D'~t

# Lorenz attractor
eqs = [g ~ Q[1] - 3 * i * u * exp(-2im * θ)]
de = ODESystem(eqs)

Jg = transpose(calculate_jacobian(de))

println(Jg)
=#

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
module Reasoning

using LogicCircuits
using ..Probabilistic
using ..Logistic
using ..Utils

export
    UpExpFlow,
    ExpFlowΔ,
    exp_pass_up,
    Expectation,
    ExpectationUpward,
    Moment

include("Expectation.jl")
include("ExpFlowCircuits.jl")

end

ext: jl


text:
using CompositeProblems
using StructuredProximalOperators
using Random
using LinearAlgebra
using BenchmarkTools
using ForwardDiff

function check_problem_smoothoracles(pb, x, η)

    f_x = f(pb, x)
    ∇f_x = ∇f(pb, x)
    ∇²f_x_η = ∇²f_h(pb, x, η)

    ηgradfx = dot(η, ∇f_x)
    ηHessf_xη = dot(η, ∇²f_x_η)

    @show ηgradfx, ηHessf_xη

    fgrad(t) = abs(f(pb, x + t*η) - (f_x + t * ηgradfx))
    fhess(t) = abs(f(pb, x + t*η) - (f_x + t * ηgradfx + 0.5 * t^2 * ηHessf_xη))

    return compare_curves(fgrad, fhess)
end

function main()
    println("Least square test")
    n, m = 20, 15
    pb = get_random_qualifiedlasso(n, m, 0.5)

    Random.seed!(1234)
    x = rand(n)
    Random.seed!(4321)
    η = rand(n)

    res = check_problem_smoothoracles(pb, x, η)

    println("\nLogistic test")
    n, m = 200, 150
    pb = get_random_logit(n=n, m=m, sparsity=0.5)

    Random.seed!(1234)
    x = rand(n)
    Random.seed!(4321)
    η = rand(n)

    @show f(pb, x)
    @btime f($pb, $x)

    println()
    println()

    ∇fx_AD = ForwardDiff.gradient(x->f(pb, x), x)
    ∇fx_∇f = zeros(n); CompositeProblems.∇f!(pb, ∇fx_∇f, x)

    @show norm(∇fx_AD), norm(∇fx_∇f)
    @show norm(∇fx_AD - ∇fx_∇f)

    @btime CompositeProblems.∇f!($pb, $∇fx_∇f, $x)

    println()
    println()
    ∇²f_h_base = zeros(n)
    @btime CompositeProblems.∇²f_h!($pb, $∇²f_h_base, $x, $η)

    ∇²f_h_ADhess = ForwardDiff.hessian(x->f(pb, x), x) * η

    CompositeProblems.∇²f_h!(pb, ∇²f_h_base, x, η)

    @show norm(∇²f_h_ADhess)
    @show norm(∇²f_h_base)
    @show norm(∇²f_h_ADhess - ∇²f_h_base)

    @btime CompositeProblems.∇²f_h!($pb, $∇²f_h_base, $x, $η)

    # @assert false
    @btime CompositeProblems.∇²f_h!($pb, $∇²f_h_base, $x, $η)

    res = check_problem_smoothoracles(pb, x, η)
    return display_curvescomparison(res)
end

main()


text:
using SumBassSolver

p=0.001
q=0.1
m=100
c0=10
alpha = 1
F = 10
T = 50
n = 50
Pi = 600

x1,v1 = sumBassHomogeneousSolver(p,q,m,c0,alpha,F,T,n,Pi)

pp=fill(p,n)
qq=fill(q,n)
mm=fill(m,n)
cc0 = fill(c0,n)
alphas=fill(alpha,n)
FF=fill(F,n)
TT=fill(T,n)

# The general algorithm does find the right solution but takes really long
# The local optimization finishes in a local optimum
x2,v2 = @time sumBassSolver(pp,qq,mm,cc0,alphas,FF,TT,Pi)
x3,v3 = @time sumBassSolver(pp,qq,mm,cc0,alphas,FF,TT,Pi,alg="SB")
x4,v4 = @time sumBassSolver(pp,qq,mm,cc0,alphas,FF,TT,Pi,alg="SB",localopt=true)

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
information = Dict("NAME: " => "Rajvi Srivastav", "EMAIL: " => "rajvisrivastav18@gmail.com", "LANGUAGE: " => "JULIA", "BIOSTACK: " => "Genomics", "SLACK USERNAME: " => "@Rajvi")
 for (i, j) in information
       println(i, j)
       end


text:
typealias Trajdata ListRecord{VehicleState, VehicleDef, Int}
Trajdata(timestep::Float64) = ListRecord(timestep, VehicleState, VehicleDef, Int)
Base.show(io::IO, trajdata::Trajdata) = @printf(io, "Trajdata(%d frames)", nframes(trajdata))


text:
@inline function ExaTronKernels.dbreakpt(n::Int, x::ROCDeviceArray{Float64,1}, xl::ROCDeviceArray{Float64,1},
                          xu::ROCDeviceArray{Float64,1}, w::ROCDeviceArray{Float64,1})
    zero = 0.0
    nbrpt = 0
    brptmin = zero
    brptmax = zero

    @inbounds for i=1:n
        if (x[i] < xu[i] && w[i] > zero)
            nbrpt = nbrpt + 1
            brpt = (xu[i] - x[i]) / w[i]
            if nbrpt == 1
                brptmin = brpt
                brptmax = brpt
            else
                brptmin = min(brpt,brptmin)
                brptmax = max(brpt,brptmax)
            end
        elseif (x[i] > xl[i] && w[i] < zero)
            nbrpt = nbrpt + 1
            brpt = (xl[i] - x[i]) / w[i]
            if nbrpt == 1
                brptmin = brpt
                brptmax = brpt
            else
                brptmin = min(brpt,brptmin)
                brptmax = max(brpt,brptmax)
            end
        end
    end

    # Handle the exceptional case.

    if nbrpt == 0
        brptmin = zero
        brptmax = zero
    end
    AMDGPU.sync_workgroup()

    return nbrpt,brptmin,brptmax
end


text:
Pkg.clone("https://github.com/dfdx/NaiveBayes.jl.git")

using NaiveBayes
using RDatasets

iris = dataset("datasets", "iris")

# observations in columns and variables in rows
X = array(iris[:, 1:4])'
p, n = size(X)
# by default species is a PooledDataArray,
y = [species for species in iris[:, 5]]

# how much data use for training
train_frac = 0.9
k = int(floor(train_frac * n))
idxs = randperm(n)
train = idxs[1:k]
test = idxs[k+1:end]

model = GaussianNB(unique(y), p)
fit(model, X[:, train], y[train])

accuracy = countnz(predict(model, X[:,test]).==
            y[test]) / countnz(test)

println("Accuracy: $accuracy")

X = [1 1 0 2 1;
     0 0 3 1 0;
     1 0 1 0 2]

y = [:a, :b, :b, :a, :a]

m = MultinomialNB(unique(y), 3)
fit(m, X, y)

Xtest = [0 4 1;
         2 2 0;
         1 1 1]

predict(m, Xtest)

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
module Rods

function get_fn(params::Array)
    function parseParam(p::Dict)
        return Dict(
            :radius => float(p["radius"]),
            :k => float(p["k"]),
            :center => [float(x) for x in p["center"]],
            :power  => p["exponent"],
            :direction => normalize([float(x) for x in p["direction"]])
        )
    end

    function nomalize{T<:AbstractFloat}(xs::Array{T,1})
        if norm(xs) < 1e-7
            return xs
        else
            return xs / norm(xs)
        end
    end

    function distance{T<:AbstractFloat}(point::Array{T,1}, p::Dict)
        v = point - p[:center]
        return norm(v - dot(v, p[:direction])*p[:direction])
    end

    function value{T<:AbstractFloat}(coordinates::Array{T,1}, p::Dict)
        d = distance(coordinates, p)
        if d >= p[:radius]
            return 0.0
        else
            # the center should have the highest energy
            return 0.5 * p[:k] * (d - p[:radius])^p[:power]
        end
    end

    function compute_aux{T<:AbstractFloat}(coordinates::Array{T,1}, param_list::Array, accum::AbstractFloat)
        if length(param_list) == 0
            return accum
        else
            return compute_aux(coordinates, param_list[2:end], accum + value(coordinates, param_list[1]))
        end
    end

    function fn_compute(coordinates)
        return compute_aux(coordinates, [parseParam(p) for p in params], 0.0)
    end

    return fn_compute
end

export get_fn
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
#
# bobyqa.jl --
#
# Julia interface to Mike Powell's BOBYQA method.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2015-2019, Éric Thiébaut
# <https://github.com/emmt/OptimPackNextGen.jl>.
#

module Bobyqa

export
    bobyqa,
    bobyqa!

using Printf

import
    ..AbstractContext,
    ..AbstractStatus,
    ..getreason,
    ..getstatus,
    ..grow!,
    ..iterate,
    ..restart

# The dynamic library implementing the method.
import ..libbobyqa

# Status returned by most functions of the library.
struct Status <: AbstractStatus
    _code::Cint
end

const SUCCESS              = Status( 0)
const BAD_NVARS            = Status(-1)
const BAD_NPT              = Status(-2)
const BAD_RHO_RANGE        = Status(-3)
const BAD_SCALING          = Status(-4)
const TOO_CLOSE            = Status(-5)
const ROUNDING_ERRORS      = Status(-6)
const TOO_MANY_EVALUATIONS = Status(-7)
const STEP_FAILED          = Status(-8)

# Get a textual explanation of the status returned by BOBYQA.
function getreason(status::Status)
    ptr = ccall((:bobyqa_reason, libbobyqa), Ptr{UInt8}, (Cint,), status._code)
    if ptr == C_NULL
        error("unknown BOBYQA status: ", status._code)
    end
    unsafe_string(ptr)
end

"""
# Minimizing a function of many variables with bound constraints

Mike Powell's **BOBYQA** algorithm attempts to find the variables `x` which
solve the bound constrained problem:

    min f(x)  s.t.  xl ≤ x ≤ xu

where `x` is a vector of variables that has `n ≥ 2` components `f(x)` is an
objective function, `xl` and `xu` are bounds on the variables.  The algorithm
employs quadratic approximations to the objective which interpolate the
objective function at `m` points, the value `m = 2n + 1` being recommended.
The parameter `rho` controls the size of the trust region and it is reduced
automatically from `rhobeg` to `rhoend` (such that `0 < rhoend ≤ rhobeg`).

The in-place version of the algorithm is called as:

    Bobyqa.minimize!(f, x, xl, xu, rhobeg, rhoend) -> (status, x, fx)

where `f` is the objective function, `x` is a vector with the initial and final
variables, `xl` and `xu` the lower and upper bounds for the variables, `rhobeg`
and `rhoend` are the initial and final sizes of the trust region.  The result
is a tuple of 3 values: `status` indicates whether the algorithm was
successful, `x` is the final value of the variables and `fx = f(x)` is the
objective function at `x`.  Normally, `status` should be `Bobyqa.SUCCESS`;
otherwise, `getreason(status)` yields a textual explanation about the failure.
The initial variables must be feasible, that is:

    xl ≤ x ≤ xu

must hold on entry.

The method:

    Bobyqa.minimize(f, x0, xl, xu, rhobeg, rhoend) -> (status, x, fx)

is identical to `Bobyqa.minimize!` but does not modify the vector `x0` of
initial variables.

## Precision and scaling of variables

Parameter `rhobeg` should be set to the typical size (in terms of Euclidean
norm of the change of variables) of the region to explorate and `rhoend`
should be set to the typical precision.  The proper scaling of the variables is
important for the success of the algorithm and the optional `scale` keyword
should be specified if the typical precision is not the same for all variables.
If specified, `scale` is an array of strictly nonnegative values and of same
size as the variables `x`, such that `scale[i]*rho` (with `rho` the trust
region radius) is the size of the trust region for the `i`-th variable.  If
keyword `scale` is not specified, a unit scaling for all the variables is
assumed.

An error occurs if any of the differences `xu[i] - xl[i]` is less than
`2*rhobeg*scale[i]`.

## Keywords

The following keywords are available:

* `scale` specifies the typical magnitudes of the variables.  If specified, it
  must have as many elements as `x`, all strictly positive.  If not specified,
  `scale[i] = 1` is assumed for any `i ∈ 1:n`.

* `check` (`true` by default) specifies whether to throw an exception if the
  algorithm is not fully successful.

* `verbose` (`0` by default) set the amount of printing.

* `maxeval` specifies the maximum number of calls to the objective function.
  The default setting is `maxeval = 30n` with `n = length(x)` the number of
  variables.

* `npt` specifies the number of points to use for the quadratic approximation
  of the objective function.  The default setting is the recommended value:
  `npt = 2n + 1` with `n = length(x)` the number of variables.

* `work` specifies a workspace to (re)use.  It must be a vector of double
  precision floating-point values.  If it is too small, its size is
  automatically adjusted (by calling [`resize!`](@ref)).  This keyword is
  useful to avoid any new allocation (and garbage colection) when several
  similar optimizations are to be performed.

## References

The algorithm is described in:

* M.J.D. Powell, "The BOBYQA algorithm for bound constrained optimization
  without derivatives", Technical Report NA2009/06 of the Department of Applied
  Mathematics and Theoretical Physics, Cambridge, England (2009).

"""
minimize(args...; kwds...) = optimize(args...; maximize=false, kwds...)
minimize!(args...; kwds...) = optimize!(args...; maximize=false, kwds...)
@doc @doc(minimize) minimize!

"""

    Bobyqa.maximize(f, x0, xl, xu, rhobeg, rhoend) -> (status, x, fx)
    Bobyqa.maximize!(f, x, xl, xu, rhobeg, rhoend) -> (status, x, fx)

are similar to `Bobyqa.minimize` and `Bobyqa.minimize!` respectively but
solve the bound constrained maximization problem:

    max f(x)  s.t.  xl ≤ x ≤ xu

"""
maximize(args...; kwds...) = optimize(args...; maximize=true, kwds...)
maximize!(args...; kwds...) = optimize!(args...; maximize=true, kwds...)
@doc @doc(maximize) maximize!

# `_wrklen(...)` yields the number of elements in BOBYQA workspace.
_wrklen(n::Integer, npt::Integer) = _wrklen(Int(n), Int(npt))
_wrklen(n::Int, npt::Int) = (npt + 5)*(npt + n) + div(3*n*(n + 5),2)
_wrklen(x::AbstractVector{<:AbstractFloat}, npt::Integer) =
    _wrklen(length(x), npt)
function _wrklen(x::AbstractVector{<:AbstractFloat},
                 scl::AbstractVector{<:AbstractFloat},
                 npt::Integer)
    return _wrklen(x, npt) + 3*length(scl)
end

# `_work(...)` yields a large enough workspace for BOBYQA.
_work(x::AbstractVector{<:AbstractFloat}, npt::Integer) =
    Vector{Cdouble}(undef, _wrklen(x, npt))
function _work(x::AbstractVector{<:AbstractFloat},
               scl::AbstractVector{<:AbstractFloat},
               npt::Integer)
    return Vector{Cdouble}(undef, _wrklen(x, scl, npt))
end

# Wrapper for the objective function in BOBYQA, the actual objective function
# is provided by the client data as a `jl_value_t*` pointer.
function _objfun(n::Cptrdiff_t, xptr::Ptr{Cdouble}, fptr::Ptr{Cvoid})::Cdouble
    x = unsafe_wrap(Array, xptr, n)
    f = unsafe_pointer_to_objref(fptr)
    return Cdouble(f(x))
end

# With precompilation, `__init__()` carries on initializations that must occur
# at runtime like `@cfunction` which returns a raw pointer.
const _objfun_c = Ref{Ptr{Cvoid}}()
function __init__()
    _objfun_c[] = @cfunction(_objfun, Cdouble,
                             (Cptrdiff_t, Ptr{Cdouble}, Ptr{Cvoid}))
end

"""
The methods:

    Bobyqa.optimize(fc, x0, xl, xu, rhobeg, rhoend) -> (status, x, fx)
    Bobyqa.optimize!(fc, x, xl, xu, rhobeg, rhoend) -> (status, x, fx)

are identical to `Bobyqa.minimize` and `Bobyqa.minimize!` respectively but have
an additional `maximize` keyword which is `false` by default and which
specifies whether to maximize the objective function; otherwise, the method
attempts to minimize the objective function.

"""
optimize(f::Function, x0::AbstractVector{<:Real}, args...; kwds...) =
    optimize!(f, copyto!(Array{Cdouble}(undef, length(x0)), x0),
              args...; kwds...)

function optimize!(f::Function, x::DenseVector{Cdouble},
                   xl::DenseVector{Cdouble}, xu::DenseVector{Cdouble},
                   rhobeg::Real, rhoend::Real;
                   scale::DenseVector{Cdouble}=Array{Cdouble}(undef, 0),
                   maximize::Bool = false,
                   npt::Integer = 2*length(x) + 1,
                   check::Bool = false,
                   verbose::Integer = 0,
                   maxeval::Integer = 30*length(x),
                   work::Vector{Cdouble} = _work(x, scale, npt))
    n = length(x)
    length(xl) == n || error("bad length for inferior bound")
    length(xu) == n || error("bad length for superior bound")
    nscl = length(scale)
    if nscl == 0
        sclptr = Ptr{Cdouble}(0)
    elseif nscl == n
        sclptr = pointer(scale)
    else
        error("bad number of scaling factors")
    end
    grow!(work, _wrklen(x, scale, npt))
    status = Status(ccall((:bobyqa_optimize, libbobyqa), Cint,
                          (Cptrdiff_t, Cptrdiff_t, Cint, Ptr{Cvoid}, Any,
                           Ptr{Cdouble}, Ptr{Cdouble},
                           Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble,
                           Cptrdiff_t, Cptrdiff_t, Ptr{Cdouble}),
                          n, npt, (maximize ? Cint(1) : Cint(0)),
                          _objfun_c[], f, x, xl, xu, sclptr, rhobeg, rhoend,
                          verbose, maxeval, work))
    if check && status != SUCCESS
        error(getreason(status))
    end
    return (status, x, work[1])
end

@doc @doc(optimize) optimize!

# Basic version similar to the FORTRAN version.
bobyqa(f::Function, x0::DenseVector{Cdouble}, args...; kwds...) =
    bobyqa!(f, copy(x0), args...; kwds...)

function bobyqa!(f::Function, x::DenseVector{Cdouble},
                 xl::DenseVector{Cdouble}, xu::DenseVector{Cdouble},
                 rhobeg::Real, rhoend::Real;
                 npt::Integer = 2*length(x) + 1,
                 verbose::Integer = 0,
                 maxeval::Integer = 30*length(x),
                 check::Bool = true,
                 work::Vector{Cdouble} = _work(x, npt))
    n = length(x)
    length(xl) == n || error("bad length for inferior bound")
    length(xu) == n || error("bad length for superior bound")
    grow!(work, _wrklen(x, npt))
    status = Status(ccall((:bobyqa, libbobyqa), Cint,
                          (Cptrdiff_t, Cptrdiff_t, Ptr{Cvoid}, Any,
                           Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                           Cdouble, Cdouble, Cptrdiff_t, Cptrdiff_t,
                           Ptr{Cdouble}),
                          n, npt, _objfun_c[], f, x, xl, xu,
                          rhobeg, rhoend, verbose, maxeval, work))
    if check && status != SUCCESS
        error(getreason(status))
    end
    return (status, x, work[1])
end

end # module Bobyqa


text:
"""
``MolecularGeometry`` provides methods for
measuring the geometric properties of molecules.
The following quantities are supported:

+-----------------------------------+--------------------------------+
| Quantity                          |Method                          |
+===================================+================================+
|geometric center                   |gage(GeometricCenter, array)    |
|                                   |                                |
|                                   |gage(GeometricCenter, obj)      |
+-----------------------------------+--------------------------------+

"""
module MolecularGeometry
include(joinpath("..", "deps", "Distances", "src", "Distances.jl"))

macro debug(expression)
    return :($expression)
end

export GeometricCenter, MinDistance, PrincipalAxes
export RotationAngle, ProjectedRotationAngle
export gage, nearly_zero, smaller_vector

# types
include(joinpath("type", "all.jl"))

# functions
include(joinpath("function", "all.jl"))

end


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(coins::Vector{Int}, n::Int)::Int
    combinations = vcat([1], repeat([0], n))
    for i ∈ 1:length(coins)
        for j ∈ coins[i] + 1:n + 1
            combinations[j] += combinations[j - coins[i]]
        end
    end
    return combinations[end]
end

coin_list = [1, 2, 5, 10, 20, 50, 100, 200]

compute(coin_list, 200)

@benchmark compute(coin_list, 200)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
module TestAnalysis
using CustomTest
using Base.Test
using ControlSystems

## TZERO ##
# Examples from the Emami-Naeini & Van Dooren Paper
# Example 3
A = [0 1 0 0 0 0;
     0 0 1 0 0 0;
     0 0 0 0 0 0;
     0 0 0 0 1 0;
     0 0 0 0 0 1;
     0 0 0 0 0 0]
B = [0 0 1 0 0 0;
     0 0 0 0 0 1]'
C = [1 1 0 0 0 0;
     0 0 0 1 -1 0]
D = [1 0;
     1 0]

ex_3 = ss(A, B, C, D)
@test_approx_eq tzero(ex_3) [0.3411639019140099 + 1.161541399997252im,
                             0.3411639019140099 - 1.161541399997252im,
                             0.9999999999999999 + 0.0im,
                             -0.6823278038280199 + 0.0im]
# Example 4
A = [-0.129    0.0   0.396e-1  0.25e-1    0.191e-1;
     0.329e-2  0.0  -0.779e-4  0.122e-3  -0.621;
     0.718e-1  0.0  -0.1       0.887e-3  -0.385e1;
     0.411e-1  0.0   0.0      -0.822e-1   0.0;
     0.361e-3  0.0   0.35e-4   0.426e-4  -0.743e-1]
B = [0.0  0.139e-2;
     0.0  0.359e-4;
     0.0  -0.989e-2;
     0.249e-4  0.0;
     0.0  -0.534e-5]
C = [1 0 0 0 0;
     0 1 0 0 0]
D = zeros(2, 2)
ex_4 = ss(A, B, C, D)
@test_approx_eq tzero(ex_4) [-0.06467751189940692,-0.3680512036036696]

# Example 5
s = tf("s")
ex_5 = 1/s^15
@test tzero(ex_5) == Float64[]
@test tzero(ss(ex_5)) == Float64[]

# Example 6
A = [2 -1 0;
     0 0 0;
     -1 0 0]
B = [0 0 1]'
C = [0 -1 0]
D = [0]
ex_6 = ss(A, B, C, D)
@test tzero(ex_6) == Float64[]

# Example 7
ex_7 = ss(zeros(2, 2), [0 1]', [-1 0], [0])
@test tzero(ex_7) == Float64[]

# Example 8
A = [-2 1 0 0 0 0;
     1 -2 1 0 1 -1;
     0 1 -2 1 0 0;
     0 0 1 -1 0 1;
     0 -1 0 0 0 0;
     0 1 0 -1 0 0]
B = [1 0 0 0 1 0]'
C = [0 0 0 1 0 0]
D = [0]
ex_8 = ss(A, B, C, D)
# TODO : there may be a way to improve the precision of this example.
@test_approx_eq_eps tzero(ex_8) [-1.0, -1.0] 1e-7

# Example 9
ex_9 = (s - 20)/s^15
@test_approx_eq tzero(ex_9) [20.0]
@test_approx_eq tzero(ss(ex_9)) [20.0]

# Example 11
A = [-2 -6 3 -7 6;
     0 -5 4 -4 8;
     0 2 0 2 -2;
     0 6 -3 5 -6;
     0 -2 2 -2 5]
B = [-2 -8 -3 1 -8;
     7 -5 0 5 0]'
C = [0 -1 2 -1 -1;
     1 1 1 0 -1;
     0 3 -2 3 -1]
D = [0 0;
     0 0;
     0 0]
ex_11 = ss(A, B, C, D)
@test_approx_eq tzero(ex_11) [4.0, -3.0]

# Test for multiple zeros, siso tf
sys = s*(s + 1)*(s^2 + 1)*(s - 3)/((s + 1)*(s + 4)*(s - 4))
@test_approx_eq tzero(sys) [-1.0, -im, im, 3.0, 0.0]

## POLE ##
@test_approx_eq pole(sys) [-1.0, 4.0, -4.0]
@test_approx_eq pole([sys sys]) [-1.0, 4.0, -4.0, -1.0, 4.0, -4.0]
@test_approx_eq pole(ex_11) eig(ex_11.A)[1]

## ZPKDATA ##
# Sort a complex vector by real, breaking ties with imag
sortcomplex(a) = sort!(sort(a, by=imag), alg=MergeSort, by=real)
# Compare each vector in an array of vectors
macro test_array_vecs_eps(a, b, tol)
    quote
        @test size($a) == size($b)
        for (res, sol) = zip($a, $b)
            @test_approx_eq_eps sortcomplex(res) sol $tol
        end
    end
end
H = [tf(0) tf([3, 0],[1, 1, 10]) ; tf([1, 1],[1, 5]) tf([2],[1, 6])]
G = ss(H)
sol_z = vecarray(Complex128, 2, 2, Complex128[], Complex128[0.0 + 0.0im],
        Complex128[-1.0 + 0.0im], Complex128[])
sol_p = vecarray(Complex128, 2, 2, Complex128[], Complex128[-0.5 - 3.1224989991991996im,
        -0.5 + 3.1224989991991996im],
        Complex128[-5.0 + 0.0im], Complex128[-6.0 + 0.0im])
sol_k = [0.0 3.0; 1.0 2.0]
z, p, k = zpkdata(H)
@test_array_vecs_eps z sol_z 2*eps(Complex128)
@test_array_vecs_eps p sol_p 2*eps(Complex128)
@test k == sol_k
z, p, k = zpkdata(G)
@test_array_vecs_eps z sol_z 10*eps(Complex128)
@test_array_vecs_eps p sol_p 10*eps(Complex128)
@test k == sol_k

## GAIN ## #Gain is confusing when referring to zpkdata. Test dcgain instead
@test [dcgain(H[1, 1]) dcgain(H[1, 2]); dcgain(H[2, 1]) dcgain(H[2, 2])] ≈ [0 0; 0.2 1/3]
@test [dcgain(G[1, 1]) dcgain(G[1, 2]); dcgain(G[2, 1]) dcgain(G[2, 2])] ≈ [0 0; 0.2 1/3]
@test_throws ErrorException dcgain(H)
@test_throws ErrorException dcgain(G)

## MARKOVPARAM ##
@test markovparam(G, 0) == [0.0 0.0; 1.0 0.0]
@test markovparam(G, 1) == [0.0 3.0; -4.0 2.0]
@test markovparam(G, 2) == [0.0 -3.0; 20.0 -12.0]

## DAMP ##
@test_approx_eq damp(sys)[1] [1.0, 4.0, 4.0]
@test_approx_eq damp(sys)[2] [1.0, -1.0, 1.0]
@test_approx_eq damp(ex_11)[1] [1.0, 1.0, 2.0, 2.0, 3.0]
@test_approx_eq damp(ex_11)[2] [1.0, -1.0, -1.0, 1.0, -1.0]

## DAMPREPORT ##
@test sprint(dampreport, sys) == (
"|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
"|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
"+---------------+---------------+---------------+---------------+\n"*
"|  -1.000e+00   |  1.000e+00    |  1.000e+00    |  1.000e+00    |\n"*
"|  4.000e+00    |  -1.000e+00   |  4.000e+00    |  -2.500e-01   |\n"*
"|  -4.000e+00   |  1.000e+00    |  4.000e+00    |  2.500e-01    |\n")
@test sprint(dampreport, ex_11) == (
"|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
"|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
"+---------------+---------------+---------------+---------------+\n"*
"|  -1.000e+00   |  1.000e+00    |  1.000e+00    |  1.000e+00    |\n"*
"|  1.000e+00    |  -1.000e+00   |  1.000e+00    |  -1.000e+00   |\n"*
"|  2.000e+00    |  -1.000e+00   |  2.000e+00    |  -5.000e-01   |\n"*
"|  -2.000e+00   |  1.000e+00    |  2.000e+00    |  5.000e-01    |\n"*
"|  3.000e+00    |  -1.000e+00   |  3.000e+00    |  -3.333e-01   |\n")

end


text:
fileid = parse(Int64, ARGS[1])
simregexp = Regex("postsim.*_.*$(fileid).*.bson")

using BNP_WMReg_Joint
using RCall

using BSON
using BSON: @load

using Statistics
using PDMats
using StatsBase
using LinearAlgebra

files = filter(x -> occursin(simregexp, x), readdir("./postsim"))
nfiles = length(files)

i = 1
for i in 1:nfiles
    file_now = files[i]
    println("starting file: $(file_now)")

    BSON.@load "./postsim/$(file_now)" sims mesg model init
    nsim = length(sims)

    ### Send to R
    R"rm(list=ls())"
    iter = model.state.iter
    accptr =  model.state.accpt ./ model.state.iter
    K = model.K
    H = model.H
    n = model.n
    priorinfo = "$(model.prior)"
    initinfo = "$(init)"

    Sigx_type = "$(model.Σx_type)"
    gam_type = "$(model.γ_type)"

    sims_llik = [ sims[i][:llik] for i=1:nsim ]
    sims_noccup = [ sims[i][:n_occup] for i=1:nsim ]
    sims_Scounts = permutedims(hcat([ sims[i][:Scounts] for i=1:length(sims) ]...))

    if model.γ_type in (:global, :fixed)

        sims_gam = permutedims(hcat([ sims[i][:γ] for i=1:length(sims) ]...))
        sims_gam_occ = "NA"
        sims_gam_w = "NA"
        sims_gam_glob = "NA"

        if haskey(sims[1], :lfc_on)
            sims_fc_on = permutedims(hcat([ exp.(sims[i][:lfc_on]) for i=1:length(sims) ]...))
            sims_pigam = deepcopy(model.state.π_γ)
        else
            sims_fc_on = "NA"
            sims_pigam = deepcopy(model.state.π_γ)
        end

        b = 0.05
        sims_gam_alt = [ ( sims[ii][:γ][k] && abs(sims[ii][:β_y][h,k]) > b )  for ii = 1:length(sims), h = 1:model.H, k = 1:model.K ]
        sims_gam_alt_occup = [ sims[ii][:Scounts]'sims_gam_alt[ii,:,k] / model.n  for ii = 1:length(sims), k = 1:model.K ]

    elseif model.γ_type == :local

        sims_gam = "NA"
        sims_gam_occ = permutedims(hcat([ vec(sims[i][:γ_occupied]) for i=1:length(sims) ]...))
        sims_gam_w = [ sum( exp.(sims[i][:lω][sims[i][:γ][:,k]]) ) for i=1:length(sims), k = 1:model.K ]
        sims_gam_glob = [ any( sims[i][:γ][ 1:maximum( findall( sims[i][:Scounts] .> 0 ) ) , k] ) for i=1:length(sims), k = 1:model.K ]

        sims_fc_on = "NA"
        sims_pigam = permutedims(hcat([ sims[i][:π_γ] for i=1:length(sims) ]...))

        b = 0.05
        sims_gam_alt = [ ( sims[ii][:γ][h,k] && abs(sims[ii][:β_y][h,k]) > b )  for ii = 1:length(sims), h = 1:model.H, k = 1:model.K ]
        sims_gam_alt_occup = [ sims[ii][:Scounts]'sims_gam_alt[ii,:,k] / model.n  for ii = 1:length(sims), k = 1:model.K ]

    end

    sims_w = permutedims(hcat([ exp.(sims[i][:lω]) for i=1:length(sims) ]...))
    sims_alpha = [ sims[i][:α] for i=1:nsim ]

    sims_mu = permutedims(hcat([ sims[i][:μ_y] for i=1:length(sims) ]...))
    sims_beta = permutedims(cat([ sims[i][:β_y] for i=1:length(sims) ]..., dims=3), (3,2,1))
    sims_sig2 = permutedims(hcat([ sims[i][:δ_y] for i=1:length(sims) ]...))

    if model.γ_type in (:global, :fixed)
        sims_int = [ sims[i][:μ_y][h] + sum( sims[i][:γ] .* sims[i][:β_y][h,:] .* sims[i][:μ_x][h,:]) for i=1:length(sims), h=1:H ]
    elseif model.γ_type == :local
        sims_int = [ sims[i][:μ_y][h] + sum( sims[i][:γ][h,:] .* sims[i][:β_y][h,:] .* sims[i][:μ_x][h,:]) for i=1:length(sims), h=1:H ]
    end

    @rput mesg nsim iter accptr K H n priorinfo initinfo Sigx_type gam_type sims_llik sims_noccup sims_Scounts sims_gam sims_gam_occ sims_gam_w sims_gam_glob sims_fc_on sims_pigam sims_gam_alt_occup sims_w sims_alpha sims_mu sims_beta sims_sig2 sims_int

    R"ls()"
    R"save.image(file=paste0('postsim/mcmc_', $(mesg), '.rda'))"

    cp("./postsim/$(file_now)", "./postsimB/$(file_now)")
    rm("./postsim/$(file_now)")

    println("completed file number: $(i) of $(nfiles) \n\n")
end


text:
function submodel_jacobian_tests(submodel;
    dx = rand(Aeroelasticity.number_of_states(submodel)),
    x = rand(Aeroelasticity.number_of_states(submodel)),
    y = rand(Aeroelasticity.number_of_inputs(submodel)),
    p = rand(Aeroelasticity.number_of_parameters(submodel)),
    t = rand(),
    atol = sqrt(eps()),
    norm = (x) -> LinearAlgebra.norm(x, Inf))

    # rate jacobian test
    f = (dx) -> Aeroelasticity.get_residual(submodel, dx, x, y, p, t)
    J = Aeroelasticity.get_rate_jacobian(submodel, dx, x, y, p, t)
    J_fd = ForwardDiff.jacobian(f, dx)
    @test isapprox(J, J_fd; atol, norm)

    # state jacobian test
    f = (x) -> Aeroelasticity.get_residual(submodel, dx, x, y, p, t)
    J = Aeroelasticity.get_state_jacobian(submodel, dx, x, y, p, t)
    J_fd = ForwardDiff.jacobian(f, x)
    @test isapprox(J, J_fd; atol, norm)

    # input jacobian test
    f = (y) -> Aeroelasticity.get_residual(submodel, dx, x, y, p, t)
    J = Aeroelasticity.get_input_jacobian(submodel, dx, x, y, p, t)
    J_fd = ForwardDiff.jacobian(f, y)
    @test isapprox(J, J_fd; atol, norm)

    # parameter jacobian test
    f = (p) -> Aeroelasticity.get_residual(submodel, dx, x, y, p, t)
    J = Aeroelasticity.get_parameter_jacobian(submodel, dx, x, y, p, t)
    J_fd = ForwardDiff.jacobian(f, p)
    @test isapprox(J, J_fd; atol, norm)

    # time derivative test
    f = (t) -> Aeroelasticity.get_residual(submodel, dx, x, y, p, t)
    dT = Aeroelasticity.get_time_gradient(submodel, dx, x, y, p, t)
    dT_fd = ForwardDiff.derivative(f, t)
    @test isapprox(dT, dT_fd; atol, norm)

    return nothing
end

function coupling_jacobian_tests(coupling;
    dx = rand(Aeroelasticity.number_of_states(coupling)),
    x = rand(Aeroelasticity.number_of_states(coupling)),
    p = rand(Aeroelasticity.number_of_parameters(coupling)),
    t = rand(),
    atol = sqrt(eps()),
    norm = (x)->norm(x, Inf))

    # rate jacobian test
    f = (dx) -> Aeroelasticity.get_coupling_inputs(coupling, dx, x, p, t)
    J = Aeroelasticity.get_coupling_rate_jacobian(coupling, dx, x, p, t)
    J_fd = ForwardDiff.jacobian(f, dx)
    @test isapprox(J, J_fd; atol, norm)

    # state jacobian test
    f = (x) -> Aeroelasticity.get_coupling_inputs(coupling, dx, x, p, t)
    J = Aeroelasticity.get_coupling_state_jacobian(coupling, dx, x, p, t)
    J_fd = ForwardDiff.jacobian(f, x)
    @test isapprox(J, J_fd; atol, norm)

    # parameter jacobian test
    f = (p) -> Aeroelasticity.get_coupling_inputs(coupling, dx, x, p, t)
    J = Aeroelasticity.get_coupling_parameter_jacobian(coupling, dx, x, p, t)
    J_fd = ForwardDiff.jacobian(f, p)
    @test isapprox(J, J_fd; atol, norm)

    # time derivative test
    f = (t) -> Aeroelasticity.get_coupling_inputs(coupling, dx, x, p, t)
    dT = Aeroelasticity.get_coupling_time_gradient(coupling, dx, x, p, t)
    dT_fd = ForwardDiff.derivative(f, t)
    @test isapprox(dT, dT_fd; atol, norm)

    return nothing
end

function submodel_io_tests(model; atol = sqrt(eps()), norm = (x)->norm(x, Inf))

    # test separating and recombining state variables
    x = rand(Aeroelasticity.number_of_states(model))
    states = Aeroelasticity.separate_states(model, x)
    new_x = Aeroelasticity.set_states(model; states...)
    @test isapprox(new_x, x; atol, norm)

    # test separating and recombining inputs
    y = rand(Aeroelasticity.number_of_inputs(model))
    inputs = Aeroelasticity.separate_inputs(model, y)
    new_y = Aeroelasticity.set_inputs(model; inputs...)
    @test isapprox(new_y, y; atol, norm)

    # test separating and recombining parameters
    p = rand(Aeroelasticity.number_of_parameters(model))
    parameters = Aeroelasticity.separate_parameters(model, p)
    new_p = Aeroelasticity.set_parameters(model; parameters...)
    @test isapprox(new_p, p; atol, norm)

    return nothing
end

function coupling_io_tests(coupling; atol = sqrt(eps()), norm = (x)->norm(x, Inf))

    # test separating and recombining additional parameters
    pc = rand(Aeroelasticity.number_of_additional_parameters(coupling))
    parameters = Aeroelasticity.separate_parameters(coupling, pc)
    new_pc = Aeroelasticity.set_parameters(coupling; parameters...)
    @test isapprox(new_pc, pc; atol, norm)

    return nothing
end


text:
mutable struct HTensorNode{T} <: AbstractTensor{T}
   data      :: Array{T,2}
   subshape  :: Vector{Int}
   nbasis    :: Int
   children  :: Vector{HTensorNode{T}}
   partition :: Vector{Int}
end

reshape2D(a::Array{T,1}) where {T} = reshape(a, 1, length(a))
reshape2D(a::Array{T,N}) where {T,N} = reshape(a, prod(size(a,i) for i=1:N-1), size(a,N))
coreshape(t::HTensorNode) = reshape(t.data, t.subshape..., t.nbasis)

function HTensorNode(data::Array{T,2}, children = HTensorNode{T}[]) where {T}
   nbasis = size(data,2)
   if isempty(children)
      subshape = [size(data,1)]
   else
      subshape = [child.nbasis for child in children]
      @assert size(data,1) == prod(subshape)
   end
   partition = Int[ nmode(c) for c in children ]
   HTensorNode(data, subshape, nbasis, children, partition)
end

function HTensorNode(data::Array{T,N}, children = HTensorNode{T}[]) where {T,N}
   nbasis = size(data,N)
   subshape = [size(data,i) for i=1:N-1]
   data2D = reshape2D(data)
   if !isempty(children)
      @assert length(children) == N-1
      for i=1:N-1
         @assert subshape[i] == children[i].nbasis
      end
   end
   partition = Int[ nmode(c) for c in children ]
   HTensorNode(data2D, subshape, nbasis, children, partition)
end

isleaf(t::HTensorNode)  = isempty(t.children)
nmode(t::HTensorNode)   = isleaf(t) ? length(t.subshape) : sum(nmode(c) for c in t.children)
fulldim(t::HTensorNode) = isleaf(t) ? prod(Float64, t.subshape) : prod(fulldim(c) for c in t.children)

function ncoeffs(ht::HTensorNode)
   n = length(ht.data)
   for child in ht.children
      n += ncoeffs(child)
   end
   return n
end

function shape(t::HTensorNode)
   if isleaf(t)
      return t.subshape
   else
      shp = Array{Int}(undef, sum(t.partition))
      f = 0
      for i = 1:length(t.children)
         nm = t.partition[i]
         ch = t.children[i]
         shp[f+1:f+nm] = shape(ch)
         f += nm
      end
      return shp
   end
end

function array(t::HTensorNode)
   core = coreshape(t)
   if isleaf(t)
      return core
   else
      nchildren = length(t.children)
      newshape  = Int[]
      for i in 1:nchildren
         child = t.children[i]
         charr = array(child)
         core  = tensordot(core, reshape2D(charr), i)
         append!(newshape, size(charr)[1:end-1])
      end
      reshape(core, newshape..., t.nbasis)
   end
end

function getelts(t::HTensorNode, idxs)
   core = coreshape(t)
   if isleaf(t)
      return core[idxs..., :]
   else
      f = 0
      nchildren = length(t.children)
      newshape  = Int[]
      for i in 1:nchildren
         ch_nmode = t.partition[i]
         ch_idxs  = idxs[f+1 : f+ch_nmode]
         ch_elts  = getelts(t.children[i], ch_idxs)
         core = tensordot(core, reshape2D(ch_elts), i)
         append!(newshape, size(ch_elts)[1:end-1])
         f += ch_nmode
      end
      reshape(core, newshape..., t.nbasis)
   end
end


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

function run_test(f::PiecewiseQuadratic, h::PiecewiseQuadratic)
    return (simplify(envelope(f)) ≈ h) && (simplify(envelope(reverse(f))) ≈ reverse(h))
end

function plot_test(f::PiecewiseQuadratic, h::PiecewiseQuadratic)
    clf()
    X, Y = get_plot(f)
    plot(X, Y)
    X, Y = get_plot(simplify(envelope(f)))
    plot(X, Y)
    X, Y = get_plot(h)
    plot(X, Y; linestyle=:dash)
    return legend(["f", "env(f)", "f**"])
end

@testset "Midpoint to endpoint" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 1.0, 1.0, 2.0, 2.0),
                            BoundedQuadratic(1.0, 1.0, 1.0, 2.0, 1.0),
                            BoundedQuadratic(1.0, Inf, 1.0, 2.0, 2.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(-Inf, 0.0, 1.0, 2.0, 2.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, 2.0, 2.0),
                            BoundedQuadratic(1.0, 2.0, 0.0, 6.0, -2.0),
                            BoundedQuadratic(2.0, Inf, 1.0, 2.0, 2.0)])
    @test run_test(f, h)
    @test run_test(reverse(f), reverse(h))
end

@testset "Convex functions" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 1.0, 0.0, -1.0, 1.0),
                            BoundedQuadratic(1.0, Inf, 0.0, 1.0, -1.0)])
    @test run_test(f, f)
    @test run_test(reverse(f), reverse(f))

    f = PiecewiseQuadratic([BoundedQuadratic(1.0, -2.0, 1.0)])
    @test run_test(f, f)
    @test run_test(reverse(f), reverse(f))
end

@testset "Nonconvex piecewise affine" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 1.0, 0.0, 1.0, 2.0),
                            BoundedQuadratic(1.0, 1.0, 0.0, 2.0, 0.0),
                            BoundedQuadratic(1.0, Inf, 0.0, 3.0, 0.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(-Inf, 1.0, 0.0, 1.0, 1.0),
                            BoundedQuadratic(1.0, Inf, 0.0, 3.0, -1.0)])
    @test run_test(f, h)
    @test run_test(reverse(f), reverse(h))
end

@testset "Midpoint to midpoint" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 0.0, 1.0, 1.0, 1.0),
                            BoundedQuadratic(0.0, Inf, 1.0, -1.0, 1.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(-Inf, -0.5, 1.0, 1.0, 1.0),
                            BoundedQuadratic(-0.5, 0.5, 0.0, 0.0, 0.75),
                            BoundedQuadratic(0.5, Inf, 1.0, -1.0, 1.0)])
    @test run_test(f, h)
end

@testset "Upper to lower" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 0.0, 1.0, -1.0, 1.0),
                            BoundedQuadratic(0.0, Inf, 1.0, 1.0, 1.0)])
    @test run_test(f, f)

    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, 0.0, 0.0, -1.0, 1.0),
                            BoundedQuadratic(0.0, Inf, 0.0, 1.0, 1.0)])
    @test run_test(f, f)
end

@testset "Irrelevant point" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 1.0),
                            BoundedQuadratic(0.0, Inf, 0.0, 1.0, 0.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(0.0, Inf, 0.0, 1.0, 0.0)])
    @test run_test(f, h)
end

@testset "Relevant point" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, -1.0),
                            BoundedQuadratic(0.0, Inf, 0.0, 1.0, 0.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(0.0, Inf, 0.0, 1.0, -1.0)])
    @test run_test(f, h)

    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, -1.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, 1.0, 0.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(0.0, 1.0, 0.0, 2.0, -1.0)])
    @test run_test(f, h)
end

@testset "Lower to upper" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-1.0, 0.0, 0.0, 1.0, 1.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, -1.0, 1.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(-1.0, 1.0, 0.0, 0.0, 0.0)])
    @test run_test(f, h)
end

@testset "Upper to upper" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, -1.0, 1.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, 0.0, 0.0)])
    @test run_test(f, h)
end

@testset "Two points" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 1.0),
                            BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 0.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 0.0)])
    @test run_test(f, h)

    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 0.0),
                            BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 1.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 0.0)])
    @test run_test(f, h)
end

@testset "Test case 1" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 1.0, 0.0, 0.0, 0.0),
                            BoundedQuadratic(1.0, 2.0, 0.0, 1.0, -1.0),
                            BoundedQuadratic(2.0, Inf, 1.0, -4.0, 5.0)])
    z1 = 0.8284271247461898
    z2 = 2.414213562373095
    h = PiecewiseQuadratic([BoundedQuadratic(0.0, 1.0, 0.0, 0.0, 0.0),
                            BoundedQuadratic(1.0, z2, 0.0, z1, -z1),
                            BoundedQuadratic(z2, Inf, 1.0, -4.0, 5.0)])

    @test run_test(f, h)
end

@testset "Test case 2" begin
    f = PiecewiseQuadratic([BoundedQuadratic(-2.0, -1.0, 0.0, 1.0, -1.0),
                            BoundedQuadratic(-1.0, 0.0, 0.0, 2.0, 0.0),
                            BoundedQuadratic(0.0, Inf, 0.0, 0.0, 0.0)])
    h = PiecewiseQuadratic([BoundedQuadratic(-2.0, Inf, 0.0, 0.0, -3.0)])

    @test run_test(f, h)
end

@testset "Test case 3" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.00000000000000000000, 0.02157170544841656434,
                                             11.64694718222335900748,
                                             -0.04528496679275312298,
                                             0.00031129231510998337),
                            BoundedQuadratic(0.02157170544841656434, 0.02157170544841656434,
                                             11.64694718222335900748,
                                             -0.04528496679275312298,
                                             0.00031129231510998337),
                            BoundedQuadratic(0.02157170544841656434, Inf,
                                             11.64694718222335900748,
                                             -0.04528497091670485170,
                                             0.00031129240407065537)])
    h = deepcopy(f)
    @test run_test(simplify(f), simplify(h))
end


text:
"""
  This function returns an array with the number of nodes classified on each
  mesh vertex, edge, face, and region (in 3d) for the different operator
  types and degrees.

  op_type    Description
    1        SBP Gamma, CG, 2d
    2        SBPGamma, DG, 2d
    3        SBPOmega, DG 2d
    4        SBPGamma, DG, 3d
    5        SBPOmega, DG, 3d

  Degree 1 - 4 is supported for all operators.
  An exception will be thrown if an unsupported op_type and degree
  combinatino is specified.

  Inputs:
    op_type: integer specifying which operator type
    p: degree of the operator (linear, quadratic, etc.)

  Outputs:
    nodecnts: array of length 3 in 2D or 4 in 3D specifying number of
              nodes classified on each dimension mesh entity
"""
function getNodeCounts(op_type, order)
# get the number of nodes classifies on each vert, edge, and element
# op type specifies the type of operator, order the order of accuracy
  if op_type == 1  # SBPGamma, CG, 2d
    if order == 1
      return nodecnts = [1, 0, 0]
    elseif order == 2
      return nodecnts = [1, 1, 1]
    elseif order == 3
      return nodecnts = [1, 2, 3]
    elseif order == 4
      return nodecnts = [1, 3, 6]
    else
      error("unsuported op_type/order")
    end
  elseif op_type == 2  # SBPGamma DG, 2D
    if order == 1
      return [0, 0, 3]
    elseif order == 2
      return [0, 0, 7]
    elseif order == 3
      return [0, 0, 12]
    elseif order == 4
      return [0, 0, 18]
    else
      error("unsupported op_type/order")
    end
  elseif op_type == 3  # SBPOmega DG, 2D
    if order == 1
      return [0, 0, 3]
    elseif order == 2
      return [0, 0, 6]
    elseif order == 3
      return [0, 0, 10]
    elseif order == 4
      return [0, 0, 15]
    else
      error("unsupported op_type/order")
    end
  elseif op_type == 4  # SBPGamma DG, 3D
    if order == 1
      return [0, 0, 0, 4]
    elseif order == 2
      return [0, 0, 0, 11]
    elseif order == 3
      return [0, 0, 0, 24]
    elseif order == 4
      return [0, 0, 0, 45]
    else
      error("unsupported op_type/order")
    end
  elseif op_type == 5  # SBPOmega DG, 3D
    if order == 1
      return [0, 0, 0, 4]
    elseif order == 2
      return [0, 0, 0, 10]
    elseif order == 3
      return [0, 0, 0, 20]
    elseif order == 4
      return [0, 0, 0, 38]
    else
      error("unsupported op_type/order")
    end

  end

end

function calcDofs(m, nodecnt::AbstractArray{T}) where T <: Integer
# calculate number of dofs in a rectangular domain with the same number of
# elements along both axes
# m = number of elements per side
# p = order of SBP operators

  if length(nodecnt) == 3
    n = m  # number of elements in horiztonal direction
    numV = (m+1)*(n+1)
    numEdge = m*(n+1) + n*(m+1) + m*n
    numEl = 2*m*n

  #  counts = [numV, numEdge, numEl]
    return numV*nodecnt[1] + numEdge*nodecnt[2] + numEl*nodecnt[3]
  else
    nx = m
    ny = nx
    nz = nx
    nvert = (nx+1)*(ny+1)*(nz+1)

    nedges = nx*(ny+1)*(nz+1) +  # cube dges
             ny*(nx+1)*(nz+1) +
             nz*(nx+1)*(ny+1) +
             nx*nz*(ny+1) +  # cube face edges
             ny*nz*(nx+1) +
             nx*ny*(nz+1) +
             nx*ny*nz  # cube interior

    nfaces = nx*ny*(nz+1)*2 +
             ny*nz*(nx+1)*2 +
             nx*ny*(ny+1)*2 +
             nx*ny*nz*6

    nel = nx*ny*nz*6

    return nvert*nodecnt[1] + nedges*nodecnt[2] + nfaces*nodecnt[3] + nel*nodecnt[4]
  end

  return 0
#  numDof = dot(nodecnts, counts)

#  println("for a $m x $m mesh, there are $numDof dofs")
#  return numDof
end

function calcDofs3D()

  ny = nx
  nz = nx
  nvert = (nx+1)*(ny+1)*(nz+1)

  nedges = nx*(ny+1)*(nz+1) +  # cube dges
           ny*(nx+1)*(nz+1) +
           nz*(nx+1)*(ny+1) +
           nx*nz*(ny+1) +  # cube face edges
           ny*nz*(nx+1) +
           nx*ny*(nz+1) +
           nx*ny*nz  # cube interior

  nfaces = nx*ny*(nz+1)*2 +
           ny*nz*(nx+1)*2 +
           nx*ny*(ny+1)*2 +
           nx*ny*nz*6

  nel = nx*ny*nz*6

  return nvert*nodecnt[1] + nedges*nodecnt[2] + nfaces*nodecnt[3] + nel*nodecnt[4]

end

"""
  Find the mesh with specified operator operator type and degree elements
  with the closest number of degrees of freedom to the specified number

  Inputs:
    numdof: the number of dofs to try to match
    op_type: operator type (see getNumNodes)
    p: degree of operator

  Outputs:
    mesh_size_i: number of degree p element in each direction
    dof_mesh_size_i: total number of degrees of freedom on the mesh

"""
function findMeshSize(numdof::Integer, op_type, p::Integer)
# find the size of the square mesh with elements of order p
# with the closest number of dofs to
# the p1 mesh with p1_el elements per side,
# returns the number of elements per side, total number of dofs

#  p1dofs = calcDofs(p1_el, 1)
#  println("the p1 mesh has ", p1dofs, " dofs")

  # need to calculate an upper bound on the number of meshes to search
  # use the number of elements per side of a p1 mesh
  println("searching for p$p mesh with $numdof dofs")

  nodecnt = getNodeCounts(op_type, p)

  p1_el = round(Int, ceil(sqrt(numdof))) + 1 + 1

  mesh_size_i = 0  # declare in this scope
  dof_mesh_size_i = 0
    mesh_size = 1
    mesh_size_upper = 0
    for i=1:p1_el  # loop over mesh sizes
      dofs_i = calcDofs(i, nodecnt)
      if dofs_i > numdof  # break on first mesh larger than needed
        mesh_size_upper = i
        break
      end
    end
    # check whether the slightly larger of slightly smaller mesh is closest
    dofs_upper = calcDofs(mesh_size_upper, nodecnt)
    dofs_lower = calcDofs(mesh_size_upper - 1, nodecnt)

    if abs(dofs_upper - numdof) < abs(dofs_lower - numdof)
      mesh_size_i = mesh_size_upper
      dof_mesh_size_i = dofs_upper
    else
      mesh_size_i = mesh_size_upper - 1
      dof_mesh_size_i = dofs_lower
    end
    println("p $p mesh should have $mesh_size_i elements per side, and therefore $dof_mesh_size_i dofs")

  return mesh_size_i, dof_mesh_size_i
end

"""
  This function finds the number of element in each direction of a square
  (or cubic) mesh that has the closest number of degrees of freedom to
  a p=1 mesh with the specified number of elements

  Inputs:
    p1_el: number of element in each direction of the p=1 mesh
    op_type: operator type (see getNodeCounts)
    p: degree of operator

  Outputs:
    mesh_size_i: number of degree p element in each direction
    dof_mesh_size_i: total number of degrees of freedom on the mesh
"""
function findMeshSizeEl(p1_el, op_type, p)
# find the size of of mesh with the closest numbere of dofs to a p1 mesh
# of the specified size
# the mesh has the same number of element along the x and y axes
# returns the number of element per side and the total number of dofs

  nodecnt = getNodeCounts(op_type, 1)
  p1dofs = calcDofs(p1_el, nodecnt)
  println("the p1 mesh has ", p1dofs, " dofs")
  return findMeshSize(p1dofs, op_type, p)
end

"""
  This function calls findMeshSizeEl for a range of element sizes
  and writes the results to a file.
"""
function run_el()
  p1_mesh_sizes = collect(11:5:50)
  p = 4
  fname = "mesh_counts.dat"
  if isfile(fname)
    rm(fname)
  end
  f = open(fname, "w")
  for i=1:length(p1_mesh_sizes)
    els, dofs = findMeshSizeEl(p1_mesh_sizes[i], p)
    println(f, els)
  end

  close(f)

  return nothing
end

"""
  This function calls findMeshSize for a range of dof values and writes the
  resutls to a file.
"""
function run_dof()
  dof_spacing = collect(linspace(2751, 12865, 10))
  dof_counts = zeros(Int, length(dof_spacing))
  for i=1:length(dof_spacing)
    dof_counts[i] = round(dof_spacing[i])
  end
  p = 2
  fname = "mesh_counts.dat"
  if isfile(fname)
    rm(fname)
  end
  f = open(fname, "w")
  for i=1:length(dof_counts)
    els, dofs = findMeshSize(dof_counts[i], p)
    println(f, els)
  end

  close(f)

  return nothing
end

#run_dof()
op_type = 4
p = 3
numel = 40
numdof = 64^3
#findMeshSizeEl(numel, op_type, p)
findMeshSize(numdof, op_type, p)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ cd305c42-8c85-11eb-0745-11ba4f9549f7
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 764348e6-8c84-11eb-3487-23a16d28cc9b
md"""# Struct and Objects

At this point you know how to use functions to organize code and built-in types to organize data. The next step is to learn how to build your own types to organize both code and data. This is a big topic; it will take a few chapters to get there.
"""

# ╔═╡ dcbcdf34-8c85-11eb-178c-41b0cc73c257
md"""## Composite Types

We have used many of Julia’s built-in types; now we are going to define a new type. As an example, we will create a type called Point that represents a point in twodimensional space.

In mathematical notation, points are often written in parentheses with a comma sepa‐ rating the coordinates. For example, ``(0, 0)`` represents the origin, and ``(x, y)`` represents the point ``x`` units to the right and ``y`` units up from the origin.

There are several ways we might represent points in Julia:

• We could store the coordinates separately in two variables, `x` and `y`.
• We could store the coordinates as elements in an array or tuple.
• We could create a new type to represent points as objects.

Creating a new type is more complicated than the other options, but it has advantages that will be apparent soon.

A programmer-defined *composite type* is also called a *struct*. The `struct` definition for a point looks like this:

```julia
struct Point
	x
	y
end
```

The header indicates that the new struct is called `Point`. The body defines the attributes or fields of the struct. The `Point` struct has two fields: `x` and `y`.

A struct is like a factory for creating objects. To create a point, you call `Point` as if it were a function having as arguments the values of the fields. When `Point` is used as a function, it is called a *constructor*:

```julia
julia> p = Point(3.0, 4.0)
Point(3.0, 4.0)
```

The return value is a reference to a `Point` object, which we assign to `p`.

Creating a new object is called *instantiation*, and the object is an instance of the type.

When you print an instance, Julia tells you what type it belongs to and what the values of the attributes are.

Every object is an instance of some type, so “object” and “instance” are interchangea‐ ble. But in this chapter I use “instance” to indicate that I am talking about a programmer-defined type.

A state diagram that shows an object and its fields is called an *object diagram*; see Figure 15-1.
"""

# ╔═╡ 631d8fb6-8c86-11eb-3f26-333749163e46
Drawing(width=720, height=90) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	text(x=275, y=60, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("p ->")
	end
	rect(x=325, y=30, width=130, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=335, y=20, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("Point")
	end
	text(x=335, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("x")
	end
	line(x1=355, y1=45, x2=395, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=415, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3.0")
	end
	text(x=335, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("y")
	end
	line(x1=355, y1=65, x2=395, y2=65, stroke="black", marker_end="url(#arrow)")
	text(x=415, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("4.0")
	end
end

# ╔═╡ 42ef9f9a-8c88-11eb-1cfc-a92de41ca588
md"*Figure 15-1. Object diagram.*"

# ╔═╡ 4cd9060e-8c88-11eb-31a9-ef66b3a0364a
md"""## Structs Are Immutable

You can get the values of the fields using `.` notation:

```julia
julia> x = p.x
3.0
julia> p.y
4.0
```

The expression `p.x` means, “Go to the object `p` refers to and get the value of `x`.” In the example, we assign that value to a variable named `x`. There is no conflict between the variable `x` and the field `x`.

You can use dot notation as part of any expression. For example:

```julia
julia> distance = sqrt(p.x^2 + p.y^2)
5.0
```

Structs are, however, by default immutable; after construction the fields cannot change value:

```julia
julia> p.y = 1.0
ERROR: setfield! immutable struct of type Point cannot be changed
```

This may seem odd at first, but it has several advantages:

* It can be more efficient.
* It is not possible to violate the invariants provided by the type’s constructors (see “Constructors”).
* Code using immutable objects can be easier to reason about.
"""

# ╔═╡ 9afc564c-8c88-11eb-369a-e768ee84d8da
md"""## Mutable Structs

Where required, mutable composite types can be declared with the keyword `mutable struct`. Here is the definition of a mutable point:

```julia
mutable struct MPoint
	x
	y
end
```

You can assign values to an instance of a mutable struct using dot notation:

```julia
julia> blank = MPoint(0.0, 0.0)
MPoint(0.0, 0.0)
julia> blank.x = 3.0
3.0
julia> blank.y = 4.0
4.0
```
"""

# ╔═╡ bd8086c0-8c88-11eb-3960-c355a4ac636c
md"""## Rectangles

Sometimes it is obvious what the fields of an object should be, but other times you have to make decisions. For example, imagine you are designing a type to represent rectangles. What fields would you use to specify the location and size of a rectangle? You can ignore angle; to keep things simple, assume that the rectangle is either vertical or horizontal.

There are at least two possibilities:

* You could specify one corner of the rectangle (or the center), the width, and the height.
* You could specify two opposing corners.

At this point it is hard to say whether one is better than the other, so we’ll implement the first one, just as an example:

```julia
\"\"\"
Represents a rectangle.
fields: width, height, corner
\"\"\"
struct Rectangle
	width
	height
	corner
end
```

The docstring lists the fields: `width` and `height` are numbers; corner is a `Point` object that specifies the lower-left corner.

To represent a rectangle, you have to instantiate a `Rectangle` object:

```julia
julia> origin = MPoint(0.0, 0.0)
MPoint(0.0, 0.0)
julia> box = Rectangle(100.0, 200.0, origin)
Rectangle(100.0, 200.0, MPoint(0.0, 0.0))
```

Figure 15-2 shows the state of this object. An object that is a field of another object is embedded. Because the corner attribute refers to a mutable object, the latter is drawn outside the `Rectangle` object.
"""

# ╔═╡ 1d636f6e-8c89-11eb-2762-9b68204ed3eb
Drawing(width=720, height=120) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	text(x=140, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("box ->")
	end
	rect(x=200, y=30, width=200, height=70, fill="rgb(242, 242, 242)", stroke="black")
	text(x=210, y=20, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("Rectangle")
	end
	text(x=270, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("width")
	end
	line(x1=280, y1=45, x2=320, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=340, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("100.0")
	end
	text(x=270, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("height")
	end
	line(x1=280, y1=65, x2=320, y2=65, stroke="black", marker_end="url(#arrow)")
	text(x=340, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("200.0")
	end
	text(x=270, y=90, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600, text_anchor="end") do
		str("corner")
	end
	line(x1=280, y1=85, x2=440, y2=85, stroke="black", marker_end="url(#arrow)")
	rect(x=450, y=60, width=130, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=460, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("MPoint")
	end
	text(x=460, y=80, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("x")
	end
	line(x1=480, y1=75, x2=520, y2=75, stroke="black", marker_end="url(#arrow)")
	text(x=540, y=80, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("0.0")
	end
	text(x=460, y=100, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("y")
	end
	line(x1=480, y1=95, x2=520, y2=95, stroke="black", marker_end="url(#arrow)")
	text(x=540, y=100, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("0.0")
	end
end

# ╔═╡ 8b59f764-8c8b-11eb-0a18-fd33f5539326
md"*Figure 15-2. Object diagram.*"

# ╔═╡ 938ee1a6-8c8b-11eb-2c51-690096526be0
md"""## Instances as Arguments

You can pass an instance as an argument in the usual way. For example:

```julia
function printpoint(p)
	println("($(p.x), $(p.y))")
end
```

`printpoint` takes a `Point` as an argument and displays it in mathematical notation. To invoke it, you can pass `p` as an argument:

```julia
julia> printpoint(blank)
(3.0, 4.0)
```
"""

# ╔═╡ c0d707c4-8c8b-11eb-2784-7b62508b0a19
md"""#### Exercise 15-1
Write a function called `distancebetweenpoints` that takes two points as arguments and returns the distance between them.
"""

# ╔═╡ d55b88b6-8c8b-11eb-27f3-aba322cb7c63
md"""If a mutable struct object is passed to a function as an argument, the function can modify the fields of the object. For example, `movepoint!` takes a mutable `Point` object and two numbers, `dx` and `dy`, and adds the numbers to, respectively, the `x` and the ``y attribute of the `Point`:

```julia
function movepoint!(p, dx, dy)
	p.x += dx
	p.y += dy
	nothing
end
```

Here is an example that demonstrates the effect:

```julia
julia> origin = MPoint(0.0, 0.0)
MPoint(0.0, 0.0)
julia> movepoint!(origin, 1.0, 2.0)

julia> origin
MPoint(1.0, 2.0)
```

Inside the function, `p` is an alias for `origin`, so when the function modifies `p`, `origin` changes.

Passing an immutable Point object to movepoint! causes an error:

```julia
julia> movepoint!(p, 1.0, 2.0)
ERROR: setfield! immutable struct of type Point cannot be changed
```

You can, however, modify the value of a mutable attribute of an immutable object. For example, moverectangle! has as arguments a `Rectangle` object and two numbers, `dx` and `dy`, and uses `movepoint!` to move the corner of the rectangle:

```julia
function moverectangle!(rect, dx, dy)
	movepoint!(rect.corner, dx, dy)
end
```

Now `p` in `movepoint!` is an alias for `rect.corner`, so when `p` is modified, `rect.corner` changes also:

```julia
julia> box
Rectangle(100.0, 200.0, MPoint(0.0, 0.0))
julia> moverectangle!(box, 1.0, 2.0)

julia> box
Rectangle(100.0, 200.0, MPoint(1.0, 2.0))
```

!!! danger
    You cannot reassign a mutable attribute of an immutable object:

    ```julia
    julia> box.corner = MPoint(1.0, 2.0)
    ERROR: setfield! immutable struct of type Rectangle cannot be changed
    ```
"""

# ╔═╡ 711693b6-8c8c-11eb-2d09-7b64f34fac84
md"""## Instances as Return Values

Functions can return instances. For example, `findcenter` takes a `Rectangle` as an argument and returns a `Point` that contains the coordinates of the center of the rectangle:

```julia
function findcenter(rect)
	Point(rect.corner.x + rect.width / 2, rect.corner.y + rect.height / 2)
end
```

The expression `rect.corner.x` means, “Go to the object `rect` refers to and select the field named `corner`; then go to that object and select the field named `x`.”

Here is an example that passes `box` as an argument and assigns the resulting `Point` to center:

```julia
julia> center = findcenter(box)
Point(51.0, 102.0)
```
"""

# ╔═╡ b3197742-8c8c-11eb-146a-2594cc0932d4
md"""## Copying

Aliasing can make a program difficult to read because changes in one place might have unexpected effects in another place. It is hard to keep track of all the variables that might refer to a given object.

Copying an object is often an alternative to aliasing. Julia provides a function called `deepcopy` that performs a *deep copy* and can duplicate any object, including the contents of any embedded objects:

```julia
julia> p1 = MPoint(3.0, 4.0)
MPoint(3.0, 4.0)
julia> p2 = deepcopy(p1)
MPoint(3.0, 4.0)
julia> p1 ≡ p2
false
julia> p1 == p2
false
```

The `≡` operator indicates that `p1` and `p2` are not the same object, which is what we expected. But you might have expected `==` to yield `true` because these points contain the same data. In that case, you will be disappointed to learn that for mutable objects, the default behavior of the `==` operator is the same as the `=== operator; it checks object identity, not object equivalence (see “Objects and Values”). That’s because for mutable composite types, Julia doesn’t know what should be considered equivalent—at least, not yet.
"""

# ╔═╡ 2762ab0a-8c8d-11eb-280d-bb3cc9a654ed
md"""#### Exercise 15-2

Create a `Point` instance, make a copy of it, and check the equivalence and the egality of the two objects. The result may surprise you, but it explains why aliasing is a non-issue for an immutable object.
"""

# ╔═╡ 567d07fa-8c8d-11eb-3912-a1ae98ba7085
md"""## Debugging

When you start working with objects, you are likely to encounter some new exceptions. If you try to access a field that doesn’t exist, you get:

```julia
julia> p = Point(3.0, 4.0)
Point(3.0, 4.0)
julia> p.z = 1.0
ERROR: type Point has no field z Stacktrace:
     [1] setproperty!(::Point, ::Symbol, ::Float64) at ./sysimg.jl:19
     [2] top-level scope at none:0
```

If you are not sure what type an object is, you can ask:

```julia
julia> typeof(p)
Point
```

You can also use `isa` to check whether an object is an instance of a type:

```julia
julia> p isa Point
true
```

If you are not sure whether an object has a particular attribute, you can use the built-in function `fieldnames`:

```julia
julia> fieldnames(Point)
(:x, :y)
```

or the function `isdefined`:

```julia
julia> isdefined(p, :x)
true
julia> isdefined(p, :z)
false
```

The first argument can be any object; the second argument is a symbol, `:`, followed by the name of the field.
"""

# ╔═╡ ae9d5610-8c8d-11eb-1d78-2d02e842a857
md"""## Glossary

*struct*:
A user-defined type consisting of a collection of named fields. Also called a composite type.

*attribute*:
One of the named values associated with an object. Also called a field.

*constructor*:
A function with the same name as a type that creates instances of the type.

*instantiate*:
To create a new object.

*instance*:
An object that belongs to a certain type.

*object diagram*:
A diagram that shows objects, their fields, and the values of the fields.

*embedded object*:
An object that is stored as a field of another object.

*deep copy*:
To copy the contents of an object as well as any embedded objects, and any objects embedded in them, and so on; implemented by the deepcopy function.
"""

# ╔═╡ 28c26714-8c8e-11eb-2966-a9f47946b1d0
md"""## Exercises

#### Exercise 15-3

1. Write a definition for a type named `Circle` with fields `center` and `radius`, where center is a `Point` object and `radius` is a number.
2. Instantiate a `Circle` object that represents a circle with its center at ``(150, 100)`` and radius ``75``.
3. Write a function named `pointincircle` that takes a `Circle` object and a `Point` object and returns `true` if the point lies in or on the boundary of the circle.
4. Write a function named `rectincircle` that takes a `Circle` object and a `Rectangle` object and returns `true` if the rectangle lies entirely in or on the boundary of the circle.
5. Write a function named `rectcircleoverlap` that takes a `Circle` object and a `Rectangle` object and returns `true` if any of the corners of the rectangle fall inside the circle. Or, as a more challenging version, return `true` if any part of the rectangle falls inside the circle.
"""

# ╔═╡ 8d10564a-8c8e-11eb-159f-63a1fbcb3aef
md"""#### Exercise 15-4

1. Write a function called `drawrect` that takes a `Turtle` object and a `Rectangle` object and uses the turtle to draw the rectangle. See Chapter 4 for examples using `Turtle` objects.
2. Write a function called `drawcircle` that takes a `Turtle` object and a `Circle` object and draws the circle.
"""

# ╔═╡ Cell order:
# ╟─cd305c42-8c85-11eb-0745-11ba4f9549f7
# ╟─764348e6-8c84-11eb-3487-23a16d28cc9b
# ╟─dcbcdf34-8c85-11eb-178c-41b0cc73c257
# ╟─631d8fb6-8c86-11eb-3f26-333749163e46
# ╟─42ef9f9a-8c88-11eb-1cfc-a92de41ca588
# ╟─4cd9060e-8c88-11eb-31a9-ef66b3a0364a
# ╟─9afc564c-8c88-11eb-369a-e768ee84d8da
# ╟─bd8086c0-8c88-11eb-3960-c355a4ac636c
# ╟─1d636f6e-8c89-11eb-2762-9b68204ed3eb
# ╟─8b59f764-8c8b-11eb-0a18-fd33f5539326
# ╟─938ee1a6-8c8b-11eb-2c51-690096526be0
# ╟─c0d707c4-8c8b-11eb-2784-7b62508b0a19
# ╟─d55b88b6-8c8b-11eb-27f3-aba322cb7c63
# ╟─711693b6-8c8c-11eb-2d09-7b64f34fac84
# ╟─b3197742-8c8c-11eb-146a-2594cc0932d4
# ╟─2762ab0a-8c8d-11eb-280d-bb3cc9a654ed
# ╟─567d07fa-8c8d-11eb-3912-a1ae98ba7085
# ╟─ae9d5610-8c8d-11eb-1d78-2d02e842a857
# ╟─28c26714-8c8e-11eb-2966-a9f47946b1d0
# ╟─8d10564a-8c8e-11eb-159f-63a1fbcb3aef


text:
export channel_match

channel_match( S::GphysData, j::Integer, fs::Float64) = (j == 0 ? 0 : S.fs[j] == fs ? j : 0)
function channel_match( S::GphysData,
                        j::Integer,
                        fs::Float64,
                        gain::AbstractFloat,
                        loc::InstrumentPosition,
                        resp::InstrumentResponse,
                        units::String
                      )
  (j == 0) && (return 0)
  return (
    try
      @assert(S.fs[j] == fs)
      @assert(S.gain[j] == gain)
      @assert(S.units[j] == units)
      @assert(loc == S.loc[j])
      @assert(resp == S.resp[j])
      j
    catch err
      0
    end
  )
end
function channel_match( S::GphysData, C::GphysChannel )
  j = findid(S, C.id)
  j = channel_match( S, j, C.fs, C.gain, C.loc, C.resp, C.units)
  return j
end


text:
using Test
using Forecast

@testset "pacf" begin
    x = rand(100);
    @test_throws AssertionError Forecast.pacf(x; lag = 100)

    res = pacf(x; lag = 20);
    @test res isa CCF
    @test res.call == "pacf(x; type=\"stepwise-real\", lag=20, alpha=(0.95, 0.99))"
    @test res.type == "pacf_stepwise-real"
    @test res.auto
    @test res.lag == 20
    @test res.N == 100
    @test res.alpha == (0.95, 0.99)
    @test res.ci[1] ≈ 0.22335862551831387
    @test res.ci[2] ≈ 0.293542992294061
end


text:
using CUDA

# Defining Matrix dimensions
row = 9
col = 8

# Declaring Matrix
pic = Array{Float32}(rand(row,col))

function scalePic(pic::Matrix{Float32}, row::Int, col::Int)
    # Allocate and transfer matrix to GPU Memory
    d_pic = CuArray(pic)

    # Kernel Invocation
    dimBlock = (2, 2) # Threads in a block
    dimGrid = (Int(ceil(col/dimBlock[1])), Int(ceil(row/dimBlock[2]))) # Blocks in a grid

    @cuda threads=dimBlock blocks=dimGrid scalePicKernel!(d_pic, row, col)

    # Transfer results back to CPU
    pic = Array(d_pic)

    return pic
end

function scalePicKernel!(pic, row, col)
    rowIdx = blockDim().y * (blockIdx().y-1) + threadIdx().y
    colIdx = blockDim().x * (blockIdx().x-1) + threadIdx().x

    if (rowIdx ≤ row) && (colIdx ≤ col)
        pic[rowIdx,colIdx] *= 2
    end

    return nothing
end

# Calling functions
pic_scaled = scalePic(pic, row, col)

@assert pic*2 == pic_scaled


text:
using LanguageServer
using CSTParser
using Test
Range = LanguageServer.Range

@testset "LanguageServer" begin

include("test_document.jl")
include("test_communication.jl")
include("test_hover.jl")

end


text:
using Test
using LinearAlgebra
using Rotations
using StaticArrays
using InteractiveUtils: subtypes
using Quaternions
import Unitful

import Random

# Check that there are no ambiguities beyond those present in StaticArrays
ramb = detect_ambiguities(Rotations, Base, Core)
samb = detect_ambiguities(StaticArrays, Base, Core)
@test isempty(setdiff(ramb, samb))

include("util_tests.jl")
include("2d.jl")
include("rotation_tests.jl")
include("derivative_tests.jl")
include("principal_value_tests.jl")
include("unitquat.jl")
include("rodrigues_params.jl")
include("quatmaps.jl")
include("rotation_error.jl")
include("distribution_tests.jl")
include("eigen.jl")
include("nearest_rotation.jl")
include("rotation_generator.jl")
include("logexp.jl")
include("deprecated.jl")

include(joinpath(@__DIR__, "..", "perf", "runbenchmarks.jl"))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SPIRV_LLVM_Translator_unified_jll
using Base
using Base: UUID
include(path) = Base.include(@__MODULE__, path)
include(joinpath("..", ".pkg", "platform_augmentation.jl"))
import JLLWrappers

JLLWrappers.@generate_main_file_header("SPIRV_LLVM_Translator_unified")
JLLWrappers.@generate_main_file("SPIRV_LLVM_Translator_unified", UUID("85f0d8ed-5b39-5caa-b1ae-7472de402361"))
end  # module SPIRV_LLVM_Translator_unified_jll


text:
"""
    convert2image(array) -> Array{Gray}

Convert the given MNIST horizontal-major tensor (or feature matrix)
to a vertical-major `Colorant` array. The values are also color
corrected according to the website's description, which means that
the digits are black on a white background.

```julia
julia> MNIST.convert2image(MNIST.traintensor()) # full training dataset
28×28×60000 Array{Gray{N0f8},3}:
[...]

julia> MNIST.convert2image(MNIST.traintensor(1)) # first training image
28×28 Array{Gray{N0f8},2}:
[...]
```
"""
function convert2image(array::AbstractArray{T}) where {T<:Number}
    nlast = size(array)[end]
    array = reshape(array, 28, 28, :)
    array = permutedims(array, (2, 1, 3))
    if size(array)[end] == 1 && nlast != 1
        array = dropdims(array, dims=3)
    end
    if any(x -> x > 1, array) # simple check if x in [0,1]
        img = _colorview(Gray, array ./ T(255))
    else
        img = _colorview(Gray, array)
    end
    img
end


text:
# This file was generated, do not modify it. # hide
rms(y[test], ŷ(rows=test))

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


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

module WignerChains

export wigner3j_chain

using WignerSymbols

################## W^3j_k ####################################################

# wigner3j_chain(ell, L, M): This is the frontend to calculate the W^3j_k. It
# would be phantastic to dispatch on the length of the arrays `ell`, `L`, and
# `M`.
function wigner3j_chain(ell, L, M)
    @assert length(ell) == length(L) == length(M)
    k = length(ell)
    if k == 1
        return calc_w3j_1(ell, L, M)
    elseif k == 2
        return calc_w3j_2(ell, L, M)
    elseif k == 3
        return calc_w3j_3(ell, L, M)
    elseif k == 4
        return calc_w3j_4(ell, L, M)
    elseif k == 5
        return calc_w3j_5(ell, L, M)
    else
        return calc_w3j_k(ell, L, M)
    end
end

# triangle(L1, L2, L3): Returns true if the three angular momenta satisfy the
# triangle condition.
function triangle(L1, L2, L3)
    return abs(L1 - L2) <= L3 <= L1 + L2
end

# all_triangles(ell, L): Return true if all triangle conditions are met.
function all_triangles(ell, L)
    for i=1:length(ell)
        im1 = (i == 1) ? length(ell) : i-1
        triangle(ell[im1], ell[i], L[i]) || return false
    end
    return true
end

# This is the general implementation. An improvement would be to perform cyclic
# permutations until the smallest ℓᵢ is in `ell[1]`.
function calc_w3j_k(ell, L, M)
    T = Float64
    k = length(ell)

    all_triangles(ell, L) || return T(0)

    # calculate
    w3jk = T(0)
    mm = fill(0, k)
    skip_m1 = false
    for m1=-ell[1]:ell[1]
        # set m's
        mm[1] = m1
        for i=2:k
            mm[i] = mm[i-1] - M[i]
            if abs(mm[i]) > ell[i]
                skip_m1 = true
                break
            end
        end
        if skip_m1
            skip_m1 = false
            continue
        end
        if -mm[end] + mm[1] + M[1] != 0
            continue
        end

        #@show m1,ell,mm,L,M

        # calculate term
        w = wigner3j(T, ell[end], ell[1], L[1], -mm[end], mm[1], M[1])
        #@show w
        for i=2:k
            w *= wigner3j(T, ell[i-1], ell[i], L[i], -mm[i-1], mm[i], M[i])
            #@show w
        end
        w *= (-1)^sum(mm)
        #@show w

        # sum it up
        w3jk += w
    end
    return w3jk
end

# k = 1
function calc_w3j_1(ell, L, M)
    T = Float64
    if L[1] != 0 || M[1] != 0
        return T(0)
    end
    return (-1)^ell[1] * √T(2*ell[1] + 1)
end

# k = 2
function calc_w3j_2(ell, L, M)
    T = Float64
    if L[1] != L[2] || M[1] != - M[2]
        return T(0)
    end
    if !triangle(L[1], ell[1], ell[2])
        return T(0)
    end
    return (-1)^(M[1]) / T(2*L[1] + 1)
end

# k = 2, simple
function calc_w3j_2_simple(ell, L, M)
    T = Float64
    l1 = ell[1]
    l2 = ell[2]
    L1 = L[1]
    L2 = L[2]
    M1 = M[1]
    M2 = M[2]
    s4 = T(0)
    for m1=-l1:l1, m2=-l2:l2
        p = (-1)^(m1 + m2)
        w1 = wigner3j(T, l2, l1, L1, -m2, m1, M1)
        w2 = wigner3j(T, l1, l2, L2, -m1, m2, M2)
        s4 += p * w1 * w2
    end
    return s4
end

# k = 3
function calc_w3j_3(ell, L, M)
    T = Float64
    triangle(L[1], L[2], L[3]) || return T(0)
    triangle(L[3], ell[2], ell[3]) || return T(0)
    triangle(L[1], ell[3], ell[1]) || return T(0)
    triangle(L[2], ell[2], ell[1]) || return T(0)
    w3jm = wigner3j(T, L[1], L[2], L[3], -M[1], -M[2], -M[3])
    w6j = wigner6j(T, L[1], L[2], L[3], ell[2], ell[3], ell[1])
    p = (-1)^(sum(ell) + sum(L))
    return p * w3jm * w6j
end

# k = 4
function calc_w3j_4(ell, L, M)
    T = Float64
    triangle(L[1], ell[1], ell[4]) || return T(0)
    triangle(L[2], ell[2], ell[1]) || return T(0)
    triangle(L[3], ell[3], ell[2]) || return T(0)
    triangle(L[4], ell[4], ell[3]) || return T(0)
    M5 = M[1] + M[4]
    (-M5 == M[2] + M[3]) || return T(0)
    L5min = max(abs(L[1] - L[4]), abs(L[2] - L[3]), abs(ell[1] - ell[3]), abs(M5))
    L5max = min(L[1] + L[4], L[2] + L[3], ell[1] + ell[3])
    w = T(0)
    for L5=L5min:L5max
        #@show L[1], L5, L[4], M[1], -M5, M[4]
        w3a = wigner3j(T, L[1], L5, L[4], M[1], -M5, M[4])
        w3b = wigner3j(T, L[2], L5, L[3], M[2], M5, M[3])
        w6a = wigner6j(T, L[1], L5, L[4], ell[3], ell[4], ell[1])
        w6b = wigner6j(T, L[2], L5, L[3], ell[3], ell[2], ell[1])
        w += (-1)^(L5-M5) * (2*L5 + 1) * w3a * w3b * w6a * w6b
    end
    return (-1)^(L[2] + L[3] + ell[2] + ell[4]) * w
end

# k = 5
function calc_w3j_5(ell, L, M)
    T = Float64
    triangle(L[1], ell[1], ell[5]) || return T(0)
    triangle(L[2], ell[2], ell[1]) || return T(0)
    triangle(L[3], ell[3], ell[2]) || return T(0)
    triangle(L[4], ell[4], ell[3]) || return T(0)
    triangle(L[5], ell[5], ell[4]) || return T(0)
    M6 = - (M[1] + M[2])
    M7 = - (M[3] + M[4])
    (M6 == M[5] - M7) || return T(0)
    L6min = max(abs(L[1] - L[2]), abs(ell[2] - ell[5]), abs(M6))
    L6max = min(L[1] + L[2], ell[2] + ell[5])
    w = T(0)
    for L6=L6min:L6max
        L7min = max(abs(L[3] - L[4]), abs(L[5] - L6), abs(ell[4] - ell[2]), abs(M7))
        L7max = min(L[3] + L[4], L[5] + L6, ell[4] + ell[2])
        for L7=L7min:L7max
            w3a = wigner3j(T, L[1], L[2], L6, M[1], M[2], M6)
            w3b = wigner3j(T, L6, L[5], L7, -M6, M[5], -M7)
            w3c = wigner3j(T, L7, L[3], L[4], M7, M[3], M[4])
            w6a = wigner6j(T, L[1], L[2], L6, ell[2], ell[5], ell[1])
            w6b = wigner6j(T, L6, L[5], L7, ell[4], ell[2], ell[5])
            w6c = wigner6j(T, L7, L[3], L[4], ell[3], ell[4], ell[2])
            p = (-1)^(L6-M6+L7-M7) * (2*L6+1) * (2*L7+1)
            w += p * w3a * w3b * w3c * w6a * w6b * w6c
        end
    end
    return (-1)^(ell[1] + ell[2] + ell[3] + L[5]) * w
end

end

# vim: set sw=4 et sts=4 :


text:
using Test
using StringEncodings

# Test round-trip to Unicode formats
for s in ("", "\0", "a", "café crème",
          "a"^(StringEncodings.BUFSIZE-1) * "€ with an incomplete codepoint between two input buffer fills",
          "a string € チャネルパートナーの選択",
          "a string \0€ チャネルパ\0ー\0トナーの選択 with embedded and trailing nuls\0"),
    enc in (enc"UTF-8", enc"UTF-16", enc"UTF-16LE", enc"UTF-16BE", enc"UTF-32")
    a = encode(s, enc)
    @test decode(a, enc) == s
    @test decode(encode(s, enc), enc) == s
end

# Test a few non-Unicode encodings
for (s, enc) in (("noël", "ISO-8859-1"),
                 ("noël €", "ISO-8859-15", "CP1252"),
                 ("Код Обмена Информацией, 8 бит", "KOI8-R"),
                 ("国家标准", "GB18030"))
    @test decode(encode(s, enc), enc) == s
end

# Test that attempt to close stream in the middle of incomplete sequence throws
let s = "a string チャネルパートナーの選択", a = unsafe_wrap(Vector{UInt8}, s)
    # First, correct version
    p = StringEncoder(IOBuffer(), "UTF-16LE")
    write(p, a)
    close(p)
    # Test that closed pipe behaves correctly
    @test_throws ArgumentError write(p, 'a')

    b = IOBuffer()
    p = StringEncoder(b, "UTF-16LE")
    @test string(p) == "StringEncoder{UTF-8, UTF-16LE}($(string(b)))"
    write(p, a[1:10])
    @test_throws IncompleteSequenceError close(p)
    # Test that closed pipe behaves correctly even after an error
    @test_throws ArgumentError write(p, 'a')

    # This time, call show
    p = StringEncoder(IOBuffer(), "UTF-16LE")
    write(p, a[1:10])
    try
        close(p)
    catch err
        @test isa(err, IncompleteSequenceError)
        io = IOBuffer()
        showerror(io, err)
        @test String(take!(io)) == "Incomplete byte sequence at end of input"
    end

    b = IOBuffer(encode(s, "UTF-16LE")[1:19])
    p = StringDecoder(b, "UTF-16LE")
    @test string(p) == "StringDecoder{UTF-16LE, UTF-8}($(string(b)))"
    @test read(p, String) == s[1:9]
    @test_throws IncompleteSequenceError close(p)
    # Test that closed pipe behaves correctly even after an error
    @test eof(p)
    @test_throws ArgumentError read(p, UInt8)

    # Test stateful encoding, which output some bytes on final reset
    # with strings containing different scripts
    x = encode(s, "ISO-2022-JP")
    @test decode(x, "ISO-2022-JP") == s

    p = StringDecoder(b, "ISO-2022-JP", "UTF-8")
    b = IOBuffer(x)
    # Test that closed pipe behaves correctly
    close(p)
    @test eof(p)
    @test_throws ArgumentError read(p, UInt8)
    close(p)
end

@test_throws InvalidSequenceError encode("qwertyé€", "ASCII")
try
    encode("qwertyé€", "ASCII")
catch err
     io = IOBuffer()
     showerror(io, err)
     @test String(take!(io)) ==
        "Byte sequence 0xc3a9e282ac is invalid in source encoding or cannot be represented in target encoding"
end

@test_throws InvalidSequenceError decode(Vector{UInt8}("qwertyé€"), "ASCII")
try
    decode(Vector{UInt8}("qwertyé€"), "ASCII")
catch err
     io = IOBuffer()
     showerror(io, err)
     @test String(take!(io)) ==
         "Byte sequence 0xc3a9e282ac is invalid in source encoding or cannot be represented in target encoding"
end

let x = encode("ÄÆä", "ISO-8859-1")
    @test_throws InvalidSequenceError decode(x, "UTF-8")
    try
        decode(x, "UTF-8")
    catch err
         io = IOBuffer()
         showerror(io, err)
         @test String(take!(io)) ==
             "Byte sequence 0xc4c6e4 is invalid in source encoding or cannot be represented in target encoding"
    end
end

mktemp() do p, io
    s = "café crème"
    write(io, encode(s, "CP1252"))
    close(io)
    @test read(p, String, enc"CP1252") == s
end

@test_throws InvalidEncodingError p = StringEncoder(IOBuffer(), "nonexistent_encoding")
@test_throws InvalidEncodingError p = StringDecoder(IOBuffer(), "nonexistent_encoding")

try
    p = StringEncoder(IOBuffer(), "nonexistent_encoding", "absurd_encoding")
catch err
    @test isa(err, InvalidEncodingError)
    io = IOBuffer()
    showerror(io, err)
    @test String(take!(io)) ==
        "Conversion from absurd_encoding to nonexistent_encoding not supported by iconv implementation, check that specified encodings are correct"
end
try
    p = StringDecoder(IOBuffer(), "nonexistent_encoding", "absurd_encoding")
catch err
    @test isa(err, InvalidEncodingError)
    io = IOBuffer()
    showerror(io, err)
    @test String(take!(io)) ==
        "Conversion from nonexistent_encoding to absurd_encoding not supported by iconv implementation, check that specified encodings are correct"
end

mktemp() do path, io
    s = "a string \0チャネルパ\0ー\0トナーの選択 with embedded and trailing nuls\0\nand a second line"
    close(io)
    open(path, enc"ISO-2022-JP", "w") do io
        @test iswritable(io) && !isreadable(io)
        write(io, s)
    end

    nc = ncodeunits(first(s, 90))

    @test String(read(path, enc"ISO-2022-JP")) == s
    @test String(open(io->read(io, enc"ISO-2022-JP"), path)) == s
    @test String(open(io->read(io), path, enc"ISO-2022-JP")) == s

    @test String(read(path, 1000, enc"ISO-2022-JP")) == s
    @test String(open(io->read(io, 1000, enc"ISO-2022-JP"), path)) == s
    @test String(open(io->read(io, 1000), path, enc"ISO-2022-JP")) == s

    @test String(read(path, nc, enc"ISO-2022-JP")) == first(s, nc)
    @test String(open(io->read(io, nc, enc"ISO-2022-JP"), path)) == first(s, nc)
    @test String(open(io->read(io, nc), path, enc"ISO-2022-JP")) == first(s, nc)

    @test read(path, String, enc"ISO-2022-JP") == s
    @test open(io->read(io, String, enc"ISO-2022-JP"), path) == s
    @test open(io->read(io, String), path, enc"ISO-2022-JP") == s

    b = zeros(UInt8, nc)
    @test open(io->read!(io, b), path, enc"ISO-2022-JP") === b
    @test String(b) == first(s, 90)

    b = zeros(UInt8, nc)
    @test open(io->readbytes!(io, b), path, enc"ISO-2022-JP") == ncodeunits(s)
    @test String(b) == s

    b = zeros(UInt8, 1000)
    @test open(io->readbytes!(io, b), path, enc"ISO-2022-JP") == ncodeunits(s)
    @test length(b) == 1000
    @test String(b[1:ncodeunits(s)]) == s

    b = UInt8[]
    @test open(io->readbytes!(io, b), path, enc"ISO-2022-JP") == 0
    @test length(b) == 0

    b = zeros(UInt8, nc)
    @test open(io->readbytes!(io, b, nc), path, enc"ISO-2022-JP") == nc
    @test String(b) == first(s, 90)

    b = zeros(UInt8, 1000)
    @test open(io->readbytes!(io, b, nc), path, enc"ISO-2022-JP") == nc
    @test length(b) == 1000
    @test String(b[1:nc]) == first(s, 90)

    b = UInt8[]
    @test open(io->readbytes!(io, b, nc), path, enc"ISO-2022-JP") == nc
    @test String(b) == first(s, 90)

    b = UInt8[]
    open(io->while !eof(io); push!(b, read(io, UInt8)) end, path, enc"ISO-2022-JP")
    @test String(b) == s

    @test readuntil(path, enc"ISO-2022-JP", '\0') == "a string "
    @test open(io->readuntil(io, enc"ISO-2022-JP", '\0'), path) == "a string "
    @test open(io->readuntil(io, enc"ISO-2022-JP", '\0', keep=true), path) == "a string \0"
    @test readuntil(path, enc"ISO-2022-JP", "チャ") == "a string \0"
    @test open(io->readuntil(io, enc"ISO-2022-JP", "チャ"), path) == "a string \0"
    @test open(io->readuntil(io, enc"ISO-2022-JP", "チャ", keep=true), path) == "a string \0チャ"

    @test readline(path, enc"ISO-2022-JP") == split(s, '\n')[1]
    @test readline(path, enc"ISO-2022-JP", keep=true) == split(s, '\n', )[1] * '\n'
    @test open(readline, path, enc"ISO-2022-JP") == split(s, '\n')[1]

    a = readlines(path, enc"ISO-2022-JP")
    b = open(readlines, path, enc"ISO-2022-JP")
    c = collect(eachline(path, enc"ISO-2022-JP"))
    d = open(io->collect(eachline(io, enc"ISO-2022-JP")), path)
    @test a[1] == b[1] == c[1] == d[1] == split(s, '\n')[1]
    @test a[2] == b[2] == c[2] == d[2] == split(s, '\n')[2]

    a = readlines(path, enc"ISO-2022-JP", keep=true)
    b = open(io->readlines(io, keep=true), path, enc"ISO-2022-JP")
    c = collect(eachline(path, enc"ISO-2022-JP", keep=true))
    d = open(io->collect(eachline(io, enc"ISO-2022-JP", keep=true)), path)
    @test a[1] == b[1] == c[1] == d[1] == split(s, '\n')[1] * '\n'
    @test a[2] == b[2] == c[2] == d[2] == split(s, '\n')[2]

    # Test alternative syntaxes for open()
    open(path, enc"ISO-2022-JP", "r") do io
        @test isreadable(io) && !iswritable(io)
        @test read(io, String) == s
    end
    open(path, enc"ISO-2022-JP", true, false, false, false, false) do io
        @test isreadable(io) && !iswritable(io)
        @test read(io, String) == s
    end
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", "r+")
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", "w+")
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", "a+")
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", true, true, false, false, false)
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", true, false, false, false, true)
end

@testset "Bytes written only when closing with Shift JIS (issue #49)" begin
    for enc in ("SHIFT_JIS", "SHIFT_JISX0213")
        @test encode("ク", enc) == [0x83, 0x4e]

        b = IOBuffer()
        s = StringEncoder(b, enc)
        write(s, "ク")
        close(s)
        @test take!(b) == [0x83, 0x4e]
    end
end

## Test encodings support
b = IOBuffer()
show(b, enc"UTF-8")
@test String(take!(b)) == "UTF-8 string encoding"
@test string(enc"UTF-8") == "UTF-8"

encodings_list = encodings()
@test "ASCII" in encodings_list
@test "UTF-8" in encodings_list

nothing


text:
using ImageCore, FixedPointNumbers, Colors, ColorVectorSpace
using Test

@testset "map" begin
    @testset "clamp01" begin
        @test clamp01(0.1) === 0.1
        @test clamp01(-0.1) === 0.0
        @test clamp01(0.9) === 0.9
        @test clamp01(1.1) === 1.0
        @test clamp01(Inf) === 1.0
        @test clamp01(-Inf) === 0.0
        @test isnan(clamp01(NaN))
        @test clamp01(N0f8(0.1)) === N0f8(0.1)
        @test clamp01(N4f12(1.2)) === N4f12(1.0)
        @test clamp01(Gray(-0.2)) === Gray(0.0)
        @test clamp01(Gray(0.7)) === Gray(0.7)
        @test clamp01(Gray(N4f12(1.2))) === Gray(N4f12(1.0))
        @test clamp01(RGB(0.2,-0.2,1.8)) === RGB(0.2,0.0,1.0)
        A = [-1.2,0.4,800.3]
        f = takemap(clamp01, A)
        fA = f.(A)
        @test eltype(fA) == Float64
        @test fA == [0, 0.4, 1]
        f = takemap(clamp01, N0f8, A)
        fA = f.(A)
        @test eltype(fA) == N0f8
        @test fA == [N0f8(0), N0f8(0.4), N0f8(1)]

        A = [-1.2,0.4,800.3]
        clamp01!(A)
        @test eltype(A) == Float64
        @test A == [0, 0.4, 1]
    end

    @testset "clamp01nan" begin
        @test clamp01nan(0.1) === 0.1
        @test clamp01nan(-0.1) === 0.0
        @test clamp01nan(0.9) === 0.9
        @test clamp01nan(1.1) === 1.0
        @test clamp01nan(Inf) === 1.0
        @test clamp01nan(-Inf) === 0.0
        @test clamp01nan(NaN) === 0.0
        @test clamp01nan(N0f8(0.1)) === N0f8(0.1)
        @test clamp01nan(N4f12(1.2)) === N4f12(1.0)
        @test clamp01nan(Gray(-0.2)) === Gray(0.0)
        @test clamp01nan(Gray(0.7)) === Gray(0.7)
        @test clamp01nan(Gray(NaN32)) === Gray(0.0f0)
        @test clamp01nan(Gray(N4f12(1.2))) === Gray(N4f12(1.0))
        @test clamp01nan(RGB(0.2,-0.2,1.8)) === RGB(0.2,0.0,1.0)
        @test clamp01nan(RGB(0.2,NaN,1.8)) === RGB(0.2,0.0,1.0)
        A = [-1.2,0.4,-Inf,NaN,Inf,800.3]
        f = takemap(clamp01nan, A)
        fA = f.(A)
        @test eltype(fA) == Float64
        @test fA == [0, 0.4, 0, 0, 1, 1]
        f = takemap(clamp01nan, N0f8, A)
        fA = f.(A)
        @test eltype(fA) == N0f8
        @test fA == [N0f8(0), N0f8(0.4), N0f8(0), N0f8(0), N0f8(1), N0f8(1)]

        A = [-1.2,0.4,-Inf,NaN,Inf,800.3]
        clamp01nan!(A)
        @test eltype(A) == Float64
        @test A == [0, 0.4, 0, 0, 1, 1]
    end

    @testset "scaleminmax" begin
        A = [0, 1, 100, 1000, 2000, -7]
        target = map(x->clamp(x, 0, 1), A/1000)
        for (f, tgt) in ((scaleminmax(0, 1000), target),
                          (scaleminmax(0, 1000.0), target),
                          (scaleminmax(N0f8, 0, 1000), N0f8.(target)),
                          (scaleminmax(N0f8, 0, 1000.0), N0f8.(target)),
                          (scaleminmax(Gray, 0, 1000), Gray{Float64}.(target)),
                          (scaleminmax(Gray{N0f8}, 0, 1000.0), Gray{N0f8}.(target)))
            fA = @inferred(map(f, A))
            @test fA == tgt
            @test eltype(fA) == eltype(tgt)
        end
        B = A.+10
        f = scaleminmax(10, 1010)
        @test f.(B) == target
        A = [0, 1, 100, 1000]
        target = A/1000
        for (f, tgt) in ((takemap(scaleminmax, A), target),
                          (takemap(scaleminmax, N0f8, A), N0f8.(target)),
                          (takemap(scaleminmax, Gray{N0f8}, A), Gray{N0f8}.(target)))
            fA = @inferred(map(f, A))
            @test fA == tgt
            @test eltype(fA) == eltype(tgt)
        end
        A = [Gray(-0.1),Gray(0.1)]
        f = scaleminmax(Gray, -0.1, 0.1)
        @test f.(A) == [Gray(0.0),Gray(1.0)]
        A = reinterpretc(RGB, [0.0 128.0; 255.0 0.0; 0.0 0.0])
        f = scaleminmax(RGB, 0, 255)
        @test f.(A) == [RGB(0,1.0,0), RGB(128/255,0,0)]
        f = scaleminmax(RGB{N0f8}, 0, 255)
        @test f.(A) == [RGB(0,1,0), RGB{N0f8}(128/255,0,0)]
        f = takemap(scaleminmax, A)
        @test f.(A) == [RGB(0,1.0,0), RGB(128/255,0,0)]
        f = takemap(scaleminmax, RGB{N0f8}, A)
        @test f.(A) == [RGB(0,1,0), RGB{N0f8}(128/255,0,0)]
    end

    @testset "scalesigned" begin
        A = [-100,1000]
        target = A/1000
        for (f, tgt) in ((scalesigned(1000), target),
                         (scalesigned(-1000, 0, 1000), target),
                         (scalesigned(-1000, 0, 1000.0), target),
                         (takemap(scalesigned, A), target))
            fA = f.(A)
            @test fA == tgt
            @test eltype(fA) == eltype(tgt)
        end
    end

    @testset "colorsigned" begin
        g, w, m = colorant"green1", colorant"white", colorant"magenta"
        for f in (colorsigned(),
                  colorsigned(g, m),
                  colorsigned(g, w, m))
            @test f(-1) == g
            @test f( 0) == w
            @test f( 1) == m
            @test f(-0.5) ≈ mapc(N0f8, 0.5g+0.5w)
            @test f( 0.5) ≈ mapc(N0f8, 0.5w+0.5m)
            @test f(-0.25) ≈ mapc(N0f8, 0.25g+0.75w)
            @test f( 0.75) ≈ mapc(N0f8, 0.75m+0.25w)
            @test f(Gray(0.75)) ≈ mapc(N0f8, 0.75m+0.25w)
        end
        g, w, m = RGBA(0.,1.,0.,1.), RGBA(1.,1.,1.,0.), RGBA(1.,0.,1.,1.)
        f = colorsigned(g, w, m)
        @test f(-1) == g
        @test f( 0) == w
        @test f( 1) == m
        @test f(-0.5) ≈ RGBA(0.5,1.0,0.5,0.5)
        @test f(-0.5) ≈ 0.5g+0.5w
        @test f( 0.5) ≈ 0.5w+0.5m
        @test f(-0.25) ≈ 0.25g+0.75w
        @test f( 0.75) ≈ 0.75m+0.25w
        w = RGBA{Float64}(colorant"white")
        f = colorsigned(g, m)
        @test f(-1) == g
        @test f( 0) == w
        @test f( 1) == m
        @test f(-0.5) ≈ 0.5g+0.5w
        @test f( 0.5) ≈ 0.5w+0.5m
        @test f(-0.25) ≈ 0.25g+0.75w
        @test f( 0.75) ≈ 0.75m+0.25w
    end
end

nothing


text:
#Tests des packages
include("../ALAS/include.jl")

using MPCCmod
using OutputRelaxationmod
using Penalty
using DDirection
using LineSearch
using AlgoSetmod

using CUTEst
using NLPModels

using BenchmarkProfiles

problems_unconstrained = CUTEst.select(contype="unc")
problems_linconstrained = CUTEst.select(contype="linear",min_var=1, max_var=10)
problems_boundconstrained = CUTEst.select(contype="bounds")
problems_quadconstrained = CUTEst.select(contype="quadratic")
problems_genconstrained = CUTEst.select(contype="general")
problems=CUTEst.select(min_var=1, max_var=10,min_con=0,max_con=10)
m=length(problems)

#USER Choice :
#NaN MAXLIKA

pb=problems
n=100

algo1=AlgoSetmod.AlgoSet(Penalty.Quadratic,DDirection.NwtdirectionSpectral,LineSearch.ArmijoWolfe)
algo2=AlgoSetmod.AlgoSet(Penalty.Quadratic,DDirection.CGHZ,LineSearch.ArmijoWolfeHZ)
algo3=AlgoSetmod.AlgoSet(Penalty.Lagrangian,DDirection.NwtdirectionSpectral,LineSearch.ArmijoWolfe)

#....

T=zeros(n,3)

@printf("Start of a benchmark test from CUTEst with %i problems.\n",n)

for i=1:min(n,m)
 nlp = CUTEstModel(pb[i])
 @printf("CUTEst %i: %s \n",i,nlp.meta.name)

 #resolution avec ALAS
 exemple_nlp=MPCCmod.MPCC(nlp,algo1)
 xb,fb,orb,nb_eval = MPCCsolve.solve(exemple_nlp)
 @printf("%s ",orb.solve_message)
 println(nb_eval)
@show orb.solve_message=="Success" sum(nb_eval)
 T[i,1]=orb.solve_message=="Success"?sum(nb_eval):Inf
@show T[i,1]

 exemple_nlp=MPCCmod.MPCC(nlp,algo2)
 xb,fb,orb,nb_eval = MPCCsolve.solve(exemple_nlp)
 @printf("%s ",orb.solve_message)
 println(nb_eval)
 T[i,2]=orb.solve_message=="Success"?sum(nb_eval):Inf

 exemple_nlp=MPCCmod.MPCC(nlp,algo3)
 xb,fb,orb,nb_eval = MPCCsolve.solve(exemple_nlp)
 @printf("%s ",orb.solve_message)
 println(nb_eval)
 T[i,3]=orb.solve_message=="Success"?sum(nb_eval):Inf

 finalize(nlp)
end
T1=copy(T)
performance_profile(T,["Solver 1","Sovler 2","Sovler 3"], title="Celebrity Death Match")