text:
module PortHamiltonianBenchmarkSystems

using LinearAlgebra, SparseArrays

include("IOFormats.jl")
include("Downloads.jl")

"""
   `gugercin_pH_msd_chain(; n_cells=50, m=2, c_i=1.0, m_i=4.0, k_i=4.0)`

This function returns the port Hamiltonian mass-spring-damper system described in
S. Gugercin et al.:
      Structure-preserving tangential interpolation for model reduction of
      port-Hamiltonian systems
# Arguments
- `n_cells`: The number of masses. The system dimension is `2n_cells`
- `c_i`: The amount of damping
- `m_i`: The weight of the masses
- `k_i`: The stiffness of the springs
# Outputs
Matrices: ``J, R, Q, B``. The resulting transfer function is ``H(s) = B^\\mathsf{T} Q  (sI-(J-R)Q)^{-1}B``.
"""
function gugercin_pH_msd_chain(;
    n_cells=50::Int,
    m=2::Int,
    c_i=1.0,
    m_i=4.0,
    k_i=4.0
  )
  n=2*n_cells;
  # B is initialized as dense matrix. Since all results of transfer function
  # computations will lead to dense results.
  B=zeros(n,m);
  [B[2*i,i]=1.0 for i in 1:m]
  J=spzeros(n,n);
  [J[i,i+1]=1.0 for i in 1:2:(n-1)];
  J=J-J'
  # Set constants.
  R=spzeros(n,n);
  [R[i,i]=c_i for i in 2:2:n]
  Q=spzeros(n,n);
  Q[1,1]=k_i
  [Q[i,i]=2*k_i for i in 3:2:(n-1)]
  [Q[i,i]=1/m_i for i in 2:2:n]
  [Q[i,i+2]=-k_i for i in 1:2:(n-2)]
  [Q[i+2,i]=-k_i for i in 1:2:(n-2)]
  return J, R, Q, B
end

"""
`poro_elasticity_model(;
    n = 980,
    rho = 1e-3,
    alpha = 0.79,
    M = 1/7.80e3,
    kappanu = 633.33,
    eta = 1e-4,
    force_download = false
  )`

This function returns a port-Hamiltonian model of linear poroelasticity in a
bounded Lipschitz domain as described in
    Altmann, Mehrmann, Unger: Port-Hamiltonian Formulations of Poroelastic
    Network Models
# Arguments
- `n`: System dimension (can only be either: 320, 980, or 1805). Default = 980.
- `rho`: density. Default = `1e-3`.
- `alpha`: Biot-Willis fluid-solid coupling coefficient. Default = 0.79.
- `bm`: Biot-Modulus. Default = `1/7.8e3`.
- `kappanu`: Quotient kappa/Nu, where kappa denotes the permeability and nu denotes the fluid viscosity. Default = 633.33.
- `eta`: artificial damping coefficient. Default = `1e-4`.
"""
function poro_elasticity_model(;
    n = 980,
    rho = 1e-3,
    alpha = 0.79,
    bm = 1/7.80e3,
    kappanu = 633.33,
    eta = 1e-4,
    force_download = false
  )
  Y, D, M, K, Bp, Bf, A = load_poro_raw_data(
    n=n,
    force_download = force_download
  )
  Y = rho*sparse(Y)
  D = alpha*sparse(D)
  M = 1/bm .* sparse(M)
  K = kappanu*sparse(K)
  A = sparse(A)
  Bp = Bp'
  Bf = Bf'
  n = size(A, 1);
  m = size(M, 1);
  E = [Y spzeros(n,n+m); spzeros(n,n) A spzeros(n,m); spzeros(m,n+n) M];
  J = [spzeros(n,n) -A D';A spzeros(n,n+m); -D spzeros(m,n+m)];
  R = [spzeros(n,2*n+m); spzeros(n,2*n+m); spzeros(m,2*n) K] + eta*I
  B = [zeros(n,1); Bf; Bp];
  return E, J, R, B
end

function load_poro_raw_data(;
    n = 980,
    force_download = false
  )
  filename = "poro-n$n.mat"
  url = "https://zenodo.org/record/5702554/files/poro-n$n.mat?download=1"
  if n == 980
    md5_hash = hex(0x2961a189be7049ffe2d476b18cb1f678)
  elseif n == 320
    md5_hash = hex(0x97afe8c34f0e9a56bbe86d0a51b7b626)
  elseif n == 1805
    md5_hash = hex(0xc61f6687da9cd26cbf2d880d7d3a9ac9)
  else
    throw(ArgumentError("Model size is either 320, 980, or 1805"))
  end
  download_system_data_if_required(
    filename,
    url,
    md5_hash,
    force_download = force_download
  )
  dd = loadMAT(get_filepath(filename))
  return dd["Y"], dd["D"], dd["M"], dd["K"], dd["Bp"], dd["Bf"], dd["A"]
end

include("RCLLadders.jl")

export gugercin_pH_msd_chain, poro_elasticity_model

end # module


text:
module CUSOLVER

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcusolver, @allowscalar, assertscalar, unsafe_free!, @retry_reclaim

using ..CUBLAS: cublasFillMode_t, cublasOperation_t, cublasSideMode_t, cublasDiagType_t
using ..CUSPARSE: cusparseMatDescr_t

using CEnum

using Memoize

using DataStructures

# core library
include("libcusolver_common.jl")
include("error.jl")
include("libcusolver.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

# thread cache for task-local library handles
const thread_dense_handles = Vector{Union{Nothing,cusolverDnHandle_t}}()
const thread_sparse_handles = Vector{Union{Nothing,cusolverSpHandle_t}}()

# cache for created, but unused handles
const old_dense_handles = DefaultDict{CuContext,Vector{cusolverDnHandle_t}}(()->cusolverDnHandle_t[])
const old_sparse_handles = DefaultDict{CuContext,Vector{cusolverSpHandle_t}}(()->cusolverSpHandle_t[])

function dense_handle()
    CUDA.detect_state_changes()
    tid = Threads.threadid()
    if @inbounds thread_dense_handles[tid] === nothing
        ctx = context()
        thread_dense_handles[tid] = get!(task_local_storage(), (:CUSOLVER, :dense, ctx)) do
            handle = if isempty(old_dense_handles[ctx])
                cusolverDnCreate()
            else
                pop!(old_dense_handles[ctx])
            end

            finalizer(current_task()) do task
                push!(old_dense_handles[ctx], handle)
            end
            # TODO: cusolverDnDestroy to preserve memory, or at exit?

            cusolverDnSetStream(handle, stream())

            handle
        end
    end
    something(@inbounds thread_dense_handles[tid])
end

function sparse_handle()
    CUDA.detect_state_changes()
    tid = Threads.threadid()
    if @inbounds thread_sparse_handles[tid] === nothing
        ctx = context()
        thread_sparse_handles[tid] = get!(task_local_storage(), (:CUSOLVER, :sparse, ctx)) do
            handle = if isempty(old_sparse_handles[ctx])
                cusolverSpCreate()
            else
                pop!(old_sparse_handles[ctx])
            end

            finalizer(current_task()) do task
                push!(old_sparse_handles[ctx], handle)
            end
            # TODO: cusolverSpDestroy to preserve memory, or at exit?

            cusolverSpSetStream(handle, stream())

            handle
        end
    end
    something(@inbounds thread_sparse_handles[tid])
end

@inline function set_stream(stream::CuStream)
    ctx = context()
    tls = task_local_storage()
    dense_handle = get(tls, (:CUSOLVER, :dense, ctx), nothing)
    if dense_handle !== nothing
        cusolverDnSetStream(dense_handle, stream)
    end
    sparse_handle = get(tls, (:CUSOLVER, :sparse, ctx), nothing)
    if sparse_handle !== nothing
        cusolverSpSetStream(sparse_handle, stream)
    end
    return
end

function __init__()
    resize!(thread_dense_handles, Threads.nthreads())
    fill!(thread_dense_handles, nothing)

    resize!(thread_sparse_handles, Threads.nthreads())
    fill!(thread_sparse_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_dense_handles[tid] = nothing
        thread_sparse_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_dense_handles[tid] = nothing
        thread_sparse_handles[tid] = nothing
    end
end

end

ext: jl


text:
function test_fixed_window(
    values, op; op_bang=nothing, approximate_equality::Bool=false
)
    # For later comparison, to ensure that we don't mutate any of these
    values_copy = deepcopy(values)

    T = typeof(first(values))

    function makestate(window)
        return if isnothing(op_bang)
            FixedWindowAssociativeOp{T,op}(window)
        else
            FixedWindowAssociativeOp{T,op,op_bang}(window)
        end
    end

    @test_throws ArgumentError makestate(0)

    for window in 1:(2 * length(values))
        state = makestate(window)
        for (i, value) in enumerate(values)
            @test update_state!(state, value) == state

            received = window_value(state)
            # This is the expected range of the input values over which we are reducing.
            i_start = max(1, i - window + 1)
            i_stop = i
            expected = reduce(op, values[i_start:i_stop])
            if approximate_equality
                @test received ≈ expected
            else
                @test received == expected
            end

            if i < window
                @test !window_full(state)
                @test window_size(state) == i
            else
                @test window_full(state)
                @test window_size(state) == window
            end
        end
    end

    @test values == values_copy
end

@testset "FixedWindowAssociativeOp" begin
    @testset "trivial" begin
        T = Int64
        op = +
        window = 2
        state = FixedWindowAssociativeOp{T,op}(window)

        @test update_state!(state, 3) == state
        @test !window_full(state)
        @test window_size(state) == 1
        @test window_value(state) == 3

        @test update_state!(state, 4) == state
        @test window_full(state)
        @test window_size(state) == window
        @test window_value(state) == 7

        @test update_state!(state, 4) == state
        @test window_full(state)
        @test window_size(state) == window
        @test window_value(state) == 8

        @test update_state!(state, -5) == state
        @test window_full(state)
        @test window_size(state) == window
        @test window_value(state) == -1
    end

    @testset "integral" begin
        for op in (+, *, max, min)
            test_fixed_window(1:20, op)
        end
    end

    @testset "random float" begin
        values = rand(20)
        for op in (+, *, max, min)
            test_fixed_window(values, op; approximate_equality=true)
        end
    end

    @testset "random matrix multiply" begin
        # This isn't commutative.
        values = [rand(-5:5, (2, 2)) for _ in 1:20]
        test_fixed_window(values, *)
    end

    @testset "set union" begin
        values = [Set([x]) for x in 1:20]
        test_fixed_window(values, union)
    end

    @testset "mutating" begin
        values = [Data(1, rand(-5:5, (2, 2))) for _ in 1:20]
        test_fixed_window(values, merge)
        test_fixed_window(values, merge; op_bang=merge!)
    end
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using Random
export BeamSearch

struct BeamSearch <: LocalSearchAlgorithm
    ssize::Int32  # sample size
    bsize::Int32  # beam size

    BeamSearch() = new(1, 1)
    BeamSearch(a::Integer, b::Integer) = new(a, b)
    BeamSearch(other::BeamSearch) =  new(other.ssize, other.bsize)
end

# const BeamType = typeof((objID=Int32(0), dist=0.0))
### local search algorithm
function beam_init(bs::BeamSearch, index::SearchGraph, dist::Function, q, res::KnnResult, navigation_state, hints)
    beam = KnnResult(bs.bsize)
    n = length(index.db)
    # range = 1:n
    # @inbounds for i in 1:bs.bsize
    @inbounds for objID in hints
        S = get(navigation_state, objID, UNKNOWN)
        # S = navigation_state[objID]
        if S == UNKNOWN
            navigation_state[objID] = VISITED
            d = dist(q, index.db[objID])
            push!(beam, objID, d) && push!(res, objID, d)
        end
    end

    beam
end

"""
Tries to reach the set of nearest neighbors specified in `res` for `q`.
- `bs`: the parameters of `BeamSearch`
- `index`: the local search index
- `q`: the query
- `res`: The result object, it stores the results and also specifies the kind of query
"""
function search(bs::BeamSearch, index::SearchGraph, dist::Function, q, res::KnnResult, navigation_state, hints=EMPTY_INT_VECTOR)
    n = length(index.db)
    n == 0 && return res
    if length(hints) == 0
        hints = rand(1:n, bs.bsize)
    end

    beam = beam_init(bs, index, dist, q, res, navigation_state, hints)
    prev_score = typemax(Float64)

    @inbounds while abs(prev_score - last(beam).dist) > 0.0  # prepared to allow early stopping
        prev_score = last(beam).dist

        for prev in beam
            cov = last(beam).dist
            S = get(navigation_state, prev.objID, UNKNOWN)
            # S = navigation_state[prev.objID]
            S == EXPLORED && continue
            navigation_state[prev.objID] = EXPLORED
            for childID in index.links[prev.objID]
                S = get(navigation_state, childID, UNKNOWN)
                # S = navigation_state[childID]
                if S == UNKNOWN
                    navigation_state[childID] = VISITED
                    d = dist(q, index.db[childID])

                    if d <= cov
                        push!(beam, childID, d) && push!(res, childID, d)
                    end
                end
            end
        end
    end

    res
end

function opt_expand_neighborhood(fun, gsearch::BeamSearch, n::Integer, iter::Integer, probes::Integer)
    #f_(w) = ceil(Int, w * (rand() - 0.5))
    #f(x, w) = max(1, x + #f_(w))

    # g(x) = max(1, x + ceil(Int, (rand()-0.5) * log2(n)))
    logn = ceil(Int, log(2, n+1))
    probes = probes == 0 ? logn : probes
    f(x) = max(1, x + rand(-logn:logn))
    for i in 1:probes
        BeamSearch(f(gsearch.ssize), f(gsearch.bsize)) |> fun
    end
    ### f(x, w) = max(1, x + w)
    ### w = 1
    ### while w <= logn  ## log log n
    ###    BeamSearch(f(gsearch.ssize,  w), gsearch.bsize) |> fun
    ###    BeamSearch(f(gsearch.ssize,  -w), gsearch.bsize) |> fun
    ###    BeamSearch(gsearch.ssize, f(gsearch.bsize, w)) |> fun
    ###    BeamSearch(gsearch.ssize, f(gsearch.bsize, -w)) |> fun
    ###    w += w
    ### end
end


text:
@testset "linj_00001.tif" begin
filename = joinpath(testbase, "linj_00001.tif")
    # $ ScanImageTiffReader image shape linj_00001.tif
    # Shape: 512 x 512 x 10 @ i16
    tsize = ScanImageTiffReader.open(filename) do io
        size(io)
    end
    @test tsize == (512, 512, 10)
    ttype = ScanImageTiffReader.open(filename) do io
        pxtype(io)
    end
    @test ttype == Int16

    # $ ScanImageTiffReader image bytes linj_00001.tif
    # 5 MB
    dat = ScanImageTiffReader.open(filename) do io
        data(io)
    end
    @test abs(sizeof(dat)/2^20 - 5.0) < 1e-5
    # @test ScanImageTiffReader.open(filename, length) == size(dat, 3) # TODO: check this assumption

    # $ ScanImageTiffReader descriptions --frame 0 linj_00001.tif
    # "frameNumbers = 1"
    # "acquisitionNumbers = 1"
    # "frameNumberAcquisition = 1"
    # "frameTimestamps_sec = 0.000000"
    # "acqTriggerTimestamps_sec = "
    # "nextFileMarkerTimestamps_sec = "
    # "endOfAcquisition =  0"
    # "endOfAcquisitionMode = 0"
    # "dcOverVoltage = 0"
    # "epoch = [2016 6 4 13 51 7.8046]"
    desc = ScanImageTiffReader.open(filename) do io
        split(description(io, 1), "\n")
    end
    @test desc[1] == "frameNumbers = 1"
    @test desc[2] == "acquisitionNumbers = 1"
    @test desc[3] == "frameNumberAcquisition = 1"
    @test desc[4] == "frameTimestamps_sec = 0.000000"
    @test desc[5] == "acqTriggerTimestamps_sec = "
    @test desc[6] == "nextFileMarkerTimestamps_sec = "
    @test desc[7] == "endOfAcquisition =  0"
    @test desc[8] == "endOfAcquisitionMode = 0"
    @test desc[9] == "dcOverVoltage = 0"
    @test desc[10] == "epoch = [2016 6 4 13 51 7.8046]"

    # $ ScanImageTiffReader metadata linj_00001.tif
    # {
    #   "SI": {
    #     "TIFF_FORMAT_VERSION": 3,
    #     "VERSION_MAJOR": "2015",
    #     "VERSION_MINOR": "4",
    #     "acqState": "grab",
    #     "acqsPerLoop": 1,
    #     "extTrigEnable": 0,
    #     "hBeams": {
    #       "beamCalibratedStatus": [0,0,0],
    # [...]
    md = ScanImageTiffReader.open(filename) do io
        JSON.parse(metadata(io))
    end
    # only descend two levels...
    @test md["SI"]["imagingSystem"] == "LinScanner"
    @test md["SI"]["acqState"] == "grab"
    @test md["SI"]["loopAcqInterval"] == 10
    @test md["SI"]["VERSION_MINOR"] == "4"
    @test md["SI"]["extTrigEnable"] == 0
    @test md["SI"]["TIFF_FORMAT_VERSION"] == 3
    @test md["SI"]["acqsPerLoop"] == 1
    @test md["SI"]["VERSION_MAJOR"] == "2015"
    @test md["SI"]["objectiveResolution"] == 16
    @test md["RoiGroups"]["photostimRoiGroups"] == nothing
end


text:
#============================================================================================
			Time Evolve
============================================================================================#

export euler, Solution, solveAdvection

#export euler
function euler(rhs::Function, dt::T,
               s0::AdvectionState{T})::AdvectionState{T} where {T}
    r0 = rhs(s0)
    s1 = s0 + dt * r0
    AdvectionState{T}(s0.time + dt, s1.u)
end

struct Solution{T}
    dx::T
    dt::T
    states::Vector{AdvectionState{T}}
end

function solveAdvection(tmax::T, nlines::Int, lambda::T)::Solution{T} where {T}
    dx = 1 / nlines
    dt = lambda * dx
    nsteps = round(Int, tmax / dt)
    sol = Solution{T}(dx, dt, AdvectionState{T}[])
    s = initialGaussian(nlines)
    push!(sol.states, s)
    for step in 1:nsteps
        s = euler(rhsAdvection, dt, s)
        push!(sol.states, s)
    end
    sol
end


text:
function cudnnCreate()
    handle = Ref{cudnnHandle_t}()
    cudnnCreate(handle)
    return handle[]
end

function cudnnGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  cudnnGetProperty(property, value_ref)
  value_ref[]
end

version() = VersionNumber(cudnnGetProperty(CUDAapi.MAJOR_VERSION),
                          cudnnGetProperty(CUDAapi.MINOR_VERSION),
                          cudnnGetProperty(CUDAapi.PATCH_LEVEL))

function cuda_version()
  ver = cudnnGetCudartVersion()
  major, ver = divrem(ver, 1000)
  minor, patch = divrem(ver, 10)

  VersionNumber(major, minor, patch)
end


text:
function check_tnep_status(sol)
    for (idx,val) in sol["ne_branch"]
        @test val["built"] == 0.0 || val["built"] == 1.0
    end
end

@testset "test soc tnep" begin
    @testset "3-bus case" begin
        result = run_tnep("../test/data/case3_tnep.m", SOCWRPowerModel, pajarito_solver; setting = Dict("output" => Dict("branch_flows" => true)))

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 2; atol = 1e-2)
    end
    @testset "5-bus rts case" begin
        result = run_tnep("../test/data/case5_tnep.m", SOCWRPowerModel, pajarito_solver)

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 1; atol = 1e-2)
    end
end

@testset "test dc tnep" begin
    @testset "3-bus case" begin
        result = run_tnep("../test/data/case3_tnep.m", DCPPowerModel, pajarito_solver)

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 2; atol = 1e-2)
    end

    #=
    # skip this one becouse it is breaking Julia package tests
    @testset "5-bus case" begin
        result = run_tnep("../test/data/case5_tnep.m", DCPPowerModel, pajarito_solver)

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 1; atol = 1e-2)
    end
    =#
