text:
abstract type ByteStreamTestStep end

abstract type ByteStreamExpectation <: ByteStreamTestStep end

abstract type ByteStreamAction <: ByteStreamTestStep end

struct ByteStreamExpectationViolation <: Exception
    msg::String
end

function property(::Type{ByteStreamExpectationViolation}, property_name, expected, actual)
    "The ByteStream should have had " * property_name * " equal to " *
    string(expected) * " but instead it was " * string(actual)
end

struct EndInput <: ByteStreamAction end

description(ei::EndInput) = "end input"

execute(::EndInput, bs::ByteStream)::Nothing = JLSponge.end_input!(bs)

mutable struct Write <: ByteStreamAction
    data::String
    bytes_written::Union{UInt, Nothing}
    Write(data::String) = new(data, nothing)
end

with_bytes_written!(wt::Write, bytes_written) = (wt.bytes_written = bytes_written; wt)

description(wt::Write) = "write \"" + wt.data + "\" to the stream"

function execute(wt::Write, bs::ByteStream)::Nothing
    bytes_written = JLSponge.write!(bs, wt.data)
    if (wt.bytes_written !== nothing && bytes_written != wt.bytes_written)
        throw(property(ByteStreamExpectationViolation, "bytes_written", wt.bytes_written, bytes_written))
    end
    nothing
end

struct Pop <: ByteStreamAction
    len::UInt
end

description(pop::Pop)::String = "pop " * string(pop.len)

execute(pop::Pop, bs::ByteStream)::Nothing = JLSponge.pop_output!(bs, pop.len)

struct InputEnded <: ByteStreamExpectation
    input_ended::Bool
end

description(inp_ed::InputEnded) = "input_ended: " + string(inp_ed.input_ended)

function execute(inp_ed::InputEnded, bs::ByteStream)::Nothing
    input_ended = JLSponge.input_ended(bs)
    (input_ended == inp_ed.input_ended) ||
        throw(property(ByteStreamExpectationViolation, "input_ended", inp_ed.input_ended, input_ended))
    nothing
end

struct BufferEmpty <: ByteStreamExpectation
    buffer_empty::Bool
end

description(bfe::BufferEmpty) = "buffer_empty: " * string(bfe.buffer_empty);

function execute(bfe::BufferEmpty, bs::ByteStream)::Nothing
    buffer_empty = JLSponge.buffer_empty(bs)
    (bfe.buffer_empty == buffer_empty) ||
        throw(property(ByteStreamExpectationViolation, "buffer_empty", bfe.buffer_empty, buffer_empty))
    nothing
end

struct Eof <: ByteStreamExpectation
    eof::Bool
end

description(e::Eof) = "eof: " * string(e.eof)

function execute(e::Eof, bs::ByteStream)::Nothing
    e1 = JLSponge.eof(bs)
    (e.eof == e1) ||
        throw(property(ByteStreamExpectationViolation, "eof", e.eof, e1))
    nothing
end

struct BufferSize <: ByteStreamExpectation
    buffer_size::UInt
end

description(bfsize::BufferSize) = "buffer_size: " * string(bfsize.buffer_size)

function execute(bfsize::BufferSize, bs::ByteStream)::Nothing
    buffer_size = JLSponge.buffer_size(bs)
    (bfsize.buffer_size == buffer_size) ||
        throw(property(ByteStreamExpectationViolation, "buffer_size", bfsize.buffer_size, buffer_size))
    nothing
end

struct RemainingCapacity <: ByteStreamExpectation
    remaining_capacity::UInt
end

description(rcap::RemainingCapacity) = "remaining_capacity: " * string(rcap.remaining_capacity)

function execute(rcap::RemainingCapacity, bs::ByteStream)::Nothing
    remaining_capacity = JLSponge.remain_cap(bs)
    (rcap.remaining_capacity == remaining_capacity) ||
        throw(property(ByteStreamExpectationViolation, "remaining_capacity", rcap.remaining_capacity, remaining_capacity))
    nothing
end

struct BytesWritten <: ByteStreamExpectation
    bytes_written::UInt
end

description(bwrt::BytesWritten) = "bytes_written: " * string(bwrt.bytes_written)

function execute(bwrt::BytesWritten, bs::ByteStream)::Nothing
    bytes_written = bs.bytes_written_count
    (bwrt.bytes_written == bytes_written) ||
        throw(property(ByteStreamExpectationViolation, "bytes_written", bwrt.bytes_written, bytes_written))
    nothing
end

struct BytesRead <: ByteStreamExpectation
    bytes_read::UInt
end

description(brd::BytesRead) = "bytes_read: " * string(brd.bytes_read)

function execute(brd::BytesRead, bs::ByteStream)::Nothing
    bytes_read = bs.bytes_read_count
    (brd.bytes_read == bytes_read) ||
        throw(property(ByteStreamExpectationViolation, "bytes_read", brd.bytes_read, bytes_read))
    nothing
end

struct Peek <: ByteStreamExpectation
    output::String
end

description(pk::Peek) = "\"" * pk.output * "\" at the front of the stream";

function execute(pk::Peek, bs::ByteStream)::Nothing
    output = JLSponge.peek_out(bs, pk.output |> length)
    (pk.output == output) ||
        throw(ByteStreamExpectationViolation("Expected \"" * pk.output * "\" at the front of the stream, but found \"" * output * "\""))
    nothing
end

struct ByteStreamTestHarness
    test_name::String
    byte_stream::ByteStream
    steps_executed::Vector{String}
    ByteStreamTestHarness(test_name, capacity) =
        new(test_name, ByteStream(capacity), ["Initialized with (capacity=$capacity)"])
end

function execute(bstest::ByteStreamTestHarness, step::ByteStreamTestStep)
    try
        execute(step, bstest.byte_stream)
        push!(bstest.steps_executed, step |> string)
    catch e
        print(stderr, "Test Failure on expectation:\n\t $(string(step))")
        if isa(e, ByteStreamExpectationViolation)
            print(stderr, "\n\nFailure message:\n\t$(e.msg)")
            new_exception = "The test \"" * bstest.test_name * "\" failed"
        else
            print(stderr, "\n\nException:\n\t$(e)")
            new_exception = "The test \"" * bstest.test_name * "\" caused your implementation to throw an exception!"
        end
        print(stderr, "\n\nList of steps that executed successfully:")
        for s in bstest.steps_executed
            print(stderr, "\n\t", s)
        end
        print(stderr, "\n\n")
        flush(stderr)
        throw(ByteStreamExpectationViolation(new_exception))
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module Stepwise

    using DataFrames
    using StatsBase
    using ..DataMod, ..Utils

    export stepwise

    # The code for `step` and `stepwise` below is inspired by the code by Bogumił Kamiński, found at
    # https://stackoverflow.com/questions/49794476/backward-elimination-forward-selection-in-multilinear-regression-in-julia

    function multibic(data::Data, rhs::Vector{Symbol})::Tuple{Float64, Vector{RegressionModel}}
        models = get_models(data, rhs)::Vector{RegressionModel}
        return (mean(bic.(models))::Float64, models)
    end
    function multibic(models::Vector{RegressionModel})::Float64
        return mean(bic.(models))
    end

    function step(data::Data, rhs::Vector{Symbol}, forward::Bool,
                  prevbest::Float64, prevmodel::Vector{RegressionModel},
                  verbose::Bool, debug::Bool)::Tuple{Vector{Symbol}, Bool, Float64, Vector{RegressionModel}}
        options = forward ? setdiff(data.xs, rhs) : rhs
        if isempty(options)
            return (rhs, false, prevbest, prevmodel)
        end
        bestbic = prevbest
        bestrhs = rhs
        bestmodel = prevmodel
        improved = false
        for opt in options
            thisrhs = forward ? [rhs; opt] : setdiff(rhs, opt)
            thisbic, thismodel = multibic(data, thisrhs)
            if debug
                print(opt); print(" - "); println(thisbic)
            end
            if thisbic < bestbic
                bestbic = thisbic
                bestrhs = thisrhs
                bestmodel = thismodel
                improved = true
            end
        end
        if verbose
            if improved
                println("$(forward ? "forward" : "backward") step - new rhs: $bestrhs\nwith mean BIC $bestbic")
            else
                println("$(forward ? "forward" : "backward") step could not improve BIC")
            end
        end
        return (bestrhs, improved, bestbic, bestmodel)
    end

    function stepwise(data::Data; verbose::Bool = true, debug::Bool = false)::Vector{RegressionModel}
        rhs = Symbol[]
        model = get_models(data, Symbol[])
        bestbic = multibic(model)
        verbose && println("BIC empty model: $bestbic")
        while true
            verbose && println("step")
            rhs, improvedfwd, bestbic, model = step(data, rhs, true, bestbic, model, verbose, debug)
            rhs, improvedbwd, bestbic, model = step(data, rhs, false, bestbic, model, verbose, debug)
            if !improvedfwd && !improvedbwd
                sort!(rhs)
                verbose && println("final rhs: $rhs")
                return model
            end
        end
    end

end

using .Stepwise
export stepwise

ext: jl

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
"""
    Quadratic(trans::AffineMap, length, ay, by, cy, az, bz, cz)

Construct a `Quadratic` object. The `AffineMap` brings the quadratic into a
frame in which the `x` axis is approximately aligned with the quadratic. In this frame,
the `x` values of the quadratic lie between zero and `length`, and the `y` and `z`
components are parameterized by `x` as:

y = ay + by*x + cy*x^2
z = az + bz*x + cz*x^2
"""
struct Quadratic{Tfm <: AffineTransformation, T <: Real}
    transform::Tfm # transformation from global frame to quadratic's frame
    length::T  # replace with xmin, xmax?
    ay::T
    by::T
    cy::T
    az::T
    bz::T
    cz::T
end

"""
    Quadratic(p_start, p_mid, p_end)

Takes three coordinates for the start, somewhere in the middle, and the end of a quadratic
and constructs the appropriate `Quadratic` object.
"""
function Quadratic(p_start::StaticVector{3,T}, p_mid::StaticVector{3,T}, p_end::StaticVector{3,T}) where {T <: Real}
    dir = SVector(p_end[1] - p_start[1], p_end[2] - p_start[2])
    len = norm(dir) # horizontal length
    @static if VERSION < v"0.7"
        θ = atan2(dir[2], dir[1]) # rotation in x-y plane
    else
        θ = atan(dir[2], dir[1]) # rotation in x-y plane
    end
    tform = inv(AffineMap(RotZ(θ), p_start))

    p2_mid = tform(p_mid)
    p2_end = tform(p_end)

    x1 = p2_mid[1]
    y1 = p2_mid[2]
    z1 = p2_mid[3]

    x2 = p2_end[1]
    y2 = p2_end[2]
    z2 = p2_end[3]

    by = (x1*y2/x2 - x2*y1/x1)/(x1 - x2)
    bz = (x1*z2/x2 - x2*z1/x1)/(x1 - x2)

    cy = (y2/x2 - y1/x1)/(x2 - x1)
    cz = (z2/x2 - z1/x1)/(x2 - x1)

    ay = zero(cy)
    az = zero(cz)

    return Quadratic(tform, len, ay, by, cy, az, bz, cz)
end

# mop up other types
function Quadratic(p1::AbstractVector{T}, p2::AbstractVector{T}, p3::AbstractVector{T}) where {T <: Real}
    Quadratic(SVector{3,T}(p1), SVector{3,T}(p2), SVector{3,T}(p3))
end

eltype(::Quadratic{T}) where {T} = SVector{3,T}
eltype(::Type{<:Quadratic{T}}) where {T} = SVector{3,T}

# Use getindex to get points along the line
function getindex(q::Quadratic, x)
    inv(q.transform)(SVector(x, q.ay+q.by*x+q.cy*x*x, q.az+q.bz*x+q.cz*x*x))
end
length(q::Quadratic) = q.length  # matches getindex and makes sense...
if VERSION < v"0.7"
    endof(q::Quadratic) = length(q) # e.g. q[end - 1] is the point 1 metre from top of quadratic
else
    lastindex(q::Quadratic) = length(q) # e.g. q[end - 1] is the point 1 metre from top of quadratic
end

getindex(q::Quadratic, v::AbstractVector) = map(x -> q[x], v)

# Transformations (which compose into AffineMap's)
# Note: Julia's dispatch system makes it difficult to support all `Transformations`.
function (trans::Translation{V})(q::Quadratic) where {V}
    Quadratic(q.transform ∘ inv(trans), q.length, q.ay, q.by, q.cy, q.az, q.bz, q.cz)
end

function (trans::LinearMap{M})(q::Quadratic) where {M}
    Quadratic(q.transform ∘ inv(trans), q.length, q.ay, q.by, q.cy, q.az, q.bz, q.cz)
