text:
# Copyright (c) 2013: Iain Dunning, Miles Lubin, and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

mutable struct IpoptProblem
    ipopt_problem::Ptr{Cvoid}   # Reference to the internal data structure
    n::Int                      # Num vars
    m::Int                      # Num cons
    x::Vector{Float64}          # Starting and final solution
    g::Vector{Float64}          # Final constraint values
    mult_g::Vector{Float64}     # lagrange multipliers on constraints
    mult_x_L::Vector{Float64}   # lagrange multipliers on lower bounds
    mult_x_U::Vector{Float64}   # lagrange multipliers on upper bounds
    obj_val::Float64            # Final objective
    status::Cint                # Final status
    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h::Union{Function,Nothing}
    intermediate::Union{Function,Nothing}
end

function _Eval_F_CB(
    n::Cint,
    x_ptr::Ptr{Float64},
    x_new::Cint,
    obj_value::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    x = unsafe_wrap(Array, x_ptr, Int(n))
    if x_new == Cint(1)
        prob.x .= x
    end
    new_obj = convert(Float64, prob.eval_f(x))::Float64
    unsafe_store!(obj_value, new_obj)
    return Cint(1)
end

function _Eval_Grad_F_CB(
    n::Cint,
    x_ptr::Ptr{Float64},
    # A Bool indicating if `x` is a new point. We don't make use of this.
    ::Cint,
    grad_f::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    new_grad_f = unsafe_wrap(Array, grad_f, Int(n))
    x = unsafe_wrap(Array, x_ptr, Int(n))
    prob.eval_grad_f(x, new_grad_f)
    return Cint(1)
end

function _Eval_G_CB(
    n::Cint,
    x_ptr::Ptr{Float64},
    x_new::Cint,
    m::Cint,
    g_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    new_g = unsafe_wrap(Array, g_ptr, Int(m))
    x = unsafe_wrap(Array, x_ptr, Int(n))
    if x_new == Cint(1)
        prob.x .= x
    end
    prob.eval_g(x, new_g)
    return Cint(1)
end

function _Eval_Jac_G_CB(
    n::Cint,
    x_ptr::Ptr{Float64},
    ::Cint,
    ::Cint,
    nele_jac::Cint,
    iRow::Ptr{Cint},
    jCol::Ptr{Cint},
    values_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    x = unsafe_wrap(Array, x_ptr, Int(n))
    rows = unsafe_wrap(Array, iRow, Int(nele_jac))
    cols = unsafe_wrap(Array, jCol, Int(nele_jac))
    if values_ptr == C_NULL
        prob.eval_jac_g(x, rows, cols, nothing)
    else
        values = unsafe_wrap(Array, values_ptr, Int(nele_jac))
        prob.eval_jac_g(x, rows, cols, values)
    end
    return Cint(1)
end

function _Eval_H_CB(
    n::Cint,
    x_ptr::Ptr{Float64},
    ::Cint,
    obj_factor::Float64,
    m::Cint,
    lambda_ptr::Ptr{Float64},
    ::Cint,
    nele_hess::Cint,
    iRow::Ptr{Cint},
    jCol::Ptr{Cint},
    values_ptr::Ptr{Float64},
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    if prob.eval_h === nothing
        # No hessian. Return FALSE for failure.
        return Cint(0)
    end
    x = unsafe_wrap(Array, x_ptr, Int(n))
    lambda = unsafe_wrap(Array, lambda_ptr, Int(m))
    rows = unsafe_wrap(Array, iRow, Int(nele_hess))
    cols = unsafe_wrap(Array, jCol, Int(nele_hess))
    if values_ptr == C_NULL
        prob.eval_h(x, rows, cols, obj_factor, lambda, nothing)
    else
        values = unsafe_wrap(Array, values_ptr, Int(nele_hess))
        prob.eval_h(x, rows, cols, obj_factor, lambda, values)
    end
    return Cint(1)  # Return TRUE for success.
end

function _Intermediate_CB(
    alg_mod::Cint,
    iter_count::Cint,
    obj_value::Float64,
    inf_pr::Float64,
    inf_du::Float64,
    mu::Float64,
    d_norm::Float64,
    regularization_size::Float64,
    alpha_du::Float64,
    alpha_pr::Float64,
    ls_trials::Cint,
    user_data::Ptr{Cvoid},
)
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    ret = prob.intermediate(
        alg_mod,
        iter_count,
        obj_value,
        inf_pr,
        inf_du,
        mu,
        d_norm,
        regularization_size,
        alpha_du,
        alpha_pr,
        ls_trials,
    )
    # Return TRUE if the optimization should continue.
    return ret ? Cint(1) : Cint(0)
end

function CreateIpoptProblem(
    n::Int,
    x_L::Vector{Float64},
    x_U::Vector{Float64},
    m::Int,
    g_L::Vector{Float64},
    g_U::Vector{Float64},
    nele_jac::Int,
    nele_hess::Int,
    eval_f,
    eval_g,
    eval_grad_f,
    eval_jac_g,
    eval_h,
)
    @assert n == length(x_L) == length(x_U)
    @assert m == length(g_L) == length(g_U)
    eval_f_cb = @cfunction(
        _Eval_F_CB,
        Cint,
        (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cvoid}),
    )
    eval_g_cb = @cfunction(
        _Eval_G_CB,
        Cint,
        (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Cvoid}),
    )
    eval_grad_f_cb = @cfunction(
        _Eval_Grad_F_CB,
        Cint,
        (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cvoid}),
    )
    eval_jac_g_cb = @cfunction(
        _Eval_Jac_G_CB,
        Cint,
        (
            Cint,
            Ptr{Float64},
            Cint,
            Cint,
            Cint,
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Float64},
            Ptr{Cvoid},
        ),
    )
    eval_h_cb = @cfunction(
        _Eval_H_CB,
        Cint,
        (
            Cint,
            Ptr{Float64},
            Cint,
            Float64,
            Cint,
            Ptr{Float64},
            Cint,
            Cint,
            Ptr{Cint},
            Ptr{Cint},
            Ptr{Float64},
            Ptr{Cvoid},
        ),
    )
    ipopt_problem = ccall(
        (:CreateIpoptProblem, libipopt),
        Ptr{Cvoid},
        (
            Cint,
            Ptr{Float64},
            Ptr{Float64},
            Cint,
            Ptr{Float64},
            Ptr{Float64},
            Cint,
            Cint,
            Cint,
            Ptr{Cvoid},
            Ptr{Cvoid},
            Ptr{Cvoid},
            Ptr{Cvoid},
            Ptr{Cvoid},
        ),
        n,
        x_L,
        x_U,
        m,
        g_L,
        g_U,
        nele_jac,
        nele_hess,
        1,  # 1 = Fortran style indexing
        eval_f_cb,
        eval_g_cb,
        eval_grad_f_cb,
        eval_jac_g_cb,
        eval_h_cb,
    )
    if ipopt_problem == C_NULL
        if n == 0
            error(
                "IPOPT: Failed to construct problem because there are 0 " *
                "variables. If you intended to construct an empty problem, " *
                "one work-around is to add a variable fixed to 0.",
            )
        else
            error("IPOPT: Failed to construct problem for some unknown reason.")
        end
    end
    prob = IpoptProblem(
        ipopt_problem,
        n,
        m,
        zeros(Float64, n),
        zeros(Float64, m),
        zeros(Float64, m),
        zeros(Float64, n),
        zeros(Float64, n),
        0.0,
        0,
        eval_f,
        eval_g,
        eval_grad_f,
        eval_jac_g,
        eval_h,
        nothing,
    )
    finalizer(FreeIpoptProblem, prob)
    return prob
end

function FreeIpoptProblem(prob::IpoptProblem)
    ccall(
        (:FreeIpoptProblem, libipopt),
        Cvoid,
        (Ptr{Cvoid},),
        prob.ipopt_problem,
    )
    return
end

function AddIpoptStrOption(prob::IpoptProblem, keyword::String, value::String)
    if !(isascii(keyword) && isascii(value))
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:AddIpoptStrOption, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}),
        prob.ipopt_problem,
        keyword,
        value,
    )
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
    return
end

function AddIpoptNumOption(prob::IpoptProblem, keyword::String, value::Float64)
    if !isascii(keyword)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:AddIpoptNumOption, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Float64),
        prob.ipopt_problem,
        keyword,
        value,
    )
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
    return
end

function AddIpoptIntOption(prob::IpoptProblem, keyword::String, value::Integer)
    if !isascii(keyword)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:AddIpoptIntOption, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint),
        prob.ipopt_problem,
        keyword,
        value,
    )
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
    return
end

function OpenIpoptOutputFile(
    prob::IpoptProblem,
    file_name::String,
    print_level::Int,
)
    if !isascii(file_name)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall(
        (:OpenIpoptOutputFile, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint),
        prob.ipopt_problem,
        file_name,
        print_level,
    )
    if ret == 0
        error("IPOPT: Couldn't open output file.")
    end
    return
end

function SetIpoptProblemScaling(
    prob::IpoptProblem,
    obj_scaling::Float64,
    x_scaling::Union{Ptr{Cvoid},Vector{Float64}},
    g_scaling::Union{Ptr{Cvoid},Vector{Float64}},
)
    ret = ccall(
        (:SetIpoptProblemScaling, libipopt),
        Cint,
        (Ptr{Cvoid}, Float64, Ptr{Float64}, Ptr{Float64}),
        prob.ipopt_problem,
        obj_scaling,
        x_scaling,
        g_scaling,
    )
    if ret == 0
        error("IPOPT: Error setting problem scaling.")
    end
    return
end

function SetIntermediateCallback(prob::IpoptProblem, intermediate::Function)
    intermediate_cb = @cfunction(
        _Intermediate_CB,
        Cint,
        (
            Cint,
            Cint,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Cint,
            Ptr{Cvoid},
        ),
    )
    ret = ccall(
        (:SetIntermediateCallback, libipopt),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}),
        prob.ipopt_problem,
        intermediate_cb,
    )
    if ret == 0
        error("IPOPT: Something went wrong setting the intermediate callback.")
    end
    prob.intermediate = intermediate
    return
end

function IpoptSolve(prob::IpoptProblem)
    final_objval = Ref{Cdouble}(0.0)
    ret = ccall(
        (:IpoptSolve, libipopt),
        Cint,
        (
            Ptr{Cvoid},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Float64},
            Ptr{Cvoid},
        ),
        prob.ipopt_problem,
        prob.x,
        prob.g,
        final_objval,
        prob.mult_g,
        prob.mult_x_L,
        prob.mult_x_U,
        pointer_from_objref(prob),
    )
    prob.obj_val = final_objval[]
    prob.status = ret
    return prob.status
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec ${JULIA:-julia} --startup-file=no -i \
    -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" \
    "$@"
=#

function start_repl(;
        interactive::Bool = true,
        quiet::Bool = true,
        banner::Bool = false,
        history_file::Bool = true,
        color_set::Bool = false,
        )
    if !(stdout isa Base.TTY)
        error("stdout is not a TTY")
    elseif !(stdin isa Base.TTY)
        error("stdin is not a TTY")
    end
    was_interactive = Base.is_interactive
    try
        # Required for Pkg.__init__ to setup the REPL mode:
        Base.eval(:(is_interactive = $interactive))

        @info "Starting REPL..."
        Base.run_main_repl(interactive, quiet, banner, history_file, color_set)
    finally
        Base.eval(:(is_interactive = $was_interactive))
    end
end

csv_path, = ARGS

let io = open(csv_path, "w")
    ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), io.handle)
    try
        @info "Loading IPython..."
        using IPython
        start_repl()
    finally
        ccall(:jl_dump_compiles, Nothing, (Ptr{Nothing},), C_NULL)
        close(io)
    end
end
exit()

ext: jl


text:
"""
JuMP interaction
"""

"""
Pass JuMP model setup to Ipopt - a nonlinear solver

See Ipopt documentation
"""
function initialize_JuMP_model(linear_solver;print_level=5,tol=1e-3,max_iter=1000,
  constr_viol_tol=1e-3,obj_scaling_factor=1.0)

# model = Model(() -> Ipopt.Optimizer(print_level=print_level,tol=tol,max_iter=max_iter,
#   acceptable_constr_viol_tol=acceptable_constr_viol_tol,linear_solver=linear_solver))

model = Model(Ipopt.Optimizer)
JuMP.set_optimizer_attributes(model,"print_level"=>print_level,"tol"=>tol,
  "max_iter"=>max_iter,"constr_viol_tol"=>constr_viol_tol,"linear_solver"=>linear_solver,
  "obj_scaling_factor"=>obj_scaling_factor)
return model

end
"""

Let JuMP know about variables

"""
function define_variables!(n_t_segs,n_data_pts,model,start_point_frac;upper_bound=410.0,lower_bound=273.0,T0=385.0,autodiff=true)
  if isa(upper_bound,Number)
  T = @variable(model, [i=1:n_t_segs],base_name="T",lower_bound=lower_bound,
    upper_bound=map(constrain_upper,upper_bound,i,n_t_segs),
    start=map(constraint_func,i,maximum(upper_bound),maximum(lower_bound),start_point_frac[i]))
  elseif isa(upper_bound,Array)
    T = @variable(model, [i=1:n_t_segs],base_name="T",lower_bound=lower_bound[i],
      upper_bound=map(constrain_upper,upper_bound[i],i,n_t_segs),
      start=map(constraint_func,i,maximum(upper_bound),maximum(lower_bound),start_point_frac[i]))
  end

  @NLparameter(model, set_dev[i = 1:n_data_pts] == 0.0 * i)
return T,set_dev
end

"""
Let JuMP know about constraints
"""
function rdaam_define_constraints(model,density,set_dev,set_val,L,times,T,U238,U235,Th232,U238_V,U235_V,Th232_V;
   n_iter=20.0,E_L=122.3*1e3,c0=0.39528,c1=0.01073,c2=-65.12969,c3=-7.91715,
   alpha=0.04672,rmr0=0.79,logD0_a2=log10(exp(9.733)),E_trap=34*1e3,omega=1e-22,
   psi=1e-13,L_dist=8.1*1e-4,eta_q=0.91,R=Rjoules)
