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

module Sys
@doc """
Provide methods for retrieving information about hardware and the operating system.
""" Sys

export BINDIR,
       CPU_CORES,
       CPU_NAME,
       WORD_SIZE,
       ARCH,
       MACHINE,
       KERNEL,
       JIT,
       cpu_info,
       cpu_summary,
       uptime,
       loadavg,
       free_memory,
       total_memory,
       isapple,
       isbsd,
       islinux,
       isunix,
       iswindows,
       isexecutable,
       which

import ..Base: show

global BINDIR = ccall(:jl_get_julia_bindir, Any, ())::String
"""
    Sys.BINDIR

A string containing the full path to the directory containing the `julia` executable.
"""
:BINDIR

# helper to avoid triggering precompile warnings

global CPU_CORES
"""
    Sys.CPU_CORES

The number of logical CPU cores available in the system.

See the Hwloc.jl package for extended information, including number of physical cores.
"""
:CPU_CORES

"""
    Sys.ARCH

A symbol representing the architecture of the build configuration.
"""
const ARCH = ccall(:jl_get_ARCH, Any, ())

"""
    Sys.KERNEL

A symbol representing the name of the operating system, as returned by `uname` of the build configuration.
"""
const KERNEL = ccall(:jl_get_UNAME, Any, ())

"""
    Sys.MACHINE

A string containing the build triple.
"""
const MACHINE = Base.MACHINE

"""
    Sys.WORD_SIZE

Standard word size on the current machine, in bits.
"""
const WORD_SIZE = Core.sizeof(Int) * 8

function __init__()
    env_cores = get(ENV, "JULIA_CPU_CORES", "")
    global CPU_CORES = if !isempty(env_cores)
        env_cores = tryparse(Int, env_cores)
        if !(env_cores isa Int && env_cores > 0)
            Core.print(Core.stderr, "WARNING: couldn't parse `JULIA_CPU_CORES` environment variable. Defaulting Sys.CPU_CORES to 1.\n")
            env_cores = 1
        end
        env_cores
    else
        Int(ccall(:jl_cpu_cores, Int32, ()))
    end
    global SC_CLK_TCK = ccall(:jl_SC_CLK_TCK, Clong, ())
    global CPU_NAME = ccall(:jl_get_cpu_name, Ref{String}, ())
    global JIT = ccall(:jl_get_JIT, Ref{String}, ())
    global BINDIR = ccall(:jl_get_julia_bindir, Any, ())::String
    nothing
end

mutable struct UV_cpu_info_t
    model::Ptr{UInt8}
    speed::Int32
    cpu_times!user::UInt64
    cpu_times!nice::UInt64
    cpu_times!sys::UInt64
    cpu_times!idle::UInt64
    cpu_times!irq::UInt64
end
mutable struct CPUinfo
    model::String
    speed::Int32
    cpu_times!user::UInt64
    cpu_times!nice::UInt64
    cpu_times!sys::UInt64
    cpu_times!idle::UInt64
    cpu_times!irq::UInt64
    CPUinfo(model,speed,u,n,s,id,ir)=new(model,speed,u,n,s,id,ir)
end
CPUinfo(info::UV_cpu_info_t) = CPUinfo(unsafe_string(info.model), info.speed,
    info.cpu_times!user, info.cpu_times!nice, info.cpu_times!sys,
    info.cpu_times!idle, info.cpu_times!irq)

function _show_cpuinfo(io::IO, info::Sys.CPUinfo, header::Bool=true, prefix::AbstractString="    ")
    tck = SC_CLK_TCK
    if header
        println(io, info.model, ": ")
        print(io, " "^length(prefix))
        println(io, "    ", lpad("speed", 5), "    ", lpad("user", 9), "    ", lpad("nice", 9), "    ",
                lpad("sys", 9), "    ", lpad("idle", 9), "    ", lpad("irq", 9))
    end
    print(io, prefix)
    unit = tck > 0 ? " s  " : "    "
    tc = max(tck, 1)
    d(i, unit=unit) = lpad(string(round(Int64,i)), 9) * unit
    print(io,
          lpad(string(info.speed), 5), " MHz  ",
          d(info.cpu_times!user / tc), d(info.cpu_times!nice / tc), d(info.cpu_times!sys / tc),
          d(info.cpu_times!idle / tc), d(info.cpu_times!irq / tc, tck > 0 ? " s" : "  "))
    if tck <= 0
        print(io, "ticks")
    end
end

show(io::IO, info::CPUinfo) = _show_cpuinfo(io, info, true, "    ")

function _cpu_summary(io::IO, cpu::AbstractVector{CPUinfo}, i, j)
    if j-i < 9
        header = true
        for x = i:j
            header || println(io)
            _show_cpuinfo(io, cpu[x], header, "#$(x-i+1) ")
            header = false
        end
    else
        summary = CPUinfo(cpu[i].model,0,0,0,0,0,0)
        count = j - i + 1
        for x = i:j
            summary.speed += cpu[i].speed
            summary.cpu_times!user += cpu[x].cpu_times!user
            summary.cpu_times!nice += cpu[x].cpu_times!nice
            summary.cpu_times!sys += cpu[x].cpu_times!sys
            summary.cpu_times!idle += cpu[x].cpu_times!idle
            summary.cpu_times!irq += cpu[x].cpu_times!irq
        end
        summary.speed = div(summary.speed,count)
        _show_cpuinfo(io, summary, true, "#1-$(count) ")
    end
    println(io)
end

function cpu_summary(io::IO=stdout, cpu::AbstractVector{CPUinfo} = cpu_info())
    model = cpu[1].model
    first = 1
    for i = 2:length(cpu)
        if model != cpu[i].model
            _cpu_summary(io, cpu, first, i-1)
            first = i
        end
    end
    _cpu_summary(io, cpu, first, length(cpu))
end

function cpu_info()
    UVcpus = Ref{Ptr{UV_cpu_info_t}}()
    count = Ref{Int32}()
    Base.uv_error("uv_cpu_info",ccall(:uv_cpu_info, Int32, (Ptr{Ptr{UV_cpu_info_t}}, Ptr{Int32}), UVcpus, count))
    cpus = Vector{CPUinfo}(undef, count[])
    for i = 1:length(cpus)
        cpus[i] = CPUinfo(unsafe_load(UVcpus[], i))
    end
    ccall(:uv_free_cpu_info, Cvoid, (Ptr{UV_cpu_info_t}, Int32), UVcpus[], count[])
    return cpus
end

"""
    Sys.uptime()

Gets the current system uptime in seconds.
"""
function uptime()
    uptime_ = Ref{Float64}()
    Base.uv_error("uv_uptime",ccall(:uv_uptime, Int32, (Ptr{Float64},), uptime_))
    return uptime_[]
end

"""
    Sys.loadavg()

Get the load average. See: https://en.wikipedia.org/wiki/Load_(computing).
"""
function loadavg()
    loadavg_ = Vector{Float64}(undef, 3)
    ccall(:uv_loadavg, Cvoid, (Ptr{Float64},), loadavg_)
    return loadavg_
end

free_memory() = ccall(:uv_get_free_memory, UInt64, ())
total_memory() = ccall(:uv_get_total_memory, UInt64, ())

"""
    Sys.get_process_title()

Get the process title. On some systems, will always return an empty string.
"""
function get_process_title()
    buf = Vector{UInt8}(undef, 512)
    err = ccall(:uv_get_process_title, Cint, (Ptr{UInt8}, Cint), buf, 512)
    Base.uv_error("get_process_title", err)
    return unsafe_string(pointer(buf))
end

"""
    Sys.set_process_title(title::AbstractString)

Set the process title. No-op on some operating systems.
"""
function set_process_title(title::AbstractString)
    err = ccall(:uv_set_process_title, Cint, (Cstring,), title)
    Base.uv_error("set_process_title", err)
end

"""
    Sys.maxrss()

Get the maximum resident set size utilized in bytes.
See also:
    - man page of getrusage(2) on Linux and FreeBSD.
    - windows api `GetProcessMemoryInfo`
"""
maxrss() = ccall(:jl_maxrss, Csize_t, ())

"""
    Sys.isunix([os])

Predicate for testing if the OS provides a Unix-like interface.
See documentation in [Handling Operating System Variation](@ref).
"""
function isunix(os::Symbol)
    if iswindows(os)
        return false
    elseif islinux(os) || isbsd(os)
        return true
    else
        throw(ArgumentError("unknown operating system \"$os\""))
    end
end

"""
    Sys.islinux([os])

Predicate for testing if the OS is a derivative of Linux.
See documentation in [Handling Operating System Variation](@ref).
"""
islinux(os::Symbol) = (os == :Linux)

"""
    Sys.isbsd([os])

Predicate for testing if the OS is a derivative of BSD.
See documentation in [Handling Operating System Variation](@ref).

!!! note
    The Darwin kernel descends from BSD, which means that `Sys.isbsd()` is
    `true` on macOS systems. To exclude macOS from a predicate, use
    `Sys.isbsd() && !Sys.isapple()`.
"""
isbsd(os::Symbol) = (os == :FreeBSD || os == :OpenBSD || os == :NetBSD || os == :DragonFly || os == :Darwin || os == :Apple)

"""
    Sys.iswindows([os])

Predicate for testing if the OS is a derivative of Microsoft Windows NT.
See documentation in [Handling Operating System Variation](@ref).
"""
iswindows(os::Symbol) = (os == :Windows || os == :NT)

"""
    Sys.isapple([os])

Predicate for testing if the OS is a derivative of Apple Macintosh OS X or Darwin.
See documentation in [Handling Operating System Variation](@ref).
"""
isapple(os::Symbol) = (os == :Apple || os == :Darwin)

for f in (:isunix, :islinux, :isbsd, :isapple, :iswindows)
    @eval $f() = $(getfield(@__MODULE__, f)(KERNEL))
end

if iswindows()
    function windows_version()
        verinfo = ccall(:GetVersion, UInt32, ())
        VersionNumber(verinfo & 0xFF, (verinfo >> 8) & 0xFF, verinfo >> 16)
    end
else
    windows_version() = v"0.0"
end

"""
    Sys.windows_version()

Return the version number for the Windows NT Kernel as a `VersionNumber`,
i.e. `v"major.minor.build"`, or `v"0.0.0"` if this is not running on Windows.
"""
windows_version

const WINDOWS_VISTA_VER = v"6.0"

"""
    Sys.isexecutable(path::String)

Return `true` if the given `path` has executable permissions.
"""
function isexecutable(path::String)
    if iswindows()
        return isfile(path)
    else
        # We use `access()` and `X_OK` to determine if a given path is
        # executable by the current user.  `X_OK` comes from `unistd.h`.
        X_OK = 0x01
        ccall(:access, Cint, (Ptr{UInt8}, Cint), path, X_OK) == 0
    end
end
isexecutable(path::AbstractString) = isexecutable(String(path))

"""
    Sys.which(program_name::String)

Given a program name, search the current `PATH` to find the first binary with
the proper executable permissions that can be run, and return an absolute
path. Raise `ErrorException` if no such program is available.  If a path with
a directory in it is passed in for `program_name`, tests that exact path
for executable permissions only (with `.exe` and `.com` extensions added on
Windows platforms); no searching of `PATH` is performed.
"""
function which(program_name::String)
    # Build a list of program names that we're going to try
    program_names = String[]
    base_pname = basename(program_name)
    if iswindows()
        # If the file already has an extension, try that name first
        if !isempty(splitext(base_pname)[2])
            push!(program_names, base_pname)
        end

        # But also try appending .exe and .com`
        for pe in (".exe", ".com")
            push!(program_names, string(base_pname, pe))
        end
    else
        # On non-windows, we just always search for what we've been given
        push!(program_names, base_pname)
    end

    path_dirs = String[]
    program_dirname = dirname(program_name)
    # If we've been given a path that has a directory name in it, then we
    # check to see if that path exists.  Otherwise, we search the PATH.
    if isempty(program_dirname)
        # If we have been given just a program name (not a relative or absolute
        # path) then we should search `PATH` for it here:
        pathsep = iswindows() ? ';' : ':'
        path_dirs = abspath.(split(get(ENV, "PATH", ""), pathsep))

        # On windows we always check the current directory as well
        if iswindows()
            pushfirst!(path_dirs, pwd())
        end
    else
        push!(path_dirs, abspath(program_dirname))
    end

    # Here we combine our directories with our program names, searching for the
    # first match among all combinations.
    for path_dir in path_dirs
        for pname in program_names
            program_path = joinpath(path_dir, pname)
            # If we find something that matches our name and we can execute
            if isexecutable(program_path)
                return realpath(program_path)
            end
        end
    end

    # If we couldn't find anything, complain
    error("$program_name not found")