end

function (trans::AffineMap{M,V})(q::Quadratic) where {M,V}
    Quadratic(q.transform ∘ inv(trans), q.length, q.ay, q.by, q.cy, q.az, q.bz, q.cz)
end


text:
# Miscellaneous and alternative algorithms to compare with the primary
# implementations. Not exported.
using Polynomials

struct Poehn1985 <: NeXLAlgorithm end

"""
    jumpratio(ashell::AtomicSubShell, ::Type{Poehn1985})

An implement of jump ratios attributed to
  * Poehn, Wernisch, Hanke (1985) X-ray Spectrom 14(3):120, 1985

Compares reasonably over available range.
"""
function jumpratio(ashell::AtomicSubShell, ::Type{Poehn1985})
    zz = ashell.z
    if (zz in 12:49) && (ashell.subshell.index == 1)
        return Polynomial( [17.54, -0.6608, 0.01427, -1.1e-4 ])(Float64(zz))
    elseif zz in 40:82
        if ashell.subshell.index == 4
            return Polynomial( [20.03, -0.7732, 0.01159, -5.835e-5])(Float64(zz))
        elseif ashell.subshell.index == 3
            return 1.41
        elseif ashell.subshell.index == 2
            return 1.16
        end
    end
    error("Unsupported edge: $(ashell)")
end

"""
    klinewidths(elm::Element)

Linewidth of the K shell according to Bambynek'1974 errata to Bambynek 1972.
Shown to be reliable for Z>36 or so.
"""
klinewidths(elm::Element) = 1.73e-6 * z(elm)^3.93

struct Burhop1965 <: NeXLAlgorithm end
struct Sogut2002 <: NeXLAlgorithm end
struct Krause1979 <: NeXLAlgorithm end
struct Kahoul2012 <: NeXLAlgorithm end
struct Reed1975ω <: NeXLAlgorithm end

let krause1979_data
    function loadKrause1979()
        nm(v) = ismissing(v) ? 0.0 : v
        for _ in length(krause1979_data)+1:4
            push!(krause1979_data, ( 0.0, 0.0, 0.0, 0.0 ))
        end
        for row in CSV.File(joinpath(dirname(pathof(@__MODULE__)), "..", "data", "krause1979.csv"))
            # "Z","Sy","wK","wL1","wL2","wL3","f1","f1_2","f1_3","fp1_3","f2_3"
            push!(krause1979_data, (nm(row.wK), nm(row.wL1), nm(row.wL2), nm(row.wL3) ))
        end
    end
    krause1979_data = NTuple{4, Float64}[]

"""
    fluorescenceyield(ass::AtomicSubShell, ::Type{Krause1979})

An alternative for K and L-line yields.  Agrees well with others.
"""
    global function fluorescenceyield(ass::AtomicSubShell, ::Type{Krause1979})
        isempty(krause1979_data) && loadKrause1979()
        if ass.subshell.index<=4
            krause1979_data[ass.z][ass.subshell.index]
        else
            error("Krause 1979 only supports K through L3 sub-shells.")
        end
    end
end

"""
    fluorescenceyield(ass::AtomicSubShell, ::Type{Sogut2002})

An alternative for M-line yields.  Not alway reasonable.
"""
function fluorescenceyield(ass::AtomicSubShell, ::Type{Sogut2002})
    if (ass.subshell.index in 5:9) && (ass.z in 20:90)
        if ass.subshell.index==5 # Seems reasonable
            if ass.z <= 49
                return Polynomial([ 0.0005, -0.00004, 6.3425e-7])(Float64(ass.z))
            elseif ass.z <= 79
                return Polynomial([ -0.00254, 0.00006])(Float64(ass.z))
            else
                return Polynomial([ -0.01485, 0.00022 ])(Float64(ass.z))
            end
        elseif ass.subshell.index==6 # Ok for Z<=56
            if ass.z <= 36
                return Polynomial([ 0.00002, -5.66510e-24])(Float64(ass.z))
            elseif ass.z <= 56
                return Polynomial([ -0.00221, 0.00006])(Float64(ass.z))
            else
                # Seems wrong!
                return Polynomial([ 0.07255, -0.00225, 0.00002])(Float64(ass.z))
            end
        elseif ass.subshell.index==7 # Seems low
            if ass.z >= 35
                if ass.z <= 58
                    return Polynomial([ -0.00336, 0.00024, -6.2005e-6, 5.6949e-8 ])(Float64(ass.z))
                elseif ass.z <= 72
                    return Polynomial([ -0.00249, 0.00006, -7.085e-20 ])(Float64(ass.z))
                else
                    return Polynomial([ -0.0226, 0.00034 ])(Float64(ass.z))
                end
            end
        elseif ass.subshell.index==8 # Seems reasonable
            if ass.z >= 32
                if ass.z <= 60
                    return Polynomial([ 0.0027, -3.801710e-21 ])(Float64(ass.z))
                else
                    return Polynomial([ 0.27987, -0.00872, 0.00007] )(Float64(ass.z))
                end
            end
        elseif ass.subshell.index==9 # Seems reasonable above 60
            return ass.z >= 60 ? Polynomial([ -0.08517, 0.00144 ])(Float64(ass.z)) : 0.0
        end
        return 0.0
    end
    error("Sogut 2002 only implements M-lines from Ca to Th.")
end

"""
    fluorescenceyield(ass::AtomicSubShell, ::Type{Burhop1965})

An approximate expression for the K-shell fluorescence yield due to
E.H.S Burhop, J. Phys. Radium, 16, 625 (1965). Seems reasonable.
"""
function fluorescenceyield(ass::AtomicSubShell, ::Type{Burhop1965})
    if ass.subshell.index==1
        d = Polynomial( [-0.044, 0.0346, 0.0, -1.35e-6])(Float64(ass.z))
        return d^4 / (1.0 + d^4)
    else
        error("Burhop 1965 only implements the K-shell.")
    end
end

fluorescenceyield(z::Int) = fluorescenceyield(z, Burhop1965)

"""
    Kahoul 2012 expression for the K-shell fluorescence yield
"""
function fluorescenceyield(ass::AtomicSubShell, ::Type{Kahoul2012} )
    if ass.subshell.index==1
        p = if ass.z in 11:20
            Polynomial([ 1.601637, -0.298757, 0.022417, -4.92221e-4 ])
        elseif ass.z in 21:50
            Polynomial([ -0.17511, 0.05177, -6.4788e-4, 6.25318e-6 ])
        elseif ass.z in 51:99
            Polynomial([ 0.59013, 0.02214, -2.90802e-5 ])
        end
        p4=p(Float64(ass.z))^4
        return p4/(1.0+p4)
    else
        error("Kahoul 2012 only implements K-shell fluorescence yields.")
    end
end

function fluorescenceyield(ass::AtomicSubShell,::Type{Reed1975ω})
    if ass.subshell.index==1
        return Float64(ass.z)^4/(1.0e6+Float64(ass.z)^4)
    else
        error("Reed 1975 only implements K-shell fluorescence yields.")
    end
end


text:
#!/usr/bin/julia

# Copyright (C) 2021 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Twenty, Trench Map."

include("Utils.jl")

using .Utils