end

@testset "test dc-losses tnep" begin
    @testset "3-bus case" begin
        result = run_tnep("../test/data/case3_tnep.m", DCPLLPowerModel, pajarito_solver)

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 2; atol = 1e-2)
    end

    @testset "5-bus case" begin
        result = run_tnep("../test/data/case5_tnep.m", DCPLLPowerModel, pajarito_solver)

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 1; atol = 1e-2)
    end
end

if (Pkg.installed("AmplNLWriter") != nothing && Pkg.installed("CoinOptServices") != nothing)

    @testset "test ac tnep" begin
        @testset "5-bus case" begin
            result = run_tnep("../test/data/case5_tnep.m", ACPPowerModel, BonminNLSolver(["bonmin.bb_log_level=0", "bonmin.nlp_log_level=0"]))

            check_tnep_status(result["solution"])

            @test result["status"] == :LocalOptimal
            @test isapprox(result["objective"], 1; atol = 1e-2)
        end
    end

end

@testset "test tnep branch flow output" begin
    @testset "3-bus case" begin
        result = run_tnep("../test/data/case3_tnep.m", SOCWRPowerModel, pajarito_solver; setting = Dict("output" => Dict("branch_flows" => true)))

        check_tnep_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 2; atol = 1e-2)

        branches = result["solution"]["branch"]
        ne_branches = result["solution"]["ne_branch"]
        flow_keys = ["pf","qf","pt","qt"]

        for fk in flow_keys
            @test !isnan(branches["1"][fk])
            @test !isnan(ne_branches["1"][fk])
            @test !isnan(ne_branches["2"][fk])
            @test !isnan(ne_branches["3"][fk])
        end
    end
end


text:
# # A demonstration of controlling models competing in parallel

# Controlling competing iterative models in parallel is useful in
# optimization. For example, one might try several different
# optimization strategies in parallel, and divert resources to the
# strategy or strategies that are performing best.

# In this toy illustration, several anteaters of varying health
# compete for food from a common supply - a stream of random "ant" or
# inedible "seed" items. An anteater's `health` attribute is the
# probability he is able to successfully catch an ant when the
# source makes it available.

# If an anteater's score (the number of ants consumed) falls too far
# behind that of the leader, he stops consuming (his "training" stops
# early). An anteater also stops if he reaches the specified goal.

# Anteaters forage for food in parallel using Julia multi-threading.

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# **Julia version** is assumed to be 1.6.*

using IterationControl
using .Threads
using Random
Random.seed!(123)
using Plots
pyplot(size = (600, 300*(sqrt(5) - 1)))

# ## Create a channel to supply food

N = 1_000_000
food = Channel{String}() do ch
    foreach(i -> put!(ch, rand(["ant", "seed"])), 1:N)
end

[take!(food) for i in 1:3]

# ## Define anteaters and how to "train" them

mutable struct Anteater
    health::Float64        # 0.0 = dead; 1.0 = excellent health
    n::Int                 # number of ants consumed
    Anteater(health) = new(health, 0)
end

function eat!(eater::Anteater)
    take!(food) == "ant" && rand() <= eater.health && (eater.n += 1)
    return eater
end

# "Training" for `n` iterations means take `n` items from the `food`
#  source and try to catch and eat each item if it is an "ant":

function IterationControl.train!(eater::Anteater, n)
    foreach(i->eat!(eater), 1:n)
    return eater
end

# ## Create a collection of competitors

eaters = [Anteater(0.7), Anteater(0.8), Anteater(0.9), Anteater(1.0)]

#-

n_eaters = length(eaters);

# For recording progress:
history = Float64[];

# ## Define the controls that will define the rules of the competition

const LEAD_TRIGGERING_DEFEAT = 5
const ANTS_TO_WIN = 100

scores(eaters) = [e.n for e in eaters]

is_too_far_behind(eater) =
    maximum(scores(eaters)) - eater.n > LEAD_TRIGGERING_DEFEAT

has_won(eater) = eater.n >= ANTS_TO_WIN

update_history(::Any) =  append!(history, scores(eaters))

controls = [Step(1),
            Callback(is_too_far_behind, stop_if_true=true),
            Callback(has_won, stop_if_true=true),
            Callback(update_history),
            NumberLimit(1000)]

# ## Run the competition

@sync for e in eaters
    Threads.@spawn IterationControl.train!(e, controls...)
end

# ## Plot the results

n_events = div(length(history), n_eaters)
history = reshape(history, n_eaters, n_events)'

plot(history,
     title="Ant eater competition",
     ylab="num of ants consumed",
     xlab="num food items taken from source")

# Generally, only the healthiest anteater is still drawing on food
# items at the end of the race (the use of multi-threading precludes
# reproducibility here).

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=true) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src