n_data_constraints=length(U238)
kappa = 1.04-rmr0

for j in eachindex(U238)
    if (j<n_data_constraints) || (n_data_constraints==1)
      my_constr = @NLconstraint(model, [i = j:j],(rdaam_forward_diffusion(
                  alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist,psi,omega,E_trap,
                  R,E_L,logD0_a2,n_iter,U238[i],U238_V[i],U235[i],U235_V[i],
                  Th232[i],Th232_V[i],L,times...,T...
                  )-set_val[i]+set_dev[i])^2<=(set_val[i]*0.005)^2)

    elseif (j==n_data_constraints)
    my_constr2 = @NLconstraint(model, [i = j:j],(rdaam_forward_diffusion(
                   alpha,c0,c1,c2,c3,rmr0,eta_q,L_dist,psi,omega,E_trap,
                   R,E_L,logD0_a2,n_iter,U238[i],U238_V[i],U235[i],U235_V[i],Th232[i],Th232_V[i],L,
                   times...,T...)-set_val[i]+set_dev[i])==(set_val[i]*0.0))
    end
    return my_constr
end

return nothing
end

"""
```julia
constraint_func(x,upper,lower,init_scaling_young,init_scaling_old)
```
"""
function constraint_func(ii,upper,lower,scaling_frac)

#randval = rand(1)[1]
x = (upper-lower)*scaling_frac+lower

return x
end

function constrain_upper(x,i,max_i)
x_orig = copy(x)
if i==max_i
  x=273.0
else
  x=x_orig
end

return x
end


text:
export gaussian, xavier, xavier_uniform, xavier_normal, bilinear
using Random: randn, rand

"""

    gaussian(a...; mean=0.0, std=0.01)

Return a Gaussian array with a given mean and standard deviation.  The
`a` arguments are passed to `randn`.

"""
function gaussian(a...; mean=0.0, std=0.01)
    r = randn(a...)
    T = eltype(r)
    r .* T(std) .+ T(mean)
end

"""

    xavier_uniform(a...; gain=1)
    xavier(a...; gain=1)

Return uniform random weights in the range `± gain * sqrt(6 / (fanin + fanout))`.  The `a`
arguments are passed to `rand` to specify type and dimensions.  See ([Glorot and Bengio
2010](http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf)) or the [PyTorch
docs](https://pytorch.org/docs/stable/nn.init.html#torch.nn.init.xavier_uniform_) for a
description.  The function implements equation (16) of the referenced paper. Also known as
Glorot initialization. The function `xavier` is an alias for `xavier_uniform`. See also
`xavier_normal`.

"""
xavier, xavier_uniform

function xavier_uniform(a...; gain=1)
    w = rand(a...)
    if ndims(w) == 1
        fanout = 1
        fanin = length(w)
    elseif ndims(w) == 2
        fanout = size(w,1)
        fanin = size(w,2)
    else
        # if a is (3,3,16,8), then there are 16 input channels and 8 output channels
        # fanin = 3*3*16 = (3*3*16*8) ÷ 8
        # fanout = 3*3*8 = (3*3*16*8) ÷ 16
        fanin = div(length(w),  a[end])
        fanout = div(length(w), a[end-1])
    end
    s = convert(eltype(w), gain*sqrt(6 / (fanin + fanout)))
    return 2s .* w .- s
end

const xavier = xavier_uniform

"""

    xavier_normal(a...; gain=1)

Return normal distributed random weights with mean 0 and std `gain * sqrt(2 / (fanin +
fanout))`.  The `a` arguments are passed to `rand`.  See ([Glorot and Bengio
2010](http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf)) and [PyTorch
docs](https://pytorch.org/docs/stable/nn.init.html#torch.nn.init.xavier_normal_) for a
description. Also known as Glorot initialization. See also `xavier_uniform`.

"""
function xavier_normal(a...; gain=1)
    w = randn(a...)
    if ndims(w) == 1
        fanout = 1
        fanin = length(w)
    elseif ndims(w) == 2
        fanout = size(w,1)
        fanin = size(w,2)
    else
        # if a is (3,3,16,8), then there are 16 input channels and 8 output channels
        # fanin = 3*3*16 = (3*3*16*8) ÷ 8
        # fanout = 3*3*8 = (3*3*16*8) ÷ 16
        fanin = div(length(w),  a[end])
        fanout = div(length(w), a[end-1])
    end
    s = convert(eltype(w), gain*sqrt(2 / (fanin + fanout)))
    return s .* w
end

"""

Bilinear interpolation filter weights; used for initializing deconvolution layers.

Adapted from https://github.com/shelhamer/fcn.berkeleyvision.org/blob/master/surgery.py#L33

Arguments:

`T` : Data Type

`fw`: Width upscale factor

`fh`: Height upscale factor

`IN`: Number of input filters

`ON`: Number of output filters

Example usage:

w = bilinear(Float32,2,2,128,128)

"""
function bilinear(T,fw,fh,IN,ON)
    @assert fw == fh "Filter must be square"
    @assert IN == ON "Number of input and output filters must be equal"
    f=fw; N=IN;

    sz = 2*f-f%2
    center = f-0.5
    if sz%2 == 1
        center = f-1
    end
    r = range(0,stop=sz-1,length=sz)
    c = range(0,stop=sz-1,length=sz)'

    kernel = (1 .- abs.(r .- center) ./ f) .* (1 .- abs.(c .- center) ./ f)
    w = zeros(T,sz,sz,N,N);
    for i=1:N
        w[:,:,i,i] = kernel
    end
    return w
end


text:
using Documenter, StanBase

makedocs(
    modules = [StanBase],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "StanBase.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/goedman/StanBase.jl.git",
)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    readgeotable(args; coordnames=(:x,:y,:z), kwargs)

Read data from disk using `CSV.File`, optionally specifying
the columns `coordnames` with spatial coordinates.

The arguments `args` and keyword arguments `kwargs` are
forwarded to the `CSV.File` function, please check their
documentation for more details.
"""
readgeotable(args...; coordnames=(:x,:y,:z), kwargs...) =
  georef(DataFrame!(CSV.File(args...; kwargs...)), coordnames)

"""
    spheredir(θ, φ)

Returns the 3D direction given polar angle `θ` and
azimuthal angle `φ` in degrees according to the ISO
convention.
"""
function spheredir(theta, phi)
  θ, φ = deg2rad(theta), deg2rad(phi)
  SVector(sin(θ)*cos(φ), sin(θ)*sin(φ), cos(θ))
end


text:
using JuMP, GLPK
m = Model(with_optimizer(GLPK.Optimizer))

@variable(m, 0 <= x <= 2 )
@variable(m, 0 <= y <= 30 )

@objective(m, Max, 5x + 3*y )

@constraint(m, 1x + 5y <= 3.0 )

JuMP.optimize!(m)
println("Objective value: ", JuMP.objective_value(m))
println("x = ", JuMP.value(x))
println("y = ", JuMP.value(y))

max_stars_repo_name: zzsza/TIL

max_issues_repo_name: zzsza/TIL

max_forks_repo_name: zzsza/TIL


text:
__precompile__(false)

module FrustumManager

using MathManager

using LinearAlgebra

"""
TODO
"""
mutable struct Plane3D
	mPoint::Vec3f
	mNormal::Vec3f
	d::Float32

  """
  TODO
  """
  Plane3D() = new(Vec3f(),Vec3f(),0)

  """
  TODO
  """
  function Plane3D(mPoint::Vec3f, mNormal::Vec3f)
    mNormal = normalize(mNormal)
    d= -(dot(mNormal, mPoint))
    new(mPoint,mNormal,d)
  end

  """
  TODO
  """
  function Plane3D(lv1::Vec3f, lv2::Vec3f, lv3::Vec3f)
    mNormal = normalize(cross(lv3 - lv2, lv1 - lv2))
    mPoint = (lv1 + lv3) / 2 #lv2
    d= -(dot(mNormal, mPoint))
    new(mPoint,mNormal,d)
  end

  """
  TODO
  """
  function Plane3D(a::Float32, b::Float32, c::Float32, d::Float32)
    mPoint = Vec3f(0,0,0)

    # Set the normal vector
    mNormal = Vec3f(a, b, c)

    # Compute the length of the vector
    lLength = length(mNormal)

    # Normalize the vector
    mNormal = Vec3f(a / lLength, b / lLength, c / lLength)

    # And divide d by the length as well
    d = d / lLength
    new(mPoint,mNormal,d)
  end
end

export Plane3D

"""
TODO
"""
function Plane3D_Rect(lbn::Vec3f, rtf::Vec3f)
  v = Vec3f(rtf.x,rtf.y,lbn.z)
  mNormal = normalize(cross(rtf - v, lbn - v))
  mPoint = (lbn + rtf) / 2
  d= -(dot(mNormal, mPoint))
  Plane3D(mPoint,mNormal,d)
end

export Plane3D_Rect

"""
TODO
"""
GetPointDistance(this::Plane3D, lPoint::Vec3f) = this.d + dot(this.mNormal, lPoint)

export GetPointDistance

"""
TODO
"""
mutable struct Frustum
  planes::Dict{Symbol,Plane3D} #Array{Plane3D,1}
  box::Dict{Symbol,Plane3D}
  pos::Dict{Symbol,Vec3f}

  nearDistance::Float32
  farDistance::Float32
  nearWidth::Float32
  nearHeight::Float32
  farWidth::Float32
  farHeight::Float32
  ratio::Float32
  angle::Float32
  tang::Float32
  centerRadius::Float32

  nearTopLeft::Vec3f
  nearTopRight::Vec3f
  nearBottomLeft::Vec3f
  nearBottomRight::Vec3f
  farTopLeft::Vec3f
  farTopRight::Vec3f
  farBottomLeft::Vec3f
  farBottomRight::Vec3f

  nearfarTopLeft::Vec3f
  nearfarTopRight::Vec3f
  nearfarBottomLeft::Vec3f
  nearfarBottomRight::Vec3f

  Frustum() = new(
  Dict(:FRUSTUM_TOP=>Plane3D(),:FRUSTUM_BOTTOM=>Plane3D(),:FRUSTUM_LEFT=>Plane3D(),:FRUSTUM_RIGHT=>Plane3D(),:FRUSTUM_NEAR=>Plane3D(),:FRUSTUM_FAR=>Plane3D()),
  Dict(:BOX_TOP=>Plane3D(),:BOX_BOTTOM=>Plane3D(),:BOX_LEFT=>Plane3D(),:BOX_RIGHT=>Plane3D()),
  Dict(:CENTER=>Vec3f(),:CAMERA=>Vec3f(),:TARGET=>Vec3f()),
  0,0,0,0,0,0,0,0,0,0,
  Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f(),Vec3f()) #Array{Plane3D,1}(6)
end

export Frustum

getNearPosition(fstm::Frustum) = Vec3f((fstm.nearBottomLeft + fstm.nearTopRight) / 2)
getFarPosition(fstm::Frustum) = Vec3f((fstm.farBottomLeft + fstm.farTopRight) / 2)
getPosition(fstm::Frustum) = Vec3f((fstm.nearBottomLeft + fstm.farTopRight) / 2)

export getNearPosition
export getFarPosition
export getPosition

"""
TODO
"""
function getBox(this::Frustum)
  [
    this.nearfarTopLeft, this.farTopLeft, this.farTopRight,
    this.nearfarTopLeft, this.farTopRight, this.farTopLeft,

    this.nearfarTopLeft, this.nearfarTopRight, this.farTopRight,
    this.nearfarTopLeft, this.farTopRight, this.nearfarTopRight,

    this.nearfarBottomLeft, this.farBottomLeft, this.farBottomRight,
    this.nearfarBottomLeft, this.farBottomRight, this.farBottomLeft,

    this.nearfarBottomLeft, this.nearfarBottomRight, this.farBottomRight,
    this.nearfarBottomLeft, this.farBottomRight, this.nearfarBottomRight,

    this.nearfarTopLeft, this.farTopLeft, this.farBottomLeft,
    this.nearfarTopLeft, this.farBottomLeft, this.farTopLeft,

    this.nearfarTopLeft, this.nearfarBottomLeft, this.farBottomLeft,
    this.nearfarTopLeft, this.farBottomLeft, this.nearfarBottomLeft,

    this.nearfarTopRight, this.farTopRight, this.farBottomRight,
    this.nearfarTopRight, this.farBottomRight, this.farTopRight,

    this.nearfarTopRight, this.nearfarBottomRight, this.farBottomRight,
    this.nearfarTopRight, this.farBottomRight, this.nearfarBottomRight,

    this.farTopLeft, this.farTopRight, this.farBottomLeft,
    this.farTopLeft, this.farBottomLeft, this.farTopRight,

    this.farTopRight, this.farBottomRight, this.farBottomLeft,
    this.farTopRight, this.farBottomLeft, this.farBottomRight,

    this.farTopLeft, this.farTopRight, this.farBottomLeft,
    this.farTopLeft, this.farBottomLeft, this.farTopRight,

    this.nearfarTopRight, this.nearfarBottomRight, this.nearfarBottomLeft,
    this.nearfarTopRight, this.nearfarBottomLeft, this.nearfarBottomRight,
  ]
end

export getBox

"""
TODO
"""
function getVertices(this::Frustum)
  [
    this.nearTopLeft, this.farTopLeft, this.farTopRight,
    this.nearTopLeft, this.farTopRight, this.farTopLeft,

    this.nearTopLeft, this.nearTopRight, this.farTopRight,
    this.nearTopLeft, this.farTopRight, this.nearTopRight,

    this.nearBottomLeft, this.farBottomLeft, this.farBottomRight,
    this.nearBottomLeft, this.farBottomRight, this.farBottomLeft,

    this.nearBottomLeft, this.nearBottomRight, this.farBottomRight,
    this.nearBottomLeft, this.farBottomRight, this.nearBottomRight,

    this.nearTopLeft, this.farTopLeft, this.farBottomLeft,
    this.nearTopLeft, this.farBottomLeft, this.farTopLeft,

    this.nearTopLeft, this.nearBottomLeft, this.farBottomLeft,
    this.nearTopLeft, this.farBottomLeft, this.nearBottomLeft,

    this.nearTopRight, this.farTopRight, this.farBottomRight,
    this.nearTopRight, this.farBottomRight, this.farTopRight,

    this.nearTopRight, this.nearBottomRight, this.farBottomRight,
    this.nearTopRight, this.farBottomRight, this.nearBottomRight,

    this.farTopLeft, this.farTopRight, this.farBottomLeft,
    this.farTopLeft, this.farBottomLeft, this.farTopRight,

    this.farTopRight, this.farBottomRight, this.farBottomLeft,
    this.farTopRight, this.farBottomLeft, this.farBottomRight,
  ]
end

export getVertices

"""
TODO
"""
function SetFrustum(this::Frustum, angle::Float32, ratio::Float32, nearD::Float32, farD::Float32)
  this.ratio = ratio
  this.angle = angle
  this.nearDistance = nearD
  this.farDistance = farD

  this.tang = Float32(tan(deg2rad(this.angle) * 0.5))
  this.nearHeight = this.nearDistance * this.tang
  this.nearWidth = this.nearHeight * this.ratio
  this.farHeight = this.farDistance  * this.tang
  this.farWidth = this.farHeight * this.ratio
end

export SetFrustum

"""
TODO
"""
function SetCamera(this::Frustum, pos::Vec3f, target::Vec3f, up::Vec3f; near=0, far=0)
  pos = -pos
  target = -target

	Z = pos - target
	Z = normalize(Z)

	X = cross(up, Z)
	X = normalize(X)

	Y = cross(Z, X)

	nc = pos - Z * (near==0 ? this.nearDistance : near)
	fc = pos - Z * (far==0 ? this.farDistance : far)

	this.nearTopLeft = nc + Y * this.nearHeight - X * this.nearWidth
	this.nearTopRight = nc + Y * this.nearHeight + X * this.nearWidth
	this.nearBottomLeft = nc - Y * this.nearHeight - X * this.nearWidth
	this.nearBottomRight = nc - Y * this.nearHeight + X * this.nearWidth

	this.farTopLeft = fc + Y * this.farHeight - X * this.farWidth
	this.farTopRight = fc + Y * this.farHeight + X * this.farWidth
	this.farBottomLeft = fc - Y * this.farHeight - X * this.farWidth
	this.farBottomRight = fc - Y * this.farHeight + X * this.farWidth

	this.nearfarTopLeft = nc + Y * this.farHeight - X * this.farWidth
	this.nearfarTopRight = nc + Y * this.farHeight + X * this.farWidth
	this.nearfarBottomLeft = nc - Y * this.farHeight - X * this.farWidth
	this.nearfarBottomRight = nc - Y * this.farHeight + X * this.farWidth

	this.planes[:FRUSTUM_TOP] = Plane3D(this.nearTopRight, this.nearTopLeft, this.farTopLeft)
	this.planes[:FRUSTUM_BOTTOM] = Plane3D(this.nearBottomLeft, this.nearBottomRight, this.farBottomRight)
	this.planes[:FRUSTUM_LEFT] = Plane3D(this.nearTopLeft, this.nearBottomLeft, this.farBottomLeft)
	this.planes[:FRUSTUM_RIGHT] = Plane3D(this.nearBottomRight, this.nearTopRight, this.farBottomRight)
	this.planes[:FRUSTUM_NEAR] = Plane3D(this.nearTopLeft, this.nearTopRight, this.nearBottomRight)
	this.planes[:FRUSTUM_FAR] = Plane3D(this.farTopRight, this.farTopLeft, this.farBottomLeft)

  this.box[:BOX_TOP] = Plane3D(this.nearfarTopRight, this.nearfarTopLeft, this.farTopLeft) #Plane3D_Rect(this.nearfarTopRight, this.farTopLeft)
  this.box[:BOX_BOTTOM] = Plane3D(this.nearfarBottomLeft, this.nearfarBottomRight, this.farBottomRight) #Plane3D_Rect(this.nearfarBottomLeft, this.farBottomRight)
  this.box[:BOX_LEFT] = Plane3D(this.nearfarTopLeft, this.nearfarBottomLeft, this.farBottomLeft) #Plane3D_Rect(this.nearfarTopLeft, this.farBottomLeft)
  this.box[:BOX_RIGHT] = Plane3D(this.nearfarBottomRight, this.nearfarTopRight, this.farBottomRight) #Plane3D_Rect(this.nearfarBottomRight, this.farBottomRight)
  this.box[:BOX_NEAR] = Plane3D(this.nearfarTopLeft, this.nearfarTopRight, this.nearfarBottomRight)
  this.box[:BOX_FAR] = this.planes[:FRUSTUM_FAR]

  this.pos[:CENTER] = Vec3f((pos + target) / 2)
  this.pos[:CAMERA] = pos
  this.pos[:TARGET] = target
  this.centerRadius = 0
end

export SetCamera

"""
TODO
"""
function checkPoint2(this::Frustum, pos::Vec3f)
  result = :FRUSTUM_INSIDE
  distances = Dict{Symbol,Float32}()

  for (k,plane) in this.planes
    if (distances[k] = GetPointDistance(plane, pos)) < 0 result = :FRUSTUM_OUTSIDE end
  end

  for (k,plane) in this.box
    distances[k] = GetPointDistance(plane, pos)
  end

  for (k,v) in this.pos
    origin = this.pos[k]
    distances[k] = euclidean(origin,pos)
    distances[Symbol(k,:_X)] = origin.x - pos.x
    distances[Symbol(k,:_Y)] = origin.y - pos.y
    distances[Symbol(k,:_Z)] = origin.z - pos.z
  end

  X = distances[:BOX_LEFT] #distances[:CAMERA_X]
  Y = distances[:CAMERA_Y]
  Z = distances[:CAMERA_Z]

  #X,Y = normalize([X,Y])

  distances[:X] = X
  distances[:Y] = Y
  distances[:Z] = Z

  (result, distances)
end

export checkPoint2

"""
TODO
"""
function checkSphere2(this::Frustum, pos::Vec3f, radius::Number)
  result = :FRUSTUM_INSIDE
  (_, distances) = checkPoint(this, pos)

	for (k,distance) in distances
    if !haskey(this.planes, k) continue end
		if distance < -radius result = :FRUSTUM_OUTSIDE
		elseif distance < radius && result != :FRUSTUM_OUTSIDE result = :FRUSTUM_INTERSECT
		end
	end

	(result, distances)
end

export checkSphere2

"""
TODO
"""
function checkPoint(this::Frustum, pos::Vec3f)
  result = :FRUSTUM_INSIDE
  distances = Dict{Symbol,Float32}()

  for (k,plane) in this.planes
    distance = distances[k] = GetPointDistance(plane, pos)
    if distance < 0 result = :FRUSTUM_OUTSIDE
    elseif distance == 0 result = :FRUSTUM_INTERSECT
    end
  end

  (result, distances)
end

export checkPoint

"""
TODO
"""
function checkSphere(this::Frustum, pos::Vec3f, radius::Number)
  result = :FRUSTUM_INSIDE
  distances = Dict{Symbol,Float32}()

  for (k,plane) in this.planes
    distance = distances[k] = GetPointDistance(plane, pos)
		if distance < -radius result = :FRUSTUM_OUTSIDE
		elseif distance <= radius && result != :FRUSTUM_OUTSIDE result = :FRUSTUM_INTERSECT
		end
	end

	(result, distances)
end

export checkSphere

"""
TODO
"""
function checkCube(this::Frustum, center::Vec3f, size::Vec3f)
  result = :FRUSTUM_INSIDE
	distance = 0f0
  distances = Dict{Symbol,Dict{Symbol,Float32}}()

  for (k,plane) in this.planes
    # Reset counters for corners in and out
    out = 0
    in = 0

    distance = Dict{Symbol,Float32}(:LDB=>0f0,:RDB=>0f0,:LDF=>0f0,:RDF=>0f0,:LUB=>0f0,:RUB=>0f0,:LUF=>0f0,:RUF=>0f0)

    if (distance[:LDB] = GetPointDistance(plane, center + Vec3f(-size.x, -size.y, -size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:RDB] = GetPointDistance(plane, center + Vec3f(size.x, -size.y, -size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:LDF] = GetPointDistance(plane, center + Vec3f(-size.x, -size.y, size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:RDF] = GetPointDistance(plane, center + Vec3f(size.x, -size.y, size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:LUB] = GetPointDistance(plane, center + Vec3f(-size.x, size.y, -size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:RUB]= GetPointDistance(plane, center + Vec3f(size.x, size.y, -size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:LUF] = GetPointDistance(plane, center + Vec3f(-size.x, size.y, size.z))) < 0 out+=1
    else in+=1
    end

    if (distance[:RUF] = GetPointDistance(plane, center + Vec3f(size.x, size.y, size.z))) < 0 out+=1
    else in+=1
    end

    distances[k] = distance

    # If all corners are out
    if in <= 0 result = :FRUSTUM_OUTSIDE
    # If some corners are out and others are in
    elseif out > 0 && result != :FRUSTUM_OUTSIDE result = :FRUSTUM_INTERSECT
    end

  end

  (result, distances)
end

export checkCube

end #FrustumManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
#-------- Inverses of DPR1

function inv(A::SymDPR1{T},i::Integer,τ::Vector{Float64}=[1e3,10.0*length(A.D)]) where T
    # COMPUTES: inverse of a shifted SymDPR1 matrix A=diagm(A.D)+A.r*A.u*A.u',
    # inv(A-A.D[i]*I) which is a SymArrow.
    # Uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # τ=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used
    n=length(A.D)
    D=Array{T}(undef,n-1)
    z=Array{T}(undef,n-1)
    wz=one(T)/A.u[i]
    σ=A.D[i]
    for k=1:i-1
        D[k]=one(T)/(A.D[k]-σ)
        z[k]=-A.u[k]*D[k]*wz
    end
    for k=i+1:n
        D[k-1]=one(T)/(A.D[k]-σ)
        z[k-1]=-A.u[k]*D[k-1]*wz
    end
    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    for k=1:i-1
        D[k]>0.0 ? P+=A.u[k]^2*D[k] : Q+=A.u[k]^2*D[k]
    end
    for k=i:n-1
        D[k]>0.0 ? P+=A.u[k+1]^2*D[k] : Q+=A.u[k+1]^2*D[k]
    end
    A.r>0 ? P=P+one(T)/A.r : Q=Q+one(T)/A.r
    Kb=(P-Q)/abs(P+Q)
    # compute Kz
    Kz=(sum(abs,A.u)-abs(A.u[i]))*abs(wz)
    # Kz=max(abs,A.z)*abs(wz)
    if Kb<τ[1] ||  Kz<τ[2]
        b=(P+Q)*wz*wz
    else  # recompute in Double or BigFoat
        if Kz<1.0/eps()^2
            Type=Double
            Qout=1
        else
        # Example of a matrix where BigFloat is neeed, courtesy of Stan Eisenstat, is:
        # A=SymDPR1([1+3*eps(), 1-3*eps(), 0, -(1-eps()), -(1+eps())],[1,1,eps()^4, 3,3],1/16)
            Type=BigFloat
            Qout=50
        end
        σ₁=map(Type,A.D[i])
        Dd=[[Type(A.D[k])-σ₁ for k=1:i-1];
            [Type(A.D[k])-σ₁ for k=i+1:length(A.D)]]
        wzd=Type(A.u[i])

        Pd,Qd=map(Type,(0.0,0.0))
        for k=1:i-1
            convert(Float64,Dd[k])>0.0 ? Pd+=Type(A.u[k])^2/Dd[k] :
            Qd+=Type(A.u[k])^2/Dd[k]
        end
        for k=i+1:n
            convert(Float64,Dd[k-1])>0.0 ? Pd+=Type(A.u[k])^2/Dd[k-1] :
            Qd+=Type(A.u[k])^2/Dd[k-1]
            # @show P,Q
        end
        A.r > 0 ?   Pd+=Type(1.0)/Type(A.r)  : Qd+=Type(1.0)/Type(A.r)
        bd=(Pd+Qd)/(wzd*wzd)
        b=convert(Float64,bd)
    end
    # return this
    SymArrow(D,z,b,i),Kb,Kz,Qout
end # inv

function inv!(B::SymArrow{T},A::SymDPR1{T},i::Integer,τ::Vector{Float64}=[1e3,10.0*length(A.D)]) where T
    # In-place version!
    # COMPUTES: inverse of a shifted SymDPR1 matrix A=diagm(A.D)+A.r*A.u*A.u',
    # inv(A-A.D[i]*I) which is a SymArrow.
    # Uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # τ=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used
    n=length(A.D)
    # D=Array{T}(undef,n-1)
    # z=Array{T}(undef,n-1)
    wz=one(T)/A.u[i]
    σ=A.D[i]
    for k=1:i-1
        B.D[k]=one(T)/(A.D[k]-σ)
        B.z[k]=-A.u[k]*B.D[k]*wz
    end
    for k=i+1:n
        B.D[k-1]=one(T)/(A.D[k]-σ)
        B.z[k-1]=-A.u[k]*B.D[k-1]*wz
    end

    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    for k=1:i-1
        B.D[k]>0.0 ? P+=A.u[k]^2*B.D[k] : Q+=A.u[k]^2*B.D[k]
    end
    for k=i:n-1
        B.D[k]>0.0 ? P+=A.u[k+1]^2*B.D[k] : Q+=A.u[k+1]^2*B.D[k]
    end
    A.r>0 ? P+=one(T)/A.r : Q+=one(T)/A.r
    Kb=(P-Q)/abs(P+Q)

    # compute Kz
    Kz=(sum(abs,A.u)-abs(A.u[i]))*abs(wz)

    # Kz=max(abs,A.z)*abs(wz)
    if Kb<τ[1] ||  Kz<τ[2]
        b=(P+Q)*wz*wz
        B.a=b
        B.i=i
    else  # recompute in Double or BigFoat
        if Kz<1.0/eps()^2
            Type=Double
            Qout=1
        else
        # Example of a matrix where BigFloat is neeed, courtesy of Stan Eisenstat, is:
        # A=SymDPR1([1+3*eps(), 1-3*eps(), 0, -(1-eps()), -(1+eps())],[1,1,eps()^4, 3,3],1/16)
            Type=BigFloat
            Qout=50
        end
        σ₁=map(Type,A.D[i])
        Dd=[[Type(A.D[k])-σ₁ for k=1:i-1];
            [Type(A.D[k])-σ₁ for k=i+1:length(A.D)]]
        wzd=Type(A.u[i])

        Pd,Qd=map(Type,(0.0,0.0))
        for k=1:i-1
            convert(Float64,Dd[k])>0.0 ? Pd+=Type(A.u[k])^2/Dd[k] :
            Qd+=Type(A.u[k])^2/Dd[k]
        end
        for k=i+1:n
            convert(Float64,Dd[k-1])>0.0 ? Pd+=Type(A.u[k])^2/Dd[k-1] :
            Qd+=Type(A.u[k])^2/Dd[k-1]
            # @show P,Q
        end
        A.r > 0 ?   Pd+=Type(1.0)/Type(A.r)  : Qd+=Type(1.0)/Type(A.r)
        bd=(Pd+Qd)/(wzd*wzd)
        b=convert(Float64,bd)
        B.a=b
        B.i=i
    end
    # SymArrow(D,z,b,i),Kb,Kz,Qout
    return Kb,Kz,Qout
end # inv!

function inv(A::SymDPR1{T}, σ::Float64, τ::Float64=1.0e3) where T
    # COMPUTES: inverse of the σed SymDPR1 A = diagm(A.D)+A.r*A.u*A.u',
    # inv(A-σ*I) = D + ρ*u*u', σ!=A.D[i], which is again a SymDPR1
    # uses DoubleDouble to compute A.r accurately, if needed.
    # τ is tolerance, usually 1e3,  0.0 forces Double, 1e50 would never use it
    # RETURNS: SymDPR1(D,u,ρ), Kρ, Qout
    # Kρ - condition Kρ, Qout = 1 / 0 - Double was / was not used
    n=length(A.D)
    D=Array{T}(undef,n)
    u=Array{T}(undef,n)
    for k=1:n
        D[k]=one(T)/(A.D[k]-σ)
        u[k]=A.u[k]*D[k]
    end
    # compute gamma and Kgamma
    #--- compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0
    for k=1:n
        D[k]>0.0 ? P+=A.u[k]^2*D[k] : Q+=A.u[k]^2*D[k]
    end
    A.r>0 ? P+=one(T)/A.r : Q+=one(T)/A.r
    # Condition of ρ
    Kρ=(P-Q)/abs(P+Q)
    if Kρ < τ
        ρ=-one(T)/(P+Q)
    else  # recompute in Double
        Qout=1
        Pd,Qd=map(Double,(0.0,0.0))
        σ₁=Double(σ)
        for k=1:n
            D[k]>0.0 ? Pd+=Double(A.u[k])^2/(Double(A.D[k])-σ₁) : Qd+=Double(A.u[k])^2/(Double(A.D[k])-σ₁)
        end
        A.r > 0 ?   Pd+=Double(1.0)/Double(A.r)  : Qd+=Double(1.0)/Double(A.r)
        Kρ=Float64((Pd-Qd)/abs(Pd+Qd))
        if Kρ<τ/eps()
            r=Double(1.0)/(Pd+Qd)
            ρ=-(r.hi+r.lo)
        else
        # Here we need quadruple working precision. We are using BigFloat.
        # Example of a matrix where this is neeed, courtesy of Stan Eisenstat, is:
        # A=SymDPR1([1+3*eps(), 1-3*eps(), -(1-eps()), -(1+eps())],[1,1,3,3.0],1/16)
            Qout=100
            Pd,Qd=map(BigFloat,(0.0,0.0))
            σ₁=BigFloat(σ)
            for k=1:n
                D[k]>0.0 ? Pd+=BigFloat(A.u[k])^2/(BigFloat(A.D[k])-σ₁) : Qd+=BigFloat(A.u[k])^2/(BigFloat(A.D[k])-σ₁)
            end
            A.r > 0 ?   Pd+=BigFloat(1.0)/BigFloat(A.r)  : Qd+=BigFloat(1.0)/BigFloat(A.r)
            Kρ=Float64((Pd-Qd)/abs(Pd+Qd))
            r=BigFloat(1.0)/(Pd+Qd)
            ρ=-Float64(r)
        end
    end
    # returns the following
    SymDPR1(D,u,ρ), Kρ, Qout
end # inv

function inv(A::SymDPR1{T}, σ::Double) where T
    # COMPUTES: inverse of the σed SymDPR1 A = diagm(A.D)+A.r*A.u*A.u',
    # inv(A-σ*I) = D + ρ*u*u', σ!=A.D[i], which is again a SymDPR1
    # here σ is Double so it uses Double to compute everything
    # RETURNS: SymDPR1(D,u,ρ), Qout
    # Qout = 1 on exit meaning Double was used
    n=length(A.D)
    D=Array(Double,n)
    u=Array(Double,n)
    for k=1:n
        D[k]=Double(A.D[k])-σ
    end
    u=map(Double,A.u)
    oned=Double(1.0)
    zerod=Double(0.0)
    for k=1:n
        u[k]=u[k]/D[k]
        D[k]=oned/D[k]
    end
    # compute ρ
    # compute the sum in a plain loop
    P,Q=zerod,zerod
    Qout=1
    for k=1:n
        D[k].hi > 0.0 ? P=P+Double(A.u[k])*u[k] : Q=Q+Double(A.u[k])*u[k]
    end
    A.r > 0 ?   P+=Double(1.0)/Double(A.r)  : Q+=Double(1.0)/Double(A.r)
    r=oned/(P+Q)
    ρ=-(r.hi+r.lo)
    # returns the following
    # SymDPR1(T[x.hi+x.lo for x=D],T[x.hi+x.lo for x=u],ρ), Qout
    D1=Array{T}(n)
    u1=Array{T}(n)
    for k=1:n
        D1[k]=D[k].hi+D[k].lo
        u1[k]=u[k].hi+u[k].lo
    end
    SymDPR1(D1,u1,ρ), Qout
end # inv

function  eigen( A::SymDPR1{T}, Ainv::SymArrow{T}, k::Integer,
    τ::Vector{Float64}=[1e3,10.0*length(A.D),1e3,1e3,1e3]) where T
    # COMPUTES: k-th eigenpair of an ordered irreducible SymDPR1
    # A = Diagonal(A.D)+A.r*A.u*A.u', A.r > 0
    # τ=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: λ, v, Info
    # where
    # λ - k-th eigenvalue in descending order
    # v - λ's normalized eigenvector
    # Info = Sind, Kb, Kz, Kν, Kρ, Qout
    # Sind = shift index i for the k-th eigenvalue
    # Kb, Kz, Kν, Kρ - condition numbers
    # Qout = 1 / 0 - Double was / was not used
    # Set the dimension
    n = length(A.D)
    # Set all conditions initially to zero
    Kb,Kz,Kν,Kρ=0.0,0.0,0.0,0.0
    Qout=0
    v=zeros(T,n)
    # Kz is former kappa_nu
    # Determine the shift σ, the shift index i, and whether λ
    # is on the left or the right side of the nearest pole
    # Exterior eigenvalues (k = 1 or k = n):
    if k == 1
        σ,i,side = A.D[1],1,'R'
    else
        # Interior eigenvalues (k in (2,...n-1) ):
        Dtemp = A.D.-A.D[k]
        middle = Dtemp[k-1]/2.0
        Fmiddle = 1.0+A.r*sum(A.u.^2 ./(Dtemp.-middle))
        σ,i,side = Fmiddle > 0.0 ? (A.D[k],k,'R') : (A.D[k-1],k-1,'L')
    end

    # Compute the inverse of the shifted matrix, A_i^(-1), Kb and Kz
    # Ainv,Kb,Kz,Qout = inv(A,i)
    Kb,Kz,Qout = inv!(Ainv,A,i)
    # Compute the eigenvalue of the inverse σed matrix
    ν = bisect( Ainv,side )

    #  ν=fastzero([invD1; 0; invD2], [w1;wz;w2], b, side); # implement later
    #  [ν-ν₁]/ν, [ν-νeva]/ν, pause  # and compare
    if abs(ν)==Inf
        # this is nonstandard
        # Deflation in dpr1eig (ν=Inf)
        v[i]=1.0
        λ=sigma
    else
        # standard case, full computation
        # ν₁ is the F- or 1-norm of the inverse of the σed matrix
        # ν₁0=maximum([sum(abs,Ainv.z)+abs(Ainv.a), maximum(abs.(Ainv.D)+abs.(Ainv.z))])
        ν₁=0.0
        for k=1:n-1
            ν₁=max(ν₁,abs(Ainv.D[k])+abs(Ainv.z[k]))
        end
        ν₁=max(sum(abs,Ainv.z)+abs(Ainv.a), ν₁)
        Kν=ν₁/abs(ν)
        while Kν>τ[3]
            # Remedies according to Remark 3 - we σ between original
            # eigenvalues and compute DPR1 matrix
            # 1/ν₁+sigma, 1/ν+sigma
            # println("Remedy 3 ")
            ν = side=='R' ? abs(ν) : -abs(ν)
            ν₁=-sign(ν)*ν₁
            σ₁=(ν₁+ν)/(2.0*ν*ν₁)+σ
            if count(isequal(σ₁),A.D)>0 # we came back with a pole
                # recompute sigmav more accurately according with Dekker
                σᵥ=(Double(ν₁)+Double(ν))/(Double(2.0)*Double(ν)*Double(ν₁))+Double(σ)
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv,Qout₁=inv(A,σᵥ) # Ainv is Float64
                ν₁=bisect(Ainv,side)
                μ₁ = 1.0/ν₁
                D0=map(A.D,Double)-σᵥ
                D1=D0.hi+D0.lo
                if count(isequal(μ₁),D1)>0
                    ind=findall(isequal(μ₁),D1)
                    v=zeros(T,n)
                    v[ind]=1.0
                else
                    v.=A.u./(D1.-μ₁)
                end
                # σ the eigenvalue back in Double
                lam = Double(1.0)/Double(ν₁)+σᵥ

                # Return this
                λ = lam.hi+lam.lo
                normalize!(v)
                κ=Info(i,Kb,Kz,Kν,Kρ,Qout)
                return λ,v,κ
            else
                # Compute the inverse of the Shifted DPR1 (DPR1)
                Ainv₁, Kρ,Qout₁=inv(A,σ₁,τ[4]) # Ainv is Float64
                # Compute the eigenvalue by bisect for DPR1
	            # Note: instead of bisection could use dlaed4 (need a wrapper) but
	            # it is not faster. There norm(u)==1
                ν= bisect(Ainv₁,side)
                # Check the code below for DPR1
                # if side=='R' && A.D[1]>0.0 && ν<0.0
                #     ν=bisect(Ainv,'L')
                # end
                μ=1.0/ν
                ν₁=maximum(abs,Ainv₁.D)+abs(Ainv₁.r)*dot(Ainv₁.u,Ainv₁.u)
                Kν=ν₁/abs(ν)
                σ=σ₁
                # standard v
                # v=A.u./((A.D-σ₁)-μ₁)
                # Return this - σ the eigenvalue back and normalize the vector
                # λ, v = μ₁+σ₁, v/norm(v)
            end
            Qout=Qout+2*Qout₁
        end
        # Accuracy is fine, compute the eigenvector
        μ = 1.0/ν
        # v=[ A.u./((A.D-σ)-μ)]
        for l=1:n
            v[l]= A.u[l]/((A.D[l]-σ)-μ)
        end
        λ = μ + σ
        normalize!(v)

        # Remedy according to Remark 1 - we recompute the the eigenvalue
        # near zero from the inverse of the original matrix (a DPR1 matrix).
        if (abs(A.D[i])+abs(1.0/ν))/abs(λ)>τ[5]
            if (k==1 && A.D[1]<0.0 || side=='L' && sign(A.D[i])+sign(A.D[i+1])==0 || i>1 &&
                side=='R' && sign(A.D[i])+sign(A.D[i-1])==0)
                # println("Remedy 1 ")
                # Compute the inverse of the original arrowhead (DPR1)
                Ainv₁,Kρ,Qout₁ = inv(A,0.0,τ[4]) # Ainv is Float64
                Qout=Qout+4*Qout₁
                if abs(Ainv₁.r)==Inf
                    λ=0.0
                else
                    # Here we do not need bisection. We compute the Rayleigh
                    # quotient by using already computed vectors which is
                    # componentwise accurate
                    ν₁=sum(v.^2 .*Ainv₁.D)+Ainv₁.r*sum(v.*Ainv₁.u)^2;
                    λ=1.0/ν₁
                end
            end
        end
    end
    # Return this
    κ=Info(i,Kb,Kz,Kν,Kρ,Qout)
    λ,v,κ
end # eigen (k)

function eigen(A::SymDPR1{T}, τ::Vector{Float64}=[1e3,10.0*length(A.D),1e3,1e3,1e3]) where T
    # COMPUTES: all eigenvalues and eigenvectors of a real symmetric SymDPR1
    # A = Diagonal(A.D)+A.r*A.u*A.u'
    # τ = [tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3] or similar
    # RETURNS: Eigen(Λ,U), κ
    # where
    # Λ = eigenvalues in decreasing order, U = eigenvectors,
    # κ[k]=Info(Sind[k], Kb[k], Kz[k], Kν[k], Kρ[k], Qout[k]
    # Sind[k] - shift index i for the k-th eigenvalue
    # Kb, Kz, Kν, Kρ [k] - respective conditions for the k-th eigenvalue
    # Qout[k] = 1 / 0 - Double was / was not used when computing k-th eigenvalue

    n=length(A.D)
    n0=n
    # Checking if A.r > 0
    signr =  A.r > 0.0 ? 1.0 : -1.0
    # Ordering the matrix
    D=signr*A.D
    is=sortperm(D,rev=true)
    D=D[is]
    z=A.u[is]
    ρ=signr*A.r
    # Eigenvecgtor matrix and eigenvalue vector
    U=Array{T}(I,n,n)
    Λ=zeros(T,n)
    κ=[Info(0, 0.0, 0.0, 0.0, 0.0, 0) for i=1:n]

    # Quick return for 1x1, this is trivial for SymArrow, not so trivial here :)
    if n==1
        if (D==0)&&((ρ==0)|(z==0))
            Λ=0
        else
            Λ=A.D[1]+A.r*A.u[1]^2
            # Higher accuracy if needed
            KD=(abs(A.D[1])+abs(A.r)*A.u[1]^2)/abs(Λ)
            if KD>τ[1]
                Λd=Double(A.D[1])+Double(A.r)*Double(A.u[1])^2
                Λ=Λd.hi+Λd.lo
            end
            # Qout=1
            κ[1]=Info(0, 0.0, 0.0, 0.0, 0.0, 1)
        end
        return Eigen([Λ], U), κ
    end

    #  test for deflation in z
    z0=findall(iszero,z)
    zx=findall(!iszero,z)
    if isempty(zx)  # nothing to do
        Λ=A.D
        isΛ=sortperm(Λ,rev=true)
        Λ=Λ[isΛ]
        U=view(U,:,isΛ) # [:,isΛ]
        return U,Λ,Sind,Kb,Kz,Kν,Kρ,Qout
    end
    if !isempty(z0)
        Λ[z0]=D[z0]
        D=D[zx]
        z=z[zx]
        if !isempty(z)
            # return U,Λ,Sind,Kb,Kz,Kν,Kρ,Qout
            n=length(z)
        end
    end

    #  Test for deflation in D
    g=D[1:n-1]-D[2:n]
    # Can play with inexact deflation
    # g0=find(abs(g)<eps)
    # gx=find(abs(g)>=eps)
    # Only exact deflation !!
    g0=findall(iszero,g)
    gx=findall(!iszero,g)
    if !isempty(g0)
        # Deflation
        Dgx=D[gx]; zgx=z[gx]
        lg0=length(g0)
        R=Array(Tuple{LinearAlgebra.Givens{Float64},Float64},lg0)
        for l=lg0:-1:1
            R[l]=givens(z[g0[l]],z[g0[l]+1],zx[g0[l]],zx[g0[l]+1])
            z[g0[l]]=R[l][2]; z[g0[l]+1]=0.0
            # A_mul_Bc!(U,R) # multiply R'*U later
            Λ[zx[g0[l]+1]]=D[g0[l]+1]
        end
        # remains
        gx=[0;gx].+1
        nn=length(gx)
        zxx=zx[gx]
        Axx=SymDPR1(D[gx],z[gx],ρ)
        # Bxx=SymArrow(Vector{T}(undef,nn-1),Vector{T}(undef,nn-1),one(T),nn)
        Bxx=[SymArrow(Vector{T}(undef,nn-1),Vector{T}(undef,nn-1),one(T),nn) for i=1:Threads.nthreads()]
        Threads.@threads for k=1:nn
            tid=Threads.threadid()
            Λ[zxx[k]],U[zxx,zxx[k]],κ[zxx[k]]=eigen(Axx,Bxx[tid],k)
        end
        for l=1:lg0
            # U=R[l][1]'*U
            lmul!(R[l][1],U)
        end
    else
        # No deflation in D
        Ax=SymDPR1(D,z,ρ)
        Bx=[SymArrow(Vector{T}(undef,n-1),Vector{T}(undef,n-1),one(T),n) for i=1:Threads.nthreads()]
        # Bx=SymArrow(Vector{T}(undef,n-1),Vector{T}(undef,n-1),one(T),n)
        Threads.@threads for k=1:n
            tid=Threads.threadid()
            Λ[zx[k]],U[zx,zx[k]],κ[zx[k]]=eigen(Ax,Bx[tid],k)
        end
    end

    # back premutation of vectors
    isi=sortperm(is)
    # change the sign if A.r was negative
    # must sort Λ once more
    Λ=signr*Λ
    es=sortperm(Λ,rev=true)
    Λ.=Λ[es]
    # U=U[isi,es]
    # Return this
    U=view(U,isi,es)
    κ.=κ[es]
    # Return this
    return Eigen(Λ,U),κ
end # eigen (all)


text:
## spectral.jl : Gaussian random field generator based on the spectral (eigenvalue) decomposition of the covariance matrix

"""
    Spectral <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) based on a spectral (eigenvalue) decomposition of the covariance matrix.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,Spectral(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a Spectral decomposition

julia> plot(grf)
[...]

```
This is also useful when computing Gaussian random fields on a Finite Element mesh using a truncated KL expansion. Here's an example that computes the first 10 eigenfunctions on an L-shaped domain.

```jldoctest
julia> p, t = Lshape();

julia> grf = grf = GaussianRandomField(CovarianceFunction(2,Matern(0.2,1.0)),Spectral(),p,t,n=10)
Gaussian random field with 2d Matérn covariance function (λ=0.2, ν=1.0, σ=1.0, p=2.0) on a mesh with 998 points and 1861 elements, using a spectral decomposition

julia> tricontourf(p[:,1],p[:,2],grf.data.eigenfunc[:,1],triangles=t-1,cmap=get_cmap("viridis"))
[...]

```
See also: [`Cholesky`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref)
"""
struct Spectral <: GaussianRandomFieldGenerator end

# container type for eigenvalues and eigenfunctions
struct SpectralData{X,Y}
    eigenval::X
    eigenfunc::Y
end

function _GaussianRandomField(mean, cov::CovarianceFunction, method::Spectral, pts...;
                              n::Integer=0, eigensolver=EigenSolver())
    C = apply(cov,pts,pts)

    # compute eigenvalue decomposition
    if n == 0
        F = eigen(C)
        idx = sortperm(F.values,rev=true)
        Λ = F.values[idx]
        U = F.vectors[:,idx]
    else # thanks to #4
        eigenval, eigenfunc = compute(C, n, eigensolver)
        idx = sortperm(eigenval,rev=true)[1:n]
        Λ = eigenval[idx]
        U = eigenfunc[:,idx]
    end

    # if negative eigenvalues detected, remove them
    m = findfirst(x -> x < 0, Λ)
    if m != nothing
        m -= 1
        @warn begin
            "$(length(Λ) - m) negative eigenvalues ≥ $(Λ[end]) detected, Gaussian " *
            "random field will be approximated (ignoring all negative eigenvalues)"
        end

        resize!(Λ, m)
        U = U[:, 1:m]
    end

    # store eigenvalues and eigenfunctions
    Λ .= sqrt.(Λ) # note: store sqrt of eigenval for more efficient sampling
    data = SpectralData(Λ, U)

    GaussianRandomField{Spectral,typeof(cov),typeof(pts),typeof(mean),typeof(data)}(mean, cov, pts, data)
end

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{Spectral}) = length(grf.data.eigenval)

# see KarhunenLoeve.jl for sample function

Base.show(io::IO, ::Spectral) = print(io, "spectral decomposition")


text:
@testset "layer" begin
    T = Float32
    in_channel = 3
    in_channel_edge = 1
    out_channel = 5

    N = 4
    E = 4
    adj = T[0. 1. 0. 1.;
            1. 0. 1. 0.;
            0. 1. 0. 1.;
            1. 0. 1. 0.]
    fg = FeaturedGraph(adj)

    @testset "GCNConv" begin
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            gc = GCNConv(fg, in_channel=>out_channel)
            @test size(gc.weight) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)
            @test GraphSignals.adjacency_matrix(gc.fg) == adj

            Y = gc(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = gc(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(x -> sum(gc(x)), X)[1]
            @test size(g) == size(X)

            g = Zygote.gradient(model -> sum(model(X)), gc)[1]
            @test size(g.weight) == size(gc.weight)
            @test size(g.bias) == size(gc.bias)
        end

        @testset "layer without graph" begin
            gc = GCNConv(in_channel=>out_channel)
            @test size(gc.weight) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)
            @test !has_graph(gc.fg)

            fg = FeaturedGraph(adj, nf=X)
            fg_ = gc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws MethodError gc(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = gc(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(x -> sum(node_feature(gc(x))), fg)[1]
            @test size(g.nf) == size(X)

            g = Zygote.gradient(model -> sum(node_feature(model(fg))), gc)[1]
            @test size(g.weight) == size(gc.weight)
            @test size(g.bias) == size(gc.bias)
        end

        @testset "bias=false" begin
            @test length(Flux.params(GCNConv(2=>3))) == 2
            @test length(Flux.params(GCNConv(2=>3, bias=false))) == 1
        end
    end

    @testset "ChebConv" begin
        k = 6
        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))
        @testset "layer with graph" begin
            cc = ChebConv(fg, in_channel=>out_channel, k)
            @test size(cc.weight) == (out_channel, in_channel, k)
            @test size(cc.bias) == (out_channel,)
            @test GraphSignals.adjacency_matrix(cc.fg) == adj
            @test cc.k == k

            Y = cc(X)
            @test size(Y) == (out_channel, N)

            # Test with transposed features
            Y = cc(Xt)
            @test size(Y) == (out_channel, N)

            g = Zygote.gradient(x -> sum(cc(x)), X)[1]
            @test size(g) == size(X)

            g = Zygote.gradient(model -> sum(model(X)), cc)[1]
            @test size(g.weight) == size(cc.weight)
            @test size(g.bias) == size(cc.bias)
        end

        @testset "layer without graph" begin
            cc = ChebConv(in_channel=>out_channel, k)
            @test size(cc.weight) == (out_channel, in_channel, k)
            @test size(cc.bias) == (out_channel,)
            @test !has_graph(cc.fg)
            @test cc.k == k

            fg = FeaturedGraph(adj, nf=X)
            fg_ = cc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws MethodError cc(X)

            # Test with transposed features
            fgt = FeaturedGraph(adj, nf=Xt)
            fgt_ = cc(fgt)
            @test size(node_feature(fgt_)) == (out_channel, N)

            g = Zygote.gradient(x -> sum(node_feature(cc(x))), fg)[1]
            @test size(g.nf) == size(X)

            g = Zygote.gradient(model -> sum(node_feature(model(fg))), cc)[1]
            @test size(g.weight) == size(cc.weight)
            @test size(g.bias) == size(cc.bias)
        end

        @testset "bias=false" begin
            @test length(Flux.params(ChebConv(2=>3, 3))) == 2
            @test length(Flux.params(ChebConv(2=>3, 3, bias=false))) == 1
        end
    end

    # @testset "GraphConv" begin
    #     X = rand(T, in_channel, N)
    #     Xt = transpose(rand(T, N, in_channel))
    #     @testset "layer with graph" begin
    #         gc = GraphConv(fg, in_channel=>out_channel)
    #         @test adjacency_list(gc.fg) == [[2,4], [1,3], [2,4], [1,3]]
    #         @test size(gc.weight1) == (out_channel, in_channel)
    #         @test size(gc.weight2) == (out_channel, in_channel)
    #         @test size(gc.bias) == (out_channel,)

    #         Y = gc(X)
    #         @test size(Y) == (out_channel, N)

    #         # Test with transposed features
    #         Y = gc(Xt)
    #         @test size(Y) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(gc(x)), X)[1]
    #         @test size(g) == size(X)

    #         g = Zygote.gradient(model -> sum(model(X)), gc)[1]
    #         @test size(g.weight1) == size(gc.weight1)
    #         @test size(g.weight2) == size(gc.weight2)
    #         @test size(g.bias) == size(gc.bias)
    #     end

    #     @testset "layer without graph" begin
    #         gc = GraphConv(in_channel=>out_channel)
    #         @test size(gc.weight1) == (out_channel, in_channel)
    #         @test size(gc.weight2) == (out_channel, in_channel)
    #         @test size(gc.bias) == (out_channel,)

    #         fg = FeaturedGraph(adj, nf=X)
    #         fg_ = gc(fg)
    #         @test size(node_feature(fg_)) == (out_channel, N)
    #         @test_throws MethodError gc(X)

    #         # Test with transposed features
    #         fgt = FeaturedGraph(adj, nf=Xt)
    #         fgt_ = gc(fgt)
    #         @test size(node_feature(fgt_)) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(node_feature(gc(x))), fg)[1]
    #         @test size(g.nf) == size(X)

    #         g = Zygote.gradient(model -> sum(node_feature(model(fg))), gc)[1]
    #         @test size(g.weight1) == size(gc.weight1)
    #         @test size(g.weight2) == size(gc.weight2)
    #         @test size(g.bias) == size(gc.bias)
    #     end

    #     @testset "bias=false" begin
    #         @test length(Flux.params(GraphConv(2=>3))) == 3
    #         @test length(Flux.params(GraphConv(2=>3, bias=false))) == 2
    #     end
    # end

    @testset "GATConv" begin
        adj1 = [1 1 0 1;
                1 1 1 0;
                0 1 1 1;
                1 0 1 1]
        fg1 = FeaturedGraph(adj1)

        # isolated_vertex
        adj2 = [1 0 0 1;
                0 1 0 0;
                0 0 1 1;
                1 0 1 1]
        fg2 = FeaturedGraph(adj2)

        X = rand(T, in_channel, N)
        Xt = transpose(rand(T, N, in_channel))

        @testset "layer with graph" begin
            for heads = [1, 2], concat = [true, false], adj_gat in [adj1, adj2]
                fg_gat = FeaturedGraph(adj_gat)
                gat = GATConv(fg_gat, in_channel=>out_channel, heads=heads, concat=concat)

                @test size(gat.weight) == (out_channel * heads, in_channel)
                @test size(gat.bias) == (out_channel * heads,)
                @test size(gat.a) == (2*out_channel, heads)

                Y = gat(X)
                @test size(Y) == (concat ? (out_channel*heads, N) : (out_channel, N))

                # Test with transposed features
                Y = gat(Xt)
                @test size(Y) == (concat ? (out_channel*heads, N) : (out_channel, N))

                g = Zygote.gradient(x -> sum(gat(x)), X)[1]
                @test size(g) == size(X)

                g = Zygote.gradient(model -> sum(model(X)), gat)[1]
                @test size(g.weight) == size(gat.weight)
                @test size(g.bias) == size(gat.bias)
                @test size(g.a) == size(gat.a)
            end
        end

        @testset "layer without graph" begin
            for heads = [1, 2], concat = [true, false], adj_gat in [adj1, adj2]
                fg_gat = FeaturedGraph(adj_gat, nf=X)
                gat = GATConv(in_channel=>out_channel, heads=heads, concat=concat)
                @test size(gat.weight) == (out_channel * heads, in_channel)
                @test size(gat.bias) == (out_channel * heads,)
                @test size(gat.a) == (2*out_channel, heads)

                fg_ = gat(fg_gat)
                Y = node_feature(fg_)
                @test size(Y) == (concat ? (out_channel*heads, N) : (out_channel, N))
                @test_throws MethodError gat(X)

                # Test with transposed features
                fgt = FeaturedGraph(adj_gat, nf=Xt)
                fgt_ = gat(fgt)
                @test size(node_feature(fgt_)) == (concat ? (out_channel*heads, N) : (out_channel, N))

                g = Zygote.gradient(x -> sum(node_feature(gat(x))), fg_gat)[1]
                @test size(g.nf) == size(X)

                g = Zygote.gradient(model -> sum(node_feature(model(fg_gat))), gat)[1]
                @test size(g.weight) == size(gat.weight)
                @test size(g.bias) == size(gat.bias)
                @test size(g.a) == size(gat.a)
            end
        end

        @testset "bias=false" begin
            @test length(Flux.params(GATConv(2=>3))) == 3
            @test length(Flux.params(GATConv(2=>3, bias=false))) == 2
        end
    end

    # @testset "GatedGraphConv" begin
    #     num_layers = 3
    #     X = rand(T, in_channel, N)
    #     Xt = transpose(rand(T, N, in_channel))
    #     @testset "layer with graph" begin
    #         ggc = GatedGraphConv(fg, out_channel, num_layers)
    #         @test adjacency_list(ggc.fg) == [[2,4], [1,3], [2,4], [1,3]]
    #         @test size(ggc.weight) == (out_channel, out_channel, num_layers)

    #         Y = ggc(X)
    #         @test size(Y) == (out_channel, N)

    #         # Test with transposed features
    #         Y = ggc(Xt)
    #         @test size(Y) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(ggc(x)), X)[1]
    #         @test size(g) == size(X)

    #         g = Zygote.gradient(model -> sum(model(X)), ggc)[1]
    #         @test size(g.weight) == size(ggc.weight)
    #     end

    #     @testset "layer without graph" begin
    #         ggc = GatedGraphConv(out_channel, num_layers)
    #         @test size(ggc.weight) == (out_channel, out_channel, num_layers)

    #         fg = FeaturedGraph(adj, nf=X)
    #         fg_ = ggc(fg)
    #         @test size(node_feature(fg_)) == (out_channel, N)
    #         @test_throws MethodError ggc(X)

    #         # Test with transposed features
    #         fgt = FeaturedGraph(adj, nf=Xt)
    #         fgt_ = ggc(fgt)
    #         @test size(node_feature(fgt_)) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(node_feature(ggc(x))), fg)[1]
    #         @test size(g.nf) == size(X)

    #         g = Zygote.gradient(model -> sum(node_feature(model(fg))), ggc)[1]
    #         @test size(g.weight) == size(ggc.weight)
    #     end
    # end

    # @testset "EdgeConv" begin
    #     X = rand(T, in_channel, N)
    #     Xt = transpose(rand(T, N, in_channel))
    #     @testset "layer with graph" begin
    #         ec = EdgeConv(fg, Dense(2*in_channel, out_channel))
    #         @test adjacency_list(ec.fg) == [[2,4], [1,3], [2,4], [1,3]]

    #         Y = ec(X)
    #         @test size(Y) == (out_channel, N)

    #         # Test with transposed features
    #         Y = ec(Xt)
    #         @test size(Y) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(ec(x)), X)[1]
    #         @test size(g) == size(X)

    #         g = Zygote.gradient(model -> sum(model(X)), ec)[1]
    #         @test size(g.nn.weight) == size(ec.nn.weight)
    #         @test size(g.nn.bias) == size(ec.nn.bias)
    #     end

    #     @testset "layer without graph" begin
    #         ec = EdgeConv(Dense(2*in_channel, out_channel))

    #         fg = FeaturedGraph(adj, nf=X)
    #         fg_ = ec(fg)
    #         @test size(node_feature(fg_)) == (out_channel, N)
    #         @test_throws MethodError ec(X)

    #         # Test with transposed features
    #         fgt = FeaturedGraph(adj, nf=Xt)
    #         fgt_ = ec(fgt)
    #         @test size(node_feature(fgt_)) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(node_feature(ec(x))), fg)[1]
    #         @test size(g.nf) == size(X)

    #         g = Zygote.gradient(model -> sum(node_feature(model(fg))), ec)[1]
    #         @test size(g.nn.weight) == size(ec.nn.weight)
    #         @test size(g.nn.bias) == size(ec.nn.bias)
    #     end
    # end

    # @testset "GINConv" begin
    #     X = rand(Float32, in_channel, N)
    #     Xt = transpose(rand(Float32, N, in_channel))
    #     nn = Flux.Chain(Dense(in_channel, out_channel))
    #     eps = 0.001

    #     @testset "layer with graph" begin
    #         gc = GINConv(FeaturedGraph(adj), nn, eps)
    #         @test size(gc.nn.layers[1].weight) == (out_channel, in_channel)
    #         @test size(gc.nn.layers[1].bias) == (out_channel, )
    #         @test GraphSignals.adjacency_matrix(gc.fg) == adj

    #         Y = gc(FeaturedGraph(adj, nf=X))
    #         @test size(node_feature(Y)) == (out_channel, N)

    #         # Test with transposed features
    #         Y = gc(FeaturedGraph(adj, nf=Xt))
    #         @test size(node_feature(Y)) == (out_channel, N)

    #         g = Zygote.gradient(x -> sum(node_feature(gc(x))),
    #                             FeaturedGraph(adj, nf=X))[1]
    #         @test size(g.nf) == size(X)

    #         g = Zygote.gradient(model -> sum(node_feature(model(FeaturedGraph(adj, nf=X)))),
    #                             gc)[1]
    #         @test size(g.nn.layers[1].weight) == size(gc.nn.layers[1].weight)
    #         @test size(g.nn.layers[1].bias) == size(gc.nn.layers[1].bias)
    #         @test !in(:eps, Flux.trainable(gc))
    #     end
    # end

    # @testset "CGConv" begin
    #     fg = FeaturedGraph(adj)
    #     X = rand(Float32, in_channel, N)
    #     E = rand(Float32, in_channel_edge, ne(fg))
    #     Xt = transpose(rand(Float32, N, in_channel))
    #     @testset "layer with graph" begin
    #         cgc = CGConv(FeaturedGraph(adj),
    #                      (in_channel, in_channel_edge))
    #         @test size(cgc.Wf) == (in_channel, 2 * in_channel + in_channel_edge)
    #         @test size(cgc.Ws) == (in_channel, 2 * in_channel + in_channel_edge)
    #         @test size(cgc.bf) == (in_channel,)
    #         @test size(cgc.bs) == (in_channel,)

    #         Y = cgc(X, E)
    #         @test size(Y) == (in_channel, N)

    #         Yg = cgc(FeaturedGraph(adj, nf=X, ef=E))
    #         @test size(node_feature(Yg)) == (in_channel, N)
    #         @test edge_feature(Yg) == E
    #     end
    # end
end


text:
using JuMP, FactCheck
using Compat

!isdefined(:conic_solvers_with_duals) && include("solvers.jl")

const TOL = 1e-4

facts("[socduals] SOC dual") do
for solver in conic_solvers_with_duals
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)

    @variable(m, x[1:5])
    @objective(m, Max, x[1] + x[2] + x[3] + x[4] + 2*x[5])
    @constraint(m, constrcon1, norm(x[2:5]) <= x[1])
    @constraint(m, constrlin1, x[1] <= 5)
    @constraint(m, constrlin2, x[1]+2*x[2] - x[3] <= 10)
    @constraint(m, constrcon2, norm([2 3;1 1]*x[2:3]-[3;4]) <= x[5] - 2)

    solve(m)
    @fact length(getdual(constrcon1)) --> 5
    @fact length(getdual(constrcon2)) --> 3
    @fact length(m.conicconstrDuals) --> 10
    @fact length(m.linconstrDuals) --> 2

    @fact dot(getdual(constrcon2),[getvalue(x[5]) - 2;[2 3;1 1]*[getvalue(x[2]);getvalue(x[3])]-[3;4]]) --> less_than_or_equal(TOL)

end
end
end

# Test for consistency between LP duals and SOC duals

facts("[socduals] LP dual vs SOC dual (Max)") do
for lp_solver in lp_solvers
context("With LP solver $(typeof(lp_solver))") do

    m1 = Model(solver=lp_solver)
    @variable(m1, x1[1:2] >= 0)
    @objective(m1, Max, x1[1] + 2x1[2])
    @constraint(m1, c11, 3x1[1] + x1[2] <= 4)
    @constraint(m1, c12, x1[1] + 2x1[2] >= 1)
    @constraint(m1, c13, -x1[1] + x1[2] == 0.5)

    solve(m1)

    @fact getdual(c11) --> roughly(0.75, TOL)
    @fact getdual(c12) --> roughly(0.0,TOL)
    @fact getdual(c13) --> roughly(1.25,TOL)

end
end
for  conic_solver in conic_solvers_with_duals
context("With conic solver $(typeof(conic_solver))") do

    m2 = Model(solver=conic_solver)
    @variable(m2, x2[1:2] >= 0)
    @objective(m2, Max, x2[1] + 2x2[2])
    @constraint(m2, c21, 3x2[1] + x2[2] <= 4)
    @constraint(m2, c22, x2[1] + 2x2[2] >= 1)
    @constraint(m2, c23, -x2[1] + x2[2] == 0.5)
    @constraint(m2, norm(x2[1]) <= x2[2])

    solve(m2)

    @fact getdual(c21) --> roughly(0.75, TOL)
    @fact getdual(c22) --> roughly(0.0,TOL)
    @fact getdual(c23) --> roughly(1.25,TOL)

end
end
end

facts("[socduals] LP vs SOC dual (Min)") do
for lp_solver in lp_solvers
context("With LP solver $(typeof(lp_solver))") do

    m1 = Model(solver=lp_solver)
    @variable(m1, x1[1:2] >= 0)
    @objective(m1, Min, -x1[1] - 2x1[2])
    @constraint(m1, c11, 3x1[1] + x1[2] <= 4)
    @constraint(m1, c12, x1[1] + 2x1[2] >= 1)
    @constraint(m1, c13, -x1[1] + x1[2] == 0.5)

    solve(m1)

    @fact getdual(c11) --> roughly(-0.75, TOL)
    @fact getdual(c12) --> roughly(-0.0,TOL)
    @fact getdual(c13) --> roughly(-1.25,TOL)

end
end
for  conic_solver in conic_solvers_with_duals
context("With conic solver $(typeof(conic_solver))") do

    m2 = Model(solver=conic_solver)
    @variable(m2, x2[1:2] >= 0)
    @objective(m2, Min, -x2[1] - 2x2[2])
    @constraint(m2, c21, 3x2[1] + x2[2] <= 4)
    @constraint(m2, c22, x2[1] + 2x2[2] >= 1)
    @constraint(m2, c23, -x2[1] + x2[2] == 0.5)
    @constraint(m2, norm(x2[1]) <= x2[2])

    solve(m2)
    @fact getdual(c21) --> roughly(-0.75, TOL)
    @fact getdual(c22) --> roughly(-0.0,TOL)
    @fact getdual(c23) --> roughly(-1.25,TOL)

end
end
end

facts("[socduals] LP vs SOC reduced costs") do
for lp_solver in lp_solvers
context("With LP solver $(typeof(lp_solver))") do

    m1 = Model(solver=lp_solver)

    @variable(m1, x1 >= 0)
    @variable(m1, y1 >= 0)
    @constraint(m1, x1 + y1 == 1)
    @objective(m1, Max, y1)

    solve(m1)

    @fact dot([getdual(y1), getdual(x1)],[getvalue(y1); getvalue(x1)]) --> less_than_or_equal(TOL)

    @fact getvalue(x1) --> roughly(0.0,TOL)
    @fact getvalue(y1) --> roughly(1.0,TOL)

    m1 = Model(solver=lp_solver)

    @variable(m1, x1 >= 0)
    @constraint(m1, x1 <= 1)
    @objective(m1, Max, x1)

    solve(m1)

    @fact getdual(x1) --> roughly(0.0, TOL)
    @fact getvalue(x1) --> roughly(1.0,TOL)

    RHS = 2
    X_LB = 0.5
    Y_UB = 5

    m1 = Model(solver=lp_solver)

    @variable(m1, x1 >= X_LB)
    @variable(m1, y1 <= Y_UB)
    @constraint(m1, c1,  x1 + y1 == RHS)
    @objective(m1, Max, 0.8*y1 + 0.1*x1)

    solve(m1)

    @fact getvalue(x1) --> roughly(0.5, TOL)
    @fact getvalue(y1) --> roughly(1.5, TOL)

    @fact getdual(x1) --> roughly(-0.7, TOL)
    @fact getdual(y1) --> roughly(0.0, TOL)

    lp_dual_obj = getdual(c1)*RHS + getdual(x1) * X_LB + getdual(y1) * Y_UB
    @fact getobjectivevalue(m1) --> roughly(lp_dual_obj, TOL)

end
end
for  conic_solver in conic_solvers_with_duals
context("With conic solver $(typeof(conic_solver))") do

    m2 = Model(solver=conic_solver)

    @variable(m2, x2 >= 0)
    @variable(m2, y2 >= 0)
    @constraint(m2, x2 + y2 == 1)
    @objective(m2, Max, y2)
    @constraint(m2, norm(x2) <= y2)

    solve(m2)

    @fact dot([getdual(y2), getdual(x2)],[getvalue(y2); getvalue(x2)]) --> less_than_or_equal(TOL)

    @fact getvalue(x2) --> roughly(0.0,TOL)
    @fact getvalue(y2) --> roughly(1.0,TOL)

    m2 = Model(solver=conic_solver)

    @variable(m2, x2 >= 0)
    @constraint(m2, x2 <= 1)
    @objective(m2, Max, x2)
    @constraint(m2, norm(x2) <= x2)

    solve(m2)

    @fact getdual(x2) --> roughly(0.0, TOL)
    @fact getvalue(x2) --> roughly(1.0,TOL)

    RHS = 2
    X_LB = 0.5
    Y_UB = 5

    m2 = Model(solver=conic_solver)

    @variable(m2, x2 >= X_LB)
    @variable(m2, y2 <= Y_UB)
    @constraint(m2, c2, x2 + y2 == RHS)
    @objective(m2, Max, 0.8*y2 + 0.1*x2)
    @constraint(m2, norm(x2) <= y2)

    solve(m2)

    @fact getvalue(x2) --> roughly(0.5, TOL)
    @fact getvalue(y2) --> roughly(1.5, TOL)

    @fact getdual(x2) --> roughly(-0.7, TOL)
    @fact getdual(y2) --> roughly(0.0, TOL)

    conic_dual_obj = getdual(c2)*RHS + getdual(x2) * X_LB + getdual(y2) * Y_UB
    @fact getobjectivevalue(m2) --> roughly(conic_dual_obj, TOL)

end
end
end

# || b-Ax || <= a - c^T x, x\in K
# (a - c^Tx, b-Ax) \in SOC, x\in K
# ((a,b) - (c^T,A)x) \in SOC, x\in K
# if the problem is infeasible, we get a y satisfying -(a,b)^Ty > 0, (c,A^T)y \in K^*, y \in SOC^*

facts("[socduals] SOC infeasibility ray test") do
for  conic_solver in conic_solvers_with_duals
context("With conic solver $(typeof(conic_solver))") do

    m2 = Model(solver=conic_solver)

    @variable(m2, x2 >= 0)
    @constraint(m2, x2 <= 1)
    @objective(m2, Max, x2)
    @constraint(m2, c2, norm(x2) <= x2 - 1)

    status = solve(m2)

    inf_ray = getdual(c2)
    @fact status --> :Infeasible
    @fact inf_ray[1] ≥ abs(inf_ray[2]) - TOL --> true
    @fact -(-inf_ray[1]) --> greater_than_or_equal(TOL)

end
end
end


text:
export ADMMtrain1st!, euclidmap, regterm, ADMMaux

using Flux
using Flux: onehotbatch
using Juno
using Base.Iterators
using StatsBase: mean
using Zygote
using Zygote: @nograd
using LinearAlgebra
using CUDA

function euclidmap(aus, cardinality)
    if cardinality <= 0 #full DARTS discretization
        for i in 1:size(aus,1)
            if size(aus[1],1) == 8
                aus[i][1] = 0
            end
            to_mask = sortperm(aus[i])[1:length(aus[i])-1]
            aus[i][to_mask] .= 0
        end
        for i in [1,15]
            for r in 1:4
                maxes = maximum(aus[i:i+r][:], dims = 2)
                keep_rows = sortperm(maxes, rev=true)[1:2]
                for k in 1:r+1
                    if !(k in keep_rows)
                        aus[i+k-1][:] .= 0
                    else
                        aus[i+k-1][:] = aus[i+k-1][:] ./ sum(aus[i+k-1][:])
                    end
                end
                i += r+1
            end
        end
    else
        for i in 1:size(aus,1)
            to_mask = sortperm(aus[i])[1:length(aus[i])-cardinality]
            aus[i][to_mask] .= 0
        end
    end
    aus
end

function collect_αs(model::DARTSModel)
    vcat([exp.(n) for n in model.normal_αs], [exp.(n) for n in model.reduce_αs])
end

function collect_αs(model::DARTSModelBN)
    vcat([tanh.(relu.(n)) for n in model.normal_αs], [tanh.(relu.(n)) for n in model.reduce_αs])
end

"""
function collect_αs(model::DARTSModelSig)
    vcat([sigmoid.(n) for n in model.normal_αs], [sigmoid.(n) for n in model.reduce_αs])
end
"""
function regterm(m, zs, us)
    as = collect_αs(m)
    out = 0.0
    for (a, z, u) in zip(as, zs, us)
        out += sum(abs2,a - z + u)
    end
    out
end

mutable struct ADMMaux
    zs::AbstractArray
    us::AbstractArray
end

function ADMMtrain1st!(loss, model, train, val, opt_w, opt_α, zu, ρ=1f-3, losses=[0f0,0f0], epoch = 1, tmax = 50, disc = -1; cbepoch = () -> (), cbbatch = () -> ())
    zs = zu.zs
    us = zu.us
    w = all_ws_sansbn(model)
    α = all_αs(model)
    local train_loss
    local val_loss
    #admmupdate = length(train)÷epoch
    admmupdate = 10
    @show admmupdate
    ρ *= epoch
    @show ρ
    @show disc
    opt_w.os[2].t = epoch - 1
    for (i, train_batch, val_batch) in zip(1:length(train), TrainCuIterator(train), TrainCuIterator(val))
        display("grad weights")
        @time begin
            gsw = gradient(w) do
                train_loss = loss(model, train_batch...)
                return train_loss
            end
            losses[1] = train_loss
        end
        foreach(CUDA.unsafe_free!, train_batch)
        Flux.Optimise.update!(opt_w, w, gsw)
        CUDA.reclaim()
        display("grad alphas")
        @time begin
            gsα = gradient(α) do
                val_loss = loss(model, val_batch...) + ρ/2*regterm(model, zs, us)
                return val_loss
            end
            losses[2] = val_loss
        end
        foreach(CUDA.unsafe_free!, val_batch)
        Flux.Optimise.update!(opt_α, α, gsα)
        CUDA.reclaim()
        for αs in [model.normal_αs, model.reduce_αs]
            for a in αs
                a[findall(<=(0),a)] .= -Inf32
                if count(<=(0),a) == length(a)-1
                    a[findall(>(0),a)] .= Inf32
                end
            end
        end
        display(collect_αs(model))
        @time begin
            if i%admmupdate == 0
                as = collect_αs(model)
                display(as)
                zs = euclidmap(as+us, disc)
                display(zs)
                us += as - zs
                display(us)
            end
            display("admm updates")
        end
        cbbatch()
    end
    cbepoch()
    zu.zs = zs
    zu.us = us
end

max_stars_repo_name: SuReLI/DifferentiableNAS

max_issues_repo_name: SuReLI/DifferentiableNAS

max_forks_repo_name: SuReLI/DifferentiableNAS


text:
export
    VectorTextNode,
    set_text!

import ..Nodes.draw

using ..Rendering

mutable struct VectorTextNode <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    text::String

    color::Palette

    # Holds font vector lines
    mesh::Geometry.Mesh

    function VectorTextNode(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.color = Rendering.White()
        o.mesh = Geometry.Mesh()

        o
    end
end

function set_text!(node::VectorTextNode, text::String)
    node.text = text
    Nodes.set_dirty!(node, true);
end

function rebuild!(node::VectorTextNode, context::Rendering.RenderContext)
    # Use glyph properties to adjust char location.
    xpos = 0.0

    for c in node.text
        glyph = Rendering.get_glyph(context.vector_font, c)

        for vertex in glyph.vectors
            Geometry.add_vertex!(node.mesh, vertex.x + xpos, vertex.y)
        end

        xpos += context.vector_font.horizontal_offset
    end

    Geometry.build!(node.mesh)
end

function Nodes.draw(node::VectorTextNode, context::Rendering.RenderContext)
    # Transform this node's vertices using the context
    if Nodes.is_dirty(node)
        rebuild!(node, context)
        Rendering.transform!(context, node.mesh)
        Nodes.set_dirty!(node, false)
    end

    Rendering.set_draw_color(context, node.color)
    Rendering.render_lines(context, node.mesh);
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module Utils
export iftrue, next!, filtersplit, depends_on,
  normalize_mod_and_name, unique_funcname,
  flatten_blocks, MacroEnv, @MacroEnv, Dict_to_normalizedType,
  SortExpr, sortexpr,
  @doc_signature, isnodocumentationfound

using ExprParsers

next!(iter::Base.Iterators.Stateful) = iterate(iter)[1]

macro iftrue(expr)
  expr = esc(expr)
  quote
    $expr || return false
  end
end

filtersplit(f, a) = filter(f, a), filter(!f, a)

""" check whether a certain Expr builds upon given symbols """
depends_on(expr::Base.Expr, symbols) = depends_on(expr.args, symbols)
depends_on(args::Vector, symbols) = any(a -> depends_on(a, symbols), args)
depends_on(symbol::Base.Symbol, symbols) = symbol in symbols
depends_on(any, symbols) = false

# flatten out blocks
flatten_blocks(expr::Expr) = flatten_blocks(Val{expr.head}(), expr.args)
flatten_blocks(any) = any
flatten_blocks(::Val{head}, args) where head = Expr(head, flatten_blocks.(args)...)
function flatten_blocks(head::Val{:block}, args)
  args′ = flatten_blocks.(args)
  args′′ = [((a isa Expr && a.head == :block) ? a.args : [a] for a in args′)...;]
  Expr(:block, args′′...)
end

@Base.kwdef struct MacroEnv
  source::LineNumberNode
  mod::Module
end

macro MacroEnv()
  quote
    MacroEnv($(esc(:__source__)), $(esc(:__module__)))
  end
end

function Dict_to_normalizedType(d::AbstractDict)
  # TODO performance improvement possible - Symbol is called once for sorting, and once for conversion, could be combined
  ks = d |> keys |> collect |> a -> sort!(a, by=Symbol)
  rows = [Pair{Symbol(k), Symbol(d[k])} for k in ks]
  Tuple{rows...}
end

function normalize_mod_and_name(mod, name)
  parser = EP.AnyOf(EP.anysymbol, EP.NestedDot())
  normalize_mod_and_name(mod, parse_expr(parser, name))
end
normalize_mod_and_name(mod, name::Symbol) = mod, name
function normalize_mod_and_name(mod, name::EP.NestedDot_Parsed)
  mod′::Module = getproperty(mod, name.base)
  for field in name.properties[1:end-1]
    mod′ = getproperty(mod′, field)
  end
  mod′, name.properties[end]  # NestedDot.properties is known to be non-empty
end
function unique_funcname(mod, funcname)
  mod′, funcname′ = normalize_mod_and_name(mod, funcname)
  Symbol(mod′, :., funcname′)
end

include("exprparsers.jl")
include("SortExpr.jl")
include("DocsHelper.jl")
using .DocsHelper

end  # module


text:
using DelayedKalmanFilter
using DelimitedFiles
using Test

@testset "Likelihood tests" begin
    loading_path = string(pwd(),"/data/");
    protein_at_observations = readdlm(string(loading_path,"kalman_filter_test_trace_observations.csv"),',');
    model_parameters = [10000.0,5.0,log(2)/30,log(2)/90,1.0,1.0,29.0];
    measurement_variance = 10000;

    # test negative parameter check works
    negative_parameter_1 = [10000.0,5.0,log(2)/30,log(2)/90,-1.0,1.0,29.0];
    negative_parameter_2 = [-10000.0,5.0,log(2)/30,log(2)/90,1.0,1.0,29.0];

    ll = calculate_log_likelihood_at_parameter_point(protein_at_observations,negative_parameter_1,measurement_variance);
    @test ll == -Inf
    ll = calculate_log_likelihood_at_parameter_point(protein_at_observations,negative_parameter_2,measurement_variance);
    @test ll == -Inf
    # test log likelihood and derivative are correct
    ll = calculate_log_likelihood_at_parameter_point(protein_at_observations,model_parameters,measurement_variance);
    @test ll ≈ readdlm(string(loading_path,"log_likelihood_value.csv"),',')[1]

    ll, lld = calculate_log_likelihood_and_derivative_at_parameter_point(protein_at_observations,model_parameters,measurement_variance);
    @test ll ≈ readdlm(string(loading_path,"log_likelihood_value.csv"),',')[1]
    @test lld ≈ reshape(readdlm(string(loading_path,"log_likelihood_derivative_value.csv"),','),7)
end


text:
# ---
# title: 1502. Can Make Arithmetic Progression From Sequence
# id: problem1502
# author: zhwang
# date: 2022-03-22
# difficulty: Easy
# categories: Array, Sort
# link: <https://leetcode.com/problems/can-make-arithmetic-progression-from-sequence/description/>
# hidden: true
# ---
#
# Given an array of numbers `arr`. A sequence of numbers is called an arithmetic
# progression if the difference between any two consecutive elements is the
# same.
#
# Return `true` if the array can be rearranged to form an arithmetic
# progression, otherwise, return `false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,5,1]
#     Output: true
#     Explanation: We can reorder the elements as [1,3,5] or [5,3,1] with differences 2 and -2 respectively, between each consecutive elements.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,2,4]
#     Output: false
#     Explanation: There is no way to reorder the elements to obtain an arithmetic progression.
#
#
#
#
# **Constraints:**
#
#   * `2 <= arr.length <= 1000`
#   * `-10^6 <= arr[i] <= 10^6`
#
#
## @lc code=start
using LeetCode

can_make_arith(arr::Vector{Int})::Bool = length(unique!(diff(sort(arr)))) == 1

## @lc code=end

## @lc test=start
@testset "1502.can-make-arithmetic-progression-from-sequence.jl" begin
    @test can_make_arith([3,5,1])
    @test !can_make_arith([1,2,4])
end
## @lc test=end

## @lc info=start
# link: [solution to 1502]()
## @lc info=end


text:
using QuantumControlBase.QuantumPropagators: propstep!, write_to_storage!, get_from_storage!
using QuantumControlBase: evalcontrols!
using QuantumControlBase.Functionals: make_chi
using QuantumControlBase.ConditionalThreads: @threadsif
using LinearAlgebra
using Printf

"""Optimize a control problem using Krotov's method.

```julia
result = optimize_krotov(problem)
```

optimizes the given
control [`problem`](@ref QuantumControlBase.ControlProblem),
returning a [`KrotovResult`](@ref).

!!! note

    It is recommended to call [`optimize`](@ref QuantumControlBase.optimize)
    with `method=:krotov` instead of calling `optimize_krotov` directly.

Keyword arguments that control the optimization are taken from the keyword
arguments used in the instantiation of `problem`.

# Required problem keyword arguments

* `J_T`: A function `J_T(ϕ, objectives)` that evaluates the final time
  functional from a list `ϕ` of forward-propagated states and
  `problem.objectives`.

# Optional problem keyword arguments

The following keyword arguments are supported (with default values):

* `chi`: A function `chi!(χ, ϕ, objectives)` what receives a list `ϕ`
  of the forward propagates state and must set ``-χₖ=∂J_T/∂⟨ϕₖ|``. If not
  given, it will be automatically determined from `J_T` via [`make_chi`](@ref)
* `force_zygote=false`: Whether to force the use of automatic differentiation
  when calling [`make_chi`](@ref).
* `sigma=nothing`: Function that calculate the second-order contribution. If
  not given, the first-order Krotov method is used.
* `iter_start=0`: the initial iteration number
* `iter_stop=5000`: the maximum iteration number
* `prop_method`/`fw_prop_method`/`bw_prop_method`: The propagation method to
  use for each objective, see below.
* `update_hook`: A function that receives the Krotov workspace, the iteration
  number, the list of updated pulses and the list of guess pulses as
  positional arguments. The function may mutate any of its arguments. This may
  be used e.g. to apply a spectral filter to the updated pulses, or to update
  propagation workspaces inside the Krotov workspace.
* `info_hook`: A function that receives the same argumens as `update_hook`, in
  order to write information about the current iteration to the screen or to a
  file. The default `info_hook` prints a table with convergence information to
  the screen. Runs after `update_hook`. The `info_hook` function may return a
  tuple, which is stored in the list of `records` inside the
  [`KrotovResult`](@ref) object.
* `check_convergence`: a function to check whether convergence has been
  reached. Receives a [`KrotovResult`](@ref) object `result`, and should set
  `result.converged` to `true` and `result.message` to an appropriate string in
  case of convergence. Multiple convergence checks can be performed by chaining
  functions with `∘`. The convergence check is performed after any calls to
  `update_hook` and `info_hook`.
* `verbose=false`: If `true`, print information during initialization

The propagation method for the forward propagation of each objective is
determined by the first available item of the following:

* a `fw_prop_method` keyword argument
* a `prop_method` keyword argument
* a property `fw_prop_method` of the objective
* a property `prop_method` of the objective
* the value `:auto`

The propagation method for the backword propagation is determined similarly,
but with `bw_prop_method` instead of `fw_prop_method`.
"""
function optimize_krotov(problem)
    sigma = get(problem.kwargs, :sigma, nothing)
    iter_start = get(problem.kwargs, :iter_start, 0)
    update_hook! = get(problem.kwargs, :update_hook, (args...) -> nothing)
    info_hook = get(problem.kwargs, :info_hook, print_table)
    check_convergence! = get(problem.kwargs, :check_convergence, res -> res)
    # note: the default `check_convergence!` is a no-op. We still always check
    # for "Reached maximum number of iterations" in `update_result!`
    verbose = get(problem.kwargs, :verbose, false)
    skip_initial_forward_propagation =
        get(problem.kwargs, :skip_initial_forward_propagation, false)

    wrk = KrotovWrk(problem; verbose)

    ϵ⁽ⁱ⁾ = wrk.pulses0
    ϵ⁽ⁱ⁺¹⁾ = wrk.pulses1

    if skip_initial_forward_propagation
        @info "Skipping initial forward propagation"
    else
        @threadsif wrk.use_threads for (k, obj) in collect(enumerate(wrk.objectives))
            krotov_initial_fw_prop!(ϵ⁽ⁱ⁾, wrk.result.states[k], obj.initial_state, k, wrk)
        end
    end

    # TODO: if sigma, fw_storage0 = fw_storage
    update_result!(wrk, 0)
    update_hook!(wrk, 0, ϵ⁽ⁱ⁺¹⁾, ϵ⁽ⁱ⁾)
    info_tuple = info_hook(wrk, 0, ϵ⁽ⁱ⁺¹⁾, ϵ⁽ⁱ⁾)
    (info_tuple !== nothing) && push!(wrk.records, info_tuple)

    i = wrk.result.iter  # = 0, unless continuing from previous optimization
    while !wrk.result.converged
        i = i + 1
        krotov_iteration(wrk, ϵ⁽ⁱ⁾, ϵ⁽ⁱ⁺¹⁾)
        update_result!(wrk, i)
        update_hook!(wrk, i, ϵ⁽ⁱ⁺¹⁾, ϵ⁽ⁱ⁾)
        info_tuple = info_hook(wrk, i, ϵ⁽ⁱ⁺¹⁾, ϵ⁽ⁱ⁾)
        (info_tuple !== nothing) && push!(wrk.records, info_tuple)
        check_convergence!(wrk.result)
        ϵ⁽ⁱ⁾, ϵ⁽ⁱ⁺¹⁾ = ϵ⁽ⁱ⁺¹⁾, ϵ⁽ⁱ⁾
    end

    finalize_result!(ϵ⁽ⁱ⁾, wrk)

    return wrk.result

end

function krotov_initial_fw_prop!(ϵ⁽⁰⁾, ϕₖ, ϕₖⁱⁿ, k, wrk)
    Φ₀ = wrk.fw_storage[k]
    copyto!(ϕₖ, ϕₖⁱⁿ)
    (Φ₀ !== nothing) && write_to_storage!(Φ₀, 1, ϕₖⁱⁿ)
    N_T = length(wrk.result.tlist) - 1
    for n = 1:N_T
        G, dt = _fw_gen(ϵ⁽⁰⁾, k, n, wrk)
        propstep!(ϕₖ, G, dt, wrk.fw_prop_wrk[k])
        (Φ₀ !== nothing) && write_to_storage!(Φ₀, n + 1, ϕₖ)
    end
    # TODO: allow a custom propstep! routine
end

function krotov_iteration(wrk, ϵ⁽ⁱ⁾, ϵ⁽ⁱ⁺¹⁾)

    ϕ = wrk.result.states  # assumed to contain the results of forward propagation
    χ = wrk.bw_states
    J_T_func = wrk.kwargs[:J_T]
    force_zygote = get(wrk.kwargs, :get_zygote, false)
    chi! = get(wrk.kwargs, :chi, make_chi(J_T_func, wrk.objectives; force_zygote))
    N_T = length(wrk.result.tlist) - 1
    N = length(wrk.objectives)
    L = length(wrk.controls)
    X = wrk.bw_storage
    Φ = wrk.fw_storage  # TODO: pass in Φ₁, Φ₀ as parameters
    ∫gₐdt = wrk.g_a_int
    Im = imag

    # backward propagation
    chi!(χ, ϕ, wrk.objectives)
    @threadsif wrk.use_threads for k = 1:N
        write_to_storage!(X[k], N_T + 1, χ[k])
        for n = N_T:-1:1
            local (G, dt) = _bw_gen(ϵ⁽ⁱ⁾, k, n, wrk)
            propstep!(χ[k], G, dt, wrk.bw_prop_wrk[k])
            write_to_storage!(X[k], n, χ[k])
        end
    end

    # pulse update and forward propagation

    @threadsif wrk.use_threads for k = 1:N
        copyto!(ϕ[k], wrk.objectives[k].initial_state)
    end

    ∫gₐdt .= 0.0
    for n = 1:N_T  # `n` is the index for the time interval
        dt = wrk.result.tlist[n+1] - wrk.result.tlist[n]
        for k = 1:N
            get_from_storage!(χ[k], X[k], n)
        end
        ϵₙ⁽ⁱ⁺¹⁾ = [ϵ⁽ⁱ⁾[l][n] for l ∈ 1:L]  # ϵₙ⁽ⁱ⁺¹⁾ ≈ ϵₙ⁽ⁱ⁾ for non-linear controls
        # TODO: we could add a self-consistent loop here for ϵₙ⁽ⁱ⁺¹⁾
        Δuₙ′ = zeros(L)  # for step size 1
        for l = 1:L  # `l` is the index for the different controls
            ∂ϵₗ╱∂u::Function = wrk.parametrization[l].de_du_derivative
            uₗ = wrk.parametrization[l].u_of_epsilon
            for k = 1:N  # k is the index over the objectives
                ∂Hₖ╱∂ϵₗ::Union{Function,Nothing} = wrk.control_derivs[k][l]
                if !isnothing(∂Hₖ╱∂ϵₗ)
                    μₗₖₙ = (∂Hₖ╱∂ϵₗ)(ϵₙ⁽ⁱ⁺¹⁾[l])
                    if wrk.is_parametrized[l]
                        ∂ϵₗ╱∂uₗ = (∂ϵₗ╱∂u)(uₗ(ϵₙ⁽ⁱ⁺¹⁾[l]))
                        Δuₙ′[l] += ∂ϵₗ╱∂uₗ * Im(dot(χ[k], μₗₖₙ, ϕ[k]))
                    else
                        Δuₙ′[l] += Im(dot(χ[k], μₗₖₙ, ϕ[k]))
                    end
                end
            end
        end
        # TODO: second order update
        for l = 1:L
            Sₗ = wrk.update_shapes[l]
            λₐ = wrk.lambda_vals[l]
            αₗ = (Sₗ[n] / λₐ)  # Krotov step size
            Δuₗₙ = αₗ * Δuₙ′[l]
            if wrk.is_parametrized[l]
                uₗ = wrk.parametrization[l].u_of_epsilon
                ϵₗ = wrk.parametrization[l].epsilon_of_u
                ϵ⁽ⁱ⁺¹⁾[l][n] = ϵₗ(uₗ(ϵ⁽ⁱ⁾[l][n]) + Δuₗₙ)
            else
                Δϵₗₙ = Δuₗₙ
                ϵ⁽ⁱ⁺¹⁾[l][n] = ϵ⁽ⁱ⁾[l][n] + Δϵₗₙ
            end
            (∫gₐdt)[l] += αₗ * abs(Δuₙ′[l])^2 * dt
        end
        # TODO: end of self-consistent loop
        @threadsif wrk.use_threads for k = 1:N
            local (G, dt) = _fw_gen(ϵ⁽ⁱ⁺¹⁾, k, n, wrk)
            propstep!(ϕ[k], G, dt, wrk.fw_prop_wrk[k])
            write_to_storage!(Φ[k], n, ϕ[k])
        end
        # TODO: update sigma
    end  # time loop
end

# The dynamical generator for the forward propagation
function _fw_gen(ϵ, k, n, wrk)
    vals_dict = wrk.vals_dict[k]
    t = wrk.result.tlist
    for (l, control) in enumerate(wrk.controls)
        vals_dict[control] = ϵ[l][n]
    end
    dt = t[n+1] - t[n]
    evalcontrols!(wrk.G[k], wrk.objectives[k].generator, vals_dict)
    return wrk.G[k], dt
end

# The dynamical generator for the backward propagation
function _bw_gen(ϵ, k, n, wrk)
    vals_dict = wrk.vals_dict[k]
    t = wrk.result.tlist
    for (l, control) in enumerate(wrk.controls)
        vals_dict[control] = ϵ[l][n]
    end
    dt = t[n+1] - t[n]
    evalcontrols!(wrk.G[k], wrk.adjoint_objectives[k].generator, vals_dict)
    return wrk.G[k], -dt
end

function update_result!(wrk::KrotovWrk, i::Int64)
    res = wrk.result
    J_T_func = wrk.kwargs[:J_T]
    res.J_T_prev = res.J_T
    res.J_T = J_T_func(res.states, wrk.objectives)
    (i > 0) && (res.iter = i)
    if i >= res.iter_stop
        res.converged = true
        res.message = "Reached maximum number of iterations"
        # Note: other convergence checks are done in user-supplied
        # check_convergence routine
    end
    prev_time = res.end_local_time
    res.end_local_time = now()
    res.secs = Dates.toms(res.end_local_time - prev_time) / 1000.0
    # TODO: calculate τ values
end

function finalize_result!(ϵ_opt, wrk::KrotovWrk)
    res = wrk.result
    res.end_local_time = now()
    for l = 1:length(ϵ_opt)
        res.optimized_controls[l] = discretize(ϵ_opt[l], res.tlist)
    end
end

"""Print optimization progress as a table.

This functions serves as the default `info_hook` for an optimization with
Krotov's method.
"""
function print_table(wrk, iteration, args...)
    J_T = wrk.result.J_T
    g_a_int = sum(wrk.g_a_int)
    J = J_T + g_a_int
    ΔJ_T = J_T - wrk.result.J_T_prev
    ΔJ = ΔJ_T + g_a_int
    secs = wrk.result.secs

    iter_stop = "$(get(wrk.kwargs, :iter_stop, 5000))"
    widths = [max(length("$iter_stop"), 6), 11, 11, 11, 11, 11, 8]

    if iteration == 0
        header = ["iter.", "J_T", "∫gₐ(t)dt", "J", "ΔJ_T", "ΔJ", "secs"]
        for (header, w) in zip(header, widths)
            print(lpad(header, w))
        end
        print("\n")
    end

    strs = (
        "$iteration",
        @sprintf("%.2e", J_T),
        @sprintf("%.2e", g_a_int),
        @sprintf("%.2e", J),
        (iteration > 0) ? @sprintf("%.2e", ΔJ_T) : "n/a",
        (iteration > 0) ? @sprintf("%.2e", ΔJ) : "n/a",
        @sprintf("%.1f", secs),
    )
    for (str, w) in zip(strs, widths)
        print(lpad(str, w))
    end
    print("\n")
    flush(stdout)
end