"""Apply the <pattern> filter to the <image> <steps> times, return the number
of lit pixels in the final image."""
function solve(pattern :: BitArray, image :: BitArray{2}, steps :: Int) :: Int
    old_size = size(image)
    new_size = old_size .+ 4steps

    new_image = zeros(Bool, new_size)
    # Insert the old image in the new pixel grid.
    new_image[(1:old_size[1]) .+ 2steps, (1:old_size[2]) .+ 2steps] = image

    for step in 1:steps
        tmp_image = similar(new_image)

        inds = CartesianIndices(new_image)
        start, final = first(inds), last(inds)

        @views for i in inds
            start_row, start_col = max(start, i - CartesianIndex(1,1)).I
            final_row, final_col = min(final, i + CartesianIndex(1,1)).I

            neighbors = new_image[start_row:final_row, start_col:final_col]
            bit = bin2dec(vcat(eachrow(neighbors)...)) + 1

            tmp_image[i] = pattern[bit]
        end

        #= “Clean” the borders from generated “artifacts” due to the infinite
         = size of the grid, do this in even turns for optimality.
         =#
        if iseven(step)
            @. (tmp_image[:,1] = tmp_image[1,:]
              = tmp_image[:,end] = tmp_image[end,:] = 0)
        end

        new_image = copy(tmp_image)
    end

    return count(new_image)
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readlines(first(ARGS))

    pattern = BitArray(collect(data[1]) .== '#')
    image = BitArray{2}(hcat([i .== '#' for i in collect.(data[3:end])]...)')

    println(solve(pattern, image, 2))
    println(solve(pattern, image, 50))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
using Statistics

"""
    psislw(lw, [wcpp, wtrunc])

Compute the Pareto smoothed importance sampling (PSIS).

# Arguments
* `lw::AbstractArray`: Array of size n x m containing m sets of n log weights. It is also possible to provide one dimensional array of length n.
* `wcpp::Real`: Percentage of samples used for GPD fit estimate (default is 20).
* `wtrunc::Float64`: Positive parameter for truncating very large weights to n^wtrunc. Providing False or 0 disables truncation. Default values is 3/4.

# Returns
* `lw_out::AbstractArray`: Smoothed log weights
* `kss::AbstractArray`: Pareto tail indices
"""
function psislw(lw::AbstractArray, wcpp::Int64=20, wtrunc::Float64=3/4)

    lw_out = copy(lw)

    if ~(1 <= ndims(lw_out) <= 2)
        throw(DimensionMismatch("Argument `lw` must be 1 or 2 dimensional."))
    end
    if size(lw_out,1) <= 1
        error("More than one log-weight needed.")
    end
    return _psislw(lw_out, wcpp, wtrunc)
end

function _psislw(lw_out::Array{Float64}, wcpp::Int64, wtrunc::Float64)

    if ndims(lw_out) == 2
        n, m = size(lw_out)
    elseif ndims(lw_out) == 1
        n = length(lw_out)
        m = 1
    end
    kss = zeros(Float64, m)

    # precalculate constants
    cutoffmin = log(eps(0.5))
    logn = log(n)

    # loop over sets of log weights
    for i in 1:m
        x = view(lw_out, :, i)
        x .-= maximum(x)
        xcutoff = max(quantile(x, 1 - wcpp/100), cutoffmin)
        expxcutoff = exp(xcutoff)
        x2 = x[x .> xcutoff]
        n2 = length(x2)
        if n2 <= 4
            # not enough tail samples for gpdfitnew
            k = Inf
        else
            # order of tail samples
            x2si = sortperm(x2)
            # fit generalized Pareto distribution to the right tail samples
            x2 .= exp.(x2) .- expxcutoff
            k, sigma = gpdfitnew(x2)
            # compute ordered statistic for the fit
            sti = collect(0.5:n2) / n2
            qq = gpinv(sti, k, sigma)
            qq .= log.(qq .+ expxcutoff)
            # place the smoothed tail into the output array
            x[x .> xcutoff] = qq[x2si]
        end
        if wtrunc > 0
            # truncate too large weights
            lwtrunc = wtrunc * logn - logn + logsumexp(x)
            x[x .> lwtrunc] .= lwtrunc
        end
        # renormalize weights
        x .-= logsumexp(x)
        # store tail index k
        kss[i] = k
    end
    # If the provided input array is one dimensional, return kss as scalar.
    if ndims(lw_out) == 1
        kss = kss[1]
    end
    return lw_out, kss
end

export
    psislw


text:
# ---
# title: 768. Max Chunks To Make Sorted II
# id: problem768
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array
# link: <https://leetcode.com/problems/max-chunks-to-make-sorted-ii/description/>
# hidden: true
# ---
#
# _This question is the same as "Max Chunks to Make Sorted" except the integers
# of the given array are not necessarily distinct, the input array could be up
# to length `2000`, and the elements could be up to `10**8`._
#
# * * *
#
# Given an array `arr` of integers ( **not necessarily distinct** ), we split
# the array into some number of "chunks" (partitions), and individually sort
# each chunk.  After concatenating them, the result equals the sorted array.
#
# What is the most number of chunks we could have made?
#
# **Example 1:**
#
#
#
#     Input: arr = [5,4,3,2,1]
#     Output: 1
#     Explanation:
#     Splitting into two or more chunks will not return the required result.
#     For example, splitting into [5, 4], [3, 2, 1] will result in [4, 5, 1, 2, 3], which isn't sorted.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [2,1,3,4,4]
#     Output: 4
#     Explanation:
#     We can split into two chunks, such as [2, 1], [3, 4, 4].
#     However, splitting into [2, 1], [3], [4], [4] is the highest number of chunks possible.
#
#
# **Note:**
#
#   * `arr` will have length in range `[1, 2000]`.
#   * `arr[i]` will be an integer in range `[0, 10**8]`.
#
#
#
#
## @lc code=start
using LeetCode

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


text:
## Use regression to compute the best-fit quadratic model (with a dense Hessian)

const scaledthresh = 100 # FIXME hardcoded

"""
    g, B = solve(Q::QmIGE)

Solve a quadratic model `Q` for the gradient `g` and Hessian `B`. See
[`build_quadratic_model`](@ref) and [`complete_quadratic_model!`](@ref) to build `Q`.
"""
function solve(Q::QmIGE{T,N}) where {T,N}
    # Importantly, this defines the storage order in Q
    gB = LowerTriangular(Q.coefs) \ Q.rhs
    g, B = Vector{T}(undef, N), Matrix{T}(undef, N, N)
    k = 0
    dimpiv = Q.dimpiv
    for i = 1:N
        di = dimpiv[i]
        g[di] = gB[k+=1]
        for j = 1:i
            dj = dimpiv[j]
            B[di,dj] = B[dj,di] = gB[k+=1]
        end
    end
    return g, B
end

# Uses available points to set up the regression problem
"""
    Q = build_quadratic_model(box, x0, scale, thresh)

Start building a quadratic model of the function, using points already stored
in the tree that contains `box`. `scale` (see [`boxscale`](@ref)) is used in part to
ensure numeric stability, and will prevent inclusion of points that are too distant along
any coordinate axis.

`Q` is not guaranteed to be complete; the tree may not contain enough points, or enough
"independent" points, to determine all the parameters of `Q`. See also [`complete_quadratic_model!`](@ref).
"""
function build_quadratic_model(box::Box{T,N}, x0, scale, thresh) where {T,N}
    Q = QmIGE{T,N}()
    c = value(box)
    xbase = position(box, x0)
    xtmp = similar(xbase)
    flag = similar(xtmp, Bool)
    nanfill = fill(T(NaN), N)
    bbs = boxbounds(box, nanfill, nanfill)
    succeeded = true
    @assert(isleaf(box))
    if !isleaf(box)
        for i = 1:3
            succeeded &= descend!(Q, c, box.children[i], x0, xbase, xtmp, flag, bbs, scale, thresh)
        end
    end
    # To increase the numeric stability of the result, we do more boxes than
    # nominally needed (see the `swap` portion of `insert!`)
    need_extra = true
    while succeeded && (Q.nzrows[] < length(Q.rhs) || need_extra) && !isroot(box)
        if Q.nzrows[] == length(Q.rhs)
            need_extra = minimum(abs.(diag(Q.coefs))) < thresh  # because of scaling this doesn't need units
        end
        cindex = box.parent_cindex
        box = box.parent
        succeeded &= all(isfinite, box.fvalues) && !box.qnconverged
        succeeded || break
        j = 1
        if j == cindex j+=1 end
        succeeded &= descend!(Q, c, box.children[j], x0, xbase, xtmp, flag, bbs, scale, thresh)
        succeeded || break
        j += 1
        if j == cindex j+=1 end
        succeeded &= descend!(Q, c, box.children[j], x0, xbase, xtmp, flag, bbs, scale, thresh)
    end
    return Q, xbase, c, succeeded
end

# Evaluates the function at new points until the quadratic model is fully specified
"""
    complete_quadratic_model!(Q::QmIGE, c, box::Box, f::Function, x0, xbase, scale, splits, lower::AbstractVector, upper::AbstractVector, thresh)

Split boxes in the tree containing `box` until all the parameters of the quadratic model can
be determined. It is allowed to fail if one or more boxes are too small to be split further,
so it is advisable to check that all diagonal coefficients in `Q` exceed `thresh` before
proceeding.

See also [`solve`](@ref).
"""
function complete_quadratic_model!(Q::QmIGE, c, box::Box{T,N}, f::Function, x0, xbase, scale, splits, lower::AbstractVector, upper::AbstractVector, thresh) where {T,N}
    @assert(all(x->x>0, Q.dimpiv))
    xtmp  = similar(x0)
    Δx    = similar(x0)
    Δxtmp = similar(x0)
    flag  = similar(x0, Bool)
    dimpiv = Q.dimpiv
    rowoffset = 0
    # Sanity check
    i = 0
    for idim = 1:N
        for irow = 0:idim
            if isassigned(Q.rowbox, i+=1)
                @assert(isleaf(Q.rowbox[i]))
                position!(xtmp, flag, Q.rowbox[i], x0)
                for j = idim+1:N
                    sd = Q.dimpiv[j]
                    @assert(xtmp[sd]==xbase[sd])
                end
            end
        end
    end
    # Process block-by-block, each block corresponding to a "new" splitting dimension
    for idim = 1:N
        all_are_set = true
        for j = 1:idim+1
            all_are_set &= abs(Q.coefs[j+rowoffset,j+rowoffset]) >= thresh
        end
        if all_are_set
            rowoffset += idim + 1
            continue
        end
        rowb = row = rowoffset + idim + 1  # within this dimension's block, go backwards
        irow = idim
        # @show Q.rowbox[rowoffset+1:row]
        # If all boxes that differed along this dimension were out-of-scale, we
        # need to split along that dimension.
        if !isassigned(Q.rowbox, rowb)
            splitdim = Q.dimpiv[idim]
            position!(xtmp, flag, box, x0)
            xcur, s = xbase[splitdim], scale[splitdim]
            p = find_parent_with_splitdim(box, splitdim)
            bb = boxbounds(p)
            @assert xtmp[splitdim] == xcur
            l = max(bb[1], xcur - s/2)
            u = min(bb[2], xcur + s/2)
            if l >= xcur
                l, xcur = xcur, (xcur + u)/2
            elseif xcur >= u
                u, xcur = xcur, (l + xcur)/2
            end
            abs(u-l) < epswidth((l, u)) && return Q #FIXME?
            split!(box, f, xtmp, splitdim, MVector3{T}(l, xcur, u), lower[splitdim], upper[splitdim], xcur, value(box))
            xcur = xbase[splitdim]
            j = 3
            j1 = 1;    if box.xvalues[j1] == xcur j = j1; j1 += 1; end
            j2 = j1+1; if box.xvalues[j2] == xcur j = j2; j2 += 1; end
            nbrbox = box.fvalues[j1] < box.fvalues[j2] ? box.children[j1] : box.children[j2]
            Q.rowbox[rowb] = nbrbox
            box = box.children[j]
        end
        while irow >= 0
            # @show row
            # if row == rowb
            #     rng = rowoffset+1:rowoffset+idim+1
            #     display(Q.coefs[rng,rng])
            #     @show thresh
            # end
            Qd = abs(Q.coefs[row,row])
            if Qd < thresh
                # @show rowb irow
                rbox = Q.rowbox[rowb]
                @assert isleaf(rbox)
                position!(xtmp, flag, rbox, x0)
                Δx .= xtmp .- xbase
                if irow > 0
                    # When irow > 0, splitting along the corresponding dimension will fill
                    # in this missing information
                    splitdim = dimpiv[irow]
                else
                    # root = get_root(box)
                    # splitprint_colored(root, box)
                    # println()
                    # splitprint_colored(root, rbox)
                    # println()
                    # @show Δx dimpiv
                    # irow == 0 corresponds to the coefficient of the gradient
                    if Δx[dimpiv[idim]] == 0
                        # There's no displacement in the gradient coordinate, so obviously
                        # we have to split it
                        splitdim = dimpiv[idim]
                    else
                        # Here it's not so obvious which coordinate we need to split.
                        # To avoid the risk of failure, we check the result of the gaussian
                        # elimination and choose our splitting dimension as the one that
                        # has largest (remaining) diagonal value.
                        Δxtmp[:] .= Δx ./ scale
                        maxdiag, splitdim = zero(T), 0
                        rng = rowoffset+1:rowoffset+idim+1
                        Q.rowzero[row] = true
                        for ii = 1:idim
                            # The precise magnitude of the step doesn't matter for these "probe" trials
                            sd = dimpiv[ii]
                            xold = Δxtmp[sd]
                            Δxtmp[sd] = Δxtmp[sd] == 0 ? 1 : 2*Δxtmp[sd] # ensure a non-zero entry no matter what
                            # @show Δxtmp
                            setrow!(Q.rowtmp, dimpiv, Δxtmp, N, sd)
                            Δxtmp[sd] = xold
                            # @show ii sd
                            rtmp, _ = incremental_gaussian_elimination!(Q.rowtmp, Q.coefs, Q.rowzero, rng; canswap=false, debug=false)
                            thisdiag = abs(Q.rowtmp[row])
                            if thisdiag > maxdiag
                                maxdiag = thisdiag
                                splitdim = sd
                            end
                        end
                        if splitdim == 0
                            rng = rowoffset+1:rowoffset+idim+1
                            display(Q.coefs[rng,rng])
                            error("no direction works")
                        end
                        # @show maxdiag splitdim
                    end
                end
                xcur = xtmp[splitdim]
                p = find_parent_with_splitdim(rbox, splitdim)
                # Check that the box is big enough to split along this dimension
                if !isroot(p)
                    bb = boxbounds(p)
                    bbeps = epswidth(bb)
                    if !(bb[2] - bb[1] > bbeps) || !(bb[2] - xtmp[splitdim] > bbeps) || !(xtmp[splitdim] - bb[1] > bbeps)
                        return Q  # fail (FIXME?)
                    end
                end
                insrows, rboxn = split_insert!(Q, rbox, p, row:rowoffset+idim+1, f, xbase, c, xtmp, splitdim, scale, splits, lower, upper)
                if Qd == 0 && (maximum(insrows) != row || Q.coefs[row,row] == 0)
                    # Debugging
                    println("Something went wrong")
                    @show numevals(f) insrows row irow rowb rowoffset idim Q.dimpiv
                    @show rbox.parent.splitdim box.parent.splitdim
                    @show splitdim boxbounds(rbox, lower, upper)
                    rng = rowoffset+1:rowoffset+idim+1
                    display(Q.coefs[rng,rng])
                    @show xcur rbox.xvalues (rbox.xvalues .- xcur)./scale[splitdim]
                end
                Qd == 0 && @assert(maximum(insrows) == row)
                if row == rowb
                    Q.rowbox[row] = rboxn  # ensure it's a leaf
                end
            end
            if isassigned(Q.rowbox, row) && isleaf(Q.rowbox[row])
                rowb = row
            end
            row -=1; irow -= 1
        end
        rowoffset += idim + 1
    end
    return Q
end

function split_insert!(Q, rbox::Box{T}, p, rng, f, xbase, c, xtmp, splitdim, scale, splits, lower, upper) where T
    debug = false
    debug && @show rng
    row = first(rng)
    xcur = xtmp[splitdim]
    if isroot(p)
        xsplit = MVector3{T}(splits[splitdim])
        @assert(xsplit[1] < xsplit[2] < xsplit[3])
        lwr, upr = lower[splitdim], upper[splitdim]
    else
        bb = boxbounds(p)
        lwr, upr = bb[1], bb[2]
        xcur, scalesd = xtmp[splitdim], scale[splitdim]
        lo, hi = max(lwr, 5*lwr/6 + xcur/6, xcur - scalesd), min(upr, xcur/6 + 5*upr/6, xcur + scalesd)
        if lo + sqrt(eps(T))*scalesd >= xcur
            lo, xcur = xcur, (xcur+hi)/2
        elseif xcur + sqrt(eps(T))*scalesd >= hi
            xcur, hi = (xcur+lo)/2, xcur
        end
        xsplit = MVector3{T}(lo, xcur, hi)
        xcur = xtmp[splitdim]
        @assert(xsplit[1] < xsplit[2] < xsplit[3])
    end
    j1 = 1
    if xsplit[j1] == xcur j1 += 1 end
    j2 = j1+1
    if xsplit[j2] == xcur j2 += 1 end
    l, h = xsplit[j1], xsplit[j2]
    Qd = abs(Q.coefs[row, row])
    if Qd != 0
        # This is a case where the existing entry is below thresh. Before proceeding, check
        # whether the new one would be an improvement---rbox might be too small along
        # splitdim to raise the value.
        xtmp[splitdim] = l
        setrow!(Q.rowtmp, Q.dimpiv, (xtmp-xbase)./scale, ndims(rbox), splitdim)
        incremental_gaussian_elimination!(Q.rowtmp, Q.coefs, Q.rowzero, rng; canswap=false, debug=false)
        debug && @show abs(Q.rowtmp[row])
        if abs(Q.rowtmp[row]) < Qd
            xtmp[splitdim] = h
            setrow!(Q.rowtmp, Q.dimpiv, (xtmp-xbase)./scale, ndims(rbox), splitdim)
            incremental_gaussian_elimination!(Q.rowtmp, Q.coefs, Q.rowzero, rng; canswap=false, debug=false)
            debug && @show abs(Q.rowtmp[row])
            if abs(Q.rowtmp[row]) < Qd
                # Neither was an improvement, so keep what we've got
                debug && println("skipping")
                return (row,-1), rbox
            end
        end
    end
    # Split rbox and insert into Q
    val = value(rbox)
    fsplit = MVector3(val,val,val)
    xtmp[splitdim] = l
    fsplit[j1] = f(xtmp)
    xtmp[splitdim] = h
    fsplit[j2] = f(xtmp)
    add_children!(rbox, splitdim, xsplit, fsplit, lwr, upr)
    if l != xbase[splitdim]
        xtmp[splitdim] = l
        debug && @show (xtmp.-xbase)./scale
        insrow1 = insert!(Q, (xtmp.-xbase)./scale, rbox.fvalues[j1]-c, splitdim, rbox.children[j1]; canswap=false, debug=debug)
    else
        debug && println("skipping l")
        insrow1 = 0
    end
    if h != xbase[splitdim]
        xtmp[splitdim] = h
        debug && @show (xtmp.-xbase)./scale
        insrow2 = insert!(Q, (xtmp.-xbase)./scale, rbox.fvalues[j2]-c, splitdim, rbox.children[j2]; canswap=false, debug=debug)
    else
        debug && println("skipping h")
        insrow2 = 0
    end
    # Return the leaf corresponding to the original
    xtmp[splitdim] = xcur
    j = 1; if j == j1 j += 1 end; if j == j2 j += 1 end
    return ((insrow1, insrow2), rbox.children[j])
end

function descend!(Q::QmIGE, c, box, x0, xbase, Δx, flag, bbs, scale, thresh, skip::Bool=false)
    box.qnconverged && return false  # don't build quadratic models that include points tagged as converged minima
    # Is this box so far away that we should prune this branch?
    boxbounds!(bbs, flag, box)
    for i = 1:ndims(box)
        xb, bb = xbase[i], bbs[i]
        min(bb[1] - xb, xb - bb[2]) > scaledthresh*scale[i] && return true
    end
    position!(Δx, flag, box, x0)
    thisx = isleaf(box) ? zero(eltype(Δx)) : Δx[box.splitdim]
    if !skip
        # Add this point to the model
        isgood = true
        leaf = find_leaf_at(box, Δx)
        for i = 1:length(Δx)
            Δx[i] = (Δx[i] - xbase[i])/scale[i]
            # To prevent distant points from numerically overwhelming the local structure of the box,
            # put in a cutoff. The concern is that large values here will dominate
            # eigenvalues that come from local data.
            isgood &= abs(Δx[i]) < scaledthresh
        end
        if isgood
            insert!(Q, Δx, value(box)-c, box.parent.splitdim, leaf)
        else
            # we have to save space for this dimension, since they are added
            # in splitdim order
            sd = box.parent.splitdim
            if sd ∉ Q.dimpiv
                Q.ndims[] += 1
                Q.dimpiv[Q.ndims[]] = sd
            end
        end
    end
    # Recursively add all children
    if !isleaf(box)
        iskip = findfirst(isequal(thisx), box.xvalues)
        for i = 1:3
            descend!(Q, c, box.children[i], x0, xbase, Δx, flag, bbs, scale, thresh, i==iskip) || return false
        end
    end
    return true
end

function Base.insert!(Q::QmIGE, Δx, Δf, splitdim::Integer, box::Box; canswap::Bool=true, debug::Bool=false)
    coefs, rhs, rowtmp = Q.coefs, Q.rhs, Q.rowtmp
    dimpiv, rowzero, ndims_old = Q.dimpiv, Q.rowzero, Q.ndims[]
    ndims = setrow!(rowtmp, dimpiv, Δx, ndims_old, splitdim)
    if ndims > ndims_old
        # Append to coefs. We always insert at the end of the
        # block corresponding to splitdim, so that we're performing
        # elimination on incoming points. That gives us a chance to
        # test for degeneracy before inserting them.
        i = ndims + (ndims*(ndims+1))÷2 # #gcoefs + #Bcoefs so far
        i == 0 && return i
        for j = 1:i
            coefs[i,j] = rowtmp[j]
        end
        rowzero[i] = false
        Q.rowbox[i] = box
        rhs[i] = Δf
        Q.ndims[] = ndims
        Q.nzrows[] += 1
        return i
    end
    ndims == 0 && return 0
    # We've seen all the nonzero dimensions in Δx previously
    # Use elimination to determine whether it provides novel information
    blockstart = ndims + ((ndims-1)*ndims)÷2    # first row for this splitting dimension
    blockend = blockstart + ndims               # last row for this splitting dimension
    finalrow, newrhs, box = incremental_gaussian_elimination!(rowtmp, coefs, rowzero, blockstart:blockend, Δf, rhs, box, Q.rowbox; canswap=canswap, debug=debug)
    if finalrow > 0
        rowzero[finalrow] = false
        for j = 1:finalrow
            coefs[finalrow,j] = rowtmp[j]
        end
        store!(rhs, newrhs, finalrow)
        store!(Q.rowbox, box, finalrow)
        Q.nzrows[] += 1
    end
    return finalrow
end

# This implements the logic of "flipped" Gaussian elimination (the zeros block on the
# upper triangle) when adding a new point. We do it in flipped form because incoming
# points tend to build a lower-triangular matrix: the tree structure
# adds dimensions one at a time, so Δx will be all-zeros for trailing dimensions
# (once permuted into the order specified by dimpiv).
function incremental_gaussian_elimination!(newrow::AbstractVector{T}, coefs::AbstractMatrix, rowzero::AbstractVector{Bool}, rng::AbstractUnitRange, newrhs=nothing, rhs=nothing, meta=nothing, metastore=nothing; canswap::Bool=true, debug::Bool=false) where T
    # @inline myapprox(x, y, rtol) = isequal(x, y) | (abs(x-y) < rtol*(abs(x) + abs(y)))
    @inline myapprox(x, y, rtol) = (x == y) | (abs(x-y) < rtol*(abs(x) + abs(y)))
    rtol = 1000*eps(T)

    debug && @show rng
    debug && println("input: ", newrow[rng])
    i = last(rng)
    @assert(first(rng) > 0)
    @inbounds while i >= first(rng)
        if newrow[i] == 0
            i -= 1
            continue
        end
        rowzero[i] && break
        if canswap && abs(newrow[i]) > abs(coefs[i,i])
            # Swap (for numeric stability)
            debug && println("swap at ", i, "($(newrow[i]) vs $(coefs[i,i])")
            for j = 1:i
                newrow[j], coefs[i,j] = coefs[i,j], newrow[j]
            end
            newrhs = swap!(rhs, newrhs, i)
            meta   = swap!(metastore, meta, i)
        end
        c = newrow[i]/coefs[i,i]
        @simd for j = 1:i-1
            sub = c * coefs[i,j] # coefsp[j,i]
            newrow[j] = ifelse(myapprox(newrow[j], sub, rtol), zero(sub), newrow[j] - sub)
        end
        newrow[i] = 0  # rather than subtracting, just set it (to avoid roundoff error)
        debug && @show i newrow[rng]
        newrhs = subtract(rhs, c, i, newrhs)
        i -= 1
    end
    if i < first(rng)
        i = 0
    end
    return i, newrhs, meta
end

function swap!(store, x, i)
    x, store[i] = store[i], x
    return x
end
swap!(::Any, ::Nothing, i) = nothing

store!(store, x, i) = store[i] = x
store!(::Any, ::Nothing, i) = nothing

subtract(rhs, c, i, newrhs) = newrhs - c*rhs[i]
subtract(rhs, c, i, ::Nothing) = nothing

function setrow!(rowtmp, dimpiv, Δx, ndims, splitdim)
    fill!(rowtmp, 0)
    k = 0
    have_splitdim = false
    for j = 1:ndims
        d = dimpiv[j]
        have_splitdim |= splitdim == d
        # The g coefs are linear in x, the B coefs quadratic
        # The implied storage order here matches `solve` below
        Δxd = Δx[d]
        if Δxd == 0
            k += j+1
            continue
        end
        rowtmp[k+=1] = Δxd  # g coef
        for i = 1:j-1
            rowtmp[k+=1] = Δxd * Δx[dimpiv[i]] # B[d,dimpiv[i]] coefficient
        end
        rowtmp[k+=1] = (Δxd * Δxd)/2 # B[d,d] coefficient
    end
    if !have_splitdim
        # We've not seen this dimension previously, so make it the
        # next one in dimpiv.
        ndims += 1
        dimpiv[ndims] = splitdim
        rowtmp[k+=1] = Δxd = Δx[splitdim]
        for i = 1:ndims-1
            rowtmp[k+=1] = Δxd * Δx[dimpiv[i]]
        end
        rowtmp[k+=1] = (Δxd * Δxd)/2
    end
    # Back up to last nonzero dimension
    while ndims > 0 && Δx[dimpiv[ndims]] == 0
        ndims -= 1
    end
    return ndims
end


text:
using PyCall

const COLUMN_NAMES = Dict{@compat(UInt64), Vector{UTF8String}}()

immutable Record{T}
    fields::T
    columns_idx::@compat(UInt64)
end
column_names(r) = COLUMN_NAMES[r.columns_idx]
Base.getindex(r::Record, idx::Integer) = r.fields[idx]

function Base.getindex(r::Record, idx::String)
    num_idx = findfirst(column_names(r), idx)
    r[num_idx]
end

Base.getindex(r::Record, idx::Symbol) = r[string(idx)]

function field_show(field)
    if isnull(field)
        "NULL"
    else
        sprint(show, get(field))
    end
end

function Base.show(io::IO, r::Record)
    columns = column_names(r)
    names = ["$column=>$(field_show(field))" for (column, field) in zip(columns, r.fields)]
    print(io, join(names, ", "))
end

function Record(pyo::PyObject)
    keys = pyo[:keys]()
    values = pyo[:values]()
    idx = hash(keys)
    COLUMN_NAMES[idx] = keys
    for i in eachindex(values)
        if values[i] == nothing
            values[i] = Nullable()  # TODO make correct type
        else
            values[i] = Nullable(values[i])
        end
    end
    Record((values...), idx)
end

immutable RecordSet
    records::Vector{Record}
    columns_idx::@compat(UInt64)
end

function Base.getindex(s::RecordSet, idx::Integer)
    s.records[idx]
end

function Base.getindex(s::RecordSet, idx)
    RecordSet(s.records[idx], s.columns_idx)
end

function RecordSet(records)
    if isempty(records)
        RecordSet(Record{()}, 0)
    else
        records = map(Record, records)
        RecordSet(records, records[1].columns_idx)
    end
end

Base.endof(rs::RecordSet) = endof(rs.records)

function Base.show(io::IO, rs::RecordSet)
    limit = min(length(rs.records), 10)
    for i=1:limit
        print(io, rs.records[i])
        if i<limit
            print(io, "\n")
        end
    end
    if length(rs.records) > limit
        println(io, "\n...")
    end
end


text:
@doc doc"""
    Sphere{N} <: Manifold

The unit sphere manifold $\mathbb S^n$ represented by $n+1$-Tuples, i.e. in by
vectors in $\mathbb R^{n+1}$ of unit length

# Constructor

    Sphere(n)

generates the $\mathbb S^{n}\subset \mathbb R^{n+1}$
"""
struct Sphere{N} <: Manifold end
Sphere(n::Int) = Sphere{n}()

@traitimpl HasMetric{Sphere,EuclideanMetric}

@doc doc"""
    manifold_dimension(S::Sphere)

Return the dimension of the manifold $\mathbb S^n$, i.e. $n$.
"""
manifold_dimension(S::Sphere{N}) where {N} = N

proj!(S::Sphere, x) = (x ./= norm(x))

project_tangent!(S::Sphere, w, x, v) = (w .= v .- dot(x, v) .* x)

@doc doc"""
    inner(S::Sphere, x, w, v)

compute the inner product of the two tangent vectors `w,v` from the tangent
plane at `x` on the sphere `S=`$\mathbb S^n$ using the restriction of the
metric from the embedding, i.e. $ (v,w)_x = v^\mathrm{T}w $.
"""
inner(S::Sphere, x, w, v) = dot(w, v)

norm(S::Sphere, x, v) = norm(v)

distance(S::Sphere, x, y) = acos(dot(x, y))

function exp!(S::Sphere, y, x, v)
    nv = norm(S, x, v)
    if nv ≈ 0.0
        y .= x
    else
        y .= cos(nv).*x .+ (sin(nv)/nv).*v
    end
    return y
end

function log!(S::Sphere, v, x, y)
    dot_xy = dot(x, y)
    θ = acos(dot_xy)
    if θ ≈ 0.0
        zero_tangent_vector!(S, v, x)
    else
        v .= (θ/sin(θ)) .* (y .- dot_xy.*x)
    end
    return v
end

injectivity_radius(S::Sphere, args...) = π

zero_tangent_vector(S::Sphere, x) = zero(x)
zero_tangent_vector!(S::Sphere, v, x) = (v .= zero(x))

"""
    uniform_sphere_distribution(S::Sphere, x)

Uniform distribution on given sphere. Generated points will be of similar
type to `x`.
"""
function uniform_distribution(S::Sphere, x)
    d = Distributions.MvNormal(zero(x), 1.0)
    return ProjectedPointDistribution(S, d, proj!, x)
end

"""
    is_manifold_point(S,x; kwargs...)

checks, whether `x` is a valid point on the [`Sphere`](@ref) `S`, i.e. is a vector
of length [`manifold_dimension`](@ref)`(S)+1` (approximately) of unit length.
The tolerance for the last test can be set using the ´kwargs...`.
"""
function is_manifold_point(S::Sphere{N},x; kwargs...) where {N}
    if length(x) != N+1
        throw(DomainError(size(x),"The point $(x) does not lie on $S, since its size is not $(N+1)."))
    end
    if !isapprox(norm(x), 1.; kwargs...)
        throw(DomainError(norm(x), "The point $(x) does not lie on the sphere $(S) since its norm is not 1."))
    end
    return true
end

"""
    is_tangent_vector(S,x,v; kwargs... )

checks whether `v` is a tangent vector to `x` on the [`Sphere`](@ref) `S`, i.e.
atfer [`is_manifold_point`](@ref)`(S,x)`, `v` has to be of same dimension as `x`
and orthogonal to `x`.
The tolerance for the last test can be set using the ´kwargs...`.
"""
function is_tangent_vector(S::Sphere{N},x,v; kwargs...) where N
    is_manifold_point(S,x)
    if length(v) != N+1
        throw(DomainError(size(v),
            "The vector $(v) is not a tangent to a point on $S since its size does not match $(N+1)."))
    end
    if !isapprox( abs(dot(x,v)), 0.; kwargs...)
        throw(DomainError(abs(dot(x,v)),
            "The vector $(v) is not a tangent vector to $(x) on $(S), since it is not orthogonal in the embedding."
        ))
    end
    return true
end

"""
    gaussian_sphere_tvector_distribution(S::Sphere, x, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `x`.
"""
function normal_tvector_distribution(S::Sphere, x, σ)
    d = Distributions.MvNormal(zero(x), σ)
    return ProjectedTVectorDistribution(S, x, d, project_tangent!, x)
end


text:
using ColorTypes, ImageView

function showImg(img)
    imshow(convert(Array{RGB24}, img), flipy=true)
    print()
end

"""
  Spiral Search 3D

  clock wise neighbour search around a center point (shell by shell search)
  first and last layer (z-axis) full 2D spiral search (whole layer is shell)
  layers in between: search outside (next shell) and already searched area (core) will be skipped
"""
function search3D_spiral(img, center)
  lastElement = true
  endOflayer = false
  endOflBox = false
  turn=false
  log=false

  MAX=reduce(*,size(img))
  xyz=[0,0,0]
  count=[0,0]
  element=0; j=0;
  xt=0; yt=0;
  layerElement = 0
  layer = 0
  edge=0
  layerCount = 0
  maxlength = 0
  layerElementCount = 0
  layerElementCountPrev =0

  X, Y, Z = center[2], center[1], center[3]

  while true
    if lastElement
      j += 1
      if j > 15 break end
      edge=0
      layer=1
      endOflayer = false
      endOflBox = false
      layerElement = 0

      dimMAX = 1 + j*2
      layerCount = dimMAX
      layerElementCountPrev = layerElementCount
      layerElementCount = dimMAX^2 # X*Y

      xyz=[0,0,j]
      xt=1; yt=-1;
      count=[1,1]
      turn=true
      lastElement = false
    end

    if !lastElement
      layerElement += 1
      element+=1

      notFirst = !(layer == 1 && layerElement == 1)
      firstLayer = layer <= 1
      lastLayer = layer >= layerCount
      max = layerElementCount-(firstLayer || lastLayer ? (lastLayer ? 1 : 0) : layerElementCountPrev)
      endOflayer = (layerElement > max)
      border = [abs(xyz[1]) >= count[1],abs(xyz[2]) >= count[2]]

      if notFirst
        switch = [turn && border[1],!turn && border[2]]
        if switch[1] xt=-xt; elseif switch[2] yt=-yt; end
        if switch[1] || switch[2] turn=!turn; end
      end

      if border[1] && border[2] edge+=1 end

      if endOflayer
        layerElement=1
        layer+=1
        edge=0
        firstLayer = layer <= 1
        lastLayer = layer >= layerCount
        lastElement = layer > layerCount
        count=[j,j]

        if lastLayer
          xt=1
          xyz=[0,0,xyz[3]]
          turn=true
          count=[1,1]
        end

        if !lastElement xyz[3]-=1; end
        if log print("\$"); if !lastElement print("\n////// |") end end
      end

      if edge > 3
        if log print("#") end
        edge=0
        xt=-xt
        count+=[1,1]
        turn = true
      end

      endOfCube = lastLayer && endOflayer
      firstLayerElement = layerElement == 1
      lastLayerElement = endOflayer || layerElement == max
      lastElement = layer > layerCount

      if notFirst; if turn xyz[1]+=xt; else xyz[2]+=yt; end; end
      if lastElement xyz=[0,0,xyz[3]]; end

      img[xyz[2]+Y,xyz[1]+X,xyz[3]+Z]=1 - clamp(((element/MAX) + (layerElement/layerElementCount)) * 0.5,0,1)

      if log print(x>=0 ? " " : "",x, "" , y>=0 ? " " : "",y,"",z>=0 ? " " : "", z," |") end
    end
  end
end

img = Array{Float32}(undef, 31, 31, 31)
hsz=(x->Int64(round(x/2))).(size(img))
fill!(img, 1)

search3D_spiral(img,hsz)

showImg(img)

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
function print_dict(results_dict, title="")
    max_len = maximum(length.(keys(results_dict)))
    line_len = max_len + 5
    if title != ""
        println(repeat("-", line_len+10))
        println(repeat(" ", 2) * title)
        println(repeat("-", line_len+10))
    end
    lines = [repeat(" ", 2) * x * repeat(" ", line_len-length(x)) * "| " * string(round(y, digits=3)) for (x, y) in results_dict]
    map(x -> println(x), lines)
end


text:
export zcom_NO
function zcom_NO(x)::Float64
    """
    Computes the z component of the center of mass of the NO projectile
    """
    z_N = x[1, 3]
    z_O = x[2, 3]
    zcom = (z_N * m_N + z_O*m_O)/(m_N + m_O)
    return zcom
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## general machinery for irrational mathematical constants

"""
    AbstractIrrational <: Real

Number type representing an exact irrational value.
"""
abstract type AbstractIrrational <: Real end

"""
    Irrational{sym} <: AbstractIrrational

Number type representing an exact irrational value denoted by the
symbol `sym`.
"""
struct Irrational{sym} <: AbstractIrrational end

show(io::IO, x::Irrational{sym}) where {sym} = print(io, sym)

function show(io::IO, ::MIME"text/plain", x::Irrational{sym}) where {sym}
    if get(io, :compact, false)
        print(io, sym)
    else
        print(io, sym, " = ", string(float(x))[1:15], "...")
    end
end

promote_rule(::Type{<:AbstractIrrational}, ::Type{Float16}) = Float16
promote_rule(::Type{<:AbstractIrrational}, ::Type{Float32}) = Float32
promote_rule(::Type{<:AbstractIrrational}, ::Type{<:AbstractIrrational}) = Float64
promote_rule(::Type{<:AbstractIrrational}, ::Type{T}) where {T<:Real} = promote_type(Float64, T)
promote_rule(::Type{S}, ::Type{T}) where {S<:AbstractIrrational,T<:Number} = promote_type(promote_type(S, real(T)), T)

AbstractFloat(x::AbstractIrrational) = Float64(x)
Float16(x::AbstractIrrational) = Float16(Float32(x))
Complex{T}(x::AbstractIrrational) where {T<:Real} = Complex{T}(T(x))

@pure function Rational{T}(x::AbstractIrrational) where T<:Integer
    o = precision(BigFloat)
    p = 256
    while true
        setprecision(BigFloat, p)
        bx = BigFloat(x)
        r = rationalize(T, bx, tol=0)
        if abs(BigFloat(r) - bx) > eps(bx)
            setprecision(BigFloat, o)
            return r
        end
        p += 32
    end
end
(::Type{Rational{BigInt}})(x::AbstractIrrational) = throw(ArgumentError("Cannot convert an AbstractIrrational to a Rational{BigInt}: use rationalize(Rational{BigInt}, x) instead"))

@pure function (t::Type{T})(x::AbstractIrrational, r::RoundingMode) where T<:Union{Float32,Float64}
    setprecision(BigFloat, 256) do
        T(BigFloat(x), r)
    end
end

float(::Type{<:AbstractIrrational}) = Float64

==(::Irrational{s}, ::Irrational{s}) where {s} = true
==(::AbstractIrrational, ::AbstractIrrational) = false

<(::Irrational{s}, ::Irrational{s}) where {s} = false
function <(x::AbstractIrrational, y::AbstractIrrational)
    Float64(x) != Float64(y) || throw(MethodError(<, (x, y)))
    return Float64(x) < Float64(y)
end

<=(::Irrational{s}, ::Irrational{s}) where {s} = true
<=(x::AbstractIrrational, y::AbstractIrrational) = x==y || x<y

# Irrationals, by definition, can't have a finite representation equal them exactly
==(x::AbstractIrrational, y::Real) = false
==(x::Real, y::AbstractIrrational) = false

# Irrational vs AbstractFloat
<(x::AbstractIrrational, y::Float64) = Float64(x,RoundUp) <= y
<(x::Float64, y::AbstractIrrational) = x <= Float64(y,RoundDown)
<(x::AbstractIrrational, y::Float32) = Float32(x,RoundUp) <= y
<(x::Float32, y::AbstractIrrational) = x <= Float32(y,RoundDown)
<(x::AbstractIrrational, y::Float16) = Float32(x,RoundUp) <= y
<(x::Float16, y::AbstractIrrational) = x <= Float32(y,RoundDown)
<(x::AbstractIrrational, y::BigFloat) = setprecision(precision(y)+32) do
    big(x) < y
end
<(x::BigFloat, y::AbstractIrrational) = setprecision(precision(x)+32) do
    x < big(y)
end

<=(x::AbstractIrrational, y::AbstractFloat) = x < y
<=(x::AbstractFloat, y::AbstractIrrational) = x < y

# Irrational vs Rational
@pure function rationalize(::Type{T}, x::AbstractIrrational; tol::Real=0) where T
    return rationalize(T, big(x), tol=tol)
end
@pure function lessrational(rx::Rational{<:Integer}, x::AbstractIrrational)
    # an @pure version of `<` for determining if the rationalization of
    # an irrational number required rounding up or down
    return rx < big(x)
end
function <(x::AbstractIrrational, y::Rational{T}) where T
    T <: Unsigned && x < 0.0 && return true
    rx = rationalize(T, x)
    if lessrational(rx, x)
        return rx < y
    else
        return rx <= y
    end
end
function <(x::Rational{T}, y::AbstractIrrational) where T
    T <: Unsigned && y < 0.0 && return false
    ry = rationalize(T, y)
    if lessrational(ry, y)
        return x <= ry
    else
        return x < ry
    end
end
<(x::AbstractIrrational, y::Rational{BigInt}) = big(x) < y
<(x::Rational{BigInt}, y::AbstractIrrational) = x < big(y)

<=(x::AbstractIrrational, y::Rational) = x < y
<=(x::Rational, y::AbstractIrrational) = x < y

isfinite(::AbstractIrrational) = true
isinteger(::AbstractIrrational) = false
iszero(::AbstractIrrational) = false
isone(::AbstractIrrational) = false

hash(x::Irrational, h::UInt) = 3*objectid(x) - h

widen(::Type{T}) where {T<:Irrational} = T

-(x::AbstractIrrational) = -Float64(x)
for op in Symbol[:+, :-, :*, :/, :^]
    @eval $op(x::AbstractIrrational, y::AbstractIrrational) = $op(Float64(x),Float64(y))
end
*(x::Bool, y::AbstractIrrational) = ifelse(x, Float64(y), 0.0)

round(x::Irrational, r::RoundingMode) = round(float(x), r)

"""
	@irrational sym val def
	@irrational(sym, val, def)

Define a new `Irrational` value, `sym`, with pre-computed `Float64` value `val`,
and arbitrary-precision definition in terms of `BigFloat`s given be the expression `def`.
"""
macro irrational(sym, val, def)
    esym = esc(sym)
    qsym = esc(Expr(:quote, sym))
    bigconvert = isa(def,Symbol) ? quote
        function Base.BigFloat(::Irrational{$qsym}, r::MPFR.MPFRRoundingMode=MPFR.ROUNDING_MODE[]; precision=precision(BigFloat))
            c = BigFloat(;precision=precision)
            ccall(($(string("mpfr_const_", def)), :libmpfr),
                  Cint, (Ref{BigFloat}, MPFR.MPFRRoundingMode), c, r)
            return c
        end
    end : quote
        function Base.BigFloat(::Irrational{$qsym}; precision=precision(BigFloat))
            setprecision(BigFloat, precision) do
                $(esc(def))
            end
        end
    end
    quote
        const $esym = Irrational{$qsym}()
        $bigconvert
        Base.Float64(::Irrational{$qsym}) = $val
        Base.Float32(::Irrational{$qsym}) = $(Float32(val))
        @assert isa(big($esym), BigFloat)
        @assert Float64($esym) == Float64(big($esym))
        @assert Float32($esym) == Float32(big($esym))
    end
end

big(x::AbstractIrrational) = BigFloat(x)
big(::Type{<:AbstractIrrational}) = BigFloat

# align along = for nice Array printing
function alignment(io::IO, x::AbstractIrrational)
    m = match(r"^(.*?)(=.*)$", sprint(show, x, context=io, sizehint=0))
    m === nothing ? (length(sprint(show, x, context=io, sizehint=0)), 0) :
    (length(m.captures[1]), length(m.captures[2]))
end

# inv
inv(x::AbstractIrrational) = 1/x

max_stars_repo_name: naveenarun/julia

max_issues_repo_name: naveenarun/julia

max_forks_repo_name: naveenarun/julia


text:
module PlayGios

export matchesupdate, pgrankingupdate, pgstatscompute

include("Atp.jl")

using Statistics
using DataFrames
using CSV
using Dates
import HTTP
using .Atp

savedir = ""

function __init__()
    savedirchange(tempdir() * "/PlayGios")
end

function savedirchange(dir)
    isdir(dir) || mkpath(dir)
    global savedir = dir
end

function matchesupdate()
    col = [:tourney_date => Date[],
                        :winner_name => String[],
                        :loser_name => String[]
                    ]
    typesdict = Dict(first.(col) .=> valtype.(last.(col)))
    fpth = "$savedir/matches.csv"
    if isfile(fpth)
        matches = CSV.read(fpth, DataFrame; types = typesdict)
    else
        matches = DataFrame(col)
        CSV.write(fpth, matches; header = true)
    end
    d0 = isempty(matches) ? Date(1968, 1, 1) : maximum(matches.tourney_date)
    baseurl_atp =
      "https://raw.githubusercontent.com/Tennismylife/TML-Database/master/"
    baseurl_qualchal =
      "https://raw.githubusercontent.com/JeffSackmann/tennis_atp/master/atp_matches_qual_chall_"
    function tourneydatevecadapt(tourney_date_vec::AbstractVector{<:Union{Missing, Date}})
        for rowid in 1:size(tourney_date_vec, 1)
            if ismissing(tourney_date_vec[rowid])
                val = findlast(.!ismissing.(tourney_date_vec[1:rowid]))
                tourney_date_vec[rowid] = isnothing(val) ?
                  Date(valyear) : tourney_date_vec[val]
            end
        end
        tourney_date_vec = convert(Vector{Date}, tourney_date_vec)
        return tourney_date_vec
    end
    for valyear in year(d0):year(today())
        println("matchesupdate ", valyear)
        url = "$baseurl_atp$valyear.csv"
        pget = HTTP.get(url)
        csv = CSV.read(IOBuffer(String(pget.body)), DataFrame;
          types = typesdict,
          dateformat = Dict(:tourney_date => dateformat"YYYYmmdd"))
        csv.tourney_date = tourneydatevecadapt(csv.tourney_date)
        idx = findall(csv.tourney_date .> d0)
        csv2wr = csv[idx, first.(col)]
        if valyear >= 1978
            url = "$baseurl_qualchal$valyear.csv"
            pget = HTTP.get(url)
            csv = CSV.read(IOBuffer(String(pget.body)), DataFrame;
              types = typesdict,
              dateformat = Dict(:tourney_date => dateformat"YYYYmmdd"))
            csv.tourney_date = tourneydatevecadapt(csv.tourney_date)
            idx = findall(csv.tourney_date .> d0)
            append!(csv2wr, csv[idx, first.(col)])
            sort!(csv2wr, :tourney_date)
        end
        CSV.write(fpth, csv2wr; append = true)
        append!(matches, csv2wr)
    end
    return matches
end

function pgrankingget()
    col = [:date => Date[],
            :pgname => String[],
            :npg => Int16[]
        ]
    typesdict = Dict(first.(col) .=> valtype.(last.(col)))
    fpth = "$savedir/pgranking.csv"
    if isfile(fpth)
        pgranking = CSV.read(fpth, DataFrame; types = typesdict)
    else
        pgranking = DataFrame(col)
        CSV.write(fpth, pgranking; header = true)
    end
    return pgranking
end

function pgrankingupdate()
    # matches.csv update and retrieval
    matches = matchesupdate()
    matches = matches[matches.tourney_date .>= Atp.atprankingstart, :]
    # pgranking retrieval
    fpth = "$savedir/pgranking.csv"
    pgranking = pgrankingget()
    # pgranking computation
    if isempty(pgranking)
        atpranking = atp_ranking_get(Atp.atprankingstart)
        pgname = atpranking.player[1]
        dateidx = findfirst(matches.winner_name .== pgname)
        isnothing(dateidx) && error("invalid pgname=$pgname")
        date = matches.tourney_date[dateidx]
        CSV.write(fpth, Tables.table([date pgname 1]); append = true)
        push!(pgranking, [date pgname 1])
    end
    pgname = pgranking.pgname[end]
    date = pgranking.date[end]
    nextpg_idx = findlast((matches.winner_name .== pgname) .&
      (matches.tourney_date .== date))
    if isnothing(nextpg_idx)
        error("Invalid matches.csv or pgranking.csv")
    end
    while true
        idx = findfirst(matches.loser_name[nextpg_idx:end] .== pgname)
        if isnothing(idx)
            playerretired =
              matches.tourney_date[nextpg_idx] < today() - Year(2) &&
              isempty(findall(matches.winner_name[nextpg_idx+1:end] .== pgname))
            if playerretired
                atpranking = atp_ranking_get(matches.tourney_date[nextpg_idx])
                pgname = atpranking.player[1]
                idx = 1 +
                  findfirst(matches.winner_name[nextpg_idx+1:end] .== pgname)
                isnothing(idx) && error("invalid pgname=$pgname")
                if matches.winner_name[nextpg_idx - 1 + idx] != pgname
                    error("You did something wrong here")
                end
            else
                break
            end
        end
        nextpg_idx = nextpg_idx - 1 + idx
        date = matches.tourney_date[nextpg_idx]
        pgname = matches.winner_name[nextpg_idx]
        npg = sum(pgranking.pgname .== pgname) + 1
        CSV.write(fpth, Tables.table([date pgname npg]); append = true)
        push!(pgranking, [date pgname npg])
        print("pgrankingupdate ", [date pgname npg], "\r")
    end
    println()
    return pgranking
end

function pgstatscompute()
    pgranking = pgrankingget()
    if isempty(pgranking)
        pgranking = pgrankingupdate()
    end
    pgranking.dateout = [pgranking.date[2:end]; today()]
    gd = groupby(pgranking, :pgname)
    pgstats = DataFrame(pgname = [], npg = [], dateinfirst = [],
                        dateoutlast = [], durationmax = [], durationmean = [])
    for gd_i in gd
        pgname = gd_i.pgname[1]
        npg = gd_i.npg[end]
        dateinfirst = gd_i.date[1]
        dateoutlast = gd_i.dateout[end]
        daydiff = getfield.(gd_i.dateout .- gd_i.date, :value)
        durationmax = maximum(daydiff)
        durationmean = mean(daydiff)
        push!(pgstats,
          [pgname npg dateinfirst dateoutlast durationmax durationmean])
    end
    sort!(pgstats, [:npg, :durationmean]; rev = true)
    CSV.write("$savedir/pgstats.csv", pgstats)
    return pgstats
end

end # module


text:
module ComplexOptInterface

import MutableArithmetics
const MA = MutableArithmetics

import MathOptInterface
const MOI = MathOptInterface

struct HermitianPositiveSemidefiniteConeTriangle <: MOI.AbstractVectorSet
    side_dimension::Int
end
function MOI.dimension(set::HermitianPositiveSemidefiniteConeTriangle)
    return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(set.side_dimension)) +
        MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(set.side_dimension - 1))