end
which(program_name::AbstractString) = which(String(program_name))

end # module Sys

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
# Interruptible demands.

# Import packages.
using Dates, TimeSeries

# Demands that can be interrupted.

"""
A simple interruptible demand.

# Type parameters
- `T <: TimeType`: timestamp
- `L`            : network location

# Fields
- `location`: network location of the demand
- `envelope`: timeseries of minimum and maximum allowable demands

```
example = InterruptibleDemand(
    "a fixed location",
    TimeArray(
        [Time(0) , Time(8)   , Time(9) , Time(17)  , Time(18)],
        [(5., 8.), (10., 12.), (8., 8.), (11., 13.), (3., 4.)]
    )
) :: InterruptibleDemand{Time,String}
```
"""
struct InterruptibleDemand{T,L} <: FlexibleDemand{T,L}
    location :: L
    envelope :: Envelope{T}
end

"""
The "envelope" of minimum and maximum allowable demands at each time pont, with a location for the demand.

# Arguments
- `demand :: InteruptibleDemand{T,L}`: the demand
"""
function envelope(demand :: InterruptibleDemand{T,L}) :: LocatedEnvelope{T,L} where L where T <: TimeType
    TimeArray(
        timestamp(demand.envelope),
        map(x -> (demand.location, x[1], x[2]), values(demand.envelope))
    )
end

ext: jl


text:
# Window functions

immutable Cwindow
    w::Cint

    function Cwindow(w)
        if w ∉ 0:5
            throw(ArgumentError("invalid window index $w: must be ∈ 0:5"))
        end
        new(w)
    end
end

function _window!(wtype::Cwindow, x::StridedVector{Cdouble}; normalize::Int=1)
    # normalize: Int
    #    0 : don't normalize
    #    1 : normalize by power
    #    2 : normalize by magnitude
    if normalize ∉ 0:2
        throw(ArgumentError("invalid normalize flag $normalize, must be ∈ 0:2"))
    end
    g = ccall((:window, libSPTK), Cdouble,
              (Cwindow, Ptr{Cdouble}, Cint, Cint), wtype, x, length(x),
              normalize)
    x
end

for (f, wtype) in [(:blackman, 0),
                   (:hamming, 1),
                   (:hanning, 2),
                   (:bartlett, 3),
                   (:trapezoid, 4),
                   (:rectangular, 5)]
    @eval begin
        function $f(n::Integer; kargs...)
            y = ones(Float64, n)
            _window!(Cwindow($wtype), y; kargs...)
        end
    end
end


text:
"""
    SPECTrecon
System matrix (forward and back-projector) for SPECT image reconstruction.
"""
module SPECTrecon

    const RealU = Number # Union{Real, Unitful.Length}

    include("foreach.jl")
    include("helper.jl")
    include("plan-rotate.jl")
    include("rotatez.jl")
    include("plan-psf.jl")
    include("psf-gauss.jl")
    include("fft_convolve.jl")
    include("spectplan.jl")
    include("project.jl")
    include("backproject.jl")

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2020, DLR Institute of System Dynamics and Control

"""
    module Modia3D.Basics

Utility constants and functions for Modia3D

# Main developers

Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Basics

export trailingPartOfName

export neps, sign_eps, radToDeg, degToRad
export getAndCheckFullLibraryPath, getEnvironmentVariable
export zeroMVector, onesMVector, nullMRotation
export ZeroMVector
export normalizeVector, BoundingBox

export readDictOfStructsFromJSON
export listKeys
export deleteItem

export linearMovement
export PositionMVector, RotationMMatrix, assertRotationMatrix

using  StaticArrays
using  LinearAlgebra
using  JSON
using  Unitful

include("environment.jl")
include("constantsAndFunctions.jl")

end


text:
# get "super node" of a node and update all nodes in the path to the supernode
# Definition: A "super node" is a node in the same connected component, with a lower index than the given node (or equal if the index is the lowest one in the connected component)
function supernode!(superNodes, i)
    if i != superNodes[i]
        superNodes[i] = supernode!(superNodes, superNodes[i])
    end
    superNodes[i]
end

# makes matrix connected
# assumes D and DFull are symmetric
function makeconnected!(D::AbstractMatrix{T},DFull::AbstractMatrix{T}) where {T<:Real}
    N = size(D,1)

    # find all connected components
    superNodes = collect(1:N) # a "super node" is another node in the same connected component, s.t. if we follow through a path of super nodes, we will always reach the same "root" vertex
    for j=2:N
        for i=1:j-1
            D[i,j] != typemax(T) && (superNodes[j] = min(superNodes[j],i)) # always choose the lower one as super node by making the higher index refer to the lower
        end
    end

    # make sure all connected components are represented by their super nodes
    for i=2:N
        superNodes[i] = superNodes[superNodes[i]] # makes sure every index up to i has the lowest index in the connected component as their supernode
    end
    nbrComponents = length(unique(superNodes))

    if nbrComponents>1
        # collect and sort edges
        M = div(N*(N-1),2)
        E    = zeros(T,M)
        from = zeros(Int,M)
        to   = zeros(Int,M)
        ind = 1
        for j=2:N
            for i=1:j-1
                E[ind]    = DFull[i,j]
                from[ind] = i
                to[ind]   = j
                ind += 1
            end
        end
        perm = sortperm(E)
        E    = E[perm]
        from = from[perm]
        to   = to[perm]

        # add edges from shortest to longest, but only if they connect two different components
        for k=1:M
            nbrComponents==1 && break # done

            i = supernode!(superNodes,from[k])
            j = supernode!(superNodes,to[k])
            i==j && continue # the nodes belong to the same component
            i,j = minmax(i,j)

            # use this edge to connect the components
            superNodes[j] = i
            D[from[k],to[k]] = D[to[k],from[k]] = E[k]
            nbrComponents -= 1
        end
    end
    D
end

# Creates distance matrix based on Isomap distances
function isomapmatrix(D::AbstractMatrix{T}, nearestNeighbours::Integer, distanceThreshold::Real) where {T<:Real}
    @assert ndims(D)==2
    @assert size(D,1)==size(D,2)
    @assert nearestNeighbours>0
    N = size(D,1)
    nearestNeighbours = min(N,nearestNeighbours+1) # add 1 to include self as neighbour

    I = similar(D) # Isomap matrix
    I[:] .= typemax(T)

    for v=1:N # for each vertex, look at outgoing edges (i.e. row), and keep all edges in the neighbourhood (i.e. below distanceThreshold, but at least nearestNeighbour edges)
        # distsSorted = sort(D[v,:])
        # k = maximum(nearestNeighbours+1, searchsortedlast(distsSorted,distanceThreshold)) # +1 to include 0-distance to the vertex itself
        # d = distsSorted[k]

        # mask = D[v,:].<=d
        # I[v,mask] = D[v,mask]

        mask = D[v,:] .<= distanceThreshold

        if count(mask)<nearestNeighbours # not enough neighbours within radius? Include at least nearestNeighbour edges.
            mask = D[v,:] .<= sort(D[v,:][:])[nearestNeighbours] # this includes ties
        end

        I[v,mask] .= D[v,mask] # copy edges
    end

    # make symmetric
    I = min.(I,I')

    # make sure graph is connected by adding shortest distance between components
    makeconnected!(I,D)

    # make the graph complete by finding shortest distances between all nodes
    floydwarshall!(I)
end

function isomap(D::AbstractMatrix{T}, p::Integer, nearestNeighbours::Integer, distanceThreshold::Real) where {T<:Real}
    I = isomapmatrix(D, nearestNeighbours, distanceThreshold)
    mds(I,p)
end

function kruskalisomap(D::AbstractMatrix{T}, p::Integer, nearestNeighbours::Integer, distanceThreshold::Real; kwargs...) where {T<:Real}
    I = isomapmatrix(D, nearestNeighbours, distanceThreshold)
    X0 = mds(I,p)
    kruskalmds(I, p, X0; kwargs...)
end


text:
import SolverTraces: safe_lerp, ilerp, sat_ilerp,
    CrayonWrap, StrWrap
using Crayons

@testset "LERP" begin
    for t = 0.0:0.1:1.0
        @test safe_lerp(0.0, 1.0, t) == t
        @test safe_lerp(0.0, 2.0, t) == 2t
    end
    @test safe_lerp(0.0, 1.0, NaN) == 0.0
    @test safe_lerp(0.0, 1.0, -Inf) == 0.0
    @test safe_lerp(0.0, 1.0, Inf) == 1.0

    @test ilerp.(0, 2, 0:0.05:1) == vcat(repeat([0], 6),
                                         repeat([1], 9),
                                         repeat([2], 6))

    @test sat_ilerp(0, 2, 1.1) == 4
    @test sat_ilerp(0, 2, -0.1) == 0
    @test sat_ilerp(0, 128, 1.1) == 255
end

@testset "Linear colorants" begin
    @testset "Continuous" begin
        lc = LinearColorant(-1.0, 1.0, ((255, 0, 0), (0, 255, 0)))
        @test lc(-1.0) == CrayonWrap(foreground=(255, 0, 0))
        @test lc(0.0) == CrayonWrap(foreground=(255, 255, 0))
        @test lc(1.0) == CrayonWrap(foreground=(0, 255, 0))
    end

    @testset "Discrete steps" begin
        lc = LinearColorant(-1.0, 1.0, [crayon"light_red", crayon"yellow", crayon"light_green"])
        @test lc(-1.0) == StrWrap(crayon"light_red")
        @test lc(-0.6) == StrWrap(crayon"light_red")
        @test lc(-0.5) == StrWrap(crayon"yellow")
        @test lc(0.0) == StrWrap(crayon"yellow")
        @test lc(0.5) == StrWrap(crayon"yellow")
        @test lc(0.6) == StrWrap(crayon"light_green")
        @test lc(1.0) == StrWrap(crayon"light_green")
    end
end


text:
using LaserTypes
using LaserTypes: pochhammer, immutable_cache
using Unitful, UnitfulAtomic
using StaticArrays
using UnPack
using Test

units = [:SI_unitful, :atomic_unitful, :SI, :atomic]
x₀ = SVector{3}(0u"μm",0u"μm",0u"μm")
t₀ = 0u"s"
x = (x₀, auconvert.(x₀), ustrip.(x₀), ustrip.(auconvert.(x₀)))
t = (t₀, auconvert.(t₀), ustrip.(t₀), ustrip.(auconvert.(t₀)))

@testset "$unit" for (unit, xᵢ, tᵢ) in zip(units, x, t)
    s = setup_laser(LaguerreGaussLaser, unit, m = 1, p = 1, profile=ConstantProfile)
    Ex, Ey, Ez = E(xᵢ, tᵢ, s)
    @test iszero(Ex)
    @test iszero(Ey)
    @test iszero(Ez)

    @unpack E₀, Nₚₘ, z_R = immutable_cache(s)
    c = s.constants.c
    Bx, By, Bz = B(xᵢ, tᵢ,s)
    @test iszero(Bx)
    @test iszero(By)
    @test Bz ≈ - factorial(s.p)/pochhammer(abs(s.m)+1,s.p)*E₀*Nₚₘ*(√2*s.w₀)/(c*z_R)
end

@testset "Default laser" begin
    s = setup_laser(LaguerreGaussLaser, :atomic; λ = 1e5, a₀ = 1, w₀ = 1e6,
        profile = GaussProfile(τ = 100, z₀ = 0.))

    @test E([0, 0, 0], 0, s)[1] == s.derived.E₀
end


text:
function block2d_dojo(mechanism::Mechanism, F;
    rtol=1e-10,
    btol=1e-10,
    undercut=1.0,
    no_progress_undercut=1.0,
    mode::Symbol=:friction)

    (mode == :friction) && (idx = 1)
    (mode == :impact) && (idx = 2)
    timestep= mechanism.timestep
    opts_grad = SolverOptions(
        rtol=rtol,
        btol=btol,
        undercut=undercut,
        no_progress_undercut=no_progress_undercut,
        verbose=false)

    initialize!(mechanism, :block2d,
        x=[0.0, 0.0],
        v=[0.0, 0.0],
        θ=0.0,
        ω=0.0)
    z0 = get_maximal_state(mechanism)
    u0 = zeros(3)
    u0[idx] += F

    ∇x, ∇u = get_minimal_gradients!(mechanism, z0, u0,
        opts=opts_grad)

    step!(mechanism, z0, u0, opts=opts_grad)

    z1 = get_next_state(mechanism)
    x0 = maximal_to_minimal(mechanism, z0)
    x1 = maximal_to_minimal(mechanism, z1)

    ∂x∂F = ∇u[idx,idx]
    Δx = x1[idx] - x0[idx]
    return Δx, ∂x∂F
end

function block2d_gradientbundle(mechanism::Mechanism, F;
    N::Int=100,
    Σ=1e-6*I,
    rtol=1e-10,
    btol=1e-10,
    undercut=1.5,
    no_progress_undercut=10.0,
    mode::Symbol=:friction)

    (mode == :friction) && (idx = 1)
    (mode == :impact) && (idx = 2)
    timestep= mechanism.timestep
    opts_grad = SolverOptions(
            rtol=rtol,
            btol=btol,
            undercut=undercut,
            no_progress_undercut=no_progress_undercut)

    initialize!(mechanism, :block2d,
        x=[0.0, 0.0],
        v=[0.0, 0.0],
        θ=0.0,
        ω=0.0)
    z0 = get_maximal_state(mechanism)
    u0 = zeros(3)
    u0[idx] += F

    GB0 = gradient_bundle_0th(mechanism, z0, u0, idx,
        opts=opts_grad,
        N=N,
        Σ=Σ)
    GB1 = gradient_bundle_1st(mechanism, z0, u0, idx,
        opts=opts_grad,
        N=N,
        Σ=Σ)
    step!(mechanism, z0, u0,
        opts=opts_grad)
    z1 = get_next_state(mechanism)
    x0 = maximal_to_minimal(mechanism, z0)
    x1 = maximal_to_minimal(mechanism, z1)

    GB0xF = GB0[idx,idx]
    GB1xF = GB1[idx,idx]
    Δx = x1[idx] - x0[idx]
    return Δx, GB0xF, GB1xF
end


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

# TODO: optimize this
function Base.string(dt::DateTime)
    y,m,d = yearmonthday(days(dt))
    h,mi,s = hour(dt),minute(dt),second(dt)
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    mm = lpad(m,2,"0")
    dd = lpad(d,2,"0")
    hh = lpad(h,2,"0")
    mii = lpad(mi,2,"0")
    ss = lpad(s,2,"0")
    ms = millisecond(dt) == 0 ? "" : string(millisecond(dt)/1000.0)[2:end]
    return "$yy-$mm-$(dd)T$hh:$mii:$ss$(ms)"
end
Base.show(io::IO,x::DateTime) = print(io,string(x))
function Base.string(dt::Date)
    y,m,d = yearmonthday(value(dt))
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    mm = lpad(m,2,"0")
    dd = lpad(d,2,"0")
    return "$yy-$mm-$dd"
end
Base.show(io::IO,x::Date) = print(io,string(x))

### Parsing
const english = Dict{String,Int}("january"=>1,"february"=>2,"march"=>3,"april"=>4,
                 "may"=>5,"june"=>6,"july"=>7,"august"=>8,"september"=>9,
                 "october"=>10,"november"=>11,"december"=>12)
const abbrenglish = Dict{String,Int}("jan"=>1,"feb"=>2,"mar"=>3,"apr"=>4,
                     "may"=>5,"jun"=>6,"jul"=>7,"aug"=>8,"sep"=>9,
                     "oct"=>10,"nov"=>11,"dec"=>12)
const MONTHTOVALUE = Dict{String,Dict{String,Int}}("english"=>english)
const MONTHTOVALUEABBR = Dict{String,Dict{String,Int}}("english"=>abbrenglish)

# Date/DateTime Parsing
abstract Slot{T<:Any}

immutable DelimitedSlot{T<:Any} <: Slot{T}
    parser::Type{T}
    letter::Char
    width::Int
    transition::Union{Regex,AbstractString}
end

immutable FixedWidthSlot{T<:Any} <: Slot{T}
    parser::Type{T}
    letter::Char
    width::Int
end

immutable DateFormat
    slots::Array{Slot,1}
    prefix::AbstractString # optional transition from the start of a string to the 1st slot
    locale::AbstractString
end

abstract DayOfWeekSlot

# Slot rules translate letters into types. Note that
# list of rules can be extended.
immutable SlotRule
    rules::Array{Type}
end
const SLOT_RULE = SlotRule(Array{Type}(256))

getindex(collection::SlotRule, key::Char) = collection.rules[Int(key)]
setindex!(collection::SlotRule, value::Type, key::Char) = collection.rules[Int(key)] = value
keys(c::SlotRule) = map(Char, filter(el -> isassigned(c.rules, el) && c.rules[el] != Void, eachindex(c.rules)))

SLOT_RULE['y'] = Year
SLOT_RULE['Y'] = Year
SLOT_RULE['m'] = Month
SLOT_RULE['u'] = Month
SLOT_RULE['U'] = Month
SLOT_RULE['e'] = DayOfWeekSlot
SLOT_RULE['E'] = DayOfWeekSlot
SLOT_RULE['d'] = Day
SLOT_RULE['H'] = Hour
SLOT_RULE['M'] = Minute
SLOT_RULE['S'] = Second
SLOT_RULE['s'] = Millisecond

duplicates(slots) = any(map(x->count(y->x.parser==y.parser,slots),slots) .> 1)

"""
    DateFormat(format::AbstractString, locale::AbstractString="english") -> DateFormat

Construct a date formatting object that can be used for parsing date strings or
formatting a date object as a string. For details on the syntax for `format` see
[`DateTime(::AbstractString, ::AbstractString)`](:ref:`parsing <man-date-parsing>`) and
[`format`](:ref:`formatting <man-date-formatting>`).
"""
function DateFormat(f::AbstractString, locale::AbstractString="english")
    slots = Slot[]
    prefix = ""
    params = ()
    last_offset = 1

    letters = join(keys(SLOT_RULE), "")
    for m in eachmatch(Regex("(?<!\\\\)([\\Q$letters\\E])\\1*"), f)
        letter = f[m.offset]
        typ = SLOT_RULE[letter]

        width = length(m.match)
        tran = replace(f[last_offset:m.offset-1], r"\\(.)", s"\1")

        if isempty(params)
            prefix = tran
        else
            slot = tran == "" ? FixedWidthSlot(params...) : DelimitedSlot(params..., tran)
            push!(slots,slot)
        end

        params = (typ,letter,width)
        last_offset = m.offset + width
    end

    if !isempty(params)
        if last_offset > endof(f)
            slot = DelimitedSlot(params..., r"(?=\s|$)")
        else
            tran = replace(f[last_offset:end], r"\\(.)", s"\1")
            if tran == ""
                slot = FixedWidthSlot(params...)
            else
                slot = DelimitedSlot(params..., tran)
            end
        end
        push!(slots,slot)
    end

    duplicates(slots) && throw(ArgumentError("Two separate periods of the same type detected"))
    return DateFormat(slots,prefix,locale)
end

const SLOTERROR = ArgumentError("Non-digit character encountered")
slotparse(slot,x,locale) = !ismatch(r"[^0-9\s]",x) ? slot.parser(x) : throw(SLOTERROR)
function slotparse(slot::Slot{Month},x,locale)
    if slot.letter == 'm'
        ismatch(r"[^0-9\s]",x) ? throw(SLOTERROR) : return Month(x)
    elseif slot.letter == 'u'
        return Month(MONTHTOVALUEABBR[locale][lowercase(x)])
    else
        return Month(MONTHTOVALUE[locale][lowercase(x)])
    end
end
slotparse(slot::Slot{Millisecond},x,locale) = !ismatch(r"[^0-9\s]",x) ? slot.parser(Base.parse(Float64,"."*x)*1000.0) : throw(SLOTERROR)
slotparse(slot::Slot{DayOfWeekSlot},x,locale) = nothing

function getslot(x,slot::DelimitedSlot,locale,cursor)
    endind = first(search(x,slot.transition,nextind(x,cursor)))
    if endind == 0 # we didn't find the next delimiter
        s = x[cursor:end]
        index = endof(x)+1
    else
        s = x[cursor:(endind-1)]
        index = nextind(x,endind)
    end
    return index, slotparse(slot,s,locale)
end
getslot(x,slot,locale,cursor) = (cursor+slot.width, slotparse(slot,x[cursor:(cursor+slot.width-1)], locale))

function parse(x::AbstractString,df::DateFormat)
    x = strip(x)
    startswith(x, df.prefix) && (x = replace(x, df.prefix, "", 1))
    isempty(x) && throw(ArgumentError("Cannot parse empty format string"))
    if isa(df.slots[1], DelimitedSlot) && first(search(x,df.slots[1].transition)) == 0
        throw(ArgumentError("Delimiter mismatch. Couldn't find first delimiter, \"$(df.slots[1].transition)\", in date string"))
    end
    periods = Period[]
    extra = Any[]  # Supports custom slot types such as TimeZone
    cursor = 1
    for slot in df.slots
        cursor, pe = getslot(x,slot,df.locale,cursor)
        pe !== nothing && (isa(pe,Period) ? push!(periods,pe) : push!(extra,pe))
        cursor > endof(x) && break
    end
    sort!(periods,rev=true,lt=periodisless)
    if isempty(extra)
        return periods
    else
        return vcat(periods, extra)
    end
end

slotformat(slot,dt,locale) = lpad(string(value(slot.parser(dt))),slot.width,"0")
function slotformat(slot::Slot{Year},dt,locale)
    s = lpad(string(value(slot.parser(dt))),slot.width,"0")
    if slot.letter == 'y'
        return s[(end-slot.width+1):end]  # Truncate the year
    else # == 'Y'
        return s
    end
end
function slotformat(slot::Slot{Month},dt,locale)
    if slot.letter == 'm'
        return lpad(month(dt),slot.width,"0")
    elseif slot.letter == 'u'
        return VALUETOMONTHABBR[locale][month(dt)]
    else
        return VALUETOMONTH[locale][month(dt)]
    end
end
function slotformat(slot::Slot{DayOfWeekSlot},dt,locale)
    if slot.letter == 'e'
        return VALUETODAYOFWEEKABBR[locale][dayofweek(dt)]
    else # == 'E'
        return VALUETODAYOFWEEK[locale][dayofweek(dt)]
    end
end
slotformat(slot::Slot{Millisecond},dt,locale) = rpad(string(millisecond(dt)/1000.0)[3:end], slot.width, "0")

function format(dt::TimeType,df::DateFormat)
    f = df.prefix
    for slot in df.slots
        f *= slotformat(slot,dt,df.locale)
        if isa(slot, DelimitedSlot)
            f *= isa(slot.transition, AbstractString) ? slot.transition : ""
        end
    end
    return f
end

# UI
const ISODateTimeFormat = DateFormat("yyyy-mm-dd\\THH:MM:SS.s")
const ISODateFormat = DateFormat("yyyy-mm-dd")
const RFC1123Format = DateFormat("e, dd u yyyy HH:MM:SS")

"""
    DateTime(dt::AbstractString, format::AbstractString; locale="english") -> DateTime

Construct a `DateTime` by parsing the `dt` date string following the pattern given in
the `format` string. The following character codes can be used to construct the `format`
string:

| Code       | Matches   | Comment                                                      |
|:-----------|:----------|:-------------------------------------------------------------|
| `y`        | 1996, 96  | Returns year of 1996, 0096                                   |
| `Y`        | 1996, 96  | Returns year of 1996, 0096. Equivalent to `y`                |
| `m`        | 1, 01     | Matches 1 or 2-digit months                                  |
| `u`        | Jan       | Matches abbreviated months according to the `locale` keyword |
| `U`        | January   | Matches full month names according to the `locale` keyword   |
| `d`        | 1, 01     | Matches 1 or 2-digit days                                    |
| `H`        | 00        | Matches hours                                                |
| `M`        | 00        | Matches minutes                                              |
| `S`        | 00        | Matches seconds                                              |
| `s`        | .500      | Matches milliseconds                                         |
| `e`        | Mon, Tues | Matches abbreviated days of the week                         |
| `E`        | Monday    | Matches full name days of the week                           |
| `yyyymmdd` | 19960101  | Matches fixed-width year, month, and day                     |

Characters not listed above are normally treated as delimiters between date and time slots.
For example a `dt` string of "1996-01-15T00:00:00.0" would have a `format` string like
"y-m-dTH:M:S.s". If you need to use a code character as a delimiter you can escape it using
backslash. The date "1995y01m" would have the format "y\\ym\\m".
"""
DateTime(dt::AbstractString,format::AbstractString;locale::AbstractString="english") = DateTime(dt,DateFormat(format,locale))

"""
    DateTime(dt::AbstractString, df::DateFormat) -> DateTime

Construct a `DateTime` by parsing the `dt` date string following the pattern given in
the [`DateFormat`](:func:`Dates.DateFormat`) object. Similar to
`DateTime(::AbstractString, ::AbstractString)` but more efficient when repeatedly parsing
similarly formatted date strings with a pre-created `DateFormat` object.
"""
DateTime(dt::AbstractString,df::DateFormat=ISODateTimeFormat) = DateTime(parse(dt,df)...)

"""
    Date(dt::AbstractString, format::AbstractString; locale="english") -> Date

Construct a `Date` object by parsing a `dt` date string following the pattern given in the
`format` string. Follows the same conventions as
`DateTime(::AbstractString, ::AbstractString)`.
"""
Date(dt::AbstractString,format::AbstractString;locale::AbstractString="english") = Date(dt,DateFormat(format,locale))

"""
    Date(dt::AbstractString, df::DateFormat) -> Date

Parse a date from a date string `dt` using a `DateFormat` object `df`.
"""
Date(dt::AbstractString,df::DateFormat=ISODateFormat) = Date(parse(dt,df)...)

"""
    format(dt::TimeType, format::AbstractString; locale="english") -> AbstractString

Construct a string by using a `TimeType` object and applying the provided `format`. The
following character codes can be used to construct the `format` string:

| Code       | Examples  | Comment                                                      |
|:-----------|:----------|:-------------------------------------------------------------|
| `y`        | 6         | Numeric year with a fixed width                              |
| `Y`        | 1996      | Numeric year with a minimum width                            |
| `m`        | 1, 12     | Numeric month with a minimum width                           |
| `u`        | Jan       | Month name shortened to 3-chars according to the `locale`    |
| `U`        | January   | Full month name according to the `locale` keyword            |
| `d`        | 1, 31     | Day of the month with a minimum width                        |
| `H`        | 0, 23     | Hour (24-hour clock) with a minimum width                    |
| `M`        | 0, 59     | Minute with a minimum width                                  |
| `S`        | 0, 59     | Second with a minimum width                                  |
| `s`        | 000, 500  | Millisecond with a minimum width of 3                        |
| `e`        | Mon, Tue  | Abbreviated days of the week                                 |
| `E`        | Monday    | Full day of week name                                        |

The number of sequential code characters indicate the width of the code. A format of
`yyyy-mm` specifies that the code `y` should have a width of four while `m` a width of two.
Codes that yield numeric digits have an associated mode: fixed-width or minimum-width.
The fixed-width mode left-pads the value with zeros when it is shorter than the specified
width and truncates the value when longer. Minimum-width mode works the same as fixed-width
except that it does not truncate values longer than the width.

When creating a `format` you can use any non-code characters as a separator. For example to
generate the string "1996-01-15T00:00:00" you could use `format`: "yyyy-mm-ddTHH:MM:SS".
Note that if you need to use a code character as a literal you can use the escape character
backslash. The string "1996y01m" can be produced with the format "yyyy\\ymm\\m".
"""
format(dt::TimeType,f::AbstractString;locale::AbstractString="english") = format(dt,DateFormat(f,locale))

# vectorized
DateTime{T<:AbstractString}(Y::AbstractArray{T},format::AbstractString;locale::AbstractString="english") = DateTime(Y,DateFormat(format,locale))
function DateTime{T<:AbstractString}(Y::AbstractArray{T},df::DateFormat=ISODateTimeFormat)
    return reshape(DateTime[DateTime(parse(y,df)...) for y in Y], size(Y))
end
Date{T<:AbstractString}(Y::AbstractArray{T},format::AbstractString;locale::AbstractString="english") = Date(Y,DateFormat(format,locale))
function Date{T<:AbstractString}(Y::AbstractArray{T},df::DateFormat=ISODateFormat)
    return reshape(Date[Date(parse(y,df)...) for y in Y], size(Y))
end

format{T<:TimeType}(Y::AbstractArray{T},format::AbstractString;locale::AbstractString="english") = Dates.format(Y,DateFormat(format,locale))
function format(Y::AbstractArray{Date},df::DateFormat=ISODateFormat)
    return reshape([Dates.format(y,df) for y in Y], size(Y))
end
function format(Y::AbstractArray{DateTime},df::DateFormat=ISODateTimeFormat)
    return reshape([Dates.format(y,df) for y in Y], size(Y))
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module CellwiseValues

using Reexport

include("Helpers.jl")
@reexport using CellwiseValues.Helpers

include("CachedArrays.jl")
@reexport using CellwiseValues.CachedArrays

include("CellValues.jl")
@reexport using CellwiseValues.CellValues

include("CellNumbers.jl")
@reexport using CellwiseValues.CellNumbers

include("CellArrays.jl")
@reexport using CellwiseValues.CellArrays

include("Maps.jl")
@reexport using CellwiseValues.Maps

include("CellMaps.jl")
@reexport using CellwiseValues.CellMaps

include("Gallery.jl")
@reexport using CellwiseValues.Gallery

include("Kernels.jl")
@reexport using CellwiseValues.Kernels

include("CellNumberApply.jl")
@reexport using CellwiseValues.CellNumberApply

include("CellArrayApply.jl")
@reexport using CellwiseValues.CellArrayApply

include("MapApply.jl")
@reexport using CellwiseValues.MapApply

include("CellMapApply.jl")
@reexport using CellwiseValues.CellMapApply

end # module


text:
using LinearAlgebra
using Statistics
include("../module2/hello.jl")

function do_some_stuff(x:: Matrix)
    return det(x)
end

function do_other_stuff(x:: Matrix)
    return x + x
end

function do_statistics_stuff(x:: Matrix)
    return mean(x)
end

println(hello("Jackson"))

display(do_some_stuff([1 2; 3 4]))
println("")
display(do_other_stuff([1 2 3 4]))
println("")
display(do_statistics_stuff([1 2 3 4]))


text:
using PFKernels
using PFKernels.PowerSystem
using PFKernels.ForwardDiff
using UnicodePlots

PFKernels.CUDA.device!(1)

AT = PFKernels.AT

# Matlab-like residual function
function residual(V, Ybus, Sbus, pv, pq)
    # form mismatch vector
    mis = V .* conj(Ybus * V) - Sbus
    # form residual vector
    F = [real(mis[pv])
        real(mis[pq])
        imag(mis[pq]) ]
    return F
end

function recursive_value(dual::ForwardDiff.Dual{Nothing, T, N}) where {T,N}
    if T == Float64
        return ForwardDiff.value(dual)
    else
        recursive_value(dual.value)
    end
end

# Driver for dual numbers
function residual!(::Type{DT}, F::Vector{Float64},
            v_m::Vector{Float64}, v_a::Vector{Float64},
            ybus_re::PFKernels.Spmat{Vector{Int64}, Vector{Float64}},
            ybus_im::PFKernels.Spmat{Vector{Int64}, Vector{Float64}},
            pinj::Vector{Float64}, qinj::Vector{Float64},
            pv::Vector{Int64}, pq::Vector{Int64},
            nbus::Int64) where DT
    npv = length(pv)
    npq = length(pq)
    adF = Vector{DT}(undef, length(F))
    adv_m = Vector{DT}(undef, length(v_m))
    adv_a = Vector{DT}(undef, length(v_a))
    adybus_re_nzval = Vector{DT}(undef, length(ybus_re.nzval))
    adybus_re_colptr = Vector{Int64}(undef, length(ybus_re.colptr))
    adybus_re_rowval = Vector{Int64}(undef, length(ybus_re.rowval))
    adybus_im_nzval = Vector{DT}(undef, length(ybus_im.nzval))
    adybus_im_colptr = Vector{Int64}(undef, length(ybus_im.colptr))
    adybus_im_rowval = Vector{Int64}(undef, length(ybus_im.rowval))
    adpinj = Vector{DT}(undef, length(pinj))
    adqinj = Vector{DT}(undef, length(qinj))
    adpv = Vector{DT}(undef, length(pv))
    adpq = Vector{DT}(undef, length(pq))

    adF .= F
    adv_m .= v_m
    adv_a .= v_a
    adybus_re_nzval .= ybus_re.nzval
    adybus_im_nzval .= ybus_im.nzval
    adpinj .= pinj
    adqinj .= qinj
    adpv .= pv
    adpq .= pq

    dF = AT(adF)
    dv_m = AT(adv_m)
    dv_a = AT(adv_a)
    dybus_re_nzval = AT(adybus_re_nzval)
    dybus_re_colptr = AT(ybus_re.colptr)
    dybus_re_rowval = AT(ybus_re.rowval)
    dybus_im_nzval = AT(adybus_im_nzval)
    dybus_im_colptr = AT(ybus_im.colptr)
    dybus_im_rowval = AT(ybus_im.rowval)
    dpinj = AT(adpinj)
    dqinj = AT(adqinj)
    dpv = AT(pv)
    dpq = AT(pq)

    PFKernels.kernel!(PFKernels.backend, dF, dv_m, dv_a,
                    dybus_re_nzval, dybus_re_colptr, dybus_re_rowval,
                    dybus_im_nzval, dybus_im_colptr, dybus_im_rowval,
                    dpinj, dqinj, dpv, dpq, nbus)
    F .= recursive_value.(dF)
end

const PS = PowerSystem

# Inject GPU kernels
PFKernels.include(joinpath(dirname(@__FILE__), "kernels.jl"))

# Get the powerflow stuff out of the way
local_case = "case14.raw"
V, Ybus, Sbus, Sload, ref, pv, pq = PFKernels.loaddata(local_case)
npv = size(pv, 1)
npq = size(pq, 1)
F = zeros(Float64, npv + 2*npq)
Fnew = zeros(Float64, npv + 2*npq)
F♯ = residual(V, Ybus, Sbus, pv, pq)
Vm = abs.(V)
Va = angle.(V)
ybus_re, ybus_im = PFKernels.Spmat{Vector{Int}, Vector{Float64}}(Ybus)
pbus = real(Sbus)
qbus = imag(Sbus)

# Let's run some differentiated kernels
t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where {N}
t2s{N} = ForwardDiff.Dual{Nothing,t1s{N}, N} where {N}
t3s{N} = ForwardDiff.Dual{Nothing,t2s{N}, N} where {N}

N = 2
F1 = similar(F)
F2 = similar(F)
F3 = similar(F)
F1 .= 0.0
residual!(t1s{N}, F1, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
F2 .= 0.0
@time residual!(t2s{N}, F2, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
F3 .= 0.0
residual!(t3s{N}, F3, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
@assert(F1 ≈ F♯)
@assert(F2 ≈ F♯)
@assert(F3 ≈ F♯)

# And now we crank up the partials
# times = Vector{Float64}(undef, 10)
const T1 = t2s{10}
const T2 = t2s{20}
ctimes = zeros(10)
ctimes[1] = @elapsed residual!(T1, F2, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
ctimes[2] = @elapsed residual!(T2, F2, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
@show ctimes
rtimes = zeros(10)
rtimes[1] = @elapsed residual!(T1, F2, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
rtimes[2] = @elapsed residual!(T2, F2, Vm, Va,
    ybus_re, ybus_im,
    pbus, qbus, pv, pq, length(V))
@show rtimes


text:
using LikelihoodProfiler, Plots

# functions with zero derivative at one of the endpoints
f1(x) = 3x[1]^4 - 6x[1]^3 + 5
f2(x) = -3(x[1]-5)^4 + 5(x[1]-5)^2 + 5
f3(x) = 3x[1]^3 - 6x[1]^2 + 5
f4(x) = 3x[1]^4 - 6x[1]^3 -  x[1]^2 +  x[1] + 5
f5(x) = 3x[1]^4 - 6x[1]^3 - 3x[1]^2 + 4x[1] + 5

# functions non-differentiable at one of the endpoints
function f6(x)
  if (1.0 <= x[1] < 5.0)
    return ((x[1]-3)^2 + 5.0)
  else
    return abs(4.5x[1] - 13.5)
  end
end

res1 = get_interval(
  [1.5],
  1,
  f1,
  :CICO_ONE_PASS;
  loss_crit = 5.,
  silent = true
)

res2 = get_interval(
  [5.0],
  1,
  f2,
  :CICO_ONE_PASS;
  loss_crit = 7.,
  silent = true
)

res3 = get_interval(
  [1.5],
  1,
  f3,
  :CICO_ONE_PASS;
  loss_crit = 5.,
  silent = true
)

res4 = get_interval(
  [1.5],
  1,
  f4,
  :CICO_ONE_PASS;
  loss_crit = 5.,
  silent = true
)

res6 = get_interval(
  [3.0],
  1,
  f6,
  :CICO_ONE_PASS;
  loss_crit = 9.,
  silent = true
)


text:
export MachOLoadDylibCmd, MachOIdDylibCmd

export dylib_name, dylib_timestamp, dylib_version, dylib_compatibility

"""
    MachODylibStub

MachO Load Commands have an annoying and nasty habit of shimming strings into
difficult-to-reach places, such as within the header itself.  Because the
strings are variable-length, it becomes difficult to go back and find the
strings after the fact; it's better to suck it into memory immediately.  This
doesn't lend itself to `StructIO`-based IO very nicely, and so we hack our way
around the problem by creating the `Stub` which contains the data which can be
read in by `StructIO` nicely, and then immediately expanding and reading in the
`LCStr` afterward within the "true" `MachOLoadDylibCmd`.
"""
@io struct MachODylibStub
    # Note this `offset` is technically a part of the `dylib` union, but we
    # don't bother to split it out into its own `struct`.
    name_offset::UInt32
    timestamp::UInt32
    current_version::UInt32
    compatibilty::UInt32
end

"""
    MachOLoadDylibCmd

The Load Command that gives information about a dylib that must be loaded for
this Mach-O file to link properly.  The API for this laod command is given as
follows:

### Creation:
  - MachOLoadCmd()

### Accessors:
  - dylib_name()
  - dylib_timestamp()
  - dylib_version()
  - dylib_compatibility()
"""
struct MachOLoadDylibCmd{H <: MachOHandle} <: MachOLoadCmd{H}
    stub::MachODylibStub
    name::String
end

"""
    MachOIdDylibCmd

MachO dylibs have an "identity", analogous to the SONAME of an ELF shared
library.  This Load Command tells the linker about the identity of this
particular shared library.  The API of this object is identical to that of the
`MachOLoadDylibCmd` object, as they are intrinsically identical, it is the
semantic meaning behind the data that changes, and so this object is all but a
typealias for `MachOLoadDylibCmd`.
"""
struct MachOIdDylibCmd{H <: MachOHandle} <: MachOLoadCmd{H}
    stub::MachODylibStub
    name::String
end

"""
    MachODylibCmd

We define some operations that work on either a `MachOLoadDylibCmd` or
`MachOIdDylibCmd` load command, and to make it all pretty and shiny we define
this union type that allows our operations to be generic.
"""
const DylibCmd = Union{MachOLoadDylibCmd, MachOIdDylibCmd}

"""
    MachOLoadCmd(oh::MachOHeader, CT::Type{MachOLoadDylibCmd}, header)
    MachOLoadCmd(oh::MachOHeader, CT::Type{MachOIdDylibCmd}, header)

Manual override to unpack a `MachOLoadDylibCmd` from a MachO file and
immediately read in the associated name string.  See `MachODylibStub` for more.
"""
function MachOLoadCmd(oh::H, CT::Type{CTT}, header) where {CTT <: DylibCmd} where {H <: MachOHandle}
    # First, unpack the stub
    stub = unpack(oh, MachODylibStub)

    # Next, read in the string, telling it to not read in more than the total
    # load command size, and not before the DylibStub object finishes.
    offset = stub.name_offset - sizeof(MachODylibStub) - sizeof(MachOLoadCmdHeader{H})
    name = unpack_lcstr(oh, offset, load_cmd_size(header))

    # Finally, return our DylibCmd
    return CT(stub, name)
end

"""
    dylib_name(cmd)

Return the name of the dylib referred to by the given Load Command.
"""
dylib_name(cmd::DylibCmd) = cmd.name

"""
    dylib_timestamp(cmd)

Return the build timestamp of the dylib referred to by the given Load Command.
"""
dylib_timestamp(cmd::DylibCmd) = cmd.stub.timestamp

"""
    dylib_version(cmd)

Return the version of the dylib referred to by the given Load Command.
"""
dylib_version(cmd::DylibCmd) = cmd.stub.current_version

"""
    dylib_compatibility(cmd)

Return the compatibility version of the dylib referred to by the given Load
Command.
"""
dylib_compatibilty(cmd::DylibCmd) = cmd.stub.compatibilty

@derefmethod dylib_name(cmd::MachOLoadCmdRef{T}) where {T <: DylibCmd}
@derefmethod dylib_timestamp(cmd::MachOLoadCmdRef{T}) where {T <: DylibCmd}
@derefmethod dylib_version(cmd::MachOLoadCmdRef{T}) where {T <: DylibCmd}
@derefmethod dylib_compatibilty(cmd::MachOLoadCmdRef{T}) where {T <: DylibCmd}

"""
    MachORPathCmd

The Load Command that holds the RPATH of this object.
"""
struct MachORPathCmd{H <: MachOHandle} <: MachOLoadCmd{H}
    rpath_offset::UInt32
    rpath::String
end

"""
    MachOLoadCmd(oh::MachOHeader, CT::Type{MachORPathCmd}, header)

Manual override to unpack a `MachORPathCmd` from a MachO file and
immediately read in the associated rpath string.
"""
function MachOLoadCmd(oh::H, CT::Type{CTT}, header) where {CTT <: MachORPathCmd} where {H <: MachOHandle}
    rpath_offset = read(oh, UInt32)

    # Next, read in the string, telling it to not read in more than the total
    # load command size, and not before the DylibStub object finishes.
    offset = rpath_offset - sizeof(UInt32) - sizeof(MachOLoadCmdHeader{H})
    name = unpack_lcstr(oh, offset, load_cmd_size(header))

    # Finally, return our RPathCmd
    return CT(rpath_offset, name)
end

rpath(cmd::MachORPathCmd) = cmd.rpath
@derefmethod rpath(cmd::MachOLoadCmdRef{H,MachORPathCmd{H}}) where {H <: MachOHandle}


text:
# Scalar identity ITensor
struct OneITensor end

(::OneITensor * A::ITensor) = A
(A::ITensor * ::OneITensor) = A

inds(::OneITensor) = ()

abstract type AbstractProjMPO end

# This is to help with generic promote_type code
# in eltype(::AbstractProjMPO)
eltype(::OneITensor) = Bool
dim(::OneITensor) = 1
isoneitensor(::OneITensor) = true
isoneitensor(::ITensor) = false

"""
    nsite(P::ProjMPO)

Retrieve the number of unprojected (open)
site indices of the ProjMPO object `P`
"""
nsite(P::AbstractProjMPO) = P.nsite

# The range of center sites
site_range(P::AbstractProjMPO) = (P.lpos + 1):(P.rpos - 1)

"""
    length(P::ProjMPO)

The length of a ProjMPO is the same as
the length of the MPO used to construct it
"""
Base.length(P::AbstractProjMPO) = length(P.H)

function lproj(P::AbstractProjMPO)::Union{ITensor,OneITensor}
  (P.lpos <= 0) && return OneITensor()
  return P.LR[P.lpos]
end

function rproj(P::AbstractProjMPO)::Union{ITensor,OneITensor}
  (P.rpos >= length(P) + 1) && return OneITensor()
  return P.LR[P.rpos]
end

"""
    product(P::ProjMPO,v::ITensor)::ITensor

    (P::ProjMPO)(v::ITensor)

Efficiently multiply the ProjMPO `P`
by an ITensor `v` in the sense that the
ProjMPO is a generalized square matrix
or linear operator and `v` is a generalized
vector in the space where it acts. The
returned ITensor will have the same indices
as `v`. The operator overload `P(v)` is
shorthand for `product(P,v)`.
"""
function product(P::AbstractProjMPO, v::ITensor)::ITensor
  itensor_map = Union{ITensor,OneITensor}[lproj(P)]
  append!(itensor_map, P.H[site_range(P)])
  push!(itensor_map, rproj(P))

  # Reverse the contraction order of the map if
  # the first tensor is a scalar (for example we
  # are at the left edge of the system)
  if dim(first(itensor_map)) == 1
    reverse!(itensor_map)
  end

  # Apply the map
  Hv = v
  for it in itensor_map
    Hv *= it
  end

  if order(Hv) != order(v)
    error(
      string(
        "The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, ",
        "this is probably due to an index mismatch.\nCommon reasons for this error: \n",
        "(1) You are trying to multiply the ProjMPO with the $(nsite(P))-site wave-function at the wrong position.\n",
        "(2) `orthognalize!` was called, changing the MPS without updating the ProjMPO.\n\n",
        "P*v inds: $(inds(Hv)) \n\n",
        "v inds: $(inds(v))",
      ),
    )
  end

  return noprime(Hv)
end

(P::AbstractProjMPO)(v::ITensor) = product(P, v)

"""
    eltype(P::ProjMPO)

Deduce the element type (such as Float64
or ComplexF64) of the tensors in the ProjMPO
`P`.
"""
function Base.eltype(P::AbstractProjMPO)::Type
  ElType = eltype(lproj(P))
  for j in site_range(P)
    ElType = promote_type(ElType, eltype(P.H[j]))
  end
  return promote_type(ElType, eltype(rproj(P)))
end

"""
    size(P::ProjMPO)

The size of a ProjMPO are its dimensions
`(d,d)` when viewed as a matrix or linear operator
acting on a space of dimension `d`.

For example, if a ProjMPO maps from a space with
indices `(a,s1,s2,b)` to the space `(a',s1',s2',b')`
then the size is `(d,d)` where
`d = dim(a)*dim(s1)*dim(s1)*dim(b)`
"""
function Base.size(P::AbstractProjMPO)::Tuple{Int,Int}
  d = 1
  for i in inds(lproj(P))
    plev(i) > 0 && (d *= dim(i))
  end
  for j in site_range(P)
    for i in inds(P.H[j])
      plev(i) > 0 && (d *= dim(i))
    end
  end
  for i in inds(rproj(P))
    plev(i) > 0 && (d *= dim(i))
  end
  return (d, d)
end

function _makeL!(P::AbstractProjMPO, psi::MPS, k::Int)::Union{ITensor,Nothing}
  # Save the last `L` that is made to help with caching
  # for DiskProjMPO
  ll = P.lpos
  if ll ≥ k
    # Special case when nothing has to be done.
    # Still need to change the position if lproj is
    # being moved backward.
    P.lpos = k
    return nothing
  end
  # Make sure ll is at least 0 for the generic logic below
  ll = max(ll, 0)
  L = lproj(P)
  while ll < k
    L = L * psi[ll + 1] * P.H[ll + 1] * dag(prime(psi[ll + 1]))
    P.LR[ll + 1] = L
    ll += 1
  end
  # Needed when moving lproj backward.
  P.lpos = k
  return L
end

function makeL!(P::AbstractProjMPO, psi::MPS, k::Int)
  _makeL!(P, psi, k)
  return P
end

function _makeR!(P::AbstractProjMPO, psi::MPS, k::Int)::Union{ITensor,Nothing}
  # Save the last `R` that is made to help with caching
  # for DiskProjMPO
  rl = P.rpos
  if rl ≤ k
    # Special case when nothing has to be done.
    # Still need to change the position if rproj is
    # being moved backward.
    P.rpos = k
    return nothing
  end
  N = length(P.H)
  # Make sure rl is no bigger than `N + 1` for the generic logic below
  rl = min(rl, N + 1)
  R = rproj(P)
  while rl > k
    R = R * psi[rl - 1] * P.H[rl - 1] * dag(prime(psi[rl - 1]))
    P.LR[rl - 1] = R
    rl -= 1
  end
  P.rpos = k
  return R
end

function makeR!(P::AbstractProjMPO, psi::MPS, k::Int)
  _makeR!(P, psi, k)
  return P
end

"""
    position!(P::ProjMPO, psi::MPS, pos::Int)

Given an MPS `psi`, shift the projection of the
MPO represented by the ProjMPO `P` such that
the set of unprojected sites begins with site `pos`.
This operation efficiently reuses previous projections
of the MPO on sites that have already been projected.
The MPS `psi` must have compatible bond indices with
the previous projected MPO tensors for this
operation to succeed.
"""
function position!(P::AbstractProjMPO, psi::MPS, pos::Int)
  makeL!(P, psi, pos - 1)
  makeR!(P, psi, pos + nsite(P))
  return P
end

"""
    noiseterm(P::ProjMPO,
              phi::ITensor,
              ortho::String)

Return a "noise term" or density matrix perturbation
ITensor as proposed in Phys. Rev. B 72, 180403 for aiding
convergence of DMRG calculations. The ITensor `phi`
is the contracted product of MPS tensors acted on by the
ProjMPO `P`, and `ortho` is a String which can take
the values `"left"` or `"right"` depending on the
sweeping direction of the DMRG calculation.
"""
function noiseterm(P::AbstractProjMPO, phi::ITensor, ortho::String)::ITensor
  if nsite(P) != 2
    error("noise term only defined for 2-site ProjMPO")
  end

  site_range_P = site_range(P)
  if ortho == "left"
    AL = P.H[first(site_range_P)]
    AL = lproj(P) * AL
    nt = AL * phi
  elseif ortho == "right"
    AR = P.H[last(site_range_P)]
    AR = AR * rproj(P)
    nt = phi * AR
  else
    error("In noiseterm, got ortho = $ortho, only supports `left` and `right`")
  end
  nt = nt * dag(noprime(nt))

  return nt
end

function checkflux(P::AbstractProjMPO)
  checkflux(P.H)
  for n in length(P.LR)
    if isassigned(P.LR, n)
      checkflux(P.LR[n])
    end
  end
  return nothing
end


text:
"""
    PiecewiseLegendrePoly <: Function

Piecewise Legendre polynomial.

Models a function on the interval ``[xmin, xmax]`` as a set of segments on the
intervals ``S[i] = [a[i], a[i+1]]``, where on each interval the function
is expanded in scaled Legendre polynomials.
"""
struct PiecewiseLegendrePoly{T} <: Function
    polyorder::Int
    xmin::T
    xmax::T

    knots::Vector{T}
    Δx::Vector{T}
    data::Matrix{T}
    symm::Int
    l::Int

    xm::Vector{T}
    inv_xs::Vector{T}
    norm::Vector{T}

    function PiecewiseLegendrePoly(
        polyorder::Integer, xmin::AbstractFloat, xmax::AbstractFloat, knots::AbstractVector,
        Δx::AbstractVector, data::AbstractMatrix, symm::Integer, l::Integer,
        xm::AbstractVector, inv_xs::AbstractVector, norm::AbstractVector,
    )
        !any(isnan, data) || error("data contains NaN")
        issorted(knots) || error("knots must be monotonically increasing")
        Δx ≈ diff(knots) || error("Δx must work with knots")
        return new{eltype(knots)}(
            polyorder, xmin, xmax, knots, Δx, data, symm, l, xm, inv_xs, norm
        )
    end
end

function PiecewiseLegendrePoly(data, p::PiecewiseLegendrePoly; symm=0)
    return PiecewiseLegendrePoly(p.polyorder, p.xmin, p.xmax, p.knots, p.Δx,
        data, symm, p.l, p.xm, p.inv_xs, p.norm)
end

function PiecewiseLegendrePoly(
    data::Matrix, knots::Vector, l::Integer; Δx=diff(knots), symm=0
)
    polyorder, nsegments = size(data)
    size(knots) == (nsegments + 1,) || error("Invalid knots array")
    xm = @views (knots[1:(end - 1)] + knots[2:end]) / 2
    inv_xs = 2 ./ Δx
    norm = sqrt.(inv_xs)

    return PiecewiseLegendrePoly(polyorder, first(knots), last(knots), knots,
        Δx, data, symm, l, xm, inv_xs, norm)
end

Base.size(::PiecewiseLegendrePoly) = ()

function Base.show(io::IO, p::PiecewiseLegendrePoly)
    return print(io, "$(typeof(p)): xmin=$(p.xmin), xmax=$(p.xmax)")
end

(poly::PiecewiseLegendrePoly)(x) = _evaluate(poly, x)

function _evaluate(poly::PiecewiseLegendrePoly, x::Number)
    i, x̃ = _split(poly, x)
    return legval(x̃, view(poly.data, :, i)) * poly.norm[i]
end

function _evaluate(poly::PiecewiseLegendrePoly, xs::AbstractVector)
    return poly.(xs)
end

"""
    overlap(poly::PiecewiseLegendrePoly, f)

Evaluate overlap integral of `poly` with arbitrary function `f`.

Given the function `f`, evaluate the integral::

    ∫ dx * f(x) * poly(x)

using adaptive Gauss-Legendre quadrature.
"""
function overlap(
    poly::PiecewiseLegendrePoly{T}, f; rtol=eps(T), return_error=false, maxevals=10^4
) where {T}
    int_result, int_error = quadgk(
        x -> poly(x) * f(x), poly.knots...; rtol, order=10, maxevals
    )
    if return_error
        return int_result, int_error
    else
        return int_result
    end
end

"""
    deriv(poly)

Get polynomial for the derivative.
"""
function deriv(poly::PiecewiseLegendrePoly)
    ddata = legder(poly.data)

    scale = poly.inv_xs
    ddata .*= transpose(scale)
    return PiecewiseLegendrePoly(ddata, poly; symm=-poly.symm)
end

"""
    roots(poly)

Find all roots of the piecewise polynomial `poly`.
"""
function roots(poly::PiecewiseLegendrePoly; tol=1e-10, alpha=2)
    grid = poly.knots
    grid = _refine_grid(grid, alpha)
    return find_all(poly, grid)
end

function check_domain(poly::PiecewiseLegendrePoly, x)
    poly.xmin ≤ x ≤ poly.xmax || throw(DomainError(x, "x is outside the domain"))
    return true
end

"""
    _split(poly, x)

Split segment.

Find segment of poly's domain that covers `x`.
"""
function _split(poly, x::Number)
    @boundscheck check_domain(poly, x)

    i = max(searchsortedlast(poly.knots, x; lt=≤), 1)
    x̃ = x - poly.xm[i]
    x̃ *= poly.inv_xs[i]
    return i, x̃
end

function scale(poly::PiecewiseLegendrePoly, factor)
    return PiecewiseLegendrePoly(
        poly.data * factor, poly.knots, poly.l; Δx=poly.Δx, symm=poly.symm
    )
end

###########################
## PiecewiseLegendrePolyVector ##
###########################

"""
    PiecewiseLegendrePolyVector{T}

Alias for `Vector{PiecewiseLegendrePoly{T}}`.
"""
const PiecewiseLegendrePolyVector{T} = Vector{PiecewiseLegendrePoly{T}}

function PiecewiseLegendrePolyVector(
    data::AbstractArray{T,3}, knots::Vector{T}; symm=zeros(Int, size(data, 3))
) where {T<:AbstractFloat}
    return [
        PiecewiseLegendrePoly(data[:, :, i], knots, i - 1; symm=symm[i]) for
        i in axes(data, 3)
    ]
end

function PiecewiseLegendrePolyVector(
    polys::PiecewiseLegendrePolyVector, knots::AbstractVector; Δx=diff(knots), symm=0
)
    length(polys) == length(symm) ||
        throw(DimensionMismatch("Sizes of polys and symm don't match"))

    polys_new = similar(polys)
    for i in eachindex(polys)
        polys_new[i] = PiecewiseLegendrePoly(
            polys[i].data, knots, polys[i].l; Δx, symm=symm[i]
        )
    end
    return polys_new
end

function PiecewiseLegendrePolyVector(
    data::AbstractArray{T,3}, polys::PiecewiseLegendrePolyVector
) where {T}
    size(data, 3) == length(polys) ||
        throw(DimensionMismatch("Sizes of data and polys don't match"))

    polys_new = copy(polys)
    for i in eachindex(polys)
        polys_new[i].data .= data[:, :, i]
    end
    return polys_new
end

(polys::PiecewiseLegendrePolyVector)(x) = map(poly -> poly(x), polys)
function (polys::PiecewiseLegendrePolyVector)(x::AbstractArray)
    return reshape(reduce(vcat, polys.(x)), (size(polys)..., size(x)...))
end

function Base.getproperty(polys::PiecewiseLegendrePolyVector, sym::Symbol)
    if sym ∈ (:xmin, :xmax, :knots, :Δx, :polyorder, :nsegments, :xm, :inv_xs, :norm)
        return getproperty(first(polys), sym)
    elseif sym == :symm
        return map(poly -> poly.symm, polys)
        # elseif sym == :data
        #     return cat((poly.data for poly in polys)...; dims=3)
    else
        error("Unknown property $sym")
    end
end

# Backward compatibility
function overlap(
    polys::PiecewiseLegendrePolyVector{T}, f; rtol=eps(T), return_error=false
) where {T}
    return overlap.(polys, f; rtol, return_error)
end

#########################
## PiecewiseLegendreFT ##
#########################

"""
    PowerModel

Model from a high-frequency series expansion::

    A(iω) == sum(A[n] / (iω)^(n+1) for n in 1:N)

where ``iω == i * π/2 * wn`` is a reduced imaginary frequency, i.e.,
``wn`` is an odd/even number for fermionic/bosonic frequencies.
"""
struct PowerModel{T<:AbstractFloat}
    moments::Vector{T}
end

const DEFAULT_GRID = [
    range(0; length=2^6)
    trunc.(Int, exp2.(range(6, 25; length=16 * (25 - 6) + 1)))
]

"""
    PiecewiseLegendreFT <: Function

Fourier transform of a piecewise Legendre polynomial.

For a given frequency index `n`, the Fourier transform of the Legendre
function is defined as:

        p̂(n) == ∫ dx exp(im * π * n * x / (xmax - xmin)) p(x)

The polynomial is continued either periodically (`freq=:even`), in which
case `n` must be even, or antiperiodically (`freq=:odd`), in which case
`n` must be odd.
"""
struct PiecewiseLegendreFT{T} <: Function
    poly::PiecewiseLegendrePoly{T}
    freq::Symbol
    ζ::Int
    n_asymp::T

    model::PowerModel{T}
end

function Base.show(io::IO, p::PiecewiseLegendreFT)
    return print(io, "$(typeof(p))")
end

function PiecewiseLegendreFT(poly, freq=:even, n_asymp=Inf)
    (poly.xmin, poly.xmax) == (-1, 1) || error("Only interval [-1, 1] is supported")
    ζ = Dict(:even => 0, :odd => 1)[freq]
    model = power_model(freq, poly)
    return PiecewiseLegendreFT(poly, freq, ζ, float(n_asymp), model)
end

const PiecewiseLegendreFTVector{T} = Vector{PiecewiseLegendreFT{T}}

function Base.getproperty(polyFTs::PiecewiseLegendreFTVector, sym::Symbol)
    if sym ∈ (:freq, :ζ, :n_asymp)
        return getproperty(first(polyFTs), sym)
    elseif sym == :poly
        return map(poly -> poly.poly, polyFTs)
    else
        error("Unknown property $sym")
    end
end

"""
    (polyFT::PiecewiseLegendreFT)(n)

Obtain Fourier transform of polynomial for given frequency index `n`.
"""
function (polyFT::Union{PiecewiseLegendreFT,PiecewiseLegendreFTVector})(n::Integer)
    n = check_reduced_matsubara(n, polyFT.ζ)

    if abs(n) < polyFT.n_asymp
        return _compute_unl_inner(polyFT.poly, n)
    else
        return giw(polyFT, n)
    end
end

(polyFT::PiecewiseLegendreFT)(n::AbstractArray) = polyFT.(n)
function (polyFTs::PiecewiseLegendreFTVector)(n::AbstractArray)
    return reshape(reduce(vcat, polyFTs.(n)), (size(polyFTs)..., size(n)...))
end

"""
    giw(polyFT, wn)

Return model Green's function for reduced frequencies
"""
function giw(polyFT, wn::Integer)
    iw = im * π / 2 * wn
    iszero(wn) && return zero(inv_iw)
    inv_iw = 1 / iw
    return inv_iw * evalpoly(inv_iw, moments(polyFT))
end

moments(polyFT::PiecewiseLegendreFT) = polyFT.model.moments
function moments(polyFTs::PiecewiseLegendreFTVector)
    n = length(first(polyFTs).model.moments)
    return [[p.model.moments[i] for p in polyFTs] for i in 1:n]
end

Base.firstindex(::PiecewiseLegendreFT) = 1

function hat(poly::PiecewiseLegendrePoly, freq; n_asymp=Inf)
    return PiecewiseLegendreFT(poly, freq, n_asymp)
end

"""
    findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=DEFAULT_GRID)

Obtain extrema of fourier-transformed polynomial.
"""
function findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=DEFAULT_GRID)
    f = _func_for_part(polyFT, part)
    x₀ = _discrete_extrema(f, grid)
    x₀ .= 2x₀ .+ polyFT.ζ

    return _symmetrize_matsubara(x₀)
end

function _func_for_part(polyFT::PiecewiseLegendreFT, part=nothing)
    if isnothing(part)
        parity = polyFT.poly.symm
        if parity == 1
            part = iszero(polyFT.ζ) ? real : imag
        elseif parity == -1
            part = iszero(polyFT.ζ) ? imag : real
        else
            error("Cannot detect parity")
        end
    end
    return n -> part(polyFT(2n + polyFT.ζ))
end

function _discrete_extrema(f::Function, xgrid)
    fx = Float64.(f.(xgrid))
    absfx = abs.(fx)

    # Forward differences: derivativesignchange[i] now means that the secant changes sign
    # fx[i+1]. This means that the extremum is STRICTLY between x[i] and
    # x[i+2]
    gx = diff(fx)
    sgx = signbit.(gx)
    derivativesignchange = sgx[1:(end - 1)] .≠ sgx[2:end]
    derivativesignchange_a = [derivativesignchange; false; false]
    derivativesignchange_b = [false; false; derivativesignchange]

    a = xgrid[derivativesignchange_a]
    b = xgrid[derivativesignchange_b]
    absf_a = absfx[derivativesignchange_a]
    absf_b = absfx[derivativesignchange_b]
    res = _bisect_discr_extremum.(f, a, b, absf_a, absf_b)

    # We consider the outer point to be extremua if there is a decrease
    # in magnitude or a sign change inwards
    sfx = signbit.(fx)
    if absfx[1] > absfx[2] || sfx[1] ≠ sfx[2]
        pushfirst!(res, first(xgrid))
    end
    if absfx[end] > absfx[end - 1] || sfx[end] ≠ sfx[end - 1]
        push!(res, last(xgrid))
    end

    return res
end

function _bisect_discr_extremum(f, a, b, absf_a, absf_b)
    d = b - a
    d ≤ 1 && return absf_a > absf_b ? a : b
    d == 2 && return a + 1

    m = (a + b) ÷ 2
    n = m + 1
    absf_m = abs(f(m))
    absf_n = abs(f(n))
    if absf_m > absf_n
        return _bisect_discr_extremum(f, a, n, absf_a, absf_n)
    else
        return _bisect_discr_extremum(f, m, b, absf_m, absf_b)
    end
end

function _symmetrize_matsubara(x₀)
    issorted(x₀) || error("set of Matsubara points not ordered")
    first(x₀) ≥ 0 || error("points must be non-negative")
    if iszero(first(x₀))
        x₀ = [-reverse(x₀); x₀[2:end]]
    else
        x₀ = [-reverse(x₀); x₀]
    end
    return x₀
end

function derivs(ppoly, x)
    res = [ppoly(x)]
    for _ in 1:(ppoly.polyorder - 1)
        ppoly = deriv(ppoly)
        push!(res, ppoly(x))
    end
    return res
end

function power_moments(stat, deriv_x1, l)
    statsign = (stat == :odd) ? -1 : 1
    mmax = length(deriv_x1)
    coeff_lm = @. ((-1)^(1:mmax) + statsign * (-1)^l) * deriv_x1
    return -statsign / √2 * coeff_lm
end

function power_model(stat, poly)
    deriv_x1 = derivs(poly, 1)
    moments = power_moments(stat, deriv_x1, poly.l)
    return PowerModel(moments)
end

"""
    check_reduced_matsubara(n[, ζ])

Checks that `n` is a reduced Matsubara frequency.

Check that the argument is a reduced Matsubara frequency, which is an
integer obtained by scaling the freqency `ω[n]` as follows:

    β / π * ω[n] == 2n + ζ

Note that this means that instead of a fermionic frequency (`ζ == 1`),
we expect an odd integer, while for a bosonic frequency (`ζ == 0`),
we expect an even one.  If `ζ` is omitted, any one is fine.
"""
check_reduced_matsubara(n::Integer) = n
function check_reduced_matsubara(n::Integer, ζ)
    n & 1 == ζ || throw(DomainError(n, "n has the wrong parity"))
    return n
end

########################
### Helper Functions ###
########################

"""
    _compute_unl_inner(poly, wn)

Compute piecewise Legendre to Matsubara transform.
"""
function _compute_unl_inner(poly::PiecewiseLegendrePoly, wn)
    t_pin = _Pqn(poly, wn)
    return dot(poly.data, transpose(t_pin))
end
function _compute_unl_inner(polys::PiecewiseLegendrePolyVector, wn)
    t_pin = _Pqn(polys, wn)
    return map(p -> dot(p.data, transpose(t_pin)), polys)
end

function _Pqn(poly, wn)
    p = transpose(range(0; length=poly.polyorder))
    wred = π / 2 * wn
    phase_wi = _phase_stable(poly, wn)
    return _get_tnl.(p, wred .* poly.Δx ./ 2) .* phase_wi ./ (√2 .* poly.norm)
end

"""
    _get_tnl(l, w)

Fourier integral of the `l`-th Legendre polynomial::

    Tₗ(ω) == ∫ dx exp(iωx) Pₗ(x)
"""
function _get_tnl(l, w)
    result = 2im^l * sphericalbesselj(l, abs(w))
    return w < 0 ? conj(result) : result
end

# Works like numpy.choices
choose(a, choices) = [choices[a[i]][i] for i in eachindex(a)]

"""
    _shift_xmid(knots, Δx)

Return midpoint relative to the nearest integer plus a shift.

Return the midpoints `xmid` of the segments, as pair `(diff, shift)`,
where shift is in `(0, 1, -1)` and `diff` is a float such that
`xmid == shift + diff` to floating point accuracy.
"""
function _shift_xmid(knots, Δx)
    Δx_half = Δx ./ 2
    xmid_m1 = cumsum(Δx) - Δx_half
    xmid_p1 = -reverse(cumsum(reverse(Δx))) + Δx_half
    xmid_0 = knots[2:end] - Δx_half

    shift = round.(Int, xmid_0)
    diff = choose(shift .+ 2, (xmid_m1, xmid_0, xmid_p1))
    return diff, shift
end

"""
    _phase_stable(poly, wn)

Phase factor for the piecewise Legendre to Matsubara transform.

Compute the following phase factor in a stable way:

    exp.(iπ/2 * wn * cumsum(poly.Δx))
"""
function _phase_stable(poly, wn)
    xmid_diff, extra_shift = _shift_xmid(poly.knots, poly.Δx)

    if wn isa Integer
        shift_arg = wn * xmid_diff
    else
        delta_wn, wn = modf(wn)
        wn = trunc(Int, wn)
        shift_arg = wn * xmid_diff
        @. shift_arg += delta_wn * (extra_shift + xmid_diff)
    end

    phase_shifted = @. cispi(shift_arg / 2)
    corr = @. im^mod(wn * (extra_shift + 1), 4)
    return corr .* phase_shifted
end


text:
import HTTP

# create and handle an HTTP request
function request(method::String, URL::String, headers::Vector{Pair{String, String}} = Pair{String, String}[]; params...)
    exception = nothing

    # retry request up to four times
    for i in 1:4
        try
            return HTTP.request(method, URL, headers; status_exception=true, params...)
        catch e
            local found_header = false
            exception = e
            # here we find the Retry-After header and sleep for the specified amount of time
            if isa(e, HTTP.StatusError) && e.response.status == 429
                for (header, value) in e.response.headers
                    if (header == "Retry-After")
                        found_header = true
                        sleep(parse(Int, value))
                        break
                    end
                end

                # if the header wasn't found, sleep for 2 seconds
                if !found_header
                    sleep(2)
                end
            else
                # we only handle HTTP 429, so for any other status, throw the error
                rethrow
            end
        end
    end

    # if we get here, we should have run out of retries
    throw(exception)
end

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
const sqrt3 = √(3)
const sqrt2pi = √(2π)

# prepare pairwise values for `σx^2 + σy^2` and `ϕx * ϕy` for all gaussians in `gmmx` and `gmmy`
function pairwise_consts(gmmx::AbstractIsotropicGMM, gmmy::AbstractIsotropicGMM)
    t = promote_type(numbertype(gmmx),numbertype(gmmy))
    pσ, pϕ = zeros(t, length(gmmx), length(gmmy)), zeros(t, length(gmmx), length(gmmy))
    for (i,gaussx) in enumerate(gmmx.gaussians)
        for (j,gaussy) in enumerate(gmmy.gaussians)
            pσ[i,j] = gaussx.σ^2 + gaussy.σ^2
            pϕ[i,j] = gaussx.ϕ * gaussy.ϕ
        end
    end
    return pσ, pϕ
end

function pairwise_consts(mgmmx::AbstractMultiGMM{N,T,K}, mgmmy::AbstractMultiGMM{N,S,K}) where {N,T,S,K}
    t = promote_type(numbertype(mgmmx),numbertype(mgmmy))
    mpσ, mpϕ = Dict{K, Matrix{t}}(), Dict{K, Matrix{t}}()
    for key in keys(mgmmx.gmms) ∩ keys(mgmmy.gmms)
        pσ, pϕ = pairwise_consts(mgmmx.gmms[key], mgmmy.gmms[key])
        push!(mpσ, Pair(key, pσ))
        push!(mpϕ, Pair(key, pϕ))
    end
    return mpσ, mpϕ
end

"""
    lowerbound, upperbound = get_bounds(x::AbstractIsotropicGaussian, y::AbstractIsotropicGaussian, rwidth, twidth, X)
    lowerbound, upperbound = get_bounds(gmmx::AbstractSingleGMM, gmmy::AbstractSingleGMM, rwidth, twidth, X)
    lowerbound, upperbound = get_bounds(mgmmx::AbstractMultiGMM, gmmy::AbstractMultiGMM, rwidth, twidth, X)

Finds the bounds for overlap between two isotropic Gaussian distributions, two isotropic GMMs, or `two sets of
labeled isotropic GMMs for a particular region in 6-dimensional rigid rotation space, defined by `rwidth`, `twidth`,
and feature vector `X`.

`X` is a feature vector containing rotation axis components (`rx`, `ry`, and `rz`), and
translation components (`tx`, `ty`, and `tz`).

`rwidth` and `twidth` represent the sizes of the rotation and translatoin cubes, respectively,
around the point defined by `X`.

See [Campbell & Peterson, 2016](https://arxiv.org/abs/1603.00150)
"""
function get_bounds(x::AbstractIsotropicGaussian, y::AbstractIsotropicGaussian, rwidth, twidth, θ::Real, R0, t0, s=x.σ^2 + y.σ^2, w=x.ϕ*y.ϕ)
    # return Inf for bounds if the rotation lies outside the π-sphere
    if θ > π
        inf = typemax(promote_type(numbertype(x),numbertype(y)))
        return inf, inf
    end

    # prepare positions and angles
    x0 = R0*x.μ
    xnorm, ynorm = norm(x.μ), norm(y.μ-t0)
    if xnorm*ynorm == 0
        cosα = one(promote_type(numbertype(x),numbertype(y)))
    else
        cosα = dot(x0, y.μ-t0)/(xnorm*ynorm)
    end
    cosβ = cos(min(sqrt3*rwidth/2, π))

    # upper bound distance at hypercube center
    ubdist = norm(x0 + t0 - y.μ)

    # lower bound distance from the nearest point on the "spherical cap"
    if cosα >= cosβ
        lbdist = max(abs(xnorm-ynorm) - sqrt3*twidth/2, 0)
    else
        lbdist = try max(√(xnorm^2 + ynorm^2 - 2*xnorm*ynorm*(cosα*cosβ+√((1-cosα^2)*(1-cosβ^2)))) - sqrt3*twidth/2, 0)  # law of cosines
        catch e     # when the argument for the square root is negative (within machine precision of 0, usually)
            0
        end
    end

    # upperbound of dot product between directional constraints (minimizes objective)
    if length(x.dirs) == 0 || length(y.dirs) == 0
        cosγ = 1.
    else
        # NOTE: Avoid list comprehension (slow), but perform more matrix multiplications
        cosγ = -1
        for xdir in x.dirs
            for ydir in y.dirs
                cosγ = max(cosγ, dot(R0*xdir, ydir))
            end
        end
    end

    if cosγ >= cosβ
        lbdot = 1.
    else
        lbdot = cosγ*cosβ + √(1-cosγ^2)*√(1-cosβ^2)
        # lbdot = cosγ*cosβ + √(1 - cosγ^2 - cosβ^2 + cosγ^2*cosβ^2)
    end

    # evaluate objective function at each distance to get upper and lower bounds
    return -overlap(lbdist^2, s, w, lbdot), -overlap(ubdist^2, s, w, cosγ)

end

get_bounds(x::AbstractGaussian, y::AbstractGaussian, rwidth, twidth, tform::AffineMap, s=x.σ^2 + y.σ^2, w=x.ϕ*y.ϕ
    ) = get_bounds(x, y, rwidth, twidth, tform.linear.theta, RotMatrix(tform.linear), tform.translation, s, w)

get_bounds(x::AbstractGaussian, y::AbstractGaussian, rwidth, twidth, X, s=x.σ^2 + y.σ^2, w=x.ϕ*y.ϕ
    ) = get_bounds(x, y, rwidth, twidth, AffineMap(X...), s, w)

function get_bounds(gmmx::AbstractSingleGMM, gmmy::AbstractSingleGMM, rwidth, twidth, θ::Real, R0, t0, pσ=nothing, pϕ=nothing)
    # prepare pairwise widths and weights, if not provided
    if isnothing(pσ) || isnothing(pϕ)
        pσ, pϕ = pairwise_consts(gmmx, gmmy)
    end

    # sum bounds for each pair of points
    lb = 0.
    ub = 0.
    for (i,x) in enumerate(gmmx.gaussians)
        for (j,y) in enumerate(gmmy.gaussians)
            lb, ub = (lb, ub) .+ get_bounds(x, y, rwidth, twidth, θ, R0, t0, pσ[i,j], pϕ[i,j])
        end
    end
    return lb, ub
end

function get_bounds(mgmmx::AbstractMultiGMM, mgmmy::AbstractMultiGMM, rwidth, twidth, θ::Real, R0, t0, mpσ=nothing, mpϕ=nothing)
    # prepare pairwise widths and weights, if not provided
    if isnothing(mpσ) || isnothing(mpϕ)
        mpσ, mpϕ = pairwise_consts(mgmmx, mgmmy)
    end

    # sum bounds for each pair of points
    lb = 0.
    ub = 0.
    for key in keys(mgmmx.gmms) ∩ keys(mgmmy.gmms)
        lb, ub = (lb, ub) .+ get_bounds(mgmmx.gmms[key], mgmmy.gmms[key], rwidth, twidth, θ, R0, t0, mpσ[key], mpϕ[key])
    end
    return lb, ub
end

get_bounds(gmmx::AbstractGMM, gmmy::AbstractGMM, rwidth, twidth, tform::AffineMap, pσ=nothing, pϕ=nothing
    ) = get_bounds(gmmx, gmmy, rwidth, twidth, tform.linear.theta, RotMatrix(tform.linear), tform.translation, pσ, pϕ)

get_bounds(gmmx::AbstractGMM, gmmy::AbstractGMM, rwidth, twidth, X, pσ=nothing, pϕ=nothing
    ) = get_bounds(gmmx, gmmy, rwidth, twidth, AffineMap(X...), pσ, pϕ)


text:
using Test
using Algorithm101
const ALGO = Algorithm101

@testset "Find max subarray" begin
    arr = [
        13, -3, -25, 20, -3, -16, -23, 18, 20, -7, 12, -5, -22, 15, -4, 7
    ]

    @test ALGO.find_max_subarray(arr, 1, 16) == (8, 11, 43)

    arr = [
        13, -3, -25, 20, -3, -16, -23, 18, 20, -7, 12, -5, -22, 15, -4, 7
    ]

    @test ALGO.iterative_find_max_subarray(arr) == (8, 11, 43)
end


text:
import Statistics

# All loss functions should be of the form
#

"""
mse_tau
Loss function for optimising a model to mean square error
y: response
y_: prediction

Note: this is the negative mean square error
"""
function mse_tau(tau_::Float64,n::Int)
    return(n*(tau_^2))
end

"""
emse_loss
Loss function for optimising a model to expected mean square error
y: response
y_: prediction
"""
function emse_tau(y,y_::Float64)
    n_sample = length(y)
    return(n_sample*(y_^2 - 2*Statistics.var(y)))
end

"""
usage: treatment_effect([3,5,4,3,6,7,5],[1,1,1,1,0,0,0])
"""
function ate(y::Array{Float64,1},w::Array{Int64,1})
    y1 = Statistics.mean(y[w .== 1])
    y0 = Statistics.mean(y[w .== 0])
    y1 - y0
end


text:
"""
    lb(x,δ)

Convienience function to check for population extinction at theshold θ.
"""
lb(x,δ) = IfElse.ifelse(x < δ, 0, 1)

"""
    micrm_system(p, δ = 1e-6; name)

Returns a microbial consumer resource system given a dictionary of parameters. δ sets the lower threshold populations size for extinction (set to 1e-6 by default).
"""
function micrm_system(p, δ = 1e-6; name)
    @assert all([:l,:ρ,:N,:M,:m,:ω,:u] .∈ Ref(collect(keys(p)))) "missing parameters"

    u = ModelingToolkit.toparam(Symbolics.variables(:u, 1:p[:N],1:p[:M]))
    m = ModelingToolkit.toparam(Symbolics.variables(:m, 1:p[:N]))
    ρ = ModelingToolkit.toparam(Symbolics.variables(:ρ, 1:p[:M]))
    ω = ModelingToolkit.toparam(Symbolics.variables(:ω, 1:p[:M]))
    l = ModelingToolkit.toparam(Symbolics.variables(:l, 1:p[:M],1:p[:M]))

    @variables t C[1:p[:N]](t) R[1:p[:M]](t)

    D = Differential(t)

    eqns = []
    #consumers
    for i = 1:p[:N]
        RHS = -C[i] * m[i]
        for α = 1:p[:M]
            RHS += C[i] * R[α] * u[i,α]
            for β = 1:p[:M]
                RHS -= C[i] * R[α] * u[i,α] * l[α,β]
            end
        end
        push!(eqns, D(C[i]) ~ RHS)
    end

    #resources
    for α = 1:p[:M]
        RHS = ρ[α] - R[α]* ω[α]
        for i = 1:p[:N]
            RHS -= C[i] * R[α] * u[i,α]
            for β = 1:p[:M]
                RHS += C[i] * R[β] * u[i,β] * l[β,α]
            end
        end
        push!(eqns, D(R[α]) ~ RHS)
    end

    #variables
    vars = vcat(u[:],m,ρ,ω,l[:])
    vals = vcat(p[:u][:],p[:m],p[:ρ],p[:ω],p[:l][:])

    sys = ODESystem(eqns, t; name, defaults = Dict(vars[i] => vals[i] for i = eachindex(vars)))
    return(sys)
end

"""
    update_parameters!(sys, p_new)

updates an `ODESystem` object with new parameters, preventing the need for (expensive) recalucuation of the symbolic represntation of the system.
"""
function update_parameters!(sys::ODESystem, p_new)
    @assert all([:l,:ρ,:N,:M,:m,:ω,:u] .∈ Ref(collect(keys(p_new)))) "missing parameters"
    @assert (p_new[:N] + p_new[:M]) ==  length(states(sys)) "system size is different"

    #overwrite defaults
    for (k,v) in p_new
        if k ∉ [:N, :M]
            for I = CartesianIndices(v)
                if length(I) == 1
                    sys.defaults[ModelingToolkit.toparam(Symbolics.variables(k,I[1]))[1]] = v[I]
                elseif length(I) == 2
                    sys.defaults[ModelingToolkit.toparam(Symbolics.variables(k,I[1],I[2]))[1]] = v[I]
                else
                    error("$k has wrong dimensions")
                end
            end
        end
    end
end

# function update_parameters!(prob::ODEProblem,sys::ODESystem, p_new)

# end


text:
# This file was generated, do not modify it.

using MLJ
using PrettyPrinting
MLJ.color_off() # hide

@load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

scitype(X.age)

pipe = @pipeline(X -> coerce(X, :age=>Continuous),
                OneHotEncoder(),
                KNNRegressor(K=3),
                target = UnivariateStandardizer());

pipe.knn_regressor.K = 2
pipe.one_hot_encoder.drop_last = true;

evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint