text:
export apply!

using Flux
using Flux: throttle, logitcrossentropy, onecold, onehotbatch
using StatsBase: mean
using Parameters
using CUDA
using Distributions
using BSON
using Dates
using Zygote
using ArgParse

function parse_commandline()
    gpumem = CUDA.totalmem(collect(CUDA.devices())[1])/(1024^3)
    if gpumem < 2.0
        batchsize_ = 4
        num_cells_ = 4
        channels_ = 4
        trainval_fraction_ = 2f-3
    elseif gpumem < 12.0
        batchsize_ = 64
        num_cells_ = 8
        channels_ = 16
        trainval_fraction_ = 1f0
    elseif gpumem < 16.0
        batchsize_ = 64
        num_cells_ = 8
        channels_ = 16
        trainval_fraction_ = 1f0
    else
        batchsize_ = 128
        num_cells_ = 8
        channels_ = 16
        trainval_fraction_ = 1f0
    end
    global batchsize_
    global num_cells_
    global channels_
    global trainval_fraction_

    s = ArgParseSettings()

    @add_arg_table! s begin
        "--epochs"
            help = "number of epochs"
            arg_type = Int
            default = 50
        "--random_seed"
            help = "random seed (-1 for no random seed)"
            arg_type = Int
            default = 32
        "--checkpoint"
            help = "how often to save checkpoints (-1 for no checkpointing)"
            arg_type = Int
            default = -1
        "--batchsize"
            help = "batchsize"
            arg_type = Int
            default = batchsize_
        "--val_split"
            help = "fraction of train/val set to use as val set"
            arg_type = Float32
            default = 5f-1
        "--trainval_fraction"
            help = "total fraction of train/val set to use"
            arg_type = Float32
            default = trainval_fraction_
        "--test_fraction"
            help = "fraction of test set to use"
            arg_type = Float32
            default = 1f0
        "--num_cells"
            help = "number of cells in supernet"
            arg_type = Int
            default = num_cells_
        "--channels"
            help = "number of initial in supernet"
            arg_type = Int
            default = channels_
        "--rho"
            help = "admm parameter"
            arg_type = Float32
            default = 1f-3
    end

    return parse_args(s)
end

function parse_commandline_eval()
    gpumem = CUDA.totalmem(collect(CUDA.devices())[1])/(1024^3)
    if gpumem < 2.0
        batchsize_ = 1
        num_cells_ = 2
        channels_ = 18
        trainval_fraction_ = 4f-4
    elseif gpumem < 12.0
        batchsize_ = 96
        num_cells_ = 20
        channels_ = 36
        trainval_fraction_ = 1f0
    elseif gpumem < 16.0
        batchsize_ = 96
        num_cells_ = 20
        channels_ = 36
        trainval_fraction_ = 1f0
    else
        batchsize_ = 128
        num_cells_ = 20
        channels_ = 36
        trainval_fraction_ = 1f0
    end
    global batchsize_
    global num_cells_
    global channels_
    global trainval_fraction_

    s = ArgParseSettings()

    @add_arg_table! s begin
        "--epochs"
            help = "number of epochs"
            arg_type = Int
            default = 600
        "--random_seed"
            help = "random seed (-1 for no random seed)"
            arg_type = Int
            default = 32
        "--checkpoint"
            help = "how often to save checkpoints (-1 for no checkpointing)"
            arg_type = Int
            default = -1
        "--batchsize"
            help = "batchsize"
            arg_type = Int
            default = batchsize_
        "--test_batchsize"
            help = "testset batchsize"
            arg_type = Int
            default = 248
        "--val_split"
            help = "fraction of train/val set to use as val set"
            arg_type = Float32
            default = 0f0
        "--trainval_fraction"
            help = "total fraction of train/val set to use"
            arg_type = Float32
            default = trainval_fraction_
        "--test_fraction"
            help = "fraction of test set to use"
            arg_type = Float32
            default = 1f0
        "--num_cells"
            help = "number of cells in supernet"
            arg_type = Int
            default = num_cells_
        "--channels"
            help = "number of initial in supernet"
            arg_type = Int
            default = channels_
        "--droppath"
            help = "droppath probability"
            arg_type = Float32
            default = 2f-1
        "--aux"
            help = "weight of auxiliary loss"
            arg_type = Float32
            default = 4f-1
    end

    return parse_args(s)
end

@with_kw struct trial_params
    epochs::Int = 50
    batchsize::Int = batchsize_
    throttle_::Int = 20
    val_split::Float32 = 0.5
    trainval_fraction::Float32 = trainval_fraction_
    test_fraction::Float32 = 1.0
    num_cells::Int = num_cells_
    channels::Int = channels_
end

function loss(m, x, y)
    out = m(x)
    loss = logitcrossentropy(squeeze(out), y)
    return loss
end

acccb() = @show(accuracy_batched(m, val))
function accuracy(m, x, y)
    out = mean(onecold(m(x), 1:10) .== onecold(y, 1:10))
end
function accuracy(m, x, y, pert)
    out = mean(onecold(m(x, pert), 1:10) .== onecold(y, 1:10))
end
function accuracy_batched(m, xy)
    CUDA.reclaim()
    GC.gc()
    score = 0f0
    count = 0
    for batch in CuIterator(xy)
        @show acc = accuracy(m, batch...)
        score += acc*length(batch)
        count += length(batch)
        CUDA.reclaim()
        GC.gc()
    end
    display(score / count)
    score / count
end
function accuracy_batched(m, xy, pert)
    CUDA.reclaim()
    GC.gc()
    score = 0f0
    count = 0
    for batch in CuIterator(xy)
        @show acc = accuracy(m, batch..., pert)
        score += acc*length(batch)
        count += length(batch)
        CUDA.reclaim()
        GC.gc()
    end
    display(score / count)
    score / count
end
function accuracy_unbatched(m, xy; pert = [])
    CUDA.reclaim()
    GC.gc()
    xy = xy | gpu
    acc = accuracy(m, xy..., pert = pert)
    foreach(CUDA.unsafe_free!, xy)
    CUDA.reclaim()
    GC.gc()
    acc
end
Base.@kwdef mutable struct histories
    normal_αs::Vector{Vector{Array{Float32, 1}}}
    reduce_αs::Vector{Vector{Array{Float32, 1}}}
    activations::Vector{Dict}
    accuracies::Vector{Float32}
end

histories() = histories([],[],[],[])

function (hist::histories)()
    push!(hist.normal_αs, deepcopy(m.normal_αs) |> cpu)
    push!(hist.reduce_αs, deepcopy(m.reduce_αs) |> cpu)
    #push!(hist.activations, m.activations |> cpu)
    #push!(hist.accuracies, accuracy_batched(m, val |> gpu))
end

Base.@kwdef mutable struct historiessm
    normal_αs_sm::Vector{Vector{Array{Float32, 1}}}
    reduce_αs_sm::Vector{Vector{Array{Float32, 1}}}
    activations::Vector{Dict}
    accuracies::Vector{Float32}
end

historiessm() = historiessm([],[],[],[])

function (hist::historiessm)()
    #@show losses
    push!(hist.normal_αs_sm, softmax.(deepcopy(m.normal_αs)) |> cpu)
    push!(hist.reduce_αs_sm, softmax.(deepcopy(m.reduce_αs)) |> cpu)
    #push!(hist.activations, copy(m.activations.currentacts) |> cpu)
    CUDA.reclaim()
    GC.gc()
end

Base.@kwdef mutable struct historiessml
    normal_αs_sm::Vector{Vector{Array{Float32, 1}}}
    reduce_αs_sm::Vector{Vector{Array{Float32, 1}}}
    activations::Vector{Dict}
    accuracies::Vector{Float32}
    train_losses::Vector{Float32}
    val_losses::Vector{Float32}
end

historiessml() = historiessml([],[],[],[],[],[])

function (hist::historiessml)()
    @show losses
    push!(hist.normal_αs_sm, softmax.(deepcopy(m.normal_αs)) |> cpu)
    push!(hist.reduce_αs_sm, softmax.(deepcopy(m.reduce_αs)) |> cpu)
    #push!(hist.activations, copy(m.activations.currentacts) |> cpu)
    push!(hist.train_losses, losses[1])
    push!(hist.val_losses, losses[2])
    CUDA.reclaim()
    GC.gc()
end

struct CbAll
    cbs
end
CbAll(cbs...) = CbAll(cbs)

function timefunc(func)
    @show func
    @time out = func()
    out
end

(cba::CbAll)() = foreach(cb -> cb(), cba.cbs)

function prepare_folder(algo::String, args::Dict)
    if "SLURM_JOB_ID" in keys(ENV)
        uniqueid = ENV["SLURM_JOB_ID"]
    else
        uniqueid = Dates.now()
    end
    if ispath("/gpfs/work/p21001/maile/dnas/models/")
        model_dir = "/gpfs/work/p21001/maile/dnas/models/"
    elseif ispath("/projets/reva/kmaile/dnas/models/")
        model_dir = "/projets/reva/kmaile/dnas/models/"
    else
        model_dir = "test/models/"
    end
    base_folder = string(model_dir, algo, "_", uniqueid)
    mkpath(base_folder)
    BSON.@save joinpath(base_folder, "args.bson") args
    base_folder
end

function save_progress()
    m_cpu = m |> cpu
    normal_αs = m_cpu.normal_αs
    reduce_αs = m_cpu.reduce_αs
    #BSON.@save joinpath(base_folder, "model.bson") m_cpu args
    BSON.@save joinpath(base_folder, "alphas.bson") normal_αs reduce_αs
    BSON.@save joinpath(base_folder, "histepoch.bson") histepoch
    BSON.@save joinpath(base_folder, "histbatch.bson") histbatch
    if args["checkpoint"] > 0 && length(histepoch.train_losses) % args["checkpoint"] == 0
        BSON.@save joinpath(base_folder, string("model", length(histepoch.train_losses), ".bson")) m_cpu args
    end
end

max_stars_repo_name: huangyxi/DifferentiableNAS

max_issues_repo_name: huangyxi/DifferentiableNAS

max_forks_repo_name: huangyxi/DifferentiableNAS


text:
"""
    PyBuffer(o, [flags=C.PyBUF_FULL_RO])

A reference to the underlying buffer of `o`, if it satisfies the buffer protocol.

Has the following properties:
- `buf`: Pointer to the data.
- `obj`: The exporting object (usually `o`).
- `len`: The length of the buffer in bytes.
- `readonly`: True if the buffer is immutable.
- `itemsize`: The size of each element.
- `format`: The struct-syntax format of the element type.
- `ndim`: The number of dimensions.
- `shape`: The length of the buffer in each dimension.
- `strides`: The strides (in bytes) of the buffer in each dimension.
- `suboffsets`: For indirect arrays. See the buffer protocol documentation.
- `isccontiguous`: True if the buffer is C-contiguous (e.g. numpy arrays).
- `isfcontiguous`: True if the buffer is Fortran-contiguous (e.g. Julia arrays).
- `eltype`: The element type.
"""
mutable struct PyBuffer
    info::Array{C.Py_buffer,0}
    function PyBuffer()
        b = new(fill(C.Py_buffer()))
        finalizer(b) do b
            if CONFIG.isinitialized
                with_gil(false) do
                    C.PyBuffer_Release(pointer(b.info))
                end
            end
        end
        b
    end
end
PyBuffer(o, flags::Integer = C.PyBUF_FULL_RO) = begin
    b = PyBuffer()
    check(C.PyObject_GetBuffer(o, pointer(b.info), flags))
    b
end
export PyBuffer

Base.getproperty(b::PyBuffer, k::Symbol) =
    if k == :buf
        b.info[].buf
    elseif k == :obj
        pyborrowedobject(b.info[].obj)
    elseif k == :len
        b.info[].len
    elseif k == :readonly
        !iszero(b.info[].readonly)
    elseif k == :itemsize
        b.info[].itemsize
    elseif k == :format
        p = b.info[].format
        p == C_NULL ? "B" : unsafe_string(p)
    elseif k == :ndim
        b.info[].ndim
    elseif k == :shape
        p = b.info[].shape
        p == C_NULL ? (fld(b.len, b.itemsize),) : ntuple(i -> unsafe_load(p, i), b.ndim)
    elseif k == :strides
        p = b.info[].strides
        p == C_NULL ? size_to_cstrides(b.itemsize, b.shape...) :
        ntuple(i -> unsafe_load(p, i), b.ndim)
    elseif k == :suboffsets
        p = b.info[].suboffsets
        p == C_NULL ? ntuple(i -> -1, b.ndim) : ntuple(i -> unsafe_load(p, i), b.ndim)
    elseif k == :isccontiguous
        b.strides == size_to_cstrides(b.itemsize, b.shape...)
    elseif k == :isfcontiguous
        b.strides == size_to_fstrides(b.itemsize, b.shape...)
    elseif k == :eltype
        pybufferformat_to_type(b.format)
    else
        getfield(b, k)
    end

ext: jl


text:
__precompile__(false)

module NMRDecompose

using NMR

include("xcorr.jl")
include("util.jl")
include("distributions.jl")

end


text:
__precompile__(true)

module DataValueArrays

using Reexport
@reexport using Base.Cartesian
@reexport using DataValues

export DataValueArray,
       DataValueVector,
       DataValueMatrix,

       # Macros

       # Methods
       dropnull,
       dropnull!,
       nullify!,
       padnull!,
       padnull

include("typedefs.jl")
include("constructors.jl")
include("primitives.jl")
include("indexing.jl")
include("broadcast.jl")
include("map.jl")
include("nullablevector.jl")
include("reduce.jl")
include("subarray.jl")

end


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type Method
    name::AbstractString
    request_type_url::AbstractString
    request_streaming::Bool
    response_type_url::AbstractString
    response_streaming::Bool
    options::Array{Option,1}
    syntax::Int32
    Method(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Method
hash(v::Method) = ProtoBuf.protohash(v)
isequal(v1::Method, v2::Method) = ProtoBuf.protoisequal(v1, v2)
==(v1::Method, v2::Method) = ProtoBuf.protoeq(v1, v2)

type Mixin
    name::AbstractString
    root::AbstractString
    Mixin(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Mixin
hash(v::Mixin) = ProtoBuf.protohash(v)
isequal(v1::Mixin, v2::Mixin) = ProtoBuf.protoisequal(v1, v2)
==(v1::Mixin, v2::Mixin) = ProtoBuf.protoeq(v1, v2)

type Api
    name::AbstractString
    methods::Array{Method,1}
    options::Array{Option,1}
    version::AbstractString
    source_context::SourceContext
    mixins::Array{Mixin,1}
    syntax::Int32
    Api(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type Api
hash(v::Api) = ProtoBuf.protohash(v)
isequal(v1::Api, v2::Api) = ProtoBuf.protoisequal(v1, v2)
==(v1::Api, v2::Api) = ProtoBuf.protoeq(v1, v2)

export Api, Method, Mixin


text:
export UpdateMasseyColley

"""
    UpdateMasseyColley

 MasseyColley update rule, e.g., see "Whos's #1", Langville and Meyer, p.25
    this is the "Colleyized Massey method", i.e., use Colley's matrix on scores

## Parameters (none)
```
"""
struct UpdateMasseyColley <: UpdateRule
# this is a batch calculation that takes no account of past ratings, and
# has no parameters
end

function update_info( rule::UpdateMasseyColley )
    info = Dict(
                :name => "MasseyColley",
                :reference => "\"Whos's #1\", Langville and Meyer, p.25",
                :computation => "simultaneous",
                :state_model => "none",
                :input => "score",
                :output => "deterministic",
                :model => "single",
                :ties => true,
                :factors => false,
                :parameters => [],
                :record => true
                )
    return info
end
update_info( ::Type{UpdateMasseyColley} ) =  update_info( UpdateMasseyColley() )

function update_ratings( rule::UpdateMasseyColley,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame)
    n = size(input_competitions,1)
    m = length( input_ratings.players )
    I = player_indexes( input_ratings.players )

    # construct Colley matrices and vectors
    C = diagm(0 => 2*ones(Int,m))
    p = zeros(Int, m)
    d = input_competitions # just an abbreviation
    point_diff = d[:,ScoreA] - d[:,ScoreB]
    for i=1:n
        p[ I[d[i,PlayerA]] ] +=  point_diff[i]
        p[ I[d[i,PlayerB]] ] -=  point_diff[i]
        C[ I[d[i,PlayerA]], I[d[i,PlayerB]] ] -= 1
        C[ I[d[i,PlayerB]], I[d[i,PlayerA]] ] -= 1
        C[ I[d[i,PlayerA]], I[d[i,PlayerA]] ] += 1
        C[ I[d[i,PlayerB]], I[d[i,PlayerB]] ] += 1
    end

    # solve Colley's equation
    r = C \ p
    ratings = Dict{String, Float64}()
    for player in input_ratings.players
        ratings[player] = r[ I[player] ]
    end

    # output ratings list
    output_ratings = RatingsList(input_ratings.players, ratings )
    return output_ratings
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie, GeometryBasics, Colors
GLMakie.activate!() # HIDE
let
    positions = vec([(i/5, j/5, k/5) for i=1:7,j=1:7,k=1:7]) # note 7 > 5 [factor in each i,j,k], whichs is misleading
    fig, ax, obj = mesh(HyperRectangle(Vec3f0(positions[1]...), Vec3f0(0.15)),
        color = RGBA(positions[1]..., 0.5),
        transparency = false)
    wireframe!(ax,HyperRectangle(Vec3f0(positions[1]...), Vec3f0(0.15)),
        linewidth = 0.1, overdraw = false)
    for i in 2:length(positions)
        mesh!(ax, HyperRectangle(Vec3f0(positions[i]...), Vec3f0(0.15)),
            color = RGBA(positions[i]..., 0.5))
        wireframe!(ax, HyperRectangle(Vec3f0(positions[i]...), Vec3f0(0.15)),
            linewidth = 0.1, overdraw = false)
    end
    fig
    save(joinpath(@__DIR__, "output", "RGBcube.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie","GeometryBasics", "Colors"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
struct VSCodeLogger <: Logging.AbstractLogger end

const logger_lock = ReentrantLock()
function Logging.handle_message(j::VSCodeLogger, level, message, _module,
    group, id, file, line; kwargs...)
    isprogress = try_process_progress(level, message, _module, group, id, file, line; kwargs...) do progress
        lock(logger_lock)
        try
            JSONRPC.send_notification(conn_endpoint[], "repl/updateProgress", progress)
            JSONRPC.flush(conn_endpoint[])
        finally
            unlock(logger_lock)
        end
    end isa Some

    if isprogress
        return nothing
    end

    previous_logger = Logging.global_logger()
    # Pass through non-progress log messages to the global logger iff the global logger would handle it:
    if (Base.invokelatest(Logging.min_enabled_level, previous_logger) <= Logging.LogLevel(level) ||
        Base.CoreLogging.env_override_minlevel(group, _module)) &&
       Base.invokelatest(Logging.shouldlog, previous_logger, level, _module, group, id)
        Logging.handle_message(previous_logger, level, message, _module,
            group, id, file, line; kwargs...)
    end
    return nothing
end

Logging.shouldlog(::VSCodeLogger, level, _module, group, id) = true

Logging.catch_exceptions(::VSCodeLogger) = true

function Logging.min_enabled_level(::VSCodeLogger)
    min(Base.invokelatest(Logging.min_enabled_level, Logging.global_logger()), Logging.LogLevel(-1))
end

const progresslogging_pkgid = Base.PkgId(
    Base.UUID("33c8b6b6-d38a-422a-b730-caa89a2f386c"),
    "ProgressLogging"
)

"""
    try_process_progress(f, args...; kwargs...) -> nothing or Some(ans)

Try to process logging record by a function `f(::Progress) -> ans`.  Arguments
`args` and `kwargs` are the ones passed to ```Logging.handle_message`.  Return
`Some(ans)` if it is a progress record.
"""
function try_process_progress(f, args...; kwargs...)
    m = get(Base.loaded_modules, progresslogging_pkgid, nothing)
    m === nothing && return nothing
    isdefined(m, :asprogress) || return nothing
    Base.invokelatest() do
        progress = m.asprogress(args...; kwargs...)
        progress === nothing && return nothing
        return Some(f(progress))
    end
end

max_stars_repo_name: xgdgsc/julia-vscode

max_issues_repo_name: xgdgsc/julia-vscode

max_forks_repo_name: xgdgsc/julia-vscode


text:
using Conda
using PyCall

function install_miplearn()
    Conda.update()
    pip = joinpath(dirname(pyimport("sys").executable), "pip")
    isfile(pip) || error("$pip: invalid path")
    run(`$pip install miplearn==0.2.0.dev13`)
end

install_miplearn()


text:
function binomial_coeffs(n::Integer)::Vector{Int}
    if n == 1
        return [1]
    end
    if n == 2
        return [1, 1]
    end
    c = [1, 1]
    for _ in 3:n
        c=[1, (view(c[1:end-1]) .+ view(c[2:end])), 1]
    end
    return c
end

function binomial_coeffs_2(n::Integer)
    if n == 1
        return [1]
    end
    if n == 2
        return [1, 1]
    end
    c = [1 , 1]
    m = 2
    for _ in 3:n
        for i in 1:(m - 1)
            c[i] = c[i] + c[i+1]
            m += 1
        end
        pushfirst!(c, 1)
    end
    return c
end


text:
module NeuronBuilder

using ModelingToolkit

const t = Num(Sym{Real}(:t))
const D = Differential(t)

include("channels.jl")
include("synapses.jl")
include("build_neuron.jl")
include("compartments.jl")

export Synapse, IonChannel, Conductance, Soma, Compartment
export NaV, CaS, CaT, H, Ka, KCa, Kdr, Leak
export Chol, Glut
export build_channel, build_synapse, build_neuron, build_group
export add_connection

end

max_stars_repo_name: wsphillips/NeuronBuilder

max_issues_repo_name: wsphillips/NeuronBuilder

max_forks_repo_name: wsphillips/NeuronBuilder


text:
set_bigfloat_precision(128)
temp = BigFloat(0)
f(x) = exp(-x) - sin(x)
steffensen(c) = c - (f(c)*f(c))/( f(c+f(c)) - f(c) )
temp2 = BigFloat(100)

witch = big"0.588532743981861077432452045702903688531271516109030533319914299511672553307351427738524061576027409562153528176982466770293849745782742957500713135"

for i = 1:10
	temp2 = temp
	temp = steffensen(temp)
	@printf("%.4e\t%.4e\n", temp , abs(temp - witch) )
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
function colormap(name::AbstractString)
	if haskey(access_dict,name)
		return access_dict[name]()
	else
		# could be color sequence name string
		return construct_from_string(name)
	end

end

max_stars_repo_name: lobingera/ColorMaps


text:
"""
struct PCurlGradMonomialBasis{...} <: AbstractArray{Monomial}

This type implements a multivariate vector-valued polynomial basis
spanning the space needed for Raviart-Thomas reference elements on simplices.
The type parameters and fields of this `struct` are not public.
This type fully implements the [`Field`](@ref) interface, with up to first order
derivatives.
"""
struct PCurlGradMonomialBasis{D,T} <: AbstractVector{Monomial}
  order::Int
  pterms::Array{CartesianIndex{D},1}
  sterms::Array{CartesianIndex{D},1}
  perms::Matrix{Int}
  function PCurlGradMonomialBasis(::Type{T},order::Int,
      pterms::Array{CartesianIndex{D},1},sterms::Array{CartesianIndex{D},1},
      perms::Matrix{Int}) where {D,T}
    new{D,T}(order,pterms,sterms,perms)
  end
end

Base.size(a::PCurlGradMonomialBasis) = (_ndofs_pgrad(a),)
# @santiagobadia : Not sure we want to create the monomial machinery
Base.getindex(a::PCurlGradMonomialBasis,i::Integer) = Monomial()
Base.IndexStyle(::PCurlGradMonomialBasis) = IndexLinear()

"""
PCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}

Returns a `PCurlGradMonomialBasis` object. `D` is the dimension
of the coordinate space and `T` is the type of the components in the vector-value.
The `order` argument has the following meaning: the divergence of the  functions
in this basis is in the P space of degree `order`.
"""
function PCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}
  @assert T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
  P_k = MonomialBasis{D}(T, order, _p_filter)
  S_k = MonomialBasis{D}(T, order, _s_filter)
  pterms = P_k.terms
  sterms = S_k.terms
  perms = _prepare_perms(D)
  PCurlGradMonomialBasis(T,order,pterms,sterms,perms)
end

"""
    num_terms(f::PCurlGradMonomialBasis{D,T}) where {D,T}
"""
function num_terms(f::PCurlGradMonomialBasis{D,T}) where {D,T}
  Int(_p_dim(f.order,D)*D + _p_dim(f.order,D-1))
end

get_order(f::PCurlGradMonomialBasis{D,T}) where {D,T} = f.order

return_type(::PCurlGradMonomialBasis{D,T}) where {D,T} = T

function return_cache(f::PCurlGradMonomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
  @assert D == length(eltype(x)) "Incorrect number of point components"
  np = length(x)
  ndof = _ndofs_pgrad(f)
  n = 1 + f.order+1
  V = VectorValue{D,T}
  r = CachedArray(zeros(V,(np,ndof)))
  v = CachedArray(zeros(V,(ndof,)))
  c = CachedArray(zeros(T,(D,n)))
  (r, v, c)
end

function evaluate!(cache,f::PCurlGradMonomialBasis{D,T},x::AbstractVector{<:Point}) where {D,T}
  r, v, c = cache
  np = length(x)
  ndof = _ndofs_pgrad(f)
  n = 1 + f.order+1
  setsize!(r,(np,ndof))
  setsize!(v,(ndof,))
  setsize!(c,(D,n))
  for i in 1:np
    @inbounds xi = x[i]
      _evaluate_nd_pcurlgrad!(v,xi,f.order+1,f.pterms,f.sterms,f.perms,c)
    for j in 1:ndof
      @inbounds r[i,j] = v[j]
    end
  end
  r.array
end

function return_cache(
  fg::FieldGradientArray{1,PCurlGradMonomialBasis{D,T}},
  x::AbstractVector{<:Point})  where {D,T}

  f = fg.fa
  @assert D == length(eltype(x)) "Incorrect number of point components"
  np = length(x)
  ndof = _ndofs_pgrad(f)
  n = 1 + f.order+1
  xi = testitem(x)
  V = VectorValue{D,T}
  G = gradient_type(V,xi)
  r = CachedArray(zeros(G,(np,ndof)))
  v = CachedArray(zeros(G,(ndof,)))
  c = CachedArray(zeros(T,(D,n)))
  g = CachedArray(zeros(T,(D,n)))
  (r, v, c, g)
end

function evaluate!(cache,
  fg::FieldGradientArray{1,PCurlGradMonomialBasis{D,T}},
  x::AbstractVector{<:Point}) where {D,T}

  f = fg.fa
  r, v, c, g = cache
  np = length(x)
  ndof = _ndofs_pgrad(f)
  n = 1 + f.order+1
  setsize!(r,(np,ndof))
  setsize!(v,(ndof,))
  setsize!(c,(D,n))
  setsize!(g,(D,n))
  V = VectorValue{D,T}
  for i in 1:np
    @inbounds xi = x[i]
    _gradient_nd_pcurlgrad!(v,xi,f.order+1,f.pterms,f.sterms,f.perms,c,g,V)
    for j in 1:ndof
      @inbounds r[i,j] = v[j]
    end
  end
  r.array
end

# Helpers

_p_filter(e,order) = (sum(e) <= order)
_s_filter(e,order) = (sum(e) == order)

function _p_dim(order,D)
  dim = 1
  for d in 1:D
    dim *= order+d
  end
  dim/factorial(D)
end

_ndofs_pgrad(f::PCurlGradMonomialBasis{D}) where D = num_terms(f)

function _evaluate_nd_pcurlgrad!(
  v::AbstractVector{V},
  x,
  order,
  pterms::Array{CartesianIndex{D},1},
  sterms::Array{CartesianIndex{D},1},
  perms::Matrix{Int},
  c::AbstractMatrix{T}) where {V,T,D}

  dim = D
  for d in 1:dim
    _evaluate_1d!(c,x,order,d)
  end

  o = one(T)
  k = 1
  m = zero(Mutable(V))
  js = eachindex(m)
  z = zero(T)

  for ci in pterms
    for j in js

      @inbounds for i in js
        m[i] = z
      end

      s = o
      @inbounds for d in 1:dim
        s *= c[d,ci[perms[d,j]]]
      end

      m[j] = s
      v[k] = m
      k += 1
    end
  end

  for ci in sterms
    @inbounds for i in js
      m[i] = z
    end
    for j in js

      s = c[j,2]
      @inbounds for d in 1:dim
        s *= c[d,ci[d]]
      end

      m[j] = s

    end
    v[k] = m
    k += 1
  end
end

function _gradient_nd_pcurlgrad!(
  v::AbstractVector{G},
  x,
  order,
  pterms::Array{CartesianIndex{D},1},
  sterms::Array{CartesianIndex{D},1},
  perms::Matrix{Int},
  c::AbstractMatrix{T},
  g::AbstractMatrix{T},
  ::Type{V}) where {G,T,D,V}

  dim = D
  for d in 1:dim
    _evaluate_1d!(c,x,order,d)
    _gradient_1d!(g,x,order,d)
  end

  z = zero(Mutable(V))
  m = zero(Mutable(G))
  js = eachindex(z)
  mjs = eachindex(m)
  o = one(T)
  zi = zero(T)
  k = 1

  for ci in pterms
    for j in js

      s = z
      for i in js
        s[i] = o
      end

      for q in 1:dim
        for d in 1:dim
          if d != q
            @inbounds s[q] *= c[d,ci[perms[d,j]]]
          else
            @inbounds s[q] *= g[d,ci[perms[d,j]]]
          end
        end
      end

      @inbounds for i in mjs
        m[i] = zi
      end

      for i in js
        @inbounds m[i,j] = s[i]
      end
      @inbounds v[k] = m
      k += 1
    end
  end

  for ci in sterms

    @inbounds for i in mjs
      m[i] = zi
    end

    for j in js

      s = z
      for i in js
        s[i] = c[j,2]
      end

      for q in 1:dim
        for d in 1:dim
          if d != q
            @inbounds s[q] *= c[d,ci[d]]
          else
            @inbounds s[q] *= g[d,ci[d]]
          end
        end
      end
      aux = o
      @inbounds for d in 1:dim
        aux *= c[d,ci[d]]
      end
      s[j] += aux

      for i in js
        @inbounds m[i,j] = s[i]
      end
    end
    @inbounds v[k] = m
    k += 1
  end
end


text:
################################################################################
#                             Euler Angle and Axis
################################################################################

export angleaxis2quat

################################################################################
#                                 Conversions
################################################################################

# Quaternions
# ==============================================================================

"""
    function angleaxis2quat(a::Number, v::Vector)

Convert the Euler angle `a` [rad] and Euler axis `v`, which must be a unit
vector, to a quaternion.

# Remarks

It is expected that the vector `v` is unitary. However, no verification is
performed inside the function. The user must handle this situation.

# Example

```julia-repl
julia> v = [1;1;1];

julia> v /= norm(v);

julia> angleaxis2quat(pi/2,v)
Quaternion{Float64}:
  + 0.7071067811865476 + 0.408248290463863.i + 0.408248290463863.j + 0.408248290463863.k
```

"""
function angleaxis2quat(a::Number, v::Vector)
    # Check the arguments.
    if length(v) > 3
        throw(ArgumentError("The provided vector for the Euler axis must have 3 elements."))
    end

    # Create the quaternion.
    Quaternion( cos(a/2), sin(a/2)*v )
end

"""
    function angleaxis2quat(angleaxis::EulerAngleAxis)

Convert a Euler angle and Euler axis `angleaxis` (see `EulerAngleAxis`) to a
quaternion.

# Remarks

It is expected that the vector `angleaxis.v` is unitary. However, no
verification is performed inside the function. The user must handle this
situation.

# Example

```julia-repl
julia> v = [1;1;1];

julia> v /= norm(v);

julia> angleaxis2quat(EulerAngleAxis(pi/2,v))
Quaternion{Float64}:
  + 0.7071067811865476 + 0.408248290463863.i + 0.408248290463863.j + 0.408248290463863.k
```

"""
angleaxis2quat(angleaxis::EulerAngleAxis) = angleaxis2quat(angleaxis.a,
                                                           angleaxis.v)


text:
module Squareness

export issquare

_num_iteration(::I) where (I <: Integer) = 10 - leading_zeros(sizeof(I) % UInt8 - 0x01)
_num_iteration(::Bool) = 0
_num_iteration(x::BigInt) = 66 - leading_zeros(sizeof(x.d) * UInt64(abs(x.size)) - 0x01)

_inttype(::Type{UInt}) = UInt
_inttype(::Type{I}) where {I<:Integer} = promote_type(I, UInt)

function issquare_with_padic_newton(x::I) where {I <: Integer}
    x == 0 && return true
    xtz = trailing_zeros(x)
    isodd(xtz) && return false
    x >>= unsigned(xtz)
    _I = _inttype(I)
    y::_I = 1
    N = _num_iteration(x)
    for _ in 1:N
        y = ((I(3) - y * y * x) * y) >> 0x1
    end
    xrt::_I = x * y
    m::_I = one(_I) << (UInt(1) << N)
    if xrt & m != 0
        xrt = -xrt
    end
    xrt &= m - one(I)
    xrt * xrt == x
end

isqrt(n::Integer) = Base.isqrt(n)
isqrt(n::I) where {I<:Union{Int64, UInt64, Int128, UInt128}} = round(I, sqrt(n))

function issquare_with_isqrt(x::Integer)
    if x < 0
        false
    else
        xrt = isqrt(x)
        xrt * xrt == x
    end
end

issquare(n::Unsigned) = issquare_with_padic_newton(n)
issquare(n::Signed) = n < 0 ? false : issquare_with_padic_newton(n % Unsigned)
issquare(b::Bool) = true
# issquare(n::BigInt) = n < 0 ? false : issquare_with_padic_newton(n)
issquare(n::BigInt) = issquare_with_isqrt(n)

end # module


text:
sym2status(sym, i) = sym in (:∈, :⟂, :∥, :≅) ? SemiFree(i) : Dependent(i)

function _add_point!(dofs, P, s)
    if haskey(dofs, P)
        dofs[P] = _update_status(dofs[P], s)
    else
        dofs[P] = s
    end
end

function _parse_statement!(dofs, intros, ex)
    if ex.head == :(:=)
        newpoint = ex.args[1]
        constr = ex.args[2]
        newstatus = sym2status(constr.args[1], length(intros.args)+1)
        _add_free_points!(dofs, constr)
        _add_point!(dofs, newpoint, newstatus)
        _parse_constraint!(intros, constr, newpoint)
    elseif ex.head == :call && (ex.args[1] in (:Points, :points))
        for P in ex.args[2:end]
            _add_point!(dofs, P, Free())
        end
    end
end

function _parse_constraint!(intros, constr, newpoint=nothing)
    constr = _quotenodify(constr)
    if !isnothing(newpoint) && constr.args[1] in (:↓, :Midpoint, :midpoint, :∩, :Circle, :circle)
         insert!(constr.args, 2, QuoteNode(newpoint))
    end
    push!(intros.args, constr)
end

function _add_free_points!(dofs, ex::Expr)
    for s in ex.args[2:end]
        _add_free_points!(dofs, s)
    end
end

_add_free_points!(dofs, s::Symbol) = _add_point!(dofs, s, Free())

function _quotenodify(ex::Expr)
    for i in 2:length(ex.args)
        ex.args[i] = _quotenodify(ex.args[i])
    end
    return ex
end

_quotenodify(s::Symbol) = QuoteNode(s)


text:
_default_table_file_name(type::AbstractString) = "$(type)_cached.tsv"
_default_jld_file_name(type::AbstractString, hash::AbstractString) = string("$(type)_$(hash).jld")

function cached_obj_exists(
	type::AbstractString, common_cache_dir::AbstractString, hash::AbstractString
)::Bool
	return isdir(common_cache_dir) &&
		isfile(joinpath(common_cache_dir, _default_jld_file_name(type, hash)))
end

function cache_obj(
	type::AbstractString,
	common_cache_dir::AbstractString,
	obj::Any,
	hash::AbstractString,
	args_string::AbstractString;
	column_separator::AbstractString = "\t",
	time_spent::Dates.Millisecond,
	use_serialize::Bool = false
)

	total_save_path = joinpath(common_cache_dir, _default_jld_file_name(type, hash))
	mkpath(dirname(total_save_path))

	cache_table_exists = isfile(joinpath(common_cache_dir, _default_table_file_name(type)))

	# cache record
	if settings.create_cache_record
		table_file = open(joinpath(common_cache_dir, _default_table_file_name(type)), "a+")

		if !cache_table_exists
			write(table_file,
				join([
					"TIMESTAMP",
					"FILE NAME",
					"COMP TIME",
					"COMMAND",
					"TYPE",
					"JULIA_VERSION",
					"\n"
				], column_separator))
		end

		write(table_file, string(
				Dates.format(Dates.now(), "dd/mm/yyyy HH:MM:SS"), column_separator,
				_default_jld_file_name(type, hash), column_separator,
				human_readable_time(time_spent), column_separator,
				args_string, column_separator,
				typeof(obj), column_separator,
				VERSION, column_separator,
				"\n"))
		close(table_file)
	end

	# log
	if settings.log
		printcheckpoint(
				settings.log_output,
				"Saving $(type) to file $(total_save_path)[.tmp]...";
				format = settings.date_format
			)
	end

	if use_serialize
		Serialization.serialize("$(total_save_path).tmp", obj)
	else
		@save "$(total_save_path).tmp" obj
	end
	mv("$(total_save_path).tmp", "$(total_save_path)"; force = true)
end

function load_cached_obj(
	type::AbstractString, common_cache_dir::AbstractString, hash::AbstractString
)
	total_load_path = joinpath(common_cache_dir, _default_jld_file_name(type, hash))

	if settings.log
		printcheckpoint(SimpleCaching.settings.log_output, "Loading $(type) from file " *
			"$(total_load_path)..."; format = settings.date_format)
	end

	obj = nothing

	# TODO use magic number check instead of try/catch
	try
		@load total_load_path obj
	catch e
		try
			obj = Serialization.deserialize(total_load_path)
		catch e
			throw_n_log("File $(total_load_path) is neither in JLD2 format nor a " *
				"Serialized object", ArgumentError)
		end
	end

	obj
end

macro _scache(type, common_cache_dir, ex)
	will_use_serialize = _use_serialize

	esc_times = escdepth(ex)
	ex = unesc_comp(ex)

	(typeof(ex) != Expr || ex.head != :call) && (throw(ArgumentError("`@scache[jld]` can " *
		"be used only with function calls: passed $(ex)")))

	# hyigene
	type = esc(type)
	common_cache_dir = esc(common_cache_dir)
	ex = esc(ex, esc_times)

	t = _convert_input(ex, true)

	rel_esc(ex) = esc(ex, esc_times+1)
	as, ks, vs, res = (_toexpr(t.args), _toexpr(t.kwargs)..., t.res)

	# TODO: make a function that interprets `res` to be inserted in `args` or `kwargs`
	return quote
		_hash = get_hash_sha256((
			args = $(rel_esc(as)),
			kwargs = merge(
				Dict{Symbol,Any}(
					zip(
						$(rel_esc(ks)),
						$(rel_esc(vs))
					)
				),
				Dict{Symbol,Any}(
					pairs($(rel_esc(res)))
				)
			)
		))

		if cached_obj_exists($(type), $(common_cache_dir), _hash)
			load_cached_obj($(type), $(common_cache_dir), _hash)
		else
			if settings.log
				printcheckpoint(
					settings.log_output, "Computing " * $(type) * "...";
					format = settings.date_format,
					start = settings.line_starter
				)
			end

			_started = _start_timer()
			_result_value = $(esc(ex))
			_finish_time = _stop_timer(_started)

			if settings.log
				printcheckpoint(
					settings.log_output,
					"Computed " * $(type) * " in " * human_readable_time_s(_finish_time) *
						" seconds (" * human_readable_time(_finish_time) * ")";
					format = settings.date_format,
					start = settings.line_starter
				)
			end

			cache_obj(
				$(type),
				$(common_cache_dir),
				_result_value,
				_hash,
				$(rel_esc(_strarg(ex, true))),
				time_spent = _finish_time,
				use_serialize = $(esc(will_use_serialize))
			)

			_result_value
		end
	end
end

const _ext_docs = """!!! note

	The file extension will be `.jld` when using both [`scache`](@ref) and
	[`scachejld`](@ref).
"""

const _same_args_docs = """Expressions in function argument will be computed as first step
so the cached file will be loaded even if the arguments are different but will evaluate to
the same result.
"""

"""
	scache(type, cache_dir, function_call)

Cache the result of `function_call` in directory `cache_dir` prefixing the saved file with
`type`.

This macro uses `Serialize` `serialize` and `deserialize` functions to save and load cached
files so it is faster and more memory efficent than [`scachejld`](@ref) macro which uses
`JLD2` which, on the other hand, is more portable between different julia versions.

$_ext_docs

## Examples
```
julia> using SimpleCaching

julia> SimpleCaching.settings.log = true
true

julia> @scache "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 09:54:08 ] Computing cute-cube...
● [ 31/12/2021 09:54:08 ] Computed cute-cube in 0.044 seconds (00:00:00)
● [ 31/12/2021 09:54:13 ] Saving cute-cube to file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld[.tmp]...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

julia> @scache "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 09:54:24 ] Loading cute-cube from file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

 shell> ls -lh cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld
 -rw-r--r--. 1 user user 1000 31 dic 09.54 cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld
```

$_same_args_docs

```
julia> @scache "cute-cube" "./" fill(0, 3, parse(Int64, "3"), 3)
● [ 31/12/2021 09:54:54 ] Loading cute-cube from file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

```
"""
macro scache(type, common_cache_dir, ex)
	global _use_serialize = true

	:(@_scache $(esc(type)) $(esc(common_cache_dir)) $(esc(ex)))
end

"""
	scachejld(type, cache_dir, function_call)

Cache the result of `function_call` in directory `cache_dir` prefixing the saved file with
`type`.

This macro uses `JLD2` `@save` and `@load` macros to save and load cached files so it is
slower and less memory efficent than [`scache`](@ref) macro which uses `serialize`
which, on the other hand, is less portable between different julia versions.

$_ext_docs

## Examples
```
julia> using SimpleCaching

julia> SimpleCaching.settings.log = true
true

julia> @scachejld "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 10:31:15 ] Computing cute-cube...
● [ 31/12/2021 10:31:15 ] Computed cute-cube in 0.042 seconds (00:00:00)
● [ 31/12/2021 10:31:19 ] Saving cute-cube to file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld[.tmp]...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

julia> @scachejld "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 10:31:23 ] Loading cute-cube from file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

```

$_same_args_docs

```
julia> @scachejld "cute-cube" "./" fill(0, 3, round(Int64, 1.5 * 2), 3)
● [ 31/12/2021 10:32:06 ] Loading cute-cube from file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

```
"""
macro scachejld(type, common_cache_dir, ex)
	global _use_serialize = false

	:(@_scache $(esc(type)) $(esc(common_cache_dir)) $(esc(ex)))
end

"""
	scache_if(condition, type, cache_dir, function_call)

Cache the result of `function_call` only if `condition` is `true`.

Note that will not be loaded the cached result even if present.

For other parameters docs see [`scache_if`](@ref).

## Examples
```
julia> using SimpleCaching

julia> SimpleCaching.settings.log = true
true

julia> @scache_if true "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 09:54:08 ] Computing cute-cube...
● [ 31/12/2021 09:54:08 ] Computed cute-cube in 0.044 seconds (00:00:00)
● [ 31/12/2021 09:54:13 ] Saving cute-cube to file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld[.tmp]...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

julia> @scache_if true "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 09:54:24 ] Loading cute-cube from file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

 shell> ls -lh cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld
 -rw-r--r--. 1 user user 1000 31 dic 09.54 cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld
```

but passing a `false` `condition` (note there is no loading log):

```
julia> @scache_if true "cute-cube" "./" fill(0, 3, 3, 3)
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0
 ```
"""
macro scache_if(condition, type, common_cache_dir, ex)
	return quote
		if $(esc(condition))
			@scache $(esc(type)) $(esc(common_cache_dir)) $(esc(ex))
		else
			$(esc(ex))
		end
	end
end

"""
	scachejld_if(condition, type, cache_dir, function_call)

Cache the result of `function_call` only if `condition` is `true`.

Note that will not be loaded the cached result even if present.

For other parameters docs see [`scachejld`](@ref).

## Examples
```
julia> using SimpleCaching

julia> SimpleCaching.settings.log = true
true

julia> @scachejld_if true "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 09:54:08 ] Computing cute-cube...
● [ 31/12/2021 09:54:08 ] Computed cute-cube in 0.044 seconds (00:00:00)
● [ 31/12/2021 09:54:13 ] Saving cute-cube to file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld[.tmp]...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

julia> @scachejld_if true "cute-cube" "./" fill(0, 3, 3, 3)
● [ 31/12/2021 09:54:24 ] Loading cute-cube from file ./cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld...
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0

 shell> ls -lh cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld
 -rw-r--r--. 1 user user 1000 31 dic 09.54 cute-cube_8ad46882688c6820fc0b59db89cfe7f6ca494e3045d7ece8acba1027c4c03c45.jld
```

but passing a `false` `condition` (note there is no loading log):

```
julia> @scachejld_if true "cute-cube" "./" fill(0, 3, 3, 3)
3×3×3 Array{Int64, 3}:
[:, :, 1] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 2] =
 0  0  0
 0  0  0
 0  0  0

[:, :, 3] =
 0  0  0
 0  0  0
 0  0  0
 ```
"""
macro scachejld_if(condition, type, common_cache_dir, ex)
	return quote
		if $(esc(condition))
			@scachejld $(esc(type)) $(esc(common_cache_dir)) $(esc(ex))
		else
			$(esc(ex))
		end
	end
end


text:
let
    p = 3
    n = BigInt(8)
    t = (1:n)//n
    Ut,Vt = egrss.generators(t,p)

    K = egrss.full(Ut,Vt)
    xt = copy(t)
    @test ( K*xt == egrss.symv(Ut,Vt,xt) )

    Wt,c = egrss.ldl(Ut,Vt)
    L = egrss.full_tril(Ut,Wt,ones(eltype(Ut),n))
    D = Diagonal(c)
    @test ( K == L*D*L' )

    b = K*xt
    z = egrss.trsv(Ut,Wt,ones(eltype(Ut),n),b,'N')
    z = z ./ c
    x = egrss.trsv(Ut,Wt,ones(eltype(Ut),n),z,'T')
    @test ( x == xt )

end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
#NetwAnalysis: Circuit element tools
#-------------------------------------------------------------------------------

#=NOTE

Assumption: use ./ & .*
(Vectors represent lists, not arrays that support matrix operations)
=#

#==Helper functions
===============================================================================#
_zeroone(v::Number) = (zero(v), one(v))
_zeroone(v::DataMD) = (zeros(v), ones(v))
_zeroone(v) = (t=eltype(v); sz=size(v); return (zeros(t, sz), ones(t, sz))) #Arrays

#==Main functions
===============================================================================#

#Create series/shunt elements:
function shunt(::NPType{:ABCD}, y::TAdmittance{T}) where T
	y = y.v; (_0, _1) = _zeroone(y)
	m = Array{T}(undef, 2, 2) #Must build array manually
		m[1,1] = _1; m[1,2] = _0
		m[2,1] =  y; m[2,2] = _1
	return Network(:ABCD, m)
end
shunt(t::NPType{:ABCD}, z::TImpedance) = shunt(t, admittance(z))
shunt(t::NPType{:ABCD}, args...; kwargs...) = throw("series: Parameter list not supported")
shunt(t::NPType{TP}, args...; kwargs...) where TP = Network(t, shunt(t, args...; kwargs...))
shunt(s::Symbol, args...; kwargs...) = shunt(NPType(s), args...; kwargs...)

function series(::NPType{:ABCD}, z::TImpedance{T}) where T
	z = z.v; (_0, _1) = _zeroone(z)
	m = Array{T}(undef, 2, 2) #Must build array manually
		m[1,1] = _1; m[1,2] = z
		m[2,1] = _0; m[2,2] = _1
	return Network(:ABCD, m)
end
series(t::NPType{:ABCD}, y::TAdmittance) = series(t, impedance(y))
series(t::NPType{:ABCD}, args...; kwargs...) = throw("series: Parameter list not supported")
series(t::NPType{TP}, args...; kwargs...) where TP = Network(t, series(t, args...; kwargs...))
series(s::Symbol, args...; kwargs...) = series(NPType(s), args...; kwargs...)

#Last line


text:
function breakpoint_hook(f, bp)
    @info "Unhandled Hook, ignored" f bp
end

function breakpoint_hook(::typeof(JuliaInterpreter.breakpoint), bp::BreakpointSignature)
    @debug "Breakpoint created" bp
    make_interpretted_mode_overdub(bp)
end

function breakpoint_hook(::typeof(JuliaInterpreter.remove), bp::BreakpointSignature)
    @debug "Breakpoint Removed" bp
    remove_interpretted_mode_overdub(bp)
end

function breakpoint_hook(::typeof(JuliaInterpreter.update_states!), bp::BreakpointSignature)
    is_enabled = bp.enabled[]
    @debug "Breakpoint update" is_enabled
    if is_enabled
        make_interpretted_mode_overdub(bp)
    else
        remove_interpretted_mode_overdub(bp)
    end
end


text:
"""
Interior point algorithms for conic models.
"""
module Solvers

using DocStringExtensions
using Printf
using LinearAlgebra
using SparseArrays
import SuiteSparse
import LinearMaps
import IterativeSolvers
using Test
import Base.convert
import Hypatia.Cones
import Hypatia.Models
import Hypatia.SparseNonSymCache
import Hypatia.SparseSymCache
import Hypatia.diag_min
import Hypatia.update_fact
import Hypatia.inv_prod
import Hypatia.free_memory
import Hypatia.int_type
import Hypatia.nonsymm_fact_copy!
import Hypatia.symm_fact_copy!
import Hypatia.posdef_fact_copy!
import Hypatia.outer_prod!

const RealOrNothing = Union{Real, Nothing}

include("point.jl")

# solver termination status codes
@enum Status begin
    NotLoaded
    Loaded
    SolveCalled
    Optimal
    PrimalInfeasible
    DualInfeasible
    IllPosed
    PrimalInconsistent
    DualInconsistent
    SlowProgress
    IterationLimit
    TimeLimit
    NumericalFailure
    UnknownStatus
end

convert(::Type{String}, status::Status) = string(status)

abstract type Stepper{T <: Real} end

abstract type SystemSolver{T <: Real} end

"""
$(TYPEDEF)

Hypatia's interior point solver type. See source for options and defaults.
"""
mutable struct Solver{T <: Real}
    # main options
    verbose::Bool
    iter_limit::Int
    time_limit::Float64
    tol_rel_opt::T
    tol_abs_opt::T
    tol_feas::T
    tol_infeas::T
    tol_illposed::T
    tol_slow::T
    preprocess::Bool
    reduce::Bool
    rescale::Bool
    init_use_indirect::Bool
    init_tol_qr::T
    stepper::Stepper{T}
    syssolver::SystemSolver{T}

    # current status of the solver object and info
    status::Status
    solve_time::Float64
    num_iters::Int

    # helpful performance metrics for major subprocedures
    time_rescale::Float64 # rescale affine data
    time_initx::Float64 # preprocess dual equalities and finding initial x point
    time_inity::Float64 # preprocess primal equalities and finding initial y point
    time_unproc::Float64 # unprocess final point
    time_loadsys::Float64 # initialize/load system solver
    time_upsys::Float64 # update LHS and factorization etc for directions solving
    time_upfact::Float64 # update inner factorization for directions solving only
    time_uprhs::Float64 # update RHSs for directions
    time_getdir::Float64 # solve for directions
    time_search::Float64 # searches for alpha parameters

    # model and preprocessed model data
    orig_model::Models.Model{T}
    model::Models.Model{T}
    x_keep_idxs::AbstractVector{Int}
    y_keep_idxs::AbstractVector{Int}
    Ap_R::UpperTriangular{T, <:AbstractMatrix{T}}
    Ap_Q::Union{UniformScaling, AbstractMatrix{T}}
    reduce_cQ1
    reduce_Rpib0
    reduce_GQ1
    reduce_Ap_R
    reduce_Ap_Q
    reduce_y_keep_idxs
    reduce_row_piv_inv

    # current iterate
    point::Point{T}
    mu::T

    # result (solution) point
    result::Point{T}

    # residuals
    x_residual::Vector{T}
    y_residual::Vector{T}
    z_residual::Vector{T}
    tau_residual::T
    x_norm_res_t::T
    y_norm_res_t::T
    z_norm_res_t::T
    x_norm_res::T
    y_norm_res::T
    z_norm_res::T

    # direction solving helpers
    res_norm_cutoff::T # desired worst-case residual norm (inf) on direction
    max_ref_steps::Int # maximum number of iterative refinement steps

    # convergence parameters
    primal_obj_t::T
    dual_obj_t::T
    primal_obj::T
    dual_obj::T
    gap::T
    x_feas::T
    y_feas::T
    z_feas::T
    tau_feas::T

    # termination condition helpers
    x_conv_tol::T
    y_conv_tol::T
    z_conv_tol::T
    # TODO count how many slow with counter field and use option for max slow count
    prev_is_slow::Bool
    prev2_is_slow::Bool
    worst_dir_res::T

    # data scaling
    used_rescaling::Bool
    b_scale::Vector{T}
    c_scale::Vector{T}
    h_scale::Vector{T}

    function Solver{T}(;
        verbose::Bool = true,
        iter_limit::Int = 1000,
        time_limit::Real = Inf,
        tol_rel_opt::RealOrNothing = nothing,
        tol_abs_opt::RealOrNothing = nothing,
        tol_feas::RealOrNothing = nothing,
        tol_infeas::RealOrNothing = nothing,
        tol_illposed::RealOrNothing = nothing,
        default_tol_power::RealOrNothing = nothing,
        default_tol_relax::RealOrNothing = nothing,
        tol_slow::Real = 1e-3,
        preprocess::Bool = true,
        reduce::Bool = true,
        rescale::Bool = true,
        init_use_indirect::Bool = false,
        init_tol_qr::Real = 1000 * eps(T),
        stepper::Stepper{T} = default_stepper(T),
        syssolver::SystemSolver{T} = default_syssolver(T),
        ) where {T <: Real}
        if isa(syssolver, QRCholSystemSolver{T})
            @assert preprocess # require preprocessing for QRCholSystemSolver # TODO only need primal eq preprocessing or reduction
        end
        if reduce
            @assert preprocess # cannot use reduction without preprocessing # TODO only need primal eq preprocessing
        end
        # @assert !(init_use_indirect && preprocess) # cannot use preprocessing and indirect methods for initial point

        if isnothing(default_tol_power)
            default_tol_power = (T <: LinearAlgebra.BlasReal ? 0.5 : 0.4)
        end
        default_tol_power = T(default_tol_power)
        default_tol_loose = eps(T) ^ default_tol_power
        default_tol_tight = eps(T) ^ (T(1.5) * default_tol_power)
        if !isnothing(default_tol_relax)
            default_tol_loose *= T(default_tol_relax)
            default_tol_tight *= T(default_tol_relax)
        end
        if isnothing(tol_rel_opt)
            tol_rel_opt = default_tol_loose
        end
        if isnothing(tol_abs_opt)
            tol_abs_opt = default_tol_tight
        end
        if isnothing(tol_feas)
            tol_feas = default_tol_loose
        end
        if isnothing(tol_infeas)
            tol_infeas = default_tol_tight
        end
        if isnothing(tol_illposed)
            tol_illposed = default_tol_tight / 100
        end
        @assert min(tol_rel_opt, tol_abs_opt, tol_feas,
            tol_infeas, tol_illposed, tol_slow) >= 0

        solver = new{T}()

        solver.verbose = verbose
        solver.iter_limit = iter_limit
        solver.time_limit = time_limit
        solver.tol_rel_opt = tol_rel_opt
        solver.tol_abs_opt = tol_abs_opt
        solver.tol_feas = tol_feas
        solver.tol_infeas = tol_infeas
        solver.tol_illposed = tol_illposed
        solver.tol_slow = tol_slow
        solver.preprocess = preprocess
        solver.reduce = reduce
        solver.rescale = rescale
        solver.init_use_indirect = init_use_indirect
        solver.init_tol_qr = init_tol_qr
        solver.stepper = stepper
        solver.syssolver = syssolver
        solver.status = NotLoaded

        return solver
    end
end

default_stepper(T) = CombinedStepper{T}()
default_syssolver(T) = QRCholDenseSystemSolver{T}()

function solve(solver::Solver{T}) where {T <: Real}
    @assert solver.status == Loaded
    solver.status = SolveCalled
    start_time = time()
    solver.num_iters = 0
    solver.solve_time = NaN

    solver.time_rescale = 0
    solver.time_initx = 0
    solver.time_inity = 0
    solver.time_unproc = 0
    solver.time_loadsys = 0
    solver.time_upsys = 0
    solver.time_upfact = 0
    solver.time_uprhs = 0
    solver.time_getdir = 0
    solver.time_search = 0

    solver.res_norm_cutoff = 0
    solver.max_ref_steps = 5

    solver.x_norm_res_t = NaN
    solver.y_norm_res_t = NaN
    solver.z_norm_res_t = NaN
    solver.x_norm_res = NaN
    solver.y_norm_res = NaN
    solver.z_norm_res = NaN

    solver.primal_obj_t = NaN
    solver.dual_obj_t = NaN
    solver.primal_obj = NaN
    solver.dual_obj = NaN
    solver.gap = NaN
    solver.x_feas = NaN
    solver.y_feas = NaN
    solver.z_feas = NaN
    solver.tau_feas = NaN

    # preprocess and find initial point
    orig_model = solver.orig_model
    result = solver.result = Point(orig_model)
    # copy original model to solver.model, which may be modified
    model = solver.model = Models.Model{T}(
        orig_model.c, orig_model.A, orig_model.b, orig_model.G, orig_model.h,
        orig_model.cones, obj_offset = orig_model.obj_offset)
    (init_z, init_s) = initialize_cone_point(solver.orig_model)

    solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver)

    if solver.reduce
        # TODO don't find point / unnecessary stuff before reduce
        solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, true)
        solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s)
    else
        solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s)
        solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, false)
    end

    if solver.status == SolveCalled
        point = solver.point = Point(model)
        point.x .= init_x
        point.y .= init_y
        point.z .= init_z
        point.s .= init_s
        point.tau[] = one(T)
        point.kap[] = one(T)
        calc_mu(solver)
        if isnan(solver.mu) || abs(one(T) - solver.mu) > sqrt(eps(T))
            @warn("initial mu is $(solver.mu) but should be 1 (this could " *
                "indicate a problem with cone barrier oracles)")
        end
        Cones.load_point.(model.cones, point.primal_views)
        Cones.load_dual_point.(model.cones, point.dual_views)

        # setup iteration helpers
        solver.x_residual = zero(model.c)
        solver.y_residual = zero(model.b)
        solver.z_residual = zero(model.h)
        solver.tau_residual = 0

        solver.x_conv_tol = inv(1 + norm(model.c, Inf))
        solver.y_conv_tol = inv(1 + norm(model.b, Inf))
        solver.z_conv_tol = inv(1 + norm(model.h, Inf))
        solver.prev_is_slow = false
        solver.prev2_is_slow = false
        solver.worst_dir_res = 0

        stepper = solver.stepper
        load(stepper, solver)
        solver.time_loadsys = @elapsed load(solver.syssolver, solver)

        solver.verbose && print_header(stepper, solver)
        flush(stdout)

        # iterate from initial point
        while true
            improv = calc_convergence_params(solver)

            if solver.verbose
                print_iteration(stepper, solver)
                flush(stdout)
            end

            check_convergence(solver) && break

            if solver.num_iters == solver.iter_limit
                solver.verbose && println("iteration limit reached; terminating")
                solver.status = IterationLimit
                break
            end

            if time() - start_time >= solver.time_limit
                solver.verbose && println("time limit reached; terminating")
                solver.status = TimeLimit
                break
            end

            if expect_improvement(stepper)
                if improv < solver.tol_slow
                    if solver.prev_is_slow && solver.prev2_is_slow
                        if solver.verbose
                            println("slow progress in consecutive " *
                                "iterations; terminating")
                        end
                        solver.status = SlowProgress
                        break
                    else
                        solver.prev2_is_slow = solver.prev_is_slow
                        solver.prev_is_slow = true
                    end
                else
                    solver.prev2_is_slow = solver.prev_is_slow
                    solver.prev_is_slow = false
                end
            end

            solver.res_norm_cutoff = T(1e-4) * max(solver.x_norm_res,
                solver.y_norm_res, solver.z_norm_res, solver.tau_feas)
            solver.worst_dir_res = 0

            step(stepper, solver) || break
            flush(stdout)
            calc_mu(solver)

            if min(point.tau[], point.kap[], solver.mu) <= 0
                @warn("numerical failure: tau is $(point.tau[]), " *
                    "kappa is $(point.kap[]), mu is $(solver.mu); terminating")
                solver.status = NumericalFailure
                break
            end

            flush(stdout)
            solver.num_iters += 1
        end

        # finalize result point
        solver.time_unproc = @elapsed postprocess(solver)
    end

    solver.solve_time = time() - start_time

    # free memory used by some system solvers
    free_memory(solver.syssolver)

    if solver.verbose
        println("\nstatus is $(solver.status) after $(solver.num_iters) " *
            "iterations and $(trunc(solver.solve_time, digits=3)) seconds\n")
    end
    flush(stdout)

    return solver
end

function calc_mu(solver::Solver{T}) where {T <: Real}
    point = solver.point
    solver.mu = (dot(point.z, point.s) + point.tau[] * point.kap[]) /
        (solver.model.nu + 1)
    return solver.mu
end

function calc_convergence_params(solver::Solver{T}) where {T <: Real}
    model = solver.model
    point = solver.point
    tau = point.tau[]

    # x_residual = -A'*y - G'*z - c*tau
    x_residual = solver.x_residual
    mul!(x_residual, model.G', point.z)
    mul!(x_residual, model.A', point.y, true, true)
    solver.x_norm_res_t = norm(x_residual, Inf)
    @. x_residual += model.c * tau
    solver.x_norm_res = norm(x_residual, Inf) / tau
    @. x_residual *= -1
    x_feas = solver.x_norm_res * solver.x_conv_tol

    # y_residual = A*x - b*tau
    y_residual = solver.y_residual
    mul!(y_residual, model.A, point.x)
    solver.y_norm_res_t = norm(y_residual, Inf)
    @. y_residual -= model.b * tau
    solver.y_norm_res = norm(y_residual, Inf) / tau
    y_feas = solver.y_norm_res * solver.y_conv_tol

    # z_residual = s + G*x - h*tau
    z_residual = solver.z_residual
    mul!(z_residual, model.G, point.x)
    @. z_residual += point.s
    solver.z_norm_res_t = norm(z_residual, Inf)
    @. z_residual -= model.h * tau
    solver.z_norm_res = norm(z_residual, Inf) / tau
    z_feas = solver.z_norm_res * solver.z_conv_tol

    # tau_residual = c'*x + b'*y + h'*z + kap
    solver.primal_obj_t = dot(model.c, point.x)
    solver.dual_obj_t = -dot(model.b, point.y) - dot(model.h, point.z)
    solver.tau_residual = solver.primal_obj_t - solver.dual_obj_t + point.kap[]
    tau_feas = abs(solver.tau_residual)

    # check improvement
    improv = zero(T)
    for (curr, prev) in ((x_feas, solver.x_feas), (y_feas, solver.y_feas),
        (z_feas, solver.z_feas), (tau_feas, solver.tau_feas))
        if isnan(prev) || isnan(curr)
            continue
        end
        improv = max(improv, (prev - curr) / (abs(prev) + eps(T)))
    end
    solver.x_feas = x_feas
    solver.y_feas = y_feas
    solver.z_feas = z_feas
    solver.tau_feas = tau_feas

    # gap
    solver.primal_obj = solver.primal_obj_t / tau + model.obj_offset
    solver.dual_obj = solver.dual_obj_t / tau + model.obj_offset
    solver.gap = dot(point.z, point.s)

    return improv
end

function check_convergence(solver::Solver{T}) where {T <: Real}
    # check convergence criteria
    # TODO nearly primal or dual infeasible or nearly optimal cases?
    tau = solver.point.tau[]
    primal_obj_t = solver.primal_obj_t
    dual_obj_t = solver.dual_obj_t

    is_feas = (max(solver.x_feas, solver.y_feas, solver.z_feas) <= solver.tol_feas)
    is_abs_opt = (solver.gap <= solver.tol_abs_opt)
    is_rel_opt = (min(solver.gap / tau, abs(primal_obj_t - dual_obj_t)) <=
        solver.tol_rel_opt * max(tau, min(abs(primal_obj_t), abs(dual_obj_t))))
    if is_feas && (is_abs_opt || is_rel_opt)
        solver.verbose && println("optimal solution found; terminating")
        solver.status = Optimal
        return true
    end

    if dual_obj_t > eps(T) && solver.x_norm_res_t <= solver.tol_infeas * dual_obj_t
        solver.verbose && println("primal infeasibility detected; terminating")
        solver.status = PrimalInfeasible
        solver.primal_obj = primal_obj_t
        solver.dual_obj = dual_obj_t
        return true
    end

    if primal_obj_t < -eps(T) && max(solver.y_norm_res_t, solver.z_norm_res_t) <=
        solver.tol_infeas * -primal_obj_t
        solver.verbose && println("dual infeasibility detected; terminating")
        solver.status = DualInfeasible
        solver.primal_obj = primal_obj_t
        solver.dual_obj = dual_obj_t
        return true
    end

    # TODO experiment with ill-posedness check
    if solver.mu <= solver.tol_illposed && tau <=
        solver.tol_illposed * min(one(T), solver.point.kap[])
        solver.verbose && println("ill-posedness detected; terminating")
        solver.status = IllPosed
        return true
    end

    return false
end

function initialize_cone_point(model::Models.Model{T}) where {T <: Real}
    init_z = zeros(T, model.q)
    init_s = zeros(T, model.q)

    for (cone, idxs) in zip(model.cones, model.cone_idxs)
        Cones.setup_data!(cone)
        primal_k = view(Cones.use_dual_barrier(cone) ? init_z : init_s, idxs)
        dual_k = view(Cones.use_dual_barrier(cone) ? init_s : init_z, idxs)
        Cones.set_initial_point!(primal_k, cone)
        Cones.load_point(cone, primal_k)
        @assert Cones.is_feas(cone)
        g = Cones.grad(cone)
        @. dual_k = -g
        Cones.load_dual_point(cone, dual_k)
        @assert Cones.is_dual_feas(cone)
    end

    return (init_z, init_s)
end

get_status(solver::Solver) = solver.status
get_solve_time(solver::Solver) = solver.solve_time
get_num_iters(solver::Solver) = solver.num_iters

get_primal_obj(solver::Solver) = solver.primal_obj
get_dual_obj(solver::Solver) = solver.dual_obj

get_s(solver::Solver) = copy(solver.result.s)
get_z(solver::Solver) = copy(solver.result.z)
get_x(solver::Solver) = copy(solver.result.x)
get_y(solver::Solver) = copy(solver.result.y)

get_tau(solver::Solver) = solver.point.tau[]
get_kappa(solver::Solver) = solver.point.kap[]
get_mu(solver::Solver) = solver.mu

function load(solver::Solver{T}, model::Models.Model{T}) where {T <: Real}
    # @assert solver.status == NotLoaded # TODO maybe want a reset function that just keeps options
    solver.orig_model = model
    solver.status = Loaded
    return solver
end

include("process.jl")

include("search.jl")

include("steppers/common.jl")

include("systemsolvers/common.jl")

# release memory used by sparse system solvers
free_memory(::SystemSolver) = nothing
free_memory(syssolver::Union{NaiveSparseSystemSolver, SymIndefSparseSystemSolver}) =
    free_memory(syssolver.fact_cache)

# verbose helpers
function print_header(stepper::Stepper, solver::Solver)
    println()
    @printf("%5s %12s %12s |%9s ", "iter", "p_obj", "d_obj", "abs_gap")
    if iszero(solver.model.p)
        @printf("%9s %9s ", "x_feas", "z_feas")
    else
        @printf("%9s %9s %9s ", "x_feas", "y_feas", "z_feas")
    end
    @printf("|%9s %9s %9s |%8s %8s ", "tau", "kap", "mu", "dir_res", "prox")

    print_header_more(stepper, solver)
    println()
    return
end
print_header_more(stepper::Stepper, solver::Solver) = nothing

function print_iteration(stepper::Stepper, solver::Solver)
    @printf("%5d %12.4e %12.4e |%9.2e ", solver.num_iters, solver.primal_obj,
        solver.dual_obj, solver.gap)
    if iszero(solver.model.p)
        @printf("%9.2e %9.2e ", solver.x_feas, solver.z_feas)
    else
        @printf("%9.2e %9.2e %9.2e ", solver.x_feas, solver.y_feas, solver.z_feas)
    end
    @printf("|%9.2e %9.2e %9.2e |", solver.point.tau[], solver.point.kap[],
        solver.mu)

    if !iszero(solver.num_iters)
        @printf("%8.1e %8.1e ", solver.worst_dir_res, stepper.searcher.prox)
        print_iteration_more(stepper, solver)
    end
    println()
    return
end
print_iteration_more(stepper::Stepper, solver::Solver) = nothing

end


text:
# This file was generated, do not modify it. # hide
nnregressor.epochs = 15

MLJ.fit!(mach, rows=train, verbosity=3)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ImageMagick_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ImageMagick")
JLLWrappers.@generate_main_file("ImageMagick", UUID("c73af94c-d91f-53ed-93a7-00f77d67a9d7"))
end  # module ImageMagick_jll


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to apply the Knitro multi-start
# procedure to solve the following simple nonlinear optimization
# problem.  This model is test problem HS15 from the Hock &
# Schittkowski collection.
#
# min   100(x1 - x0^2)^2 +(1 - x0)^2
# s.t.  x0 x1 >= 1
#       x0 + x1^2 >= 0
#       x0 <= 0.5
#
# Some solves should converge to the minimum at(0.5, 2.0),
# with final objective = 306.5, while others may converge to
# another local minimum at(-0.79212, -1.26243), with final
# objective = 360.4.
#
# The example also shows how to set a callback function to
# perform some user-defined task after each multi-start solve.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO
using Test

function example_multistart()
    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalF                                       *
    #*------------------------------------------------------------------*
    # The signature of this function matches KN_eval_callback in knitro.jl.
    # Only "obj" is set in the KNITRO.KN_eval_result structure.
    function callbackEvalF(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate nonlinear objective
        dTmp = x[2] - x[1] * x[1]
        evalResult.obj[1] = 100.0 * (dTmp * dTmp) + ((1.0 - x[1]) * (1.0 - x[1]))

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalG                                       *
    #*------------------------------------------------------------------*
    # The signature of this function matches KN_eval_callback in knitro.h.
    # Only "objGrad" is set in the KN_eval_result structure.
    function callbackEvalG!(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        # Evaluate gradient of nonlinear objective
        dTmp = x[2] - x[1] * x[1]
        evalResult.objGrad[1] = (-400.0 * dTmp * x[1]) - (2.0 * (1.0 - x[1]))
        evalResult.objGrad[2] = 200.0 * dTmp

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackEvalH                                       *
    #*------------------------------------------------------------------*
    # The signature of this function matches KN_eval_callback in knitro.h.
    # Only "hess" and "hessVec" are set in the KN_eval_result structure.
    function callbackEvalH!(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x
        # Scale objective component of hessian by sigma
        sigma = evalRequest.sigma

        # Evaluate the hessian of the nonlinear objective.
        # Note: Since the Hessian is symmetric, we only provide the
        #       nonzero elements in the upper triangle(plus diagonal).
        #       These are provided in row major ordering as specified
        #       by the setting KN_DENSE_ROWMAJOR in "KN_set_cb_hess()".
        # Note: The Hessian terms for the quadratic constraints
        #       will be added internally by Knitro to form
        #       the full Hessian of the Lagrangian.
        evalResult.hess[1] = sigma * ((-400.0 * x[2]) + (1200.0 * x[1] * x[1]) + 2.0) #(0,0)
        evalResult.hess[2] = sigma * (-400.0 * x[1]) #(0,1)
        evalResult.hess[3] = sigma * 200.0           #(1,1)

        return 0
    end

    #*------------------------------------------------------------------*
    #*     FUNCTION callbackMSProcess                                   *
    #*------------------------------------------------------------------*
    # The signature of this function matches KNITRO.KN_user_callback in knitro.h.
    # Argument "kcSub" is the context pointer for the last multi-start
    # subproblem solved inside Knitro.
    function  callbackMSProcess(kcSub, x, lambda_, userParams)
        # Print solution of the just completed multi-start solve.

        n = KNITRO.KN_get_number_vars(kc)
        println("callbackMSProcess: ")
        println("    Last solution: obj= ", KNITRO.KN_get_obj_value(kc))
        for i in 1:n
            println("                   x[$i]= ", x[i])
        end

        return 0
    end

    #*------------------------------------------------------------------*
    #*     main                                                         *
    #*------------------------------------------------------------------*
    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Initialize knitro with the problem definition.

    # Add the variables and set their bounds.
    # Note: any unset lower bounds are assumed to be
    # unbounded below and any unset upper bounds are
    # assumed to be unbounded above.
    n = 2
    KNITRO.KN_add_vars(kc, n)
    KNITRO.KN_set_var_lobnds(kc,  [-KNITRO.KN_INFINITY, -KNITRO.KN_INFINITY]) # not necessary since infinite
    KNITRO.KN_set_var_upbnds(kc,  [0.5, KNITRO.KN_INFINITY])
    # Define an initial point.  If not set, Knitro will generate one.
    KNITRO.KN_set_var_primal_init_values(kc, [-2.0, 1.0])

    # Add the constraints and set their lower bounds
    m = 2
    KNITRO.KN_add_cons(kc, m)
    KNITRO.KN_set_con_lobnds(kc,  [1.0, 0.0])

    # Both constraints are quadratic so we can directly load all the
    # structure for these constraints.

    # First load quadratic structure x0*x1 for the first constraint
    KNITRO.KN_add_con_quadratic_struct(kc, 0, 0, 1, 1.0)

    # Load structure for the second constraint.  below we add the linear
    # structure and the quadratic structure separately, though it
    # is possible to add both together in one call to
    # "KNITRO.KN_add_con_quadratic_struct()" since this api function also
    # supports adding linear terms.

    # Add linear term x0 in the second constraint
    KNITRO.KN_add_con_linear_struct(kc, 1, 0, 1.0)

    # Add quadratic term x1^2 in the second constraint
    KNITRO.KN_add_con_quadratic_struct(kc, 1, 1, 1, 1.0)

    # Add a callback function "callbackEvalF" to evaluate the nonlinear
    # (non-quadratic) objective.  Note that the linear and
    # quadratic terms in the objective could be loaded separately
    # via "KNITRO.KN_add_obj_linear_struct()" / "KNITRO.KN_add_obj_quadratic_struct()".
    # However, for simplicity, we evaluate the whole objective
    # function through the callback.
    cb = KNITRO.KN_add_objective_callback(kc, callbackEvalF)

    # Also add a callback function "callbackEvalG" to evaluate the
    # objective gradient.  If not provided, Knitro will approximate
    # the gradient using finite-differencing.  However, we recommend
    # providing callbacks to evaluate the exact gradients whenever
    # possible as this can drastically improve the performance of Knitro.
    # We specify the objective gradient in "dense" form for simplicity.
    # However for models with many constraints, it is important to specify
    # the non-zero sparsity structure of the constraint gradients
    # (i.e. Jacobian matrix) for efficiency(this is true even when using
    # finite-difference gradients).
    KNITRO.KN_set_cb_grad(kc, cb, callbackEvalG!)

    # Add a callback function "callbackEvalH" to evaluate the Hessian
    # (i.e. second derivative matrix) of the objective.  If not specified,
    # Knitro will approximate the Hessian. However, providing a callback
    # for the exact Hessian(as well as the non-zero sparsity structure)
    # can greatly improve Knitro performance and is recommended if possible.
    # Since the Hessian is symmetric, only the upper triangle is provided.
    # Again for simplicity, we specify it in dense(row major) form.
    KNITRO.KN_set_cb_hess(kc, cb, KNITRO.KN_DENSE_ROWMAJOR, callbackEvalH!)

    # Specify that the user is able to provide evaluations
    # of the hessian matrix without the objective component.
    # turned off by default but should be enabled if possible.
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_HESSIAN_NO_F, KNITRO.KN_HESSIAN_NO_F_ALLOW)

    # Set minimize or maximize(if not set, assumed minimze)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)

    # Example of how to register a callback function that performs some
    # task after each multistart solve.
    KNITRO.KN_set_ms_process_callback(kc, callbackMSProcess)
    #
    # Disable automatic scaling.
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_SCALE, KNITRO.KN_SCALE_NO)

    # Enable multi-start
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_MULTISTART, KNITRO.KN_MULTISTART_YES)

    # Perform multistart in parallel using max number of available threads
    nThreads = Sys.CPU_THREADS
    if nThreads > 1
        println("Running Knitro multistart in parallel with $nThreads threads.")
        KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_PAR_MSNUMTHREADS, 1)
    end

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.

    nStatus = KNITRO.KN_solve(kc)

    # An example of obtaining solution information.
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)
    println("Optimal objective value  = ", objSol)
    println("Optimal x(with corresponding multiplier)")
    for i in 1:n
        println("  x[$i] = ", x[i], "(lambda = ",  lambda_[m+i], ")")
    end
    println("Optimal constraint values(with corresponding multiplier)")
    c = KNITRO.KN_get_con_values(kc)
    for j in 1:m
        println("  c[$j] = ", c[j], "(lambda = ",  lambda_[j], ")")
    end
    println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
    println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example HS15 multistart" begin
        @test nStatus == 0
        @test objSol  ≈ 306.5
        @test x ≈ [0.5, 2]
    end
end

example_multistart()


text:
module API

using Scratch

using unixODBC_jll
const unixODBC_dm = unixODBC_jll.libodbc
const unixODBC_inst = unixODBC_jll.libodbcinst

@static if Sys.iswindows()
    const odbc32_dm = "odbc32"
    const odbc32_inst = "odbccp32"
    const iODBC_dm = odbc32_dm
    const iODBC_inst = odbc32_inst
else
    using iODBC_jll
    const iODBC_dm = iODBC_jll.libiodbc
    const iODBC_inst = iODBC_jll.libiodbcinst
    const odbc32_dm = unixODBC_jll.libodbc
    const odbc32_inst = unixODBC_jll.libodbcinst
end

@enum DM_TYPE unixODBC iODBC odbc32

@static if Sys.iswindows()
    const odbc_dm = Ref(odbc32)
elseif Sys.isapple()
    const odbc_dm = Ref(iODBC)
else
    const odbc_dm = Ref(unixODBC)
end

function setunixODBC(; kw...)
    odbc_dm[] = unixODBC
    setupenv(; kw...)
    return
end

function setiODBC(; kw...)
    odbc_dm[] = iODBC
    setupenv(; kw...)
    return
end

function setodbc32(; kw...)
    odbc_dm[] = odbc32
    setupenv(; kw...)
    return
end

const SQLWCHAR = Cushort
const SQLWCHAR32 = UInt32

sqlwcharsize() = odbc_dm[] == iODBC ? SQLWCHAR32 : SQLWCHAR

include("consts.jl")
include("apitypes.jl")

str(x::Vector{UInt8}, len) = String(x[1:len])
str(x::Vector{UInt16}, len) = transcode(String, view(x, 1:len))
str(x::Vector{UInt32}, len) = transcode(String, x[1:len])
str(x::Vector{UInt32}) = unsafe_string(pointer(transcode(UInt8, x)))

function swapsqlwchar(expr)
    for i = 1:length(expr.args)
        if expr.args[i] == :(Ptr{SQLWCHAR})
            expr.args[i] = :(Ptr{SQLWCHAR32})
        end
    end
    return expr
end

macro odbc(func,args,vals...)
    esc(quote
        ret = SQL_SUCCESS
        @static if Sys.iswindows() # odbc_dm[] == odbc32
            # This branch is guarded by `@static` to avoid issues on Apple Silicon
            while true
                ret = ccall( ($func, odbc32_dm), stdcall, SQLRETURN, $args, $(vals...))
                ret == SQL_STILL_EXECUTING || break
                sleep(0.001)
            end
        else
            if odbc_dm[] == iODBC
                while true
                    ret = ccall( ($func, iODBC_dm), SQLRETURN, $(swapsqlwchar(args)), $(vals...))
                    ret == SQL_STILL_EXECUTING || break
                    sleep(0.001)
                end
            elseif odbc_dm[] == unixODBC
                while true
                    ret = ccall( ($func, unixODBC_dm), SQLRETURN, $args, $(vals...))
                    ret == SQL_STILL_EXECUTING || break
                    sleep(0.001)
                end
            end
        end
        ret
    end)
end

macro odbcinst(func,args,vals...)
    esc(quote
        @static if Sys.iswindows() # odbc_dm[] == odbc32
            # This branch is guarded by `@static` to avoid issues on Apple Silicon
            ccall( ($func, odbc32_inst), stdcall, SQLRETURN, $args, $(vals...))
        else
            if odbc_dm[] == iODBC
                ccall( ($func, iODBC_inst), SQLRETURN, $(swapsqlwchar(args)), $(vals...))
            elseif odbc_dm[] == unixODBC
                ccall( ($func, unixODBC_inst), SQLRETURN, $args, $(vals...))
            end
        end
    end)
end

macro checksuccess(h, expr)
    esc(quote
        ret = $expr
        if ret == SQL_SUCCESS_WITH_INFO
            @debug diagnostics($h)
        elseif ret == SQL_ERROR || ret == SQL_INVALID_HANDLE
            error(diagnostics($h))
        end
        ret
    end)
end

const SQL_HANDLE_ENV  = Int16(1)
const SQL_HANDLE_DBC  = Int16(2)
const SQL_HANDLE_STMT = Int16(3)
const SQL_HANDLE_DESC = Int16(4)
const SQL_NULL_HANDLE = C_NULL

function SQLAllocHandle(handletype::SQLSMALLINT, parenthandle::Ptr{Cvoid}, handle::Ref{Ptr{Cvoid}})
    @odbc(:SQLAllocHandle,
        (SQLSMALLINT, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
        handletype, parenthandle, handle)
end

function SQLFreeHandle(handletype::SQLSMALLINT,handle::Ptr{Cvoid})
    @odbc(:SQLFreeHandle,
        (SQLSMALLINT, Ptr{Cvoid}),
        handletype, handle)
end

const SQL_ATTR_ODBC_VERSION = 200
const SQL_OV_ODBC3 = 3

function SQLSetEnvAttr(env_handle::Ptr{Cvoid}, attribute, value)
    @odbc(:SQLSetEnvAttr,
        (Ptr{Cvoid}, SQLINTEGER, UInt, SQLINTEGER),
        env_handle, attribute, value, 0)
end

mutable struct Handle
    type::Int16
    ptr::Ptr{Cvoid}
    function Handle(type, parent=SQL_NULL_HANDLE)
        ref = Ref{Ptr{Cvoid}}()
        @checksuccess parent SQLAllocHandle(type, parent isa Handle ? parent.ptr : parent, ref)
        ptr = ref[]
        ptr == C_NULL && error("error creating ODBC.API.Handle; null pointer encountered")
        if type == SQL_HANDLE_ENV
            @checksuccess parent SQLSetEnvAttr(ptr, SQL_ATTR_ODBC_VERSION, SQL_OV_ODBC3)
        end
        h = new(type, ptr)
        finalizer(h) do x
            if x.ptr != C_NULL
                if x.type == SQL_HANDLE_DBC
                    SQLDisconnect(x.ptr)
                end
                SQLFreeHandle(x.type, x.ptr)
                x.ptr = C_NULL
            end
        end
        return h
    end
end

gettype(h::Handle) = h.type
getptr(h::Handle) = h.ptr == C_NULL ? error("invalid odbc handle") : h.ptr
getptr(x::Ptr) = x

const ODBC_ENV = Ref{Handle}()

const ODBCCONFIGDIR = Ref{String}()
const ODBCINI = Ref{String}()
const ODBCINSTINI = Ref{String}()

function setupenv(; kw...)
    ODBCCONFIGDIR[] = @get_scratch!("odbcconfig")
    ODBCINI[] = abspath(joinpath(ODBCCONFIGDIR[], "odbc.ini"))
    ODBCINSTINI[] = abspath(joinpath(ODBCCONFIGDIR[], "odbcinst.ini"))
    if isdefined(ODBC_ENV, :x)
        finalize(ODBC_ENV[])
    end
    if !haskey(ENV, "OVERRIDE_ODBCJL_CONFIG")
        delete!(ENV, "ODBCINI")
        delete!(ENV, "ODBCSYSINI")
        delete!(ENV, "ODBCINSTINI")
        ENV["ODBCINI"] = ODBCINI[]
        if odbc_dm[] == iODBC
            ENV["ODBCINSTINI"] = ODBCINSTINI[]
        elseif odbc_dm[] == unixODBC
            ENV["ODBCSYSINI"] = realpath(ODBCCONFIGDIR[])
        end
    end
    for (k, v) in pairs(kw)
        ENV[string(k)] = v
    end
    ODBC_ENV[] = Handle(SQL_HANDLE_ENV)
    return
end

const ODBC_DBC = Ref{Ptr{Cvoid}}()

function setdebug(trace, tracefile)
    if !isdefined(ODBC_DBC, :x) || ODBC_DBC[] == C_NULL
        SQLAllocHandle(SQL_HANDLE_DBC, getptr(ODBC_ENV[]), ODBC_DBC)
    end
    if trace
        SQLSetConnectAttr(ODBC_DBC[], SQL_ATTR_TRACE, SQL_OPT_TRACE_OFF)
        if tracefile != ""
            SQLSetConnectAttr(ODBC_DBC[], SQL_ATTR_TRACEFILE, tracefile)
        end
        SQLSetConnectAttr(ODBC_DBC[], SQL_ATTR_TRACE, SQL_OPT_TRACE_ON)
    else
        SQLSetConnectAttr(ODBC_DBC[], SQL_ATTR_TRACE, SQL_OPT_TRACE_OFF)
    end
end

function __init__()
    setupenv()
    return
end

const SQL_COMMIT = 0
const SQL_ROLLBACK = 1

function SQLEndTran(dbc, type)
    @odbc(:SQLEndTran,
        (SQLSMALLINT, Ptr{Cvoid}, SQLSMALLINT),
        SQL_HANDLE_DBC, getptr(dbc), type)
end

endtran(dbc, commit::Bool) = @checksuccess dbc SQLEndTran(dbc, commit ? SQL_COMMIT : SQL_ROLLBACK)

const SQL_ATTR_TRACE = SQLINTEGER(104)
const SQL_ATTR_TRACEFILE = SQLINTEGER(105)
const SQL_OPT_TRACE_OFF = SQLUINTEGER(0)
const SQL_OPT_TRACE_ON = SQLUINTEGER(1)

const SQL_ATTR_AUTOCOMMIT = SQLINTEGER(102)
const SQL_AUTOCOMMIT_OFF = SQLUINTEGER(0)
const SQL_AUTOCOMMIT_ON = SQLUINTEGER(1)

function SQLSetConnectAttr(dbc, attr::SQLINTEGER, value::SQLUINTEGER)
    @odbc(:SQLSetConnectAttr,
        (Ptr{Cvoid}, SQLINTEGER, SQLUINTEGER, SQLINTEGER),
        getptr(dbc), attr, value, SQL_IS_UINTEGER)
end

setcommitmode(dbc, on::Bool) = @checksuccess dbc SQLSetConnectAttr(dbc, SQL_ATTR_AUTOCOMMIT, on ? SQL_AUTOCOMMIT_ON : SQL_AUTOCOMMIT_OFF)

function SQLSetConnectAttr(dbc, attr::SQLINTEGER, value::String)
    @odbc(:SQLSetConnectAttr,
        (Ptr{Cvoid}, SQLINTEGER, Ptr{UInt8}, SQLINTEGER),
        getptr(dbc), attr, value, sizeof(value))
end

const SQL_ATTR_CONNECTION_DEAD = 1209
const SQL_CD_TRUE = 1
const SQL_CD_FALSE = 0

function SQLGetConnectAttr(dbc, attr, ret)
    ref = Ref{SQLINTEGER}()
    @odbc(:SQLGetConnectAttr,
        (Ptr{Cvoid}, SQLINTEGER, Ref{SQLUINTEGER}, SQLINTEGER, Ref{SQLINTEGER}),
        getptr(dbc), attr, ret, 0, ref)
end

function getconnectattr(dbc, attr)
    ret = Ref{SQLUINTEGER}()
    SQLGetConnectAttr(dbc, attr, ret)
    return ret[]
end

const SQL_BATCH_SUPPORT = 121
const SQL_BS_SELECT_EXPLICIT = 0x00000001
const SQL_BS_ROW_COUNT_EXPLICIT = 0x00000002
const SQL_BS_SELECT_PROC = 0x00000004
const SQL_BS_ROW_COUNT_PROC = 0x00000008

const SQL_FORWARD_ONLY_CURSOR_ATTRIBUTES2 = 147
const SQL_CA2_CRC_EXACT = 0x00001000
const SQL_CA2_CRC_APPROXIMATE = 0x00002000

function getinfosqluinteger(dbc::Handle, type=121)
    ref = Ref{SQLUINTEGER}()
    len = Ref{SQLSMALLINT}()
    @checksuccess dbc @odbc(:SQLGetInfo,
        (Ptr{Cvoid}, SQLUSMALLINT, Ref{SQLUINTEGER}, SQLSMALLINT, Ref{SQLSMALLINT}),
        getptr(dbc), type, ref, 0, len)
    return ref[]
end

const SQL_IDENTIFIER_QUOTE_CHAR = 29

function getinfostring(dbc::Handle, type)
    buf = Vector{UInt8}(undef, 6)
    len = Ref{SQLSMALLINT}()
    @checksuccess dbc @odbc(:SQLGetInfo,
        (Ptr{Cvoid}, SQLUSMALLINT, Ptr{UInt8}, SQLSMALLINT, Ref{SQLSMALLINT}),
        getptr(dbc), type, buf, sizeof(buf), len)
    return str(buf, len[])
end

const SQL_ALL_TYPES = 0

function SQLGetTypeInfo(stmt)
    @odbc(:SQLGetTypeInfo,
        (Ptr{Cvoid}, SQLSMALLINT),
        getptr(stmt), SQL_ALL_TYPES)
end

function gettypes(dbc)
    stmt = Handle(SQL_HANDLE_STMT, getptr(dbc))
    SQLGetTypeInfo(stmt)
    return stmt
end

const SQL_DRIVER_COMPLETE = UInt16(1)
const SQL_DRIVER_COMPLETE_REQUIRED = UInt16(3)
const SQL_DRIVER_NOPROMPT = UInt16(0)
const SQL_DRIVER_PROMPT = UInt16(2)

function SQLDriverConnect(dbc,connstr)
    # need to call SQLDriverConnectW to ensure our application gets labelled "unicode"
    # and subsequent library calls behave properly
    c = transcode(sqlwcharsize(), connstr)
    # push!(c, sqlwcharsize()(0))
    ret = @odbc(:SQLDriverConnectW,
        (Ptr{Cvoid},Ptr{Cvoid},Ptr{SQLWCHAR},SQLSMALLINT,Ptr{SQLCHAR},SQLSMALLINT,Ptr{SQLSMALLINT},SQLUSMALLINT),
        getptr(dbc),C_NULL,c,length(c),C_NULL,0,C_NULL,SQL_DRIVER_NOPROMPT)
    if ret == SQL_ERROR
        # @warn diagnostics(dbc)
        ret = @odbc(:SQLDriverConnect,
            (Ptr{Cvoid},Ptr{Cvoid},Ptr{SQLCHAR},SQLSMALLINT,Ptr{SQLCHAR},SQLSMALLINT,Ptr{SQLSMALLINT},SQLUSMALLINT),
            getptr(dbc),C_NULL,connstr,SQL_NTS,C_NULL,0,C_NULL,SQL_DRIVER_NOPROMPT)
    end
    return ret
end

function driverconnect(connstr)
    dbc = Handle(SQL_HANDLE_DBC, ODBC_ENV[])
    @checksuccess dbc SQLDriverConnect(dbc, connstr)
    return dbc
end

connect(dsn, extraauth) = driverconnect("$dsn;$extraauth")

function SQLDisconnect(dbc::Ptr{Cvoid})
    @odbc(:SQLDisconnect,
        (Ptr{Cvoid},),
        dbc)
end

function disconnect(h::Handle)
    if h.type == SQL_HANDLE_DBC
        @checksuccess(h, SQLDisconnect(h.ptr))
        finalize(h)
    end
end

function cwstring(s::AbstractString)
    bytes = codeunits(String(s))
    0 in bytes && throw(ArgumentError("embedded NULs are not allowed in input strings: $(repr(s))"))
    return transcode(sqlwcharsize(), bytes)
end

function SQLPrepare(stmt::Ptr{Cvoid},query::AbstractString)
    q = cwstring(query)
    @odbc(:SQLPrepareW,
        (Ptr{Cvoid},Ptr{SQLWCHAR},Int16),
        stmt,q,length(q))
end

function prepare(dbc::Handle, sql)
    stmt = Handle(SQL_HANDLE_STMT, getptr(dbc))
    enableasync(stmt)
    @checksuccess stmt SQLPrepare(getptr(stmt), sql)
    return stmt
end

function SQLNumParams(stmt::Ptr{Cvoid},param_count::Ref{SQLSMALLINT})
    @odbc(:SQLNumParams,
        (Ptr{Cvoid},Ptr{SQLSMALLINT}),
        stmt,param_count)
end

function numparams(stmt::Handle)
    out = Ref{SQLSMALLINT}()
    @checksuccess stmt SQLNumParams(getptr(stmt), out)
    return out[]
end

const SQL_PARAM_INPUT = Int16(1)
const SQL_PARAM_OUTPUT = Int16(4)
const SQL_PARAM_INPUT_OUTPUT = Int16(2)

function SQLBindParameter(stmt,x::Int,iotype::Int16,ctype::Int16,sqltype::Int16,column_size::Int,decimal_digits::Int,param_value,param_size::Int,len::Ptr{SQLLEN})
    @odbc(:SQLBindParameter,
        (Ptr{Cvoid},UInt16,Int16,Int16,Int16,UInt,Int16,Ptr{Cvoid},Int,Ptr{SQLLEN}),
        getptr(stmt),x,iotype,ctype,sqltype,column_size,decimal_digits,param_value,param_size,len)
end

bindparam(stmt,x,iotype,ctype,sqltype,column_size,decimal_digits,param_value,param_size,len) =
    @checksuccess stmt SQLBindParameter(stmt,x,iotype,ctype,sqltype,column_size,decimal_digits,param_value,param_size,len)

const SQL_CLOSE = UInt16(0)

function SQLFreeStmt(stmt::Ptr{Cvoid},param::UInt16)
    @odbc(:SQLFreeStmt,
        (Ptr{Cvoid},UInt16),
        stmt, param)
end

function freestmt(stmt)
    if stmt.ptr != C_NULL
        SQLFreeStmt(getptr(stmt), SQL_CLOSE)
    end
end

function SQLExecute(stmt::Ptr{Cvoid})
    @odbc(:SQLExecute,
        (Ptr{Cvoid},),
        stmt)
end

execute(stmt::Handle) = @checksuccess stmt SQLExecute(getptr(stmt))

function SQLExecDirect(stmt::Ptr{Cvoid},query::AbstractString)
    q = cwstring(query)
    @odbc(:SQLExecDirectW,
        (Ptr{Cvoid},Ptr{SQLWCHAR},Int),
        stmt,q,length(q))
end

function execdirect(stmt::Handle, sql)
    @checksuccess stmt SQLExecDirect(getptr(stmt), sql)
    return stmt
end

function SQLNumResultCols(stmt::Ptr{Cvoid},cols::Ref{Int16})
    @odbc(:SQLNumResultCols,
        (Ptr{Cvoid},Ref{Int16}),
        stmt, cols)
end

function numcols(stmt::Handle)
    out = Ref{SQLSMALLINT}()
    @checksuccess stmt SQLNumResultCols(getptr(stmt), out)
    return out[]
end

function SQLRowCount(stmt::Ptr{Cvoid},rows::Ref{SQLLEN})
    @odbc(:SQLRowCount,
        (Ptr{Cvoid},Ref{SQLLEN}),
        stmt, rows)
end

function numrows(stmt::Handle)
    out = Ref{SQLLEN}()
    @checksuccess stmt SQLRowCount(getptr(stmt), out)
    return out[]
end

function SQLDescribeCol(stmt,i,nm::Vector,len::Vector,dt::Vector,cs::Vector,dd::Vector,nul::Vector)
    @odbc(:SQLDescribeColW,
        (Ptr{Cvoid},SQLUSMALLINT,Ptr{SQLWCHAR},SQLSMALLINT,Ptr{SQLSMALLINT},Ptr{SQLSMALLINT},Ptr{SQLULEN},Ptr{SQLSMALLINT},Ptr{SQLSMALLINT}),
        stmt,i,nm,length(nm),pointer(len, i),pointer(dt, i),pointer(cs, i),pointer(dd, i),pointer(nul, i))
end

const SQL_ATTR_ROW_ARRAY_SIZE = 27
const SQL_IS_UINTEGER = -5
const SQL_IS_INTEGER = -6

const SQL_ATTR_ASYNC_ENABLE = 4
const SQL_ASYNC_ENABLE_ON = 1

function SQLSetStmtAttr(stmt::Ptr{Cvoid},attribute,value,value_length)
    @odbc(:SQLSetStmtAttrW,
        (Ptr{Cvoid},SQLINTEGER,SQLULEN,SQLINTEGER),
        stmt,attribute,value,value_length)
end

setrowset(stmt::Handle, rowset) = @checksuccess stmt SQLSetStmtAttr(getptr(stmt), SQL_ATTR_ROW_ARRAY_SIZE, rowset, SQL_IS_UINTEGER)
enableasync(stmt::Handle) = SQLSetStmtAttr(getptr(stmt), SQL_ATTR_ASYNC_ENABLE, SQL_ASYNC_ENABLE_ON, SQL_IS_UINTEGER)

const SQL_ATTR_ROWS_FETCHED_PTR  = 26

function SQLSetStmtAttr(stmt::Ptr{Cvoid},attribute,value::Ref{SQLLEN},value_length)
    @odbc(:SQLSetStmtAttrW,
        (Ptr{Cvoid},SQLINTEGER,Ref{SQLLEN},SQLINTEGER),
        stmt,attribute,value,value_length)
end

const SQL_NTS = -3

function setrowsfetched(stmt)
    rowsfetchedref = Ref{SQLLEN}(0)
    SQLSetStmtAttr(getptr(stmt), SQL_ATTR_ROWS_FETCHED_PTR, rowsfetchedref, SQL_NTS)
    return rowsfetchedref
end

const SQL_FETCH_NEXT = Int16(1)

function SQLFetchScroll(stmt::Ptr{Cvoid},fetch_orientation::Int16,fetch_offset::Int)
    @odbc(:SQLFetchScroll,
        (Ptr{Cvoid},Int16,Int),
        stmt,fetch_orientation,fetch_offset)
end

fetch(stmt::Handle) = @checksuccess stmt SQLFetchScroll(getptr(stmt), SQL_FETCH_NEXT, 0)

function SQLBindCol(stmt::Ptr{Cvoid},x,ctype,mem,jlsize,indicator::Vector{SQLLEN})
    @odbc(:SQLBindCol,
        (Ptr{Cvoid},SQLUSMALLINT,SQLSMALLINT,Ptr{Cvoid},SQLLEN,Ptr{SQLLEN}),
        stmt,x,ctype,mem,jlsize,indicator)
end

function SQLGetData(stmt::Ptr{Cvoid},i,ctype,mem,jlsize,indicator::Vector{SQLLEN})
    @odbc(:SQLGetData,
        (Ptr{Cvoid},SQLUSMALLINT,SQLSMALLINT,Ptr{Cvoid},SQLLEN,Ptr{SQLLEN}),
        stmt,i,ctype,mem,jlsize,indicator)
end

function SQLMoreResults(stmt::Ptr{Cvoid})
    @odbc(:SQLMoreResults,
        (Ptr{Cvoid},),
        stmt)
end

moreresults(stmt::Handle) = @checksuccess stmt SQLMoreResults(getptr(stmt))

function SQLGetDiagRec(handletype, handle, i, state, native, error_msg, msg_length)
    @odbc(:SQLGetDiagRecW,
        (SQLSMALLINT,Ptr{Cvoid},SQLSMALLINT,Ptr{SQLWCHAR},Ref{SQLINTEGER},Ptr{SQLWCHAR},SQLSMALLINT,Ref{SQLSMALLINT}),
        handletype,handle,i,state,native,error_msg,length(error_msg),msg_length)
end

function diagnostics(h::Handle)
    state = Vector{sqlwcharsize()}(undef, 6)
    native = Ref{SQLINTEGER}()
    error = Vector{sqlwcharsize()}(undef, 1024)
    len = Ref{SQLSMALLINT}()
    i = 1
    io = IOBuffer()
    while SQLGetDiagRec(gettype(h), getptr(h), i, state, native, error, len) == SQL_SUCCESS
        write(io, "$(str(state, 5)): $(str(error, len[]))")
        i += 1
    end
    return String(take!(io))
end

macro checkinst(expr)
    esc(quote
        ret = $expr
        if ret == 0
            error(installererror())
        end
        ret
    end)
end

function installererror()
    buf = Vector{UInt8}(undef, 512)
    code = Ref{UInt16}()
    len = Ref{UInt16}()
    i = 1
    err = ""
    while @odbcinst(:SQLInstallerError,
        (UInt16, Ref{UInt16}, Ptr{UInt8}, UInt16, Ref{UInt16}),
        i, code, buf, sizeof(buf), len) != SQL_NO_DATA
        err *= str(buf, len[])
        i += 1
    end
    return err
end

function getdrivers()
    name = Vector{UInt8}(undef, 1024)
    desc = Vector{UInt8}(undef, 1024)
    namelen = Ref{SQLSMALLINT}()
    desclen = Ref{SQLSMALLINT}()
    drivers = Dict{String, String}()
    while @checksuccess(ODBC_ENV[], @odbc(:SQLDrivers,
        (Ptr{Cvoid}, SQLUSMALLINT, Ptr{SQLCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}, Ptr{SQLCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}),
        getptr(ODBC_ENV[]), SQL_FETCH_NEXT, name, sizeof(name), namelen, desc, sizeof(desc), desclen)) == SQL_SUCCESS
        drivers[str(name, namelen[])] = str(desc, desclen[])
    end
    return drivers
end

function adddriver(name, path; kw...)
    ex = length(kw) > 0 ? "\0" * join((string(k, "=", v) for (k, v) in kw), '\0') : ""
    driver = "$name\0Driver=$path$ex\0\0"
    out = Vector{UInt8}(undef, 1024)
    ref = Ref{UInt16}()
    usage = Ref{UInt16}()
    ret = @checkinst @odbcinst(:SQLInstallDriverEx,
        (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, UInt16, Ref{UInt16}, UInt16, Ref{UInt16}),
        driver, C_NULL, out, length(out), ref, 2, usage)
    return ret
end

function removedriver(name, removedsns)
    usage = Ref{UInt16}()
    ret = @checkinst @odbcinst(:SQLRemoveDriver,
        (Ptr{UInt8}, Bool, Ref{UInt16}),
        name, removedsns, usage)
    return ret
end

function getdsns()
    name = Vector{UInt8}(undef, 1024)
    desc = Vector{UInt8}(undef, 1024)
    namelen = Ref{SQLSMALLINT}()
    desclen = Ref{SQLSMALLINT}()
    dsns = Dict{String, String}()
    while @checksuccess(ODBC_ENV[], @odbc(:SQLDataSources,
        (Ptr{Cvoid}, SQLUSMALLINT, Ptr{SQLCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}, Ptr{SQLCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}),
        getptr(ODBC_ENV[]), SQL_FETCH_NEXT, name, sizeof(name), namelen, desc, sizeof(desc), desclen)) == SQL_SUCCESS
        dsns[str(name, namelen[])] = str(desc, desclen[])
    end
    return dsns
end

function writeinientry(file, section, key, value)
    @checkinst @odbcinst(:SQLWritePrivateProfileString,
        (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}),
        section, string(key), string(value), file)
end

function adddsn(name, driver; kw...)
    ret = @checkinst @odbcinst(:SQLWriteDSNToIni,
        (Ptr{UInt8}, Ptr{UInt8}),
        name, driver)
    for (k, v) in kw
        writeinientry("odbc.ini", name, k, v)
    end
    return
end

function removedsn(name)
    ret = @checkinst @odbcinst(:SQLRemoveDSNFromIni,
        (Ptr{UInt8},),
        name)
    return
end

end # module API