text:
"""
    gradient_hessian_unequal_rates(η, i, j, t)

Compute the gradient of the log-probability of a simple birth and death process
evaluated at the point ``η = (λ, μ)^{\\prime}``. Variable `i` is
the initial population size, `j` is the final population size, and `t` is the
elapsed time.
"""
function gradient_hessian_unequal_rates(
  η::Vector{F},
  i::F,
  j::F,
  t::F
)::Tuple{Vector{F}, Symmetric{F, Matrix{F}}} where {
  F <: AbstractFloat
}
  g1, g2, h11, h21, h22 = if η[1] > η[2]
    gradient_hessian_v1(η, i, j, t)
  else
    gradient_hessian_v2(η, i, j, t)
  end

  [g1, g2], Symmetric([h11 h21; h21 h22])
end

"""
    gradient_hessian_v1(η, i, j, t)

Compute the gradient of the log-probability of a simple birth and death process
evaluated at the point ``η = (λ, μ)^{\\prime}`` where `η[1] > η[2]`. Variable
`i` is the initial population size, `j` is the final population size, and `t`
is the elapsed time.
"""
function gradient_hessian_v1(
  η::Vector{F},
  i::F,
  j::F,
  t::F
)::Tuple{F, F, F, F, F} where {
  F <: AbstractFloat
}
  a, b = (j <= i) ? (i, j) : (j, i)

  θ = (η[2] - η[1]) * t
  ω = log(η[2] / η[1])

  ψ = η[1] - η[2]
  ρ = η[1] + η[2]
  σ = η[1] * η[2]

  ξ = -ω / ψ

  if j > 0
    u1 = zero(F)
    u2 = zero(F)
    fa = zero(F)

    if (i > 1) && (j > 1)
      q1, q2, q3 = if t <= ξ
        x = logexpm1(θ - ω) + log1mexp(θ + ω) - 2 * log1mexp(θ)
        log_hypergeometric_joint(a, b, x)
      else
        x = log(σ) + θ + 2 * (logexpm1(-θ) - log(ψ))
        log_meixner_ortho_poly_joint(a, b, x)
      end

      u1 = q2 - q1
      u2 = q3 - q2
      fa = (i - 1) * (j - 1) * exp(u2) / (i + j - 2) -
           i * j * exp(u1) / (i + j - 1)
    else
      # hypergeometric(i - 1, j - 1, z, k=0) is one because either i or j is one
      # hypergeometric(i - 2, j - 2, z, k=-1) is zero
      u1 = log(σ) - θ + 2 * (log1mexp(θ) - log(ψ))
      u2 = -Inf
      fa = -i * j * exp(u1) / (i + j - 1)
    end

    w0 = log(-θ)
    w1 = log1mexp(θ)
    w2 = log1pexp(θ)
    w3 = expm1(ω + θ)
    w4 = exp(θ)
    w5 = exp(2 * θ)
    w6 = expm1(θ)
    w7 = exp(ω + θ)
    w8 = exp(2 * ω + θ)

    y1 = w4 + 1
    y2 = (w3 * w6)^2
    y3 = exp(u1 + θ - 4 * w1)

    f1 = expm1(log(η[1] / ρ) + w0 + w2 - w1)
    f2 = expm1(log(η[2] / ρ) + w0 + w2 - w1)

    # firt-order partial derivative with respect to λ
    x1 = expm1(w0 + θ - w1) / w3
    x2 = ρ * exp(u1 + w0 + θ - 2 * w1) / (σ * t)

    r1 = -((i + j) * x1 + i * j * f1 * x2 / (i + j - 1) - j) / η[1]

    # firt-order partial derivative with respect to μ
    x1 = w7 * expm1(w0 - w1) / w3

    r2 = ((i + j) * x1 + i * j * f2 * x2 / (i + j - 1) + i) / η[2]

    # second-order partial derivative with respect to λ twice
    x1 = (1 + w4 * ((2 - θ) * expm1(log(η[2] * t) + 2 * θ) +
                    (1 - θ - (η[1] + 3 * η[2]) * t) * w4 +
                    ((1 + θ) * η[1] + η[2]) * t)) / y2
    x2 = y3 / σ
    x3 = (η[2] * w6 - η[1] * θ * y1)^2 +
         (η[2] * w6)^2 -
         2 * (w6 * y1 - θ * w4) * θ * η[1]^2
    x4 = w4 * (ψ * ρ * f1)^2 / σ

    r3 = ((i + j) * x1 +
          i * j * x2 * (x3 + x4 * fa) / (i + j - 1) -
          j) * η[1]^(-2)

    # second-order partial derivative with respect to λ and then μ
    x1 = w4 * (1 + ((1 + θ) * η[1] + η[2]) * t -
               (2 - θ + (η[1] + 3 * η[2]) * t) * w4 +
               (η[2] * t + (1 - θ) * (η[2] * t + 1)) * w5) / y2
    x2 = y3 * η[2]^(-2)
    x3 = (η[1]^2 + η[2]^2) * w6^2 +
         σ * θ^2 * (y1^2 + 2 * w4) -
         ρ^2 * θ * y1 * w6
    x4 = w7 * (ψ * ρ / η[2])^2 * f1 * f2

    r4 = -((i + j) * x1 +
           i * j * x2 * (x3 + x4 * fa) / (i + j - 1)) * η[1]^(-2)

    # second-order partial derivative with respect to μ twice
    x1 = w8 * ((2 + θ) * η[1] * t +
               w4 * (1 + θ - (3 * η[1] + η[2]) * t -
                     (2 + θ - (η[1] + (1 - θ) * η[2]) * t) * w4 + w5)) / y2
    x2 = y3 / σ
    x3 = (η[1] * w6 - η[2] * θ * y1)^2 +
         (η[1] * w6)^2 -
         2 * (w6 * y1 - θ * w4) * θ * η[2]^2
    x4 = w4 * (ψ * ρ * f2)^2 / σ

    r5 = ((i + j) * x1 +
          i * j * x2 * (x3 + x4 * fa) / (i + j - 1) -
          i) * η[2]^(-2)

    r1, r2, r3, r4, r5
  else
    w1 = exp(θ)
    w2 = expm1(θ)
    w3 = expm1(ω + θ)

    y1 = η[1] * w3 * w2
    y2 = y1^2

    # hypergeometric(i - 1, j - 1, z, k=0) and
    # hypergeometric(i - 2, j - 2, z, k=-1) are both zero
    dλ_logϕ = -(1 - w1 * (1 - θ)) / y1
    dμ_logϕ = w1 * (1 + θ - w1) / y1
    dλλ_logϕ = (1 - w1 * (2 * (1 - η[2] * t) - θ * η[1] * t -
                          w1 * (1 - 4 * η[2] * t +
                                w1 * (2 - θ) * η[2] * t))) / y2
    dλμ_logϕ = -w1 * (1 + (θ * η[1] + ρ) * t -
                      w1 * (2 * (1 + ρ * t) -
                            w1 * (1 - (θ * η[2] - ρ) * t))) / y2
    dμμ_logϕ = w1 * ((2 + θ) * η[1] * t +
                     w1 * (1 - 4 * η[1] * t -
                           w1 * (2 * (1 - η[1] * t) + θ * η[2] * t - w1))) / y2

    r1 = i * dλ_logϕ
    r2 = i * (1 / η[2] + dμ_logϕ)
    r3 = i * dλλ_logϕ
    r4 = i * dλμ_logϕ
    r5 = i * (-1 / η[2]^2 + dμμ_logϕ)

    r1, r2, r3, r4, r5
  end
end

"""
    gradient_hessian_v2(η, i, j, t)

Compute the gradient of the log-probability of a simple birth and death process
evaluated at the point ``η = (λ, μ)^{\\prime}`` where `η[2] > η[1]`. Variable
`i` is the initial population size, `j` is the final population size, and `t`
is the elapsed time.
"""
function gradient_hessian_v2(
  η::Vector{F},
  i::F,
  j::F,
  t::F
)::Tuple{F, F, F, F, F} where {
  F <: AbstractFloat
}
  a, b = (j <= i) ? (i, j) : (j, i)

  θ = (η[1] - η[2]) * t
  ω = log(η[1] / η[2])

  ψ = η[1] - η[2]
  ρ = η[1] + η[2]
  σ = η[1] * η[2]

  ξ = ω / ψ

  if j > 0
    u1 = F(0)
    u2 = F(0)
    fa = F(0)

    if (i > 1) && (j > 1)
      q1, q2, q3 = if t <= ξ
        x = logexpm1(θ - ω) + log1mexp(θ + ω) - 2 * log1mexp(θ)
        log_hypergeometric_joint(a, b, x)
      else
        x = log(2 * σ) + logexpm1(log1pexp(2 * θ) - (θ + log(F(2)))) -
            2 * log(η[2] - η[1])
        log_meixner_ortho_poly_joint(a, b, x)
      end

      u1 = q2 - q1
      u2 = q3 - q2
      fa = (i - 1) * (j - 1) * exp(u2) / (i + j - 2) -
           i * j * exp(u1) / (i + j - 1)
    else
      # hypergeometric(i - 1, j - 1, z, k=0) is one because either i or j is one
      # hypergeometric(i - 2, j - 2, z, k=-1) is zero
      u1 = log(σ) - θ + 2 * (log1mexp(θ) - log(η[2] - η[1]))
      u2 = -Inf
      fa = -i * j * exp(u1) / (i + j - 1)
    end

    w0 = log(-θ)
    w1 = log1mexp(θ)
    w2 = log1pexp(θ)
    w3 = expm1(ω + θ)
    w4 = exp(θ)
    w5 = exp(2 * θ)
    w6 = expm1(θ)
    w7 = exp(ω + θ)
    w8 = exp(2 * ω + θ)

    y1 = w4 + 1
    y2 = (w3 * w6)^2
    y3 = exp(u1 + θ - 4 * w1)

    f1 = expm1(log(η[1] / ρ) + w0 + w2 - w1)
    f2 = expm1(log(η[2] / ρ) + w0 + w2 - w1)

    # firt-order partial derivative with respect to λ
    x1 = -exp(θ + ω - log1mexp(θ + ω)) * expm1(w0 - w1)
    x2 = ρ * exp(u1 + w0 + θ - 2 * w1) / (σ * t)

    r1 = ((i + j) * x1 + i * j * f1 * x2 / (i + j - 1) + j) / η[1]

    # firt-order partial derivative with respect to μ
    x1 = expm1(w0 + θ - w1) / w3

    r2 = -((i + j) * x1 + i * j * f2 * x2 / (i + j - 1) - i) / η[2]

    # second-order partial derivative with respect to λ twice
    x1 = w8 * ((2 + θ) * η[2] * t +
               w4 * (1 + θ - (η[1] + 3 * η[2]) * t -
                     (2 + θ - ((1 - θ) * η[1] + η[2]) * t) * w4 + w5)) / y2
    x2 = y3 / σ
    x3 = (η[1] * θ * y1 - η[2] * w6)^2 +
          (η[2] * w6)^2 +
          2 * (θ * w4 - w6 * y1) * θ * η[1]^2
    x4 = w4 * (ψ * ρ * f1)^2 / σ

    r3 = ((i + j) * x1 +
           i * j * x2 * (x3 + x4 * fa) / (i + j - 1) -
           j) * η[1]^(-2)

    # second-order partial derivative with respect to λ and then μ
    x1 = w8 * ((1 + θ) + (2 + θ) * η[2] * t -
               (2 + θ + (η[1] + 3 * η[2]) * t) * w4 +
               (1 + ((1 - θ) * η[1] + η[2]) * t) * w5) / y2
    x2 = y3 * η[2]^(-2)
    x3 = (η[1]^2 + η[2]^2) * w6^2 +
         σ * θ^2 * (y1^2 + 2 * w4) -
         ρ^2 * θ * y1 * w6
    x4 = w4 * (ψ * ρ)^2 * f1 * f2 / σ

    r4 = -((i + j) * x1 +
           i * j * x2 * (x3 + x4 * fa) / (i + j - 1)) * η[1]^(-2)

    # second-order partial derivative with respect to μ twice
    x1 = (1 + w4 * ((2 - θ) * expm1(log(η[1] * t) + 2 * θ) +
                    (1 - θ - (3 * η[1] + η[2]) * t) * w4 +
                    (η[1] + (1 + θ) * η[2]) * t)) / y2
    x2 = y3 / σ
    x3 = (η[2] * θ * y1 - η[1] * w6)^2 +
         (η[1] * w6)^2 +
         2 * (θ * w4 - w6 * y1) * θ * η[2]^2
    x4 = w4 * (ψ * ρ * f2)^2 / σ

    r5 = ((i + j) * x1 +
          i * j * x2 * (x3 + x4 * fa) / (i + j - 1) -
          i) * η[2]^(-2)

    r1, r2, r3, r4, r5
  else
    w1 = exp(θ)
    w2 = expm1(θ)
    w3 = expm1(ω + θ)

    y1 = η[2] * w3 * w2
    y2 = y1^2

    # hypergeometric(i - 1, j - 1, z, k=0) and
    # hypergeometric(i - 2, j - 2, z, k=-1) are both zero
    dλ_logϕ = w1 * (1 + θ - w1) / y1
    dμ_logϕ = -(1 - w1 * (1 - θ)) / y1
    dλλ_logϕ = w1 * ((2 + θ) * η[2] * t -
                     w1 * (4 * η[2] * t - 1 +
                           w1 * (θ * η[1] * t + 2 * (1 - η[2] * t) - w1))) / y2
    dλμ_logϕ = -w1 * (1 + (θ * η[2] + ρ) * t -
                      w1 * (2 * (1 + ρ * t) -
                            w1 * (1 - (θ * η[1] - ρ) * t))) / y2
    dμμ_logϕ = (1 + w1 * (θ * η[2] * t + 2 * (η[1] * t - 1) +
                    w1 * ((1 - 4 * η[1] * t) +
                    w1 * (2 - θ) * η[1] * t))) / y2

    r1 = i * dλ_logϕ
    r2 = i * (1 / η[2] + dμ_logϕ)
    r3 = i * dλλ_logϕ
    r4 = i * dλμ_logϕ
    r5 = i * (-1 / η[2]^2 + dμμ_logϕ)

    r1, r2, r3, r4, r5
  end
end


text:
using Statistics, MultivariateStats, RDatasets, LinearAlgebra, Plots; pyplot()

data = dataset("datasets", "iris")
data = data[:,[:SepalLength,:SepalWidth,:PetalLength,:PetalWidth]]
x = convert(Array{Float64,2},data)'

model = fit(PCA, x, maxoutdim=4, pratio = 0.999)
M = projection(model)

function manualProjection(x)
    covMat = cov(x')
    ev = eigvals(covMat)
    eigOrder = sortperm(eigvals(covMat),rev=true)
    eigvecs(covMat)[:,eigOrder]
end

println("Manual vs. package: ",maximum(abs.(M-manualProjection(x))))

pcVar = principalvars(model) ./ tvar(model)
cumVar = cumsum(pcVar)
pcDat = M[:,1:2]'*x

p1 = plot(pcVar, c=:blue, label="Variance due to PC")
	plot!(1:length(cumVar), cumVar, label="Cumulative Variance", c=:red,
        xlabel="Principle component",ylabel="Variance Proportion",ylims=(0,1))
p2 = scatter(pcDat[1,:],pcDat[2,:], c=:blue, xlabel="PC 1", ylabel="PC 2",
            msw=0, legend=:none)
plot(p1, p2, size=(800,400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#
# Date created: 2022-03-10
# Author: aradclif
#
#
############################################################################################
# Why not just exp! the result of lvlogsoftmax? Consider that such an approach would perform
# the operations:
# C[i_1,…] = A[i_1,…] - B[i_1, 1, …]
# C[i_1,…] = exp(C[i_1,…])
# One can combine this in order to traverse the memory once, i.e.
# C[i_1,…] = exp(A[i_1,…] - B[i_1, 1, …])

function aminusb_expbody(N::Int, D)
    body = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> Symbol(:i_, d), N)...)
    e = Expr(:(=), c, Expr(:call, :exp, Expr(:call, :-, a, b)))
    push!(body.args, e)
    body
end

function aminusb_exp_quote(N::Int, D)
    ls = loopgen(N)
    body = aminusb_expbody(N, D)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return C
    end
end

@generated function aminusb_exp!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    aminusb_exp_quote(N, D)
end

function lvsoftmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = lvsoftmax1(A)
    else
        B = lvlse(A, dims=dims)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(B, d), N)
        C = similar(A, promote_type(T, eltype(B)))
        aminusb_exp!(C, A, B, Dᴮ′)
    end
    return C
end

lvsoftmax(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvsoftmax(A, (dims,))
lvsoftmax(A::AbstractArray{T, N}; dims=:) where {T, N} = lvsoftmax(A, dims)
lvsoftmax(A::AbstractArray{T, N}) where {T, N} = lvsoftmax1(A)
lvsoftmax(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvsoftmax1(A)

function lvsoftmax1(A::AbstractArray{T, N}) where {T, N}
    b = lvlse1(A)
    C = similar(A, Base.promote_op(exp, T))
    @turbo for i ∈ eachindex(A)
        C[i] = exp(A[i] - b)
    end
    C
end

################ threaded version

function taminusb_exp_quote(N::Int, D)
    ls = loopgen(N)
    body = aminusb_expbody(N, D)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return C
    end
end

@generated function taminusb_exp!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    taminusb_exp_quote(N, D)
end

function lvtsoftmax(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = lvtsoftmax1(A)
    else
        B = lvtlse(A, dims=dims)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : size(B, d), N)
        C = similar(A, promote_type(T, eltype(B)))
        taminusb_exp!(C, A, B, Dᴮ′)
    end
    return C
end

lvtsoftmax(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtsoftmax(A, (dims,))
lvtsoftmax(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtsoftmax(A, dims)
lvtsoftmax(A::AbstractArray{T, N}) where {T, N} = lvtsoftmax1(A)
lvtsoftmax(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtsoftmax1(A)

function lvtsoftmax1(A::AbstractArray{T, N}) where {T, N}
    b = lvtlse1(A)
    C = similar(A, Base.promote_op(exp, T))
    @tturbo for i ∈ eachindex(A)
        C[i] = exp(A[i] - b)
    end
    C
end


text:
@testset "SpecialFunctions" begin
    @testset "besselj" begin
        # Check recurrence relations
        # besselj(ν, z) = z / 2ν * (besselj(ν - 1, z) + besselj(ν + 1, z))
        rng = MersenneTwister(42)
        setprecision(Arb, 64) do
            for _ = 1:1000
                deg = rand(rng, 0:10)
                ν = Arb(20rand(rng) - 10)
                z = ArbSeries(20rand(rng, deg + 1) .- 10)

                lhs = besselj(ν, z)
                rhs = z / 2ν * (besselj(ν - 1, z) + besselj(ν + 1, z))
                @test Arblib.degree(lhs) == Arblib.degree(rhs) == deg
                @test all(Arblib.overlaps.(lhs[:], rhs[:]))
            end
        end

        z = ArbSeries([1, 2, 3, 4, 5])
        @test isequal(besselj0(z), besselj(zero(Arb), z))
        @test isequal(besselj1(z), besselj(one(Arb), z))
    end
end


text:
module Button

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export button, buttongroup

Genie.Renderer.Html.register_normal_element("q__btn", context = @__MODULE__)
Genie.Renderer.Html.register_normal_element("q__btn__group", context = @__MODULE__)

function button(
                label::String = "",
                args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...)

  wrap() do
    q__btn(args...; attributes([:label => label, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

function buttongroup(args...;
                    wrap::Function = StippleUI.DEFAULT_WRAPPER,
                    kwargs...)
  wrap() do
    q__btn__group(args...; kwargs...)
  end
end

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# 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.

using Test

# This file gets called first. It it doesn't crash, all is well.
include("issue980.jl")

for file in readdir(@__DIR__)
    if file in ["issue980.jl", "dummy.jl", "hygiene.jl", "runtests.jl"]
        continue
    elseif !endswith(file, ".jl")
        continue
    end
    @testset "$(file)" begin
        include(file)
    end
end

@testset "MOI.$(submodule)" for submodule in [
    "Nonlinear",
    "Bridges",
    "FileFormats",
    "Test",
    "Utilities",
    "Benchmarks",
]
    include("$(submodule)/$(submodule).jl")
end

# Test hygiene of @model macro
include("hygiene.jl")

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
module DerivativeTest

import Calculus

using Base.Test
using ForwardDiff

include(joinpath(dirname(@__FILE__), "utils.jl"))

########################
# test vs. Calculus.jl #
########################

x = 1

for f in NUMBER_TO_NUMBER_FUNCS
    println("  ...testing $f")
    v = f(x)
    d = ForwardDiff.derivative(f, x)
    @test_approx_eq_eps d Calculus.derivative(f, x) FINITEDIFF_ERROR

    out = DerivativeResult(zero(v))
    ForwardDiff.derivative!(out, f, x)
    @test_approx_eq ForwardDiff.value(out) v
    @test_approx_eq ForwardDiff.derivative(out) d
end

for f in NUMBER_TO_ARRAY_FUNCS
    println("  ...testing $f")
    v = f(x)
    d = ForwardDiff.derivative(f, x)
    @test_approx_eq_eps d Calculus.derivative(f, x) FINITEDIFF_ERROR

    out = similar(v)
    ForwardDiff.derivative!(out, f, x)
    @test_approx_eq out d

    out = DerivativeResult(similar(v))
    ForwardDiff.derivative!(out, f, x)
    @test_approx_eq ForwardDiff.value(out) v
    @test_approx_eq ForwardDiff.derivative(out) d
end

end # module


text:
include("SSA.jl")

function addSim(sys::System, tauVar::TauVar, n::Int64)
    ens = newEnsemble(sys, n)
    npoints = length(sys.data.dd)
    if sys.routine.ssa == :Direct
        for i in 2:npoints
            if i == 2 && sys.routine.shoot
                reShoot!(sys, ens, i-1)
            end
            _sim!(sys, ens, i)
            for k in eachindex(ens)
                for xx in ens[k].x
                    push!(ens[k].xa,xx)
                end
                push!(ens[k].ta,sys.times[i])
            end
        end
    elseif sys.routine.ssa == :Tau
        for i in 2:npoints
            if i == 2 && sys.routine.shoot
                reShoot!(sys, ens, i-1)
            end
            _sim!(sys, ens, i, tauVar)
            for k in eachindex(ens)
                for xx in ens[k].x
                    push!(ens[k].xa,xx)
                end
                push!(ens[k].ta,sys.times[i])
            end
        end
    end
    ens
end

function sim(sys::System, θ::Vector{Float64}, σ2::Vector{Float64}, n::Int64)
    sys.state.θ = θ
    sys.state.σ2 = diagm(σ2)
    tauVar = preAllocTauVar(sys)
    ens = Ensemble()
    append!(ens, addSim(sys, tauVar, n))
    ens
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
function AMPGO02()
   nlp = Model()

   @variable(nlp, x, start=0.0)

   @NLobjective(
    nlp,
    Min,
    sin(x)+sin((10.0/3.0)*x)
   )

   return nlp
end


text:
#=

=#

"""
    Group nodes into the corresponding Strongly Connected Components and
    classify each component to correctly initialize the fibration algorithm.
    From that, it provides the initial partitioning of the network together
    with the queue of pivot sets to run the fibration partitioning algorithm.

    Args:
        graph:
            'Graph' structure storing the vertices and edges' objects.
    Return:
        partition:
            List containing 'Fiber' objects. This list corresponds to the
            initial partitioning to be used for the fast fibration algorithm.
        pivot_queue:
            List containing 'Fiber' objects. This list corresponds to the queue
            of pivot sets. Each pivot set is used in one iteration of the algorithm
            to determine the next splitting.
"""
function initialize(graph::Graph)
    N = length(graph.vertices)

    scc_info = extract_strong(graph, true)
    node_labels = scc_info[1]
    unique_labels = scc_info[2]
    components = scc_info[3]

    # -- 'sccs' holds a list of 'StrongComponent' objects --
    sccs = StrongComponent[]
    for label in collect(keys(components))
        new_scc = StrongComponent()
        insert_nodes(components[label], new_scc)
        push!(sccs, new_scc)
    end

    # -- Check if each SCC receives or not input from other
    # -- components not itself --
    partition = Fiber[Fiber()]
    autopivot = Fiber[]
    for strong in sccs
        check_input(strong, graph)
        classify_strong(strong, graph)

        if strong.type == 0
            insert_nodes(strong.nodes, partition[1])
        elseif strong.type == 1
            new_fiber = Fiber()
            insert_nodes(strong.nodes, new_fiber)
            append!(partition, [new_fiber])
        else
            new_fiber = Fiber()
            insert_nodes(strong.nodes, new_fiber)
            append!(autopivot, [new_fiber])
            insert_nodes(strong.nodes, partition[1])
        end
    end

    # -- Define the pivot set queue and fiber pointer for each
    # -- node. Also push the defined classes to the queue.
    pivot_queue = Fiber[]
    fiber_index = [ -1 for j in 1:N ]
    for (index, class) in enumerate(partition)
        push!(pivot_queue, copy_fiber(class))
        for v in class.nodes
            fiber_index[v] = index
        end
        partition[index].index = index
    end
    # -- Push the isolated self-loop nodes to the queue --
    for isolated in autopivot
        push!(pivot_queue, copy_fiber(isolated))
    end
    # -- Save the graph property 'fiber_index' --
    set_vertices_properties("fiber_index", fiber_index, graph)
    return partition, pivot_queue
end

"""
    Push all the splitted classes to the pivot queue, with exception
    of the largest class.

    Args:
        new_classes:
            List of 'Fiber' objects containing all the classes that must be
            pushed to pivot set (excluding the larger one)..
        pivot_queue:
            List of 'Fiber' objects. It corresponds to the pivot set that must
            be updated with new splitted classes.
"""
function enqueue_splitted(new_classes::Array{Fiber}, pivot_queue::Array{Fiber})
    classes_size = [ fiber.number_nodes for fiber in new_classes ]
    mxval, mxind = findmax(classes_size)

    for (j, fiber) in enumerate(new_classes)
        if j!=mxind
            new_pivot = copy_fiber(fiber)
            push!(pivot_queue, new_pivot)
        end
    end
end

"""
    Returns a list of indexes of the possible input set unstable classes of the
    'graph'. Also, it returns a dictionary which maps a node to its input set. This
    function is called by the 'fast_partitioning' function.

    Args:
        graph:
            'Graph' structure storing the vertices and edges' objects.
        pivot:
            'Fiber' object representing the current pivot set which is used for the
            splitting process.
        partition:

        number_edgetype:
            Total number of edge types within the parsed 'graph'.

        eprop_name:

    Return:
        receivers:
            Array of 'Fiber' objects. Each fiber contained in this array receives at least
            one input from the parsed pivot set.
        input_sets:
            Dictionary with integer indexes and value as arrays of integers. The integer
            indexes represents the nodes of the graph, while their associated array of
            integers represents the input set with respect to the pivot set.
            Example: input_sets[2] => [0 2 1] -> Node 2 receives 2 input of type 2 and
            1 input of type 3 from the pivot set.
"""
function calculate_input_set(graph::Graph, pivot::Fiber, partition::Array{Fiber},
                             number_edgetype::Int, eprop_name::String)
    fiber_index = graph.int_vproperties["fiber_index"]
    edgetype_prop = graph.int_eproperties[eprop_name]

    receivers_classes = Int[]
    input_sets = Dict{Int, Array{Int}}()

    for pivot_v in pivot.nodes
        out_edges = graph.vertices[pivot_v].edges_source
        for edge in out_edges
            target = edge.target
            etype = edgetype_prop[edge.index]
            push!(receivers_classes, fiber_index[target])
            # In case 'target' is not included in the dict yet.
            if get(input_sets, target, -1)==-1
                input_sets[target] = [ 0 for j in 1:number_edgetype ]
            end
            input_sets[target][etype] += 1
        end
    end
    receivers_classes = collect(Int, Set(receivers_classes))
    receivers = Fiber[ copy_fiber(partition[j]) for j in receivers_classes ]
    return receivers, input_sets
end

"""
    Given the input sets of all relevant nodes of the graph, define the
    partitioning necessary to define the new input set stable partition.

    Args:
        receivers:
            Array of 'Fiber' objects. Each fiber contained in this array receives at least
            one input from the parsed pivot set. This list is provided by another function
            called 'calculate_input_set'.
        input_sets:
            Dictionary with integer indexes and value as arrays of integers. The integer
            indexes represents the nodes of the graph, while their associated array of
            integers represents the input set with respect to the pivot set.
            Example: input_sets[2] => [0 2 1] -> Node 2 receives 2 input of type 2 and
            1 input of type 3 from the pivot set. This dictionary is provided by another
            function called 'calculate_input_set'.
        number_edgetype:
            Total number of edge types within the parsed 'graph'.
    Return:
        final_splitting:
            Dictionary (key -> Integer, value -> Array of Array of Integer).
            Each key represent an index of the 'receivers' list and it holds
            the splitting of the fiber stored in the index. If the array of
            arrays contains only one array it means that the class was not
            splitted.
"""
function split_from_input_set(receivers::Array{Fiber}, input_sets::Dict{Int, Array{Int}},
                              number_edgetype::Int)

    # Default input set string for those who does not receive input from 'pivot'.
    default_str = ""
    for j in 1:number_edgetype default_str = default_str*"0" end

    final_splitting = Dict{Int, Array{Array{Int}}}()

    for j in 1:length(receivers)
        current_fiber = receivers[j]
        sub_input = Dict{String, Array{Int}}()
        for node in current_fiber.nodes
            if get(input_sets, node, -1)==-1
                if get(sub_input, default_str, -1)==-1
                    sub_input[default_str] = Int[]
                end
                push!(sub_input[default_str], node)
            else
                inputset_str = array2string(input_sets[node])
                if get(sub_input, inputset_str, -1)==-1
                    sub_input[inputset_str] = Int[]
                end
                push!(sub_input[inputset_str], node)
            end
        end
        # At this point, each key of 'sub_input' represents a new subpartition
        # of the original fiber. If there is only one key, then the fiber was not
        # divided by the current pivot set.
        final_splitting[j] = [ sub_input[key] for key in keys(sub_input) ]
    end
    return final_splitting
end

"""
    Function to be used according to the outputs of 'calculate_input_set' and 'split_from_input_set'.
    After all possible input set unstable classes are defined, all the relevant input sets are
    calculated and the proper splitting is done, this function is responsible to define the final
    partitioning for a given pivot set.

    The return is null, the main change is performed in the mutable 'partition' and 'pivot_queue'
    structures.
"""
function define_new_partition(final_splitting::Dict{Int, Array{Array{Int}}}, receivers::Array{Fiber},
                              partition::Array{Fiber}, pivot_queue::Array{Fiber}, graph::Graph)
    fiber_index = graph.int_vproperties["fiber_index"]

    for j in 1:length(receivers)
        if length(final_splitting[j])==1 continue end

        current_fiber = receivers[j]
        new_classes = Fiber[]
        nodes_to_remove = Int[]
        partition_fiber = partition[current_fiber.index]

        # Choose the first splitted class to keep in the original fiber.
        for arr in final_splitting[j][2:length(final_splitting[j])]
            new_fiber = Fiber()
            new_fiber.index = length(partition)+1
            insert_nodes(arr, new_fiber)
            for u in new_fiber.nodes
                fiber_index[u] = new_fiber.index
            end
            push!(partition, new_fiber)
            push!(new_classes, copy_fiber(new_fiber))
            append!(nodes_to_remove, arr)
        end
        delete_nodes(nodes_to_remove, partition_fiber)
        push!(new_classes, copy_fiber(partition_fiber))

        enqueue_splitted(new_classes, pivot_queue)
    end
end

"""
    Generate a input set stable partition with respect to the pivot set.

    Given a pivot set, we select all the possible unstable classes of
    the current 'partition' and define these classes as 'receivers',
    meaning classes that have input from the pivot set. After this, we
    check all edges coming out from the pivot set and calculate the input set
    for its targets.

    Using a hash table, each node receiving edges from the pivot set holds
    a list representing its input set. With these we make the proper
    splitting of the unstable classes.

    'eprop_name' {default to "edgetype"} is the edge property that should
    be used in the algorithm.

    At the end of the function, 'partition' is input-set stable with
    respect to 'pivot' and the 'pivot_queue' will be properly updated.

    Args:
        graph:
            'Graph' structure holding the vertices and edges' objects that
            defines the network.
        pivot:
            'Fiber' object representing the current pivot set that must be used
            to split unstable classes.
        partition:
            List of 'Fiber' object representing the current partitioning of the
            network. This list is modified if any class is splitted during this
            iteration.
        pivot_queue:
            List of 'Fiber' objects. It corresponds to the pivot set that must
            be updated with new splitted classes.
        n_edgetype:
            Integer value representing the number of edge types within the network.
        eprop_name:
            String representing the edge property of the network that should be used
            to query the type of each edge.
"""
function fast_partitioning(graph::Graph, pivot::Fiber, partition::Array{Fiber},
                           pivot_queue::Array{Fiber}, n_edgetype::Int, eprop_name="edgetype")
    # Define the input sets and the fiber receiving information from 'pivot'.
    receivers, input_sets = calculate_input_set(graph, pivot, partition, n_edgetype, eprop_name)
    # After defined the 'receivers' and the proper input sets, we can calculate the
    # proper splitting.
    final_splitting = split_from_input_set(receivers, input_sets, n_edgetype)
    # Define the new partitioning of the graph with respect to the pivot set.
    define_new_partition(final_splitting, receivers, partition, pivot_queue, graph)
end

"""
    Divide a directed graph into several groups according to the fibration
    partitioning of the graph.

    Args:
        graph:
            'Graph' structure where the vertices and edges are placed.
        eprop_name:
            String representing the integer edge property of the network that should
            be used to query the type of each edge.
    Return:
        partition:
            Final fibration partitioning.
"""
function fast_fibration(graph::Graph, eprop_name="edgetype")
    if !graph.is_directed
        print("Undirected network\n")
        return
    end

    edgetype_prop = graph.int_eproperties[eprop_name]
    number_edgetype = length(collect(Int, Set(edgetype_prop)))

    partition, pivot_queue = initialize(graph)
    while length(pivot_queue)>0
        pivot_set = pop!(pivot_queue)
        fast_partitioning(graph, pivot_set, partition, pivot_queue, number_edgetype)
    end
    return partition
end

"""
    Given a list of fibers, returns the count of nontrivial fibers,
    the total count of fibers and also a dictionary containing all
    the nodes for each fiber.

    Args:
        partition:
            List of 'Fiber' objects.
    Return:
        nontrivial_fibers:
            Count of fibers containing more than one node.
        total_fibers:
            Total count of fibers, including trivial and nontrivial
            fibers.
        fibers_map:
            Dictionary where the indexes are the list index of the fiber
            in the parsed 'partition' and the contained value is an array
            of integers representing the nodes inside the current fiber.
"""
function extract_fiber_groups(partition::Array{Fiber})
    nontrivial_fibers = 0
    total_fibers = 0
    fibers_map = Dict{Int, Array{Int}}()
    for j in 1:length(partition)
        fibers_map[j] = partition[j].nodes
        if partition[j].number_nodes > 1
            nontrivial_fibers += 1
        end
        total_fibers += 1
    end
    return nontrivial_fibers, total_fibers, fibers_map
end


text:
mutable struct IB # 物品栏管理器
	l::Int
	n::Vector{UInt8}
	i::Vector{Item}
end
function IB(l::Int)
	return IB(l,Vector{UInt8}(undef,(l+1)>>1),[EI for i in 1:l])
end
function get_n(ib::IB,id::Int)::UInt8
	idd=(id+1)>>1
	n=id&1==0 ? ib.n[idd]>>4 : ib.n[idd]&15
	return n+0x1
end
function set_n!(ib::IB,id::Int,x::UInt8)
	v=x-0x1
	idd=(id+1)>>1
	id&1==0 ? ib.n[idd]=ib.n[idd]&15|(v<<4) : ib.n[idd]=ib.n[idd]&240|v
end
function get(ib::IB,id::Int)
	return Pair(ib.i[id],get_n(ib,id))
end
function set!(ib::IB,id::Int,x::Pair{Item,UInt8})
	ib.i[id]=x.first
	set_n!(ib,id,x.second)
end
function remove!(ib::IB,id::Int)
	ib.i[id]=EI()
end
function clear!(ib::IB)
	for i in 1:ib.l
		ib.i[i]=EI()
	end
end
function reduce!(ib::IB,id::Int,n::Int)
	orin=get_n(ib,id)
	if n<orin
		set_n!(ib,id,orin-n)
	else
		ib.i[id]=EI()
	end
end
function give!(ib::IB,it::Item,n::Int)::Int
	last=n;st=stack(it);itid=id(it)
	for i in 1:ib.l
		if last==0
			return
		end
		iid=id(ib.i[i])
		if iid=="" # EI
			ib.i[i]=deepcopy(it)
			if last>st
				last-=st
				set_n!(ib,i,st)
			else
				set_n!(ib,i,last)
				return 0
			end
		elseif iid==itid
			orin=get_n(ib,i)
			if orin==st # 很可能出现
				continue
			elseif orin+last>st
				set_n!(ib,i,st)
				last-=(st-orin)
			else
				set_n!(ib,i,orin+last)
				return 0
			end
		end
	end
	return last
end


text:
import ScikitLearnBase: fit!

AbstractOptimiser = Flux.Optimise.AbstractOptimiser

export fit!

"""
    fit!(model::MatFacModel, D::AbstractMatrix;
         capacity=10^8,
         opt=ADAGrad(), lr=0.01, max_epochs=1000,
         abs_tol=1e-9, rel_tol=1e-6,
         verbosity=1)

Fit a MatFacModel to dataset D.

* `capacity` refers to memory capacity. It controls
  the size of minibatches during training. I.e., larger
  `capacity` means larger minibatches.
* `opt` is an optional Flux `AbstractOptimiser` object.
  Overrides the `lr` kwarg.
* `lr` is learning rate. Defaults to 0.01.
* `abs_tol` and `rel_tol` are standard convergence criteria.
* `max_epochs`: an upper bound on the number of epochs.
* `verbosity`: larger values make the method print more
  information to stdout.

We recommend loading `model` and `D` to GPU _before_ fitting:
```
model = gpu(model)
D = gpu(D)
```
"""
function fit!(model::MatFacModel, D::AbstractMatrix;
              capacity::Integer=10^8,
              max_epochs=1000, lr::Number=0.01,
              opt::Union{Nothing,AbstractOptimiser}=nothing,
              abs_tol::Number=1e-9, rel_tol::Number=1e-6,
              verbosity::Integer=1)

    #############################
    # Preparations
    #############################

    # Define a verbose print
    function vprint(a...; level=1)
        if verbosity >= level
            print(string(a...))
        end
    end

    M_D, N_D = size(D)
    K, M = size(model.X)
    N = size(model.Y,2)

    # Validate input size
    if (M != M_D)|(N != N_D)
        throw(ArgumentError, string("Incompatible sizes! Data:",
                                    size(D),
                                    "; Model: ", (M,N)
                                   )
             )
    end

    inv_MN = 1.0/(M*N)
    inv_MK = 1.0/(M*K)
    inv_NK = 1.0/(N*K)

    col_batch_size = div(capacity,M)
    row_batch_size = div(capacity,N)

    # Prep the row and column transformations
    col_layers = make_viewable(model.col_transform)
    row_layers = make_viewable(model.row_transform)

    # Define the likelihood function
    likelihood = (X,Y,
                  r_layers,
                  c_layers,
                  noise, D)-> inv_MN*invlinkloss(noise,
                                                  c_layers(
                                                   r_layers(
                                                    transpose(X)*Y
                                                   )
                                                  ),
                                               D)

    # Prep the regularizers
    col_layer_regs = make_viewable(model.col_transform_reg)
    row_layer_regs = make_viewable(model.row_transform_reg)

    col_regularizer = (layers, reg) -> inv_NK*reg(layers)
    row_regularizer = (layers, reg) -> inv_MK*reg(layers)
    X_regularizer = (X, reg) -> inv_MK*reg(X)
    Y_regularizer = (Y, reg) -> inv_NK*reg(Y)

    # Initialize some objects to store gradients
    Y_grad = zero(model.Y)
    X_grad = zero(model.X)
    col_layer_grads = fmapstructure(tozero, rec_trainable(col_layers))
    row_layer_grads = fmapstructure(tozero, rec_trainable(row_layers))
    noise_model_grads = fmapstructure(tozero, rec_trainable(model.noise_model))

    # If no optimiser is provided, initialize
    # the default (an ADAGrad optimiser)
    if opt == nothing
        opt = Flux.Optimise.ADAGrad(lr)
    end

    # Track the loss
    prev_loss = Inf
    loss = Inf

    #############################
    # Main loop
    #############################
    epoch = 1
    t_start = time()
    while epoch <= max_epochs

        vprint("Epoch ", epoch,":  ")

        loss = 0.0

        Y_grad .= 0
        col_layer_grads = fmap(tozero, col_layer_grads)
        noise_model_grads = fmap(tozero, noise_model_grads)

        ######################################
        # Iterate through the ROWS of data
        for row_batch in BatchIter(M, row_batch_size)

            X_view = view(model.X, :, row_batch)
            D_v = view(D, row_batch, :)
            row_layers_view = view(row_layers, row_batch, 1:N)
            col_layers_view = view(col_layers, row_batch, 1:N)

            # Define the likelihood for this batch
            col_likelihood = (Y, cl, noise) -> likelihood(X_view, Y,
                                                          row_layers_view, cl,
                                                          noise, D_v)

            # Accumulate the gradient
            batchloss, grads = Zygote.withgradient(col_likelihood,
                                                   model.Y,
                                                   col_layers_view,
                                                   model.noise_model)

            binop!(.+, Y_grad, grads[1])
            binop!(.+, col_layer_grads, grads[2])
            binop!(.+, noise_model_grads, grads[3])

            loss += batchloss

        end

        # Accumulate Y regularizer gradient
        regloss, Y_reg_grad = Zygote.withgradient(Y_regularizer, model.Y, model.Y_reg)
        loss += regloss
        binop!(.+, Y_grad, Y_reg_grad[1])

        # Update Y and the Y regularizer
        update!(opt, model.Y, Y_grad)
        update!(opt, model.Y_reg, Y_reg_grad[2])

        # Accumulate layer regularizer gradients
        regloss, reg_grads = Zygote.withgradient(col_regularizer, col_layers, col_layer_regs)
        loss += regloss
        binop!(.+, col_layer_grads, reg_grads[1])

        # Update layers and layer regularizers
        update!(opt, col_layers, col_layer_grads)
        update!(opt, col_layer_regs, reg_grads[2])

        # Update the noise model
        update!(opt, model.noise_model, noise_model_grads)

        X_grad .= 0
        row_layer_grads = fmap(tozero, row_layer_grads)

        ######################################
        # Iterate through the COLUMNS of data
        for col_batch in BatchIter(N, col_batch_size)

            D_v = view(D, :, col_batch)
            noise_view = view(model.noise_model, col_batch)
            Y_view = view(model.Y, :, col_batch)

            row_layers_view = view(row_layers, 1:M, col_batch)
            col_layers_view = view(col_layers, 1:M, col_batch)

            row_likelihood = (X, rl) -> likelihood(X, Y_view,
                                                   rl, col_layers_view,
                                                   noise_view, D_v)

            batchloss, g = Zygote.withgradient(row_likelihood, model.X, row_layers_view)
            binop!(.+, X_grad, g[1])
            binop!(.+, row_layer_grads, g[2])

            loss += batchloss
        end

        # Accumulate X regularizer gradients
        regloss, X_reg_grads = Zygote.withgradient(X_regularizer, model.X, model.X_reg)

        loss += regloss
        binop!(.+, X_grad, X_reg_grads[1])

        # Update X and the X regularizer
        update!(opt, model.X, X_grad)
        update!(opt, model.X_reg, X_reg_grads[2])

        # Accumulate the layer regularizer gradients
        regloss, reg_grads = Zygote.withgradient(row_regularizer, row_layers, row_layer_regs)
        loss += regloss
        binop!(.+, row_layer_grads, reg_grads[1])

        # Update the layers and regularizers
        update!(opt, row_layers, row_layer_grads)
        update!(opt, row_layer_regs, reg_grads[2])

        # Report the loss
        elapsed = time()-t_start
        vprint("Loss=",loss, " (", round(Int, elapsed), "s elapsed)\n")

        # Check termination conditions
        loss_diff = prev_loss - loss
        if loss_diff < abs_tol
            vprint("Terminated: reached abs_tol<",abs_tol, "\n"; level=0)
            break
        elseif loss_diff/loss < abs_tol
            vprint("Terminated: reached rel_tol<",rel_tol, "\n"; level=0)
            break
        else
            prev_loss = loss
            epoch += 1
        end
    end
    if epoch >= max_epochs
        vprint("Terminated: reached max_epochs=",max_epochs, "\n"; level=0)
    end

    #############################
    # Termination
    #############################

    return model
end

max_stars_repo_name: dpmerrell/BatchMatFac

max_issues_repo_name: dpmerrell/BatchMatFac

max_forks_repo_name: dpmerrell/BatchMatFac


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb02qd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    MMAX = 20
    NRHSMX = 20
    LDA = MMAX
    LDB = max( MMAX, NMAX )
    LDWORK = max(   min( MMAX, NMAX) + 3*NMAX + 1, 2*min( MMAX, NMAX) + NRHSMX )
    A = Array{Float64,2}(undef, LDA,NMAX)
    B = Array{Float64,2}(undef, LDB,NRHSMX)
    Y = Array{Float64,1}(undef, NMAX*NRHSMX)
    JPVT = Array{BlasInt,1}(undef, NMAX)
    SVAL = Array{Float64,1}(undef, 3)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    NRHS = parse(BlasInt, vs[3])
    RCOND = parse(Float64, replace(vs[4],'D'=>'E'))
    SVLMAX = parse(Float64, replace(vs[5],'D'=>'E'))
    JOB =  vs[6][1]
    INIPER =  vs[7][1]
if ( M<0 || M>MMAX )
else
if ( N<0 || N>NMAX )
else
if ( NRHS<0 || NRHS>NRHSMX )
else

    vs = String[]
    _isz,_jsz = (M,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,NRHS)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    RANK, INFO = SLICOT.mb02qd!(JOB, INIPER, M, N, NRHS, RCOND, SVLMAX, A, B, Y, JPVT, SVAL, LDWORK)
    println(io, "RANK = $RANK")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"B:")
    _nc = NRHS
    _nr = N
    show(io,"text/plain",B[1:_nr,1:_nc])
    println(io,)

end # if
end # if
end # if
end # if
    close(f)
end # run_X()


text:
@testset "453.minimum-moves-to-equal-array-elements.jl" begin
    @test min_moves_453([1,2,3]) == 3
end


text:
using SampleDistributions
using Distributions
using Base.Test

include("distributionsapi.jl")


text:
using LinearAlgebra, SparseArrays

export vec2mat, vec2matinv, vec2matadj

"""
    vec2mat(v)
Convert the vector `v` to a symmetric matrix by stacking the entries.

# Examples
```jldoctest
julia> vec2mat([1, 2, 3, 4, 5, 6])
3×3 Array{Int64,2}:
 1  2  4
 2  3  5
 4  5  6
```
"""
function vec2mat(v::SparseVector{Tv,Ti}) where {Tv, Ti}
    (indices, values) = findnz(v)
    sd = Int( 0.5*(sqrt(8*length(v)+1) - 1) )
    mat = spzeros(Tv, sd,sd)
    col = 1
    total = 0
    for i in 1:length(indices)
        idx = indices[i]
        while !(idx in total+1:total+col)
            total += col
            col += 1
        end
        mat[col, idx-total] = values[i]
        mat[idx-total, col] = values[i]
    end
    return mat
end
function vec2mat(v::Array{T,1}) where T
    sd = Int( 0.5*(sqrt(8*length(v)+1) - 1) )
    mat = zeros(T, sd,sd)
    col = 1
    total = 0
    for j = 1:sd
        mat[1:j,j] = v[total+1:total+j]
        mat[j,1:(j-1)] = v[total+1:total+j-1]
        total += j
    end
    return mat
end

"""
    vec2matadj(M)
Convert matrix `M` to a vector such that `dot(vec2matadj(M), u) == dot(M, vec2mat(u))`
for any vector `u`.

# Examples
```jldoctest
julia> M = [1 2; 2 3];

julia> vec2matadj(M)
3-element Array{Int64,1}:
 1
 4
 3

julia> using LinearAlgebra: dot

julia> dot(vec2matadj(M), [0, 1, 0]) == dot(M, vec2mat([0, 1, 0]))
true
```
"""
function vec2matadj(M::Union{SparseMatrixCSC{Tv,Ti}, Array{Tv,2}}) where {Tv, Ti}
    @assert issymmetric(M)
    sd = size(M,1)
    v = typeof(M) <: SparseMatrixCSC ?
    spzeros(Int( 0.5*sd *(sd+1) )) : zeros(Tv, Int( 0.5*sd *(sd+1) ))
    total = 0
    for j in 1:sd
        # Double the off-diagonal entries
        v[total+1:total+j-1] = 2 * M[1:(j-1),j]
        v[total+j] = M[j,j]
        total += j
    end
    return v
end

"""
    vec2matinv(M)
Convert matrix `M` to a vector such that `vec2mat(vec2matinv(M))` equals `M`.

# Examples
```jldoctest
julia> M = vec2mat([1, 2, 3, 4, 5, 6])
3×3 Array{Int64,2}:
 1  2  4
 2  3  5
 4  5  6

julia> vec2matinv(M)
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6
```
"""
function vec2matinv(M::Union{SparseMatrixCSC{Tv,Ti}, Array{Tv,2}, Symmetric{Tv,Ti}}) where {Tv, Ti}
    @assert issymmetric(M)
    sd = size(M,1)
    vecsize = sd * (sd+1) ÷ 2
    v = if M isa SparseMatrixCSC spzeros(vecsize)
    elseif M isa Array zeros(Tv, vecsize)
    else Array{Tv,1}(undef, vecsize)
    end
    total = 0
    for j in 1:sd
        v[total+1:total+j] = M[1:j,j]
        total += j
    end
    return v
end


text:
function packmultihypo(fnc::CommonConvWrapper{T}) where {T<:FunctorInferenceType}
  fnc.hypotheses != nothing ? string(fnc.hypotheses) : ""
end
function parsemultihypostr(str::AS) where {AS <: AbstractString}
  mhcat=nothing
  if length(str) > 0
    mhcat = extractdistribution(str)
  end
  return mhcat
end

## packing converters-----------------------------------------------------------
# heavy use of multiple dispatch for converting between packed and original data types during DB usage

function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: FunctorInferenceType}
  # println("convert(::Type{PackedFunctionNodeData{$P}}, d::FunctionNodeData{$T})")
  @error("convert GenericWrapParam is deprecated, use CommonConvWrapper instead.")
  mhstr = packmultihypo(d.fnc)
  return PackedFunctionNodeData(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          string(d.frommodule), convert(P, d.fnc.usrfnc!), mhstr)
end

function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: ConvolutionObject}
  mhstr = packmultihypo(d.fnc)  # this is where certainhypo error occurs
  return PackedFunctionNodeData(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          string(d.frommodule), convert(P, d.fnc.usrfnc!),
          mhstr, d.fnc.certainhypo )  # extract two values from ccw for storage -- ccw thrown away
end

## unpack converters------------------------------------------------------------

function convert(
            ::Type{IncrementalInference.GenericFunctionNodeData{IncrementalInference.CommonConvWrapper{F},Symbol}},
            d::IncrementalInference.GenericFunctionNodeData{P,String} ) where {F <: FunctorInferenceType, P <: PackedInferenceType}
  #
  # TODO store threadmodel=MutliThreaded,SingleThreaded in persistence layer
  usrfnc = convert(F, d.fnc)
  mhcat = parsemultihypostr(d.multihypo)

  # TODO -- improve prepgenericconvolution for hypotheses and certainhypo field recovery when deserializing
  # reconstitute from stored data
  ccw = prepgenericconvolution(DFG.DFGVariable[], usrfnc, multihypo=mhcat)
  ccw.certainhypo = d.certainhypo

  ret = FunctionNodeData{CommonConvWrapper{typeof(usrfnc)}}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          Symbol(d.frommodule), ccw)
  # error("what what $(ret.fnc.certainhypo)")
  return ret
end

function convert(::Type{PT}, ::T) where {PT <: PackedInferenceType, T <:FunctorInferenceType}
  getfield(T.name.module, Symbol("Packed$(T.name.name)"))
end
function convert(::Type{T}, ::PT) where {T <: FunctorInferenceType, PT <: PackedInferenceType}
  getfield(PT.name.module, Symbol(string(PT.name.name)[7:end]))
end

# function getmodule(t::T) where T
#   T.name.module
# end
# function getname(t::T) where T
#   T.name.name
# end
# function getpackedtype(typestring::AS) where {AS <: AbstractString}
#   # println("Caesar.getpackedtype($(typestring))")
#   eval(Meta.parse(typestring))() # TODO consider caching or better
# end
# function encodePackedType(topackdata::VariableNodeData)
#   @warn "'encodePackedType' Deprecated..."
#   # error("IncrementalInference.encodePackedType(::VariableNodeData): Unknown packed type encoding of $(topackdata)")
#   pack(nothing, topackdata)
# end
# function encodePackedType(topackdata::GenericFunctionNodeData{CommonConvWrapper{T}, Symbol}) where {T <: FunctorInferenceType}
#   # println("IncrementalInference.encodePackedType(::GenericFunctionNodeData{T,Symbol}): Unknown packed type encoding T=$(T) of $(topackdata)")
#   fnctype = getfnctype(topackdata)
#   fnc = getfield(getmodule(fnctype), Symbol("Packed$(getname(fnctype))"))
#   convert(PackedFunctionNodeData{fnc}, topackdata)
# end
# function encodePackedType(topackdata::GenericFunctionNodeData{T, <:AbstractString}) where {T <: PackedInferenceType}
#   error("IncrementalInference.encodePackedType(::FunctionNodeData{T, <:AbstractString}): Unknown packed type encoding T=$(T) of $(topackdata)")
#   # @show T, typeof(topackdata)
#   # @warn "Yes, its packed!"
#   # fnctype = getfnctype(topackdata)
#   # @show fnc = getfield(getmodule(fnctype), Symbol("Packed$(getname(fnctype))"))
#   # convert(PackedFunctionNodeData{T}, topackdata)
#   topackdata
# end

"""
    $(SIGNATURES)

Encode complicated function node type to related 'Packed<type>' format assuming a user supplied convert function .
"""
function convert2packedfunctionnode(fgl::G,
                                    fsym::Symbol ) where G <: AbstractDFG
  #
  # fid = fgl.fIDs[fsym]
  fnc = getfnctype(fgl, fsym)
  usrtyp = convert(PackedInferenceType, fnc)
  cfnd = convert(PackedFunctionNodeData{usrtyp}, solverData(getFactor(fgl, fsym)) )
  return cfnd, usrtyp
end

# Variables

function decodePackedType(dfg::G, packeddata::PackedVariableNodeData) where G <: AbstractDFG
  @warn "Deprecated?"
  convert(IncrementalInference.VariableNodeData, packeddata)
end
# Factors
function decodePackedType(dfg::G, packeddata::GenericFunctionNodeData{PT,<:AbstractString}) where {PT, G <: AbstractDFG}
  usrtyp = convert(FunctorInferenceType, packeddata.fnc)
  fulltype = FunctionNodeData{CommonConvWrapper{usrtyp}}
  factor = convert(fulltype, packeddata)
  return factor
end

"""
  $(SIGNATURES)
After deserializing a factor using decodePackedType, use this to
completely rebuild the factor's CCW and user data.
"""
function rebuildFactorMetadata!(dfg::G, factor::DFGFactor)::DFGFactor where G <: AbstractDFG
  # Set up the neighbor data
  neighbors = map(vId->getVariable(dfg, vId), getNeighbors(dfg, factor))
  neighborUserData = map(v->solverData(v).softtype, neighbors)

  # Rebuilding the CCW
  setDefaultFactorNode!(dfg, factor, neighbors, factor.data.fnc.usrfnc!)

  #... Copying neighbor data into the factor?
  ccw_new = solverData(factor)
  for i in 1:Threads.nthreads()
    ccw_new.fnc.cpt[i].factormetadata.variableuserdata = deepcopy(neighborUserData)
  end

  # Copying all other fields in the factor
  # TODO: Confirm whether we still need to do this?
  ## Rebuild solverData(fcnode).fncargvID, however, the list is order sensitive
  # out_neighbors does not gaurantee ordering -- i.e. why is it not being saved
  # for field in fieldnames(typeof(ccw_jld))
  #   if field != :fnc
  #     setfield!(ccw_new, field, getfield(ccw_jld, field))
  #   end
  # end

  return factor
end

"""
    $(SIGNATURES)

Make a full memory copy of the graph and encode all composite function node
types -- assuming that convert methods for 'Packed<type>' formats exist.  The same converters
are used for database persistence with CloudGraphs.jl.
"""
function encodefg(fgl::G ) where G <: AbstractDFG
  #
  fgs = deepcopy(fgl)
  # fgs.g = Graphs.incdict(Graphs.ExVertex,is_directed=false)

  # @showprogress 1 "Encoding variables..."
  for vsym in getVariableIds(fgl)
    # cpvert = deepcopy(  )
    var = getVariable(fgl, vsym)
    # addVariable!(fgs, cpvert)
  end

  # @showprogress 1 "Encoding factors..."
  for (fsym,fid) in fgs.fIDs
    data,ftyp = convert2packedfunctionnode(fgl, fsym)
    data = FunctionNodeData{ftyp}(Int[], false, false, Int[], m, gwpf)
    # newvert = ExVertex(fid,string(fsym))
    # for (key,val) in getVert(fgl,fid,api=api).attributes
    #   newvert.attributes[key] = val
    # end
    ## losing fgl.fncargvID before setdata
    # setData!(newvert, data)
    # api.addvertex!(fgs, newvert)
  end
  fgs.g.inclist = typeof(fgl.g.inclist)()

  # iterated over all edges
  # @showprogress 1 "Encoding edges..."
  for (eid, edges) in fgl.g.inclist
    fgs.g.inclist[eid] = Vector{typeof(edges[1])}()
    for ed in edges
      newed = Graphs.Edge(ed.index,
          fgs.g.vertices[ed.source.index],
          fgs.g.vertices[ed.target.index]  )
      push!(fgs.g.inclist[eid], newed)
    end
  end

  return fgs
end

# import IncrementalInference: decodefg, loadjld

veeCategorical(val::Categorical) = val.p
veeCategorical(val::Union{Nothing, Vector{Float64}}) = val

#


text:
import Dolo

errs = Dolo.lint("LAMP_2s.yaml")

for err in errs
    println(err)
end


text:
#####################
# Upward Flow circuits
#####################

abstract type UpFlowΔNode{O,F} <: DecoratorΔNode{O} end
abstract type UpFlowLeafNode{O,F} <: UpFlowΔNode{O,F} end
abstract type UpFlowInnerNode{O,F} <: UpFlowΔNode{O,F} end

struct UpFlowLiteral{O,F} <: UpFlowLeafNode{O,F}
    origin::O
    pr::F
end

struct UpFlowConstant{O,F} <: UpFlowLeafNode{O,F}
    origin::O
    pr::F
end

abstract type UpFlow⋀{O,F} <: UpFlowInnerNode{O,F} end
abstract type UpFlow⋁{O,F} <: UpFlowInnerNode{O,F} end

# use this version of UpFlow⋀ when ⋀ gates are unique to their parent ⋁ to save space (compact=true)
struct UpFlow⋀Compact{O,F} <: UpFlow⋀{O,F}
    origin::O
    children::Vector{<:UpFlowΔNode{<:O,F}}
    cached_pr_factors::Vector{F}
end

# use this version of UpFlow⋀ when ⋀ gates can have multiple parents to save computation (compact=false)
struct UpFlow⋀Cached{O,F} <: UpFlow⋀{O,F}
    origin::O
    children::Vector{<:UpFlowΔNode{<:O,F}}
    pr::F
end

# this version of UpFlow⋁ is efficient for arity-1 gates
struct UpFlow⋁Compact{O,F} <: UpFlow⋁{O,F}
    origin::O
    child::UpFlowΔNode{<:O,F}
    cached_pr_factors::Vector{F}
end

struct UpFlow⋁Cached{O,F} <: UpFlow⋁{O,F}
    origin::O
    children::Vector{<:UpFlowΔNode{<:O,F}}
    pr::F
end

const UpFlowΔ{O,F} = AbstractVector{<:UpFlowΔNode{O,F}}

#####################
# traits
#####################

@inline GateType(::Type{<:UpFlowLiteral}) = LiteralGate()
@inline GateType(::Type{<:UpFlowConstant}) = ConstantGate()
@inline GateType(::Type{<:UpFlow⋀}) = ⋀Gate()
@inline GateType(::Type{<:UpFlow⋁}) = ⋁Gate()

const HasPr = Union{UpFlow⋁Cached,UpFlow⋀Cached,UpFlowLeafNode}

#####################
# constructors and conversions
#####################

const flow_opts★ = (max_factors= 2, #only save space using compact nodes when it does not cause too much computation with large numbers of factors. 2 or 3 are good values
                     compact⋀=true,
                     compact⋁=true)

"""Construct a upward flow circuit from a given other circuit"""
function UpFlowΔ(circuit::Δ, m::Int, ::Type{El}, opts = flow_opts★)  where El
    # TODO get rid of the arguments above, they are mostly useless

    O = grapheltype(circuit) # type of node in the origin
    F = (El == Bool) ? BitVector : Vector{El}
    fmem  = () -> some_vector(El, m) # note: fmem's return type will determine type of all UpFlows in the circuit (should be El)

    cache = Dict{ΔNode, UpFlowΔNode}()
    sizehint!(cache, length(circuit)*4÷3)

    upflow_node(::LiteralGate, n::ΔNode) = UpFlowLiteral{O,F}(n, fmem())
    upflow_node(::ConstantGate, n::ΔNode) = UpFlowConstant{O,F}(n, fmem())

    upflow_node(::⋀Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        @assert length(children)>=2
        pr_fs = pr_factors(children)
        if opts.compact⋀ && length(pr_fs) <= opts.max_factors
            UpFlow⋀Compact{O,F}(n, children, pr_fs)
        else
            UpFlow⋀Cached{O,F}(n, children, fmem())
        end
    end

    upflow_node(::⋁Gate, n::ΔNode) = begin
        children = map(c -> cache[c], n.children)
        @assert length(children)>=1
        pr_fs = pr_factors(children)
        if opts.compact⋁ && length(children)==1 && length(pr_fs) <= opts.max_factors
            UpFlow⋁Compact{O,F}(n, children[1], pr_fs)
        else
            UpFlow⋁Cached{O,F}(n, children, fmem())
        end
    end

    map(circuit) do node
        fnode = upflow_node(GateType(node), node)
        cache[node] = fnode
        fnode
    end
end

#####################
# methods
#####################

@inline literal(n::UpFlowLiteral)::Lit  = literal(n.origin)
@inline constant(n::UpFlowConstant)::Bool = constant(n.origin)
@inline children(n::UpFlow⋁Compact) = [n.child]
@inline children(n::Union{UpFlow⋀,UpFlow⋁Cached}) = n.children

@inline pr(n::UpFlowΔNode) = prod_fast(pr_factors(n))
@inline pr(n::HasPr) = n.pr

@inline pr_factors(n::UpFlowΔNode) = n.cached_pr_factors
@inline pr_factors(n::HasPr) = [n.pr]
@inline function pr_factors(ns::Vector{<:UpFlowΔNode{O,F}})::Vector{F} where {O,F}
    flatmap(c -> pr_factors(c)::Vector{F}, ns, F[])
end

flow_length(circuit::UpFlowΔ) = length(pr_factors(circuit[end])[1])

#####################

# change the size of the flow vectors throughout this circuit (e.g., change minibatch size)
function resize_flows(circuit::UpFlowΔ{F}, size::Int) where {F}
    # check if resize is needed
    if size != flow_length(circuit)
        for n in circuit
            resize_pr(n,size)
        end
    end
end

@inline resize_pr(n::UpFlowΔNode, size) = () #do nothing
@inline resize_pr(n::HasPr, size) = resize!(n.pr, size)

#####################

function pass_up(circuit::UpFlowΔ{O,F}, data::PlainXData{E}) where {E <: eltype(F)} where {F,O}
    resize_flows(circuit, num_examples(data))
    for n in circuit
        pass_up_node(n, data)
    end
    @assert all(pr(circuit[end]) .≈ one(eltype(F))) "Root probabilities are expected to be 1: $(pr(circuit[end])) ≆ $(one(eltype(F)))"
end

pass_up_node(n::UpFlowΔNode, ::PlainXData) = () #do nothing

function pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}) where {E <: eltype(F)} where {F,O}
    npr = pr(n)
    if positive(n)
        @. npr = $feature_matrix(data)[:,$variable(n)]
    else
        @. npr = @fastmath $one(eltype(F)) - @view $feature_matrix(data)[:,$variable(n)] # using views here strangely seems to save a few allocations and a little time
    end
end

function pass_up_node(n::UpFlowConstant{O,F}, data::PlainXData{E}) where {E <: eltype(F)} where {F,O}
    npr = pr(n)
    if constant(n)
        @. npr = $one(eltype(F))
    else
        @. npr = $zero(eltype(F))
    end
end

function pass_up_node(n::UpFlow⋀Cached, ::PlainXData)
    assign_prod(pr(n), pr_factors(n.children))
end

function pass_up_node(n::UpFlow⋁Cached, ::PlainXData)
    # it is critical to allow fusing of broadcasts here -- this means mapping to child prs outside of the broadcast
    assign_prod(pr(n), pr_factors(n.children[1]))
    for c in n.children[2:end]
        accumulate_prod(pr(n), pr_factors(c))
    end
end


text:
function get_sensors_ring(; num_pts = 180, φ_inicial = 0, φ_final = 2π, kR = 1, θ = π / 2)
    ϕ_range = range(φ_inicial, φ_final; length = num_pts)
    θ_range = ones(length(ϕ_range)) .* θ

    x = kR .* cos.(ϕ_range) .* sin.(θ_range)
    y = kR .* sin.(ϕ_range) .* sin.(θ_range)
    z = kR .* cos.(θ_range)
    return Array(transpose([x y z]))
end

"""
    get_sensors_sphereSurface_Fibonacci(; num_pts::Integer=100, R=5k₀, θ_lims=(0,π))

Function uses the `Fibonacci Sphere Algorithm` to create `num_pts` points over a sphere shell of radius `R`. This is a Legacy code, prefer the Lebedev points.

Returns each axis (in cartesian coordinates) in row: [x y z]

size([x y z]) = (3,num_pts)
"""
function get_sensors_sphereSurface_Fibonacci(;
    num_pts::Integer = 100,
    R = 5k₀,
    θ_lims = (0, π),
)
    indices = 1:num_pts
    θ_min, θ_max = θ_lims

    θ = acos.(range(cos(θ_min), stop = cos(θ_max), length = num_pts))
    ϕ = π * (1 + sqrt(5)) .* indices

    x, y, z = R .* cos.(ϕ) .* sin.(θ), R .* sin.(ϕ) .* sin.(θ), R .* cos.(θ)
    return Array(transpose([x y z]))
end

"""
    get_sensors_sphereSurface(; R=100k₀, leb_order=125)

return xyz, w

---

Returns points in a sphere of radius R and their corresponding weights, following Lebedev quadrature.

The number of points depends on the order, check Lebedev.jl docs.

We used `leb_order = 125` to procude N=5294 points. That was enough for all inhouse tests.
Sensors will be positioned in Far Field, therefore, `R` should be "large" - which dependent on the cloud size.

**Note**: Increase this value if your system size is already ~ R=100k₀.
"""
function get_sensors_sphereSurface(; R = 100k₀, leb_order = 125)
    x, y, z, w = lebedev_by_order(leb_order)

    xyz = zeros(3, length(x))
    for i ∈ eachindex(x)
        _inSpherical = cart2sph([x[i], y[i], z[i]])

        _inSpherical[3] = R * _inSpherical[3] # change sensor radius

        xyz[:, i] = sph2cart(_inSpherical)
    end

    return xyz, w
end


text:
Z = Strang(1)
@test Matrix(Z) == reshape([2.0],(1,1))

n = rand(1:10)
Z = Strang(n)

for i in 1:n, j in 1:n
    i==j && @test Z[i,j] == 2
    abs(i-j)==1 && @test Z[i,j] == -1
    abs(i-j)>1 && @test Z[i,j] == 0
end

A = Strang(10)
u = ones(10)
@test A*u == [1.0;zeros(8);1.0]

#Matvec product
b = randn(n)
@test Z*b ≈ Matrix(Z)*b

m = rand(1:10)
A = randn(m, n)
@test A*Z ≈ A*Matrix(Z)


text:
f=[1,2]
while f[end]<4000000
    append!(f,f[end-1]+f[end])
end
pop!(f)
print(sum(filter(x->(x%2==0),f)))


text:
## Types

# Model parameters (to be estimated)
"""
    HePPCATModel{S<:Number,T<:Real}

Model parameters for probabilistic PCA with noise that is heteroscedastic across samples.
This is the return type of [`heppcat(_)`](@ref), the corresponding estimation function.
# Properties
+ `F  :: Matrix`    factor matrix (computed via `F = U*sqrt(Diagonal(λ))*Vt`)
+ `U  :: Matrix{S}` eigenvectors of factor covariance `F*F'`
+ `λ  :: Vector{T}` eigenvalues of factor covariance `F*F'` (spike eigenvalues)
+ `Vt :: Matrix{S}` (transposed) eigenvectors of `F'*F` (i.e., right singular vectors of `F`)
+ `v  :: Vector{T}` noise variances
"""
struct HePPCATModel{S<:Number,T<:Real}
    U::Matrix{S}   # eigvecs of FF' (factor/spike covariance)
    λ::Vector{T}   # eigvals of FF' (factor/spike covariance)
    Vt::Matrix{S}  # (transposed) eigvecs of F'F (i.e., right singvecs of F)
    v::Vector{T}   # block noise variances
    function HePPCATModel{S,T}(U,λ,Vt,v) where {S<:Number,T<:Real}
        size(U,2) == length(λ) || throw(DimensionMismatch("U has dimensions $(size(U)) but λ has length $(length(λ))"))
        size(Vt,1) == length(λ) || throw(DimensionMismatch("Vt has dimensions $(size(Vt)) but λ has length $(length(λ))"))
        new{S,T}(U,λ,Vt,v)
    end
end
HePPCATModel(U::Matrix{S},λ::Vector{T},Vt::Matrix{S},v::Vector{T}) where {S<:Number,T<:Real} = HePPCATModel{S,T}(U,λ,Vt,v)

"""
    HePPCATModel(U::AbstractMatrix,λ::AbstractVector,Vt::AbstractMatrix,v::AbstractVector)

Construct HePPCATModel object from factor eigenstructure and noise variances.
"""
function HePPCATModel(U::AbstractMatrix,λ::AbstractVector,Vt::AbstractMatrix,v::AbstractVector)
    S = promote_type(eltype(U),eltype(Vt))
    T = promote_type(eltype(λ),eltype(v))
    HePPCATModel(convert(Matrix{S},U),
        convert(Vector{T},λ),
        convert(Matrix{S},Vt),
        convert(Vector{T},v))
end

"""
    HePPCATModel(F::AbstractMatrix,v::AbstractVector)

Construct HePPCATModel object from factor matrix and noise variances.
"""
function HePPCATModel(F::AbstractMatrix,v::AbstractVector)
    U,s,V = svd(F)
    return HePPCATModel(U,s.^2,V',v)
end

Base.:(==)(F::HePPCATModel, G::HePPCATModel) = all(f -> getfield(F, f) == getfield(G, f), 1:nfields(F))

function Base.getproperty(M::HePPCATModel, d::Symbol)
    if d === :F
        U = getfield(M, :U)
        λ = getfield(M, :λ)
        Vt = getfield(M, :Vt)
        return U*sqrt(Diagonal(λ))*Vt
    else
        return getfield(M, d)
    end
end

Base.propertynames(M::HePPCATModel, private::Bool=false) =
    private ? (:F, fieldnames(typeof(M))...) : (:F, :U, :λ, :Vt, :v)

# Update methods
"""
    RootFinding

Root-finding (global maximization) update method.
No fields.
"""
struct RootFinding end

"""
    ExpectationMaximization

Expectation maximization update method.
No fields.
"""
struct ExpectationMaximization end

"""
    MinorizeMaximize

Minorize maximize update method.
No fields.
"""
struct MinorizeMaximize end

"""
    ProjectedGradientAscent{T}

Projected gradient ascent update method.
One field: `stepsize::T`.
"""
struct ProjectedGradientAscent{T}
    stepsize::T
end

"""
    StiefelGradientAscent{T}

Stiefel gradient ascent update method.
One field: `stepsize::T`.
"""
struct StiefelGradientAscent{T}
    stepsize::T
end

"""
    DifferenceOfConcave

Difference of concave update method.
No fields.
"""
struct DifferenceOfConcave end

"""
    QuadraticSolvableMinorizer

Minorize maximize update with quadratic solvable minorizer.
No fields.
"""
struct QuadraticSolvableMinorizer end

"""
    CubicSolvableMinorizer

Minorize maximize update with cubic solvable minorizer.
No fields.
"""
struct CubicSolvableMinorizer end

"""
    QuadraticMinorizer

Minorize maximize update using quadratic minorizer.
No fields.
"""
struct QuadraticMinorizer end

"""
    OptimalQuadraticMinorizer

Minorize maximize update using quadratic minorizer with optimized curvature.
No fields.
"""
struct OptimalQuadraticMinorizer end

# Step sizes
"""
    InverseLipschitz{T}

Inverse Lipschitz step size using bound function specified by field `bound::T`.
Default choice is `bound=LipBoundU2`.
"""
struct InverseLipschitz{T}
    bound::T
end
InverseLipschitz() = InverseLipschitz(LipBoundU2)   # default bound

"""
    ArmijoSearch{S<:Integer,T<:AbstractFloat}

Armijo line search with parameters:
+ `maxsearches :: S`  maximum number of line searches
+ `stepsize    :: T`  initial stepsize
+ `contraction :: T`  contraction factor
+ `tol         :: T`  tolerance for sufficient decrease
"""
struct ArmijoSearch{S<:Integer,T<:AbstractFloat}
    maxsearches::S  # maximum number of line searches
    stepsize::T     # initial stepsize
    contraction::T  # contraction factor
    tol::T          # tolerance for sufficient decrease
end