end

function MOI.Utilities.set_dot(x::AbstractVector, y::AbstractVector,
                               set::HermitianPositiveSemidefiniteConeTriangle)
    sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension)
    result = MOI.Utilities.set_dot(x, y, sym)
    for k in (MOI.dimension(sym) + 1):MOI.dimension(set)
        result = MA.add_mul!!(result, 2, x[k], y[k])
    end
    return result
end

include("Bridges/Bridges.jl")

end # module


text:
module test_lupin_partition_by_condition

using Test
using Lupin

@test partition(<=(3), [1, 2, 3, 4, 5]) == ([1, 2, 3], [4, 5])
@test partition(<=(5), [1, 2, 3, 4, 5]) == ([1, 2, 3, 4, 5], [])

@test partition(<=(3), 1:5) == (1:3, 4:5)
@test partition(<=(5), 1:5) == (1:5, [])

end # module test_lupin_partition_by_condition


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>2.0, :alg=>TMJets(max_steps=1_000, abs_tol=1e-12, orderT=6, orderQ=2)),
    "X0"   => Hyperrectangle(low=fill(0.99, 9), high=fill(1.01, 9)),
    "info" => Dict(
        "name" => "biomodel9d",
        "dim" => 9,
        "linear" => false,
        "hybrid" => false
    )
])


text:
import ProximalOperators: Conjugate

Conjugate(f::IndFree) = IndZero()
Conjugate(f::IndZero) = IndFree()
Conjugate(f::SqrNormL2) = SqrNormL2(1.0/f.lambda)

# TODO: Add other useful functions and calculus rules such as translation


text:
using JuMP
using Gurobi

abstract OPSolver

type GurobiExactSolver <: OPSolver
    output::Bool
    multithreaded::Bool
    time_limit::Float64
    max_license_tries::Int
end
GurobiExactSolver(;output=false,
                   multithreaded=true,
                   time_limit=Inf,
                   max_license_tries=10) = GurobiExactSolver(output,
                                                       multithreaded,
                                                       time_limit,
                                                       max_license_tries)

function first_move(solver::OPSolver, op::SimpleOP)
    return solve_op(solver, op)[2]
end

type OPSolution
    x
    u
end
OPSolution() = OPSolution(nothing,nothing)

build_path(op, s::OPSolution) = build_path(op, s.x)
first_move(op, s::OPSolution) = build_path(op, s)[2]

function build_path(op, x)
    xopt = round(Int,x)
    path = [op.start]
    current = -1
    dist_sum = 0.0
    while current != op.stop
        if current==-1
            current=op.start
        end
        current = findfirst(xopt[current,:])
        push!(path, current)
    end
    return path
end

function solve_op(solver::GurobiExactSolver, op)
    path = build_path(op, gurobi_solve(op, solver=solver))
    @assert distance(op, path) <= op.distance_limit + 1e-3*op.distance_limit
    return path
end

function gurobi_solve(op; solver=GurobiExactSolver(), initial::OPSolution=OPSolution())
    threads = solver.multithreaded ? 0: 1
    local m
    for l in 1:solver.max_license_tries
        try
            m = Model(solver=Gurobi.GurobiSolver(OutputFlag=solver.output, Threads=threads, TimeLimit=solver.time_limit))
        catch ex
            if l >= solver.max_license_tries
                println("Maximum number of tries ($l) to create Gurobi model failed.")
                sleep(10.0)
                rethrow(ex)
            else
                sleep(0.1)
                continue
            end
        end
        break
    end

    N = length(op)

    without_start = [1:op.start-1; op.start+1:N]
    without_stop = [1:op.stop-1; op.stop+1:N]
    without_both = intersect(without_start, without_stop)

    @defVar(m, x[1:N,1:N], Bin)
    @defVar(m, 2 <= u[without_start] <= N, Int)

    if !is(initial.x, nothing)
        setValue(x, initial.x)
        for keytuple in keys(u)
            key = keytuple[1]
            setValue(u[key], initial.u[key])
        end
    end

    @setObjective(m, Max, sum{ sum{op.r[i]*x[i,j], j=1:N}, i=1:N })

    @addConstraint(m, sum{x[op.start,j], j=without_start} == 1)
    @addConstraint(m, sum{x[i,op.stop], i=without_stop} == 1)
    # problem

    @addConstraint(m, connectivity[k=without_both], sum{x[i,k], i=1:N} == sum{x[k,j], j=1:N})

    @addConstraint(m, once[k=1:N], sum{x[k,j], j=1:N} <= 1)
    @addConstraint(m, sum{ sum{op.distances[i,j]*x[i,j], j=1:N}, i=1:N } <= op.distance_limit)
    @addConstraint(m, nosubtour[i=without_start,j=without_start], u[i]-u[j]+1 <= (N-1)*(1-x[i,j]))

    if op.start != op.stop
        @addConstraint(m, sum{x[op.stop,i],i=1:N} == 0)
    end

    local status
    for l in 1:solver.max_license_tries
        try
            status = JuMP.solve(m, suppress_warnings=true)
        catch ex
            if l >= solver.max_license_tries
                println("Maximum number of tries ($l) to create Gurobi model failed.")
                sleep(10.0)
                rethrow(ex)
            else
                sleep(0.1)
                continue
            end
        end
        break
    end

    if solver.time_limit < Inf
        if status != :Optimal && status != :UserLimit
            warn("Not solved to optimality:\n$op")
        end
    else
        if status != :Optimal
            warn("Not solved to optimality:\n$op")
        end
    end

    soln = OPSolution(getValue(x), getValue(u))

    if distance(op,build_path(op,soln)) > op.distance_limit
        warn("Path Length: $(distance(op,build_path(op,soln))), Limit: $(op.distance_limit)")
    end
    return soln
end

function solve_opcsp_feedback(op, d::Vector{Float64}, solver=GurobiExactSolver())
    d_belief = MVN(zeros(length(op)), copy(op.covariance))
    openset = collect(1:length(op))
    total_dist = 0.0
    choice = -1
    path = [op.start]
    apply_measurement!(d_belief, op.start, d[op.start])
    while choice != op.stop
        mean_future = SimpleOP(op.r[openset] + d_belief.mean[openset],
                          [],
                          op.distance_limit-total_dist,
                          findfirst(openset, path[end]),
                          findfirst(openset, op.stop),
                          op.distances[openset, openset])
        # mpc_soln = gurobi_solve(mean_future)
        # mpc_path = build_path(mean_future, mpc_soln)
        mpc_path = solve_op(solver, mean_future)
        choice = openset[mpc_path[2]]
        total_dist += op.distances[path[end], choice]
        push!(path, choice)
        splice!(openset, mpc_path[1])
        apply_measurement!(d_belief, choice, d[choice])
    end
    return path
end

type Cheater end

function cheat(op, d::Vector{Float64})
    complete_knowledge = SimpleOP(op.r+d, op.positions, op.distance_limit, op.start, op.stop)
    soln = gurobi_solve(complete_knowledge)
    return build_path(op, soln)
end

type OnlyUnobservableUncertainty end

function unobservable_dist(op, d::Vector{Float64})
    ucov = zeros(op.covariance)
    ud = zeros(d)
    for i in 1:length(op)
        cov = copy(op.covariance)
        mvn = MVN(zeros(d),cov)
        for j in 1:length(op)
            if j!=i
                mvn = apply_measurement(mvn, j, d[j])
            end
        end
        @assert sum(mvn.covariance) == mvn.covariance[i,i]
        ucov[i,i] = mvn.covariance[i,i]
        ud[i] = mvn.mean[i]
    end
    @assert isdiag(ucov)
    return MVN(ud, ucov)
end

function cheat_observe_others(op, d::Vector{Float64})
    udist = unobservable_dist(op, d)
    uproblem = SimpleOP(op.r+udist.mean, op.positions, op.distance_limit, op.start, op.stop)
    soln = gurobi_solve(uproblem)
    return build_path(op, soln)
end


text:
"""
Fit Garch model to financial data.
"""
module Garch

using Random, Turing, Distributed, MCMCChains, Distributions, Statistics, StatsBase, LinearAlgebra, Optim, NLSolversBase

export simulate, garch, garchZD, fitBayes, fitMLE, garchGJR

"""
    simulate(a, b, c; y₀=0.001, h₀=0.0, n=100, d=Normal(0., 1.))

Simulate Garch(1,1) process:
`rₜ = √hₜ ξₜ`
`hₜ = a + b yₜ² + c hₜ₋₁`

"""
function simulate(a, b, c; y₀ = 0.001, h₀ = 0.0, n = 100, d = Normal(0.0, 1.0))
    @assert c >= 0 && a >= 0 && b > 0 "a,c shoud be >=0, b > 0"
    h = fill(0.0, n)
    y = similar(h)
    y[1] = y₀
    ξ = rand(d, n - 1)
    for i = 2:n
        h[i] = c + a * y[i-1] * y[i-1] + b * h[i-1]
        y[i] = ξ[i-1] * √h[i]
    end
    return y, h
end

"""
    garchZD(y::Vector{Float64})

Model Zero drift Garch(1,1) process for Bayesian estimation.
"""
@model function garchZD(y)
    n = length(y)
    h = 0.0
    a ~ Uniform(0.0001, 1)
    b ~ Uniform(0.0, 1)
    y[1] ~ Normal(y[1], h)
    for i = 2:n
        h = a * y[i-1] * y[i-1] + b * h
        y[i] ~ Normal(0, √h)
    end
end

"""
    garch(y::Vector{Float64})

Model Garch(1,1) for Bayesian estimation.
"""
@model function garch(y)
    n = length(y)
    h = 0.01
    a ~ Uniform(0.0001, 1.0)
    b ~ Uniform(0.0001, 1.0)
    w ~ Uniform(0.0, 1.0)
    y[1] ~ Normal(y[1], h)
    for i = 2:n
        h = w + a * y[i-1] * y[i-1] + b * h
        y[i] ~ Normal(0, √h)
    end
end

"""
    garchGJR(y::Vector{Float64}, z::Vector{Float64})

Model GJR-Garch(1,1) for Bayesian estimation.
"""
@model function garchGJR(y, z)
    n = length(y)
    h = 0.01
    a1 ~ Uniform(0.0001, 10.0)
    a2 ~ Uniform(0.0001, 10.0)
    β ~ Uniform(0.0001, 10.0)
    a0 ~ Uniform(0.0, 10.0)
    δ ~ Uniform(0.0001, 10.0)
    y[1] ~ Normal(y[1], h)
    for i = 2:n
        j = y[i-1] < 0.0 ? 1 : 0
        h = a0 + (a1 + j * a2) * y[i-1] * y[i-1] + β * h + δ * z[i-1]
        y[i] ~ Normal(0, sqrt(h))
    end
end

"""
    fitMLE(x::Vector{Float64}; isZD=true)

Fit Garch(1,1) on `x` using Maximum Likelihood estimation.
"""
function fitMLE(x; isZD = true)
    model = isZD ? garchZD(x) : garch(x)
    out = optimize(model, MLE())
    m, s = coef(out), stderror(out)
    return (pars = m, serr = s)
end

"""
    fitBayes(x::Vector{Float64}; isZD=true, nchains=4, nsims=1000, isDistributed=false)

Fit Garch(1,1) on `x` using Bayesian MCMC.
"""
function fitBayes(x; isZD = true, nchains = 4, nsims = 1000, isDistributed = false)
    model = isZD ? garchZD(x) : garch(x)
    if isDistributed
        chn = sample(model, NUTS(0.65), MCMCDistributed(), nsims, nchains, progress = false)
    else
        chn = sample(model, sampler, nsims, progress = false)
    end
    return chn
end

end


text:
#=
Modia library with 1D heat transfer component models (inspired from Modelica Standard Library).

Developer: Martin Otter, DLR-SR
Copyright 2021, DLR Institute of System Dynamics and Control
License: MIT (expat)
=#

using ModiaLang
using Unitful

HeatPort = Model( T = potential,   # Absolute temperature
                  Q_flow = flow )  # Heat flow into the component

FixedTemperature = Model(
    T    = 293.15u"K",
    port = HeatPort,
    equations = :[port.T = T]
)

FixedHeatFlow = Model(
    Q_flow = 0.0u"W",   # Fixed heat flow into the connected component
    port   = HeatPort,
    equations = :[port.Q_flow = -Q_flow]
)

HeatCapacitor = Model(
    C = 1.0u"J/K",
    port = HeatPort,
    T = Var(init = 293.15u"K"),
    equations = :[T = port.T,
                  der(T) = port.Q_flow/C]
)

ThermalConductor = Model(
    G = 1.0u"W/K",
    port_a = HeatPort,
    port_b = HeatPort,
    equations = :[dT = port_a.T - port_b.T
                  0  = port_a.Q_flow + port_b.Q_flow
                  port_a.Q_flow = G*dT]
)

# Fully insulated rod with 1D heat transfer and port_a/port_b on left/right side

T_grad1(T,Ta,dx,i) = i == 1         ? (Ta - T[1])/(dx/2) : (T[i-1] - T[i]  )/dx
T_grad2(T,Tb,dx,i) = i == length(T) ? (T[i] - Tb)/(dx/2) : (T[i]   - T[i+1])/dx

InsulatedRod = Model(
    L       = 1.0u"m",            # Length of rod
    A       = 0.0004u"m^2",       # Rod area
    rho     = 7500.0u"kg/m^3",    # Density of rod material
    lambda  = 74.0u"W/(m*K)",     # Thermal conductivity of rod material
    c       = 450.0u"J/(kg*K)",   # Specific heat capacity of rod material
    port_a  = HeatPort,           # Heat port on left side
    port_b  = HeatPort,           # Heat port on right side
    T       = Map(init = fill(293.15u"K", 1)),  # Initial temperature and number of nodes
    equations = :[
        n   = length(T)
        dx  = L/n
        Ce  = c*rho*A*dx
        k   = lambda*A/Ce
        der(T) = k*[T_grad1(T,port_a.T,dx,i) - T_grad2(T,port_b.T,dx,i) for i in eachindex(T)]
        port_a.Q_flow = lambda*A*(port_a.T - T[1])/(dx/2)
        port_b.Q_flow = lambda*A*(port_b.T - T[n])/(dx/2)
    ]
)

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
@add_problem sdp_and_exp function sdp_and_exp_log_det_atom(
    handle_problem!,
    ::Val{test},
    atol,
    rtol,
    ::Type{T},
) where {T,test}
    x = Variable(2, 2)
    p = maximize(logdet(x), [x[1, 1] == 1, x[2, 2] == 1]; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol = 10atol rtol = rtol
        @test evaluate(logdet(x)) ≈ 0 atol = 10atol rtol = rtol
    end
end


text:
using RDatasets, Clustering, Random, LinearAlgebra, Plots; pyplot()
Random.seed!(0)

df = dataset("cluster", "xclara")
n,_ = size(df)
dataPoints = [convert(Array{Float64,1},df[i,:]) for i in 1:n]
shuffle!(dataPoints)
D = [norm(pt1 - pt2) for pt1 in dataPoints, pt2 in dataPoints]

result = hclust(D)
for K in 1:30
    clusters = cutree(result,k=K)
    println("K=$(K): ",[sum(clusters .== i) for i in 1:K])
end

cluster(ell,K) = (1:n)[cutree(result,k=K) .== ell]

C1, C2, C3 = cluster(1,30),cluster(2,30),cluster(3,30)

plt = scatter( first.(dataPoints[C1]),last.(dataPoints[C1]),c=:blue, msw=0)
      scatter!( first.(dataPoints[C2]),last.(dataPoints[C2]), c=:red, msw=0)
      scatter!( first.(dataPoints[C3]),last.(dataPoints[C3]), c=:green, msw=0)
for ell in 4:30
    clst = cluster(ell,30)
    scatter!(first.(dataPoints[clst]),last.(dataPoints[clst]),
        ms=10, c=:purple, shape=:xcross, ratio=:equal, legend=:none,
        xlabel="V1", ylabel="V2")
end
plot(plt)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
modeled after the DRCSIM AtlasCommand ROS message, but updated for the atlas hardware API.
control parameter spec:
  * q, qd, f are sensed position, velocity, torque, from AtlasJointState
  * q_d, qd_d, f_d are desired position, velocity, torque, from
    AtlasJointDesired

  The final joint command will be:
   k_q_p   * ( q_d - q )     +
   k_q_i   * 1/s * ( q_d - q ) +
   k_qd_p  * ( qd_d - qd )   +
   k_f_p   * ( f_d - f )     +
   ff_qd   * qd        +
   ff_qd_d   * qd_d        +
   ff_f_d  * f_d         +
   ff_const
"""
mutable struct atlas_command_t <: LCMType
    utime::Int64
    num_joints::Int32
    joint_names::Vector{String}
    position::Vector{Float64}
    velocity::Vector{Float64}
    effort::Vector{Float64}

    k_q_p::Vector{Float64}
    k_q_i::Vector{Float64}
    k_qd_p::Vector{Float64}
    k_f_p::Vector{Float64}
    ff_qd::Vector{Float64}  # maps to kd_position in drcsim API (there isnt an equivalent gain in the bdi api)
    ff_qd_d::Vector{Float64}
    ff_f_d::Vector{Float64}
    ff_const::Vector{Float64}

    # k_torque can be an unsigned int 8value from 0 to 255,
    # at run time, a double between 0 and 1 is obtained
    # by dividing by 255.0d.
    k_effort::Vector{UInt8}

    # max allowed controller update period in milli-seconds simulation time
    # for controller synchronization. See SynchronizationStatistics.msg for
    # tracking synchronization status.
    desired_controller_period_ms::UInt8
end

@lcmtypesetup(atlas_command_t,
    joint_names => (num_joints,),
    position => (num_joints,),
    velocity => (num_joints,),
    effort => (num_joints,),
    k_q_p => (num_joints,),
    k_q_i => (num_joints,),
    k_qd_p => (num_joints,),
    k_f_p => (num_joints,),
    ff_qd => (num_joints,),
    ff_qd_d => (num_joints,),
    ff_f_d => (num_joints,),
    ff_const => (num_joints,),
    k_effort => (num_joints,),
)


text:
#
# Structure to read mmCIF field data indexes
#
Base.@kwdef mutable struct Indexes_mmCIF_fields
  index::Int = 0
  name::Int = 0
  resname::Int = 0
  chain::Int = 0
  resnum::Int = 0
  x::Int = 0
  y::Int = 0
  z::Int = 0
  beta::Int = 0
  occup::Int = 0
end


text:
using MyExample
using Test

@testset "MyExample.jl" begin
@test my_f(2,1) == 7
@test my_f(2,3) == 13
@test my_f(1,3) == 11
end

@testset "Derivative Tests" begin
@test derivative_of_my_f(2,1) == 2
end


text:
using SafeTestsets
const LONGER_TESTS = false

const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Core"
    @time @safetestset "DynamicHMC" begin include("dynamicHMC.jl") end
    @time @safetestset "Turing" begin include("turing.jl") end
    # @time @safetestset "ABC" begin include("abc.jl") end
end

if GROUP == "Stan" || GROUP == "All"
    @time @safetestset "Stan_String" begin include("stan_string.jl") end
    @time @safetestset "Stan" begin include("stan.jl") end
end


text:
# lex indicates if its lexigraphical (i.e., x, y) or reverse (y, x)
# If in lexigraphical order the coefficient vector's entries
# corrrespond to the following basis polynomials:
# [T0(x) * T0(y), T1(x) * T0(y), T0(x) * T1(y), T2(x) * T0(y), T1(x) * T1(y), T0(x) * T2(y), ...]
# else, if not in lexigraphical order:
# [T0(x) * T0(y), T0(x) * T1(y), T1(x) * T0(y), T0(x) * T2(y), T1(x) * T1(y), T2(x) * T0(y), ...]
"""
Pre-plan an Inverse Padua Transform.
"""
struct IPaduaTransformPlan{lex,IDCTPLAN,T}
    cfsmat::Matrix{T}
    idctplan::IDCTPLAN
end

IPaduaTransformPlan(cfsmat::Matrix{T},idctplan,::Type{Val{lex}}) where {T,lex} =
    IPaduaTransformPlan{lex,typeof(idctplan),T}(cfsmat,idctplan)

"""
Pre-plan an Inverse Padua Transform.
"""
function plan_ipaduatransform!(::Type{T},N::Integer,lex) where T
    n=Int(cld(-3+sqrt(1+8N),2))
    if N ≠ div((n+1)*(n+2),2)
        error("Padua transforms can only be applied to vectors of length (n+1)*(n+2)/2.")
    end
    IPaduaTransformPlan(Array{T}(undef,n+2,n+1),FFTW.plan_r2r!(Array{T}(undef,n+2,n+1),FFTW.REDFT00),lex)
end

plan_ipaduatransform!(::Type{T},N::Integer) where {T} = plan_ipaduatransform!(T,N,Val{true})
plan_ipaduatransform!(v::AbstractVector{T},lex...) where {T} = plan_ipaduatransform!(eltype(v),length(v),lex...)

function *(P::IPaduaTransformPlan,v::AbstractVector{T}) where T
    cfsmat=trianglecfsmat(P,v)
    n,m=size(cfsmat)
    rmul!(view(cfsmat,:,2:m-1),0.5)
    rmul!(view(cfsmat,2:n-1,:),0.5)
    tensorvals=P.idctplan*cfsmat
    paduavec!(v,P,tensorvals)
end

ipaduatransform!(v::AbstractVector,lex...) = plan_ipaduatransform!(v,lex...)*v
"""
Inverse Padua Transform maps the 2D Chebyshev coefficients to the values of the interpolation polynomial at the Padua points.
"""
ipaduatransform(v::AbstractVector,lex...) = plan_ipaduatransform!(v,lex...)*copy(v)

"""
Creates ``(n+2)x(n+1)`` Chebyshev coefficient matrix from triangle coefficients.
"""
function trianglecfsmat(P::IPaduaTransformPlan{true},cfs::AbstractVector)
    N=length(cfs)
    n=Int(cld(-3+sqrt(1+8N),2))
    cfsmat=fill!(P.cfsmat,0)
    m=1
    for d=1:n+1
        @inbounds for k=1:d
            j=d-k+1
            cfsmat[k,j]=cfs[m]
            if m==N
                return cfsmat
            else
                m+=1
            end
        end
    end
    return cfsmat
end

function trianglecfsmat(P::IPaduaTransformPlan{false},cfs::AbstractVector)
    N=length(cfs)
    n=Int(cld(-3+sqrt(1+8N),2))
    cfsmat=fill!(P.cfsmat,0)
    m=1
    for d=1:n+1
        @inbounds for k=d:-1:1
            j=d-k+1
            cfsmat[k,j]=cfs[m]
            if m==N
                return cfsmat
            else
                m+=1
            end
        end
    end
    return cfsmat
end

"""
Vectorizes the function values at the Padua points.
"""
function paduavec!(v,P::IPaduaTransformPlan,padmat::Matrix)
    n=size(padmat,2)-1
    N=(n+1)*(n+2)
    if iseven(n)>0
        d=div(n+2,2)
        m=0
        @inbounds for i=1:n+1
            v[m+1:m+d]=view(padmat,1+mod(i,2):2:n+1+mod(i,2),i)
            m+=d
        end
    else
        @inbounds v[:]=view(padmat,1:2:N-1)
    end
    return v
end

"""
Pre-plan a Padua Transform.
"""
struct PaduaTransformPlan{lex,DCTPLAN,T}
    vals::Matrix{T}
    dctplan::DCTPLAN
end

PaduaTransformPlan(vals::Matrix{T},dctplan,::Type{Val{lex}}) where {T,lex} =
    PaduaTransformPlan{lex,typeof(dctplan),T}(vals,dctplan)

"""
Pre-plan a Padua Transform.
"""
function plan_paduatransform!(::Type{T},N::Integer,lex) where T
    n=Int(cld(-3+sqrt(1+8N),2))
    if N ≠ ((n+1)*(n+2))÷2
        error("Padua transforms can only be applied to vectors of length (n+1)*(n+2)/2.")
    end
    PaduaTransformPlan(Array{T}(undef,n+2,n+1),FFTW.plan_r2r!(Array{T}(undef,n+2,n+1),FFTW.REDFT00),lex)
end

plan_paduatransform!(::Type{T},N::Integer) where {T} = plan_paduatransform!(T,N,Val{true})
plan_paduatransform!(v::AbstractVector{T},lex...) where {T} = plan_paduatransform!(eltype(v),length(v),lex...)

function *(P::PaduaTransformPlan,v::AbstractVector{T}) where T
    N=length(v)
    n=Int(cld(-3+sqrt(1+8N),2))
    vals=paduavalsmat(P,v)
    tensorcfs=P.dctplan*vals
    m,l=size(tensorcfs)
    rmul!(tensorcfs,T(2)/(n*(n+1)))
    rmul!(view(tensorcfs,1,:),0.5)
    rmul!(view(tensorcfs,:,1),0.5)
    rmul!(view(tensorcfs,m,:),0.5)
    rmul!(view(tensorcfs,:,l),0.5)
    trianglecfsvec!(v,P,tensorcfs)
end

paduatransform!(v::AbstractVector,lex...) = plan_paduatransform!(v,lex...)*v
"""
Padua Transform maps from interpolant values at the Padua points to the 2D Chebyshev coefficients.
"""
paduatransform(v::AbstractVector,lex...) = plan_paduatransform!(v,lex...)*copy(v)

"""
Creates ``(n+2)x(n+1)`` matrix of interpolant values on the tensor grid at the ``(n+1)(n+2)/2`` Padua points.
"""
function paduavalsmat(P::PaduaTransformPlan,v::AbstractVector)
    N=length(v)
    n=Int(cld(-3+sqrt(1+8N),2))
    vals=fill!(P.vals,0.)
    if iseven(n)>0
        d=div(n+2,2)
        m=0
        @inbounds for i=1:n+1
            vals[1+mod(i,2):2:n+1+mod(i,2),i]=view(v,m+1:m+d)
            m+=d
        end
    else
        @inbounds vals[1:2:end]=view(v,:)
    end
    return vals
end

"""
Creates length ``(n+1)(n+2)/2`` vector from matrix of triangle Chebyshev coefficients.
"""
function trianglecfsvec!(v,P::PaduaTransformPlan{true},cfs::Matrix)
    m=size(cfs,2)
    l=1
    for d=1:m
        @inbounds for k=1:d
            j=d-k+1
            v[l]=cfs[k,j]
            l+=1
        end
    end
    return v
end

function trianglecfsvec!(v,P::PaduaTransformPlan{false},cfs::Matrix)
    m=size(cfs,2)
    l=1
    for d=1:m
        @inbounds for k=d:-1:1
            j=d-k+1
            v[l]=cfs[k,j]
            l+=1
        end
    end
    return v
end

"""
Returns coordinates of the ``(n+1)(n+2)/2`` Padua points.
"""
function paduapoints(::Type{T}, n::Integer) where T
    N=div((n+1)*(n+2),2)
    MM=Matrix{T}(undef,N,2)
    m=0
    delta=0
    NN=fld(n+2,2)
    @inbounds for k=n:-1:0
        if isodd(n)>0
            delta=mod(k,2)
        end
        @inbounds for j=NN+delta:-1:1
            m+=1
            MM[m,1]=sinpi(T(k)/n-T(0.5))
            if isodd(n-k)>0
                MM[m,2]=sinpi((2j-one(T))/(n+1)-T(0.5))
            else
                MM[m,2]=sinpi(T(2j-2)/(n+1)-T(0.5))
            end
        end
    end
    return MM
end

paduapoints(n::Integer) = paduapoints(Float64,n)


text:
using Revise
includet(joinpath("..", "init.jl"))

db_url_in = joinpath(
    "sqlite:///$(spineitemsPath)", "case_study_b2", "case_study_b2.sqlite"
)
db_url_out = joinpath(
    "sqlite:///$(spineitemsPath)", "case_study_b2_out", "case_study_b2_out.sqlite"
)

m = SpineOpt.run_spinemodel(db_url_in, db_url_out; cleanup=false, log_level=3)


text:
# Shorthand type names
typealias Vec Vector{Float64}
typealias Mat Matrix{Float64}
typealias SymMat Symmetric{Float64,Matrix{Float64}}

# Utility functions
AAt(A::Mat) = Symmetric(A*A')
AtA(A::Mat) = Symmetric(A'*A)

# sumsq
sumsq(x::Vec) = norm(x)^2

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench