text:
#Author: Joshua Ong
#Summary run robust with same number of ambualnce, different alpha values, different training sets to see when the model
#1) saturates 2) performs best

using DelimitedFiles
using AmbulanceDeployment

include("single_robust.jl")

alpha_list = [.1,.05,.01,.001,.0001]
cross_validation = [1,2,3]

scenarios = Dict{String, Dict{Int, Vector{Vector{Int}}}}()
generated_deployment = Dict{String, Dict{Int, Vector{Vector{Int}}}}()
upperbounds = Dict{String, Dict{Int, Vector{Float64}}}()
upptiming = Dict{String, Dict{Int, Vector{Float64}}}()
lowtiming = Dict{String, Dict{Int, Vector{Float64}}}()
amb_deployment = Dict{String, Dict{Int, Vector{Int}}}()

namb = 40

for j in cross_validation
    for i in alpha_list
        name = "robust" * string(i)

        model = generate_robust(i, 40, j) #alpha,namb,cross_validation

        scenarios[name] = Dict()
        generated_deployment[name] = Dict()
        upperbounds[name] = Dict()
        upptiming[name] = Dict()
        lowtiming[name] = Dict()
        scenarios[name][namb] = model.scenarios
        generated_deployment[name][namb] = model.deployment
        upperbounds[name][namb] = model.upperbounds
        upptiming[name][namb] = model.upptiming
        lowtiming[name][namb] = model.lowtiming

        filename = "cross" * string(j) * "robust" * string(i) * ".csv"
        print(filename)
        writedlm(filename,  deployment(model), ',')
    end
end

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
using SDE.Schauder
using Base.Test
srand(10)
## procudes the same as computing mu' and then performing pick_up(mu).
function pickedup_mu(y, L)
	n = 2^(L-1) #number of even element/elements in lowest level!
	N = length(y)
	mu = zeros(2n-1)

	dy = [diff(y),0] #h_t = f(y_t)(x_t+h - x_t)
	m = zeros(2n-1)
	fn = float(n)
	for t in 1:N
		for i in 1:n - 1
			fi = float(i)
			m[2i-1] = hat(y[t]*fn-fi + 1.)
			m[2i] = hat(y[t]*fn-fi +.5)
		end
		m[2n-1] = hat(y[t]*fn-fn+1)
		m = Schauder.pickup_mu!(m)
		mu = mu .+ m .* dy[t]
	end
	mu
end

L = 3
N = 100

n = 2^L -1
z = rand(N)
p = finger_pm(L,0)
pm = Schauder.permutationmatrix(p)

th1 = rand(n)*0.1 .+ linspace(0,1.,n)
th2 = copy(th1)
Schauder.pickup!(th2)

mu1 = Schauder.fe_mu(z,L, 0)
mu1b = mu1
if isdefined(:fe_mu_c)
	mu1b = Schauder.fe_mu_c(z,L, 0)
end
#mu2 = pickedup_mu(z, L)
mu2 = Schauder.pickup_mu!(Schauder.fe_mu(z,L, 0))
mu2b = Schauder.pickup_mu!(copy(mu1))

tic()
Sigma1 = Schauder.fe_Sigma_at(z,0.1, L)
toc()
tic()
Sigma1b = Schauder.fe_Sigma_dot(z,0.1, L)
toc()

Sigma2 = copy(Sigma1)
Sigma2 = Schauder.pickup_Sigma!(Sigma2)

println([th1 th2 mu1 mu1b mu2 mu2b])

@test norm(mu2-mu2b) < eps()
@test norm(mu1-mu1b) < eps()
@test norm(Sigma1-Sigma1b) < 4eps()

@test norm(th1'*mu1 - th2'*mu2) < 8eps()
@test norm(th1'*mu1 - th2'*mu2b) < 8eps()

@test norm(th1'*Sigma1*th1 - th1'*Sigma1b*th1) < 8eps()
@test norm(th1'*Sigma1*th1 - th2'*Sigma2*th2) < 8eps()

ext: jl


text:
using Distributions

# function ∑matrix(x::Vector, y::Vector)
#     z = x .* y
#     z .* z'
# end

function shapiro{T<:Float64}(x::Vector{T})
    s = sort(x)
    d = Truncated(Normal(mean(s),1),minimum(s), maximum(s))
    m = sort(rand(d, length(s)))

    V  = cov(m)
    a1 = m' * inv(V)
    a2 = sqrt(m' * inv(V) * inv(V) * m)
    a  = a1 ./ a2

    # ∑ ai * si and then ^ 2
    # so a cross product and then square it
    # a's ^ 2
    # num = (a's) ^ 2

    n   = a .* s
    nu  = sum(n)
    num = nu ^ 2

    # ∑ (xi - xbar)  ^ 2
    # so first get y = xi - xbar vector
    # y'y
    # y = s.-mean(s)
    # den = y'y

    d   = s .- mean(s)
    de  = d .^ 2
    den = sum(de)

    num/den
end


text:
export
	ConstraintSense,
	Inequality,
	Equality,
	Stage,
	State,
	Control,
	Coupled,
	Dynamical

export
	evaluate,
	jacobian

abstract type GeneralConstraint end

"Specifies whether the constraint is an equality or inequality constraint.
Valid subtypes are `Equality`, `Inequality`, and `Null`"
abstract type ConstraintSense end
"Inequality constraints of the form ``h(x) \\leq 0``"
abstract type Equality <: ConstraintSense end
"Equality constraints of the form ``g(x) = 0``"
abstract type Inequality <: ConstraintSense end
abstract type Null <: ConstraintSense end

"""
Specifies the ``bandedness'' of the constraint. This ends up being the width
of the constraint Jacobian, or the total number of input variables. This is
important to reduce the size of arrays needed to store the Jacobian, as well as
special-case the block matrix algebra.

Current subtypes:
* [`Stage`](@ref)
* [`State`](@ref) <: `Stage`
* [`Control`](@ref) <: `Stage`
* [`Coupled`](@ref)
* [`Dynamical`](@ref) <: `Coupled`
* [`CoupledState`](@ref) <: `Coupled`
* [`CoupledControl`](@ref) <: `Coupled`
* [`General`](@ref)
* [`GeneralState`](@ref) <: `General`
* [`GeneralControl`](@ref) <: `General`
"""
abstract type ConstraintType end
"Only a function of states and controls at a single knotpoint"
abstract type Stage <: ConstraintType end
"Only a function of states at a single knotpoint"
abstract type State <: Stage end
"Only a function of controls at a single knotpoint"
abstract type Control <: Stage end
"Only a function of states and controls at two adjacent knotpoints"
abstract type Coupled <: ConstraintType end
"Only a function of states and two adjacent knotpoints,
and the control at the previous knotpoint, i.e. f(x,u) - x′"
abstract type Dynamical <: Coupled end
"Only a function of states at adjacent knotpoints"
abstract type CoupledState <: Coupled end
"Only a function of controls at adjacent knotpoints"
abstract type CoupledControl <: Coupled end
"A function of all states and controls in the trajectory"
abstract type General <: ConstraintType end
"A function of all states in the trajectory"
abstract type GeneralState <: General end
"A function of all controls in the trajectory"
abstract type GeneralControl <: General end

""" $(TYPEDEF)
Abstract vector-valued constraint of size `P` for a trajectory optimization problem.
May be either inequality or equality (specified by `S<:ConstraintSense`), and be function of
single, adjacent, or all knotpoints (specified by `W<:ConstraintType`).

Interface:
Any constraint type must implement the following interface:
```julia
n = state_dim(::MyCon)
m = control_dim(::MyCon)
p = Base.length(::MyCon)
c = evaluate(::MyCon, args...)   # args determined by W
∇c = jacobian(::MyCon, args...)  # args determined by W
```

The `evaluate` and `jacobian` (identical signatures) methods should have the following signatures
* W <: State: `evaluate(::MyCon, x::SVector)`
* W <: Control: `evaluate(::MyCon, u::SVector)`
* W <: Stage: `evaluate(::MyCon, x, u)`
* W <: Dynamical: `evaluate(::MyCon, x′, x, u)`
* W <: Coupled: `evaluate(::MyCon, x′, u′ x, u)`

Or alternatively,
* W <: Stage: `evaluate(::MyCon, z::KnotPoint)`
* W <: Coupled: `evaluate(::MyCon, z′::KnotPoint, z::KnotPoint)`

The Jacobian method for [`State`](@ref) or [`Control`](@ref) is optional, since it will
	be automatically computed using ForwardDiff. Automatic differentiation
	for other types of constraints is not yet supported.

For W <: State, `control_dim(::MyCon)` doesn't need to be defined. Equivalently, for
	W <: Control, `state_dim(::MyCon)` doesn't need to be defined.

For W <: General, the more general `evaluate` and `jacobian` methods must be used
```julia
evaluate!(vals::Vector{<:AbstractVector}, ::MyCon, Z::Traj, inds=1:length(Z)-1)
jacobian!(∇c::Vector{<:AbstractMatrix}, ::MyCon, Z::Traj, inds=1:length(Z)-1)
```
These methods can be specified for any constraint, instead of the not-in-place functions
	above.
"""
abstract type AbstractConstraint{S<:ConstraintSense,W<:ConstraintType,P} <: GeneralConstraint end

# Getters
"Get type of constraint (bandedness)"
contype(::AbstractConstraint{S,W}) where {S,W} = W
"Get constraint sense (Inequality vs Equality)"
sense(::AbstractConstraint{S}) where S = S

"Returns the width of the constraint Jacobian, i.e. the total number of inputs
to the constraint"
width(con::AbstractConstraint{S,Stage}) where S = state_dim(con) + control_dim(con)
width(con::AbstractConstraint{S,State}) where S = state_dim(con)
width(con::AbstractConstraint{S,Control}) where S = control_dim(con)
width(con::AbstractConstraint{S,Coupled}) where S = 2*(state_dim(con) + control_dim(con))
width(con::AbstractConstraint{S,Dynamical}) where S = 2*state_dim(con) + control_dim(con)
width(con::AbstractConstraint{S,CoupledState}) where S = 2*state_dim(con)
width(con::AbstractConstraint{S,CoupledControl}) where S = 2*control_dim(con)
width(con::AbstractConstraint{S,<:General}) where S = Inf

width(::Type{Stage},N,M) = N+M
width(::Type{State},N,M) = N
width(::Type{Control},N,M) = M
width(::Type{Coupled},N,M) = 2N + 2M
width(::Type{Dynamical},N,M) = 2N + M
width(::Type{CoupledState},N,M) = 2N
width(::Type{CoupledControl},N,M) = 2M

"Upper bound of the constraint, as a vector, which is 0 for all constraints
(except bound constraints)"
upper_bound(con::AbstractConstraint{Inequality,W,P}) where {P,W} = @SVector zeros(P)
upper_bound(con::AbstractConstraint{Equality,W,P}) where {P,W} = @SVector zeros(P)

"Upper bound of the constraint, as a vector, which is 0 equality and -Inf for inequality
(except bound constraints)"
lower_bound(con::AbstractConstraint{Inequality,W,P}) where {P,W} = -Inf*@SVector ones(P)
lower_bound(con::AbstractConstraint{Equality,W,P}) where {P,W} = @SVector zeros(P)

"Is the constraint a bound constraint or not"
@inline is_bound(con::AbstractConstraint) = false

"Check whether the constraint is consistent with the specified state and control dimensions"
@inline check_dims(con::AbstractConstraint{S,State},n,m) where S = state_dim(con) == n
@inline check_dims(con::AbstractConstraint{S,Control},n,m) where S = control_dim(con) == m
@inline function check_dims(con::AbstractConstraint{S,W},n,m) where {S,W<:ConstraintType}
	state_dim(con) == n && control_dim(con) == m
end

"Size of state vector"
control_dims(::AbstractConstraint{S,State}) where S =
	throw(ErrorException("Cannot get control dimension from a state-only constraint"))

"Size of control vector"
state_dims(::AbstractConstraint{S,Control}) where S =
	throw(ErrorException("Cannot get state dimension from a control-only constraint"))

Base.length(::AbstractConstraint{S,W,P}) where {S,W,P} = P
Base.size(con::AbstractConstraint) = length(con), width(con)

con_label(::AbstractConstraint, i::Int) = "index $i"

############################################################################################
# 								EVALUATION METHODS 										   #
############################################################################################
"""```
evaluate!(vals::Vector{<:AbstractVector}, con::AbstractConstraint{S,W,P},
	Z, inds=1:length(Z)-1)
```
Evaluate constraints for entire trajectory. This is the most general method used to evaluate
	constraints, and should be the one used in other functions.

For W<:Stage this will loop over calls to `evaluate(con,Z[k])`

For W<:Coupled this will loop over calls to `evaluate(con,Z[k+1],Z[k])`

For W<:General,this must function must be explicitly defined. Other types may define it
	if desired.
"""
function evaluate!(vals::Vector{<:AbstractVector}, con::AbstractConstraint{P,<:Stage},
		Z::Traj, inds=1:length(Z)) where P
	for (i,k) in enumerate(inds)
		vals[i] = evaluate(con, Z[k])
	end
end

function evaluate!(vals::Vector{<:AbstractVector}, con::AbstractConstraint{P,<:Coupled},
		Z::Traj, inds=1:length(Z)-1) where P
	for (i,k) in enumerate(inds)
		vals[i] = evaluate(con, Z[k+1], Z[k])
	end
end

"""```
jacobian!(vals::Vector{<:AbstractVector}, con::AbstractConstraint{S,W,P},
	Z, inds=1:length(Z)-1)
```
Evaluate constraint Jacobians for entire trajectory. This is the most general method used to
	evaluate constraint Jacobians, and should be the one used in other functions.

For W<:Stage this will loop over calls to `jacobian(con,Z[k])`

For W<:Coupled this will loop over calls to `jacobian(con,Z[k+1],Z[k])`

For W<:General,this must function must be explicitly defined. Other types may define it
	if desired.
"""
function jacobian!(∇c::Vector{<:SizedMatrix}, con::AbstractConstraint{P,<:Stage},
		Z::Traj, inds=1:length(Z)) where P
	for (i,k) in enumerate(inds)
		jacobian!(∇c[i], con, Z[k])
	end
end

# Default methods for converting KnotPoints to states and controls
for method in [:evaluate]#, :jacobian]
	@eval begin
			@inline $(method)(con::AbstractConstraint{S,Stage},   Z::AbstractKnotPoint) where S = $(method)(con, state(Z), control(Z))
			@inline $(method)(con::AbstractConstraint{S,State},   Z::AbstractKnotPoint) where S = $(method)(con, state(Z))
			@inline $(method)(con::AbstractConstraint{S,Control}, Z::AbstractKnotPoint) where S = $(method)(con, control(Z))

			@inline $(method)(con::AbstractConstraint{S,Coupled}, Z′::AbstractKnotPoint, Z::AbstractKnotPoint) where S =
				$(method)(con, state(Z′), control(Z′), state(Z), control(Z))
			@inline $(method)(con::AbstractConstraint{S,Dynamical}, Z′::AbstractKnotPoint, Z::AbstractKnotPoint) where S =
				$(method)(con, state(Z′), state(Z), control(Z))
			@inline $(method)(con::AbstractConstraint{S,CoupledState}, Z′::AbstractKnotPoint, Z::AbstractKnotPoint) where S =
				$(method)(con, state(Z′), state(Z))
			@inline $(method)(con::AbstractConstraint{S,CoupledControl}, Z′::AbstractKnotPoint, Z::AbstractKnotPoint) where S =
				$(method)(con, control(Z′), control(Z))
	end
end

# Default method for converting KnotPoints to states and controls
@inline jacobian!(∇c, con::AbstractConstraint{P,State}, z::AbstractKnotPoint) where P =
	jacobian!(∇c, con, state(z))
@inline jacobian!(∇c, con::AbstractConstraint{P,Control}, z::AbstractKnotPoint) where P =
	jacobian!(∇c, con, control(z))

function jacobian!(∇c, con::AbstractConstraint{S,W}, x::SVector{N}) where {S,N,W<:Union{State,Control}}
	eval_c(x) = evaluate(con, x)
	∇c .= ForwardDiff.jacobian(eval_c, x)
end


text:
function ddp_solve!(prob::ProblemData;
    max_iter = 10,
    grad_tol = 1.0e-5,
    verbose = true,
	cache = false)

	println()
    (verbose && prob.m_data.obj isa StageCosts) && printstyled("Differential Dynamic Programming\n",
		color = :red, bold = true)

	# data
	p_data = prob.p_data
	m_data = prob.m_data
	s_data = prob.s_data

    # compute objective
    # s_data.obj = objective(m_data, mode = :nominal)
	objective!(s_data, m_data, mode = :nominal)
	# println("obj")
    for i = 1:max_iter
        # derivatives
        derivatives!(m_data)
		# println("derivatives")

		# backward pass
        backward_pass!(p_data, m_data)
		# println("backward pass")

        # forward pass
        forward_pass!(p_data, m_data, s_data)
		# println("forward pass")

		# cache solver data
		cache && cache!(s_data)

        # check convergence
        grad_norm = norm(s_data.gradient, Inf)
        verbose && println("     iter: $i
             cost: $(s_data.obj)
			 grad_norm: $(grad_norm)
			 c_max: $(s_data.c_max)
			 α: $(s_data.α)")
		grad_norm < grad_tol && break
        !s_data.status && break
    end
end

"""
    gradient of Lagrangian
        https://web.stanford.edu/class/ee363/lectures/lqr-lagrange.pdf
"""
function lagrangian_gradient!(s_data::SolverData, p_data::PolicyData, n, m, T)
    p = p_data.p
    Qx = p_data.Qx
    Qu = p_data.Qu

    for t = 1:T-1
        idx_x = (t == 1 ? 0 : (t - 1) * n[t-1]) .+ (1:n[t])
        s_data.gradient[idx_x] = Qx[t] - p[t]
        # NOTE: gradient wrt xT is satisfied implicitly

        idx_u = sum(n) + (t == 1 ? 0 : (t - 1) * m[t-1]) .+ (1:m[t])
        s_data.gradient[idx_u] = Qu[t]
    end
end

function lagrangian_gradient!(s_data::SolverData, p_data::PolicyData, m_data::ModelData)
	lagrangian_gradient!(s_data, p_data,
		m_data.n, m_data.m, m_data.T)
end

"""
    augmented Lagrangian solve
"""
function constrained_ddp_solve!(prob::ProblemData;
    max_iter = 10,
	max_al_iter = 5,
    grad_tol = 1.0e-5,
	con_tol = 1.0e-3,
	con_norm_type = Inf,
	ρ_init = 1.0,
	ρ_scale = 10.0,
	ρ_max = 1.0e8,
	cache = false,
    verbose = true)

	println()
	verbose && printstyled("Differential Dynamic Programming\n",
		color = :red, bold = true)

	# initial penalty
	prob.m_data.obj.ρ = ρ_init

	for i = 1:max_al_iter
		verbose && println("  al iter: $i")

		# primal minimization
		ddp_solve!(prob,
		    max_iter = max_iter,
		    grad_tol = grad_tol,
			cache = cache,
		    verbose = verbose)

		# update trajectories
		objective!(prob.s_data, prob.m_data, mode = :nominal)

		# constraint violation
		# c_max = constraint_violation(prob.m_data.obj.cons,
		# 	prob.m_data.x̄, prob.m_data.ū,
		# 	norm_type = con_norm_type)
		# verbose && println("    c_max: $(prob.s_data.c_max)\n")
		prob.s_data.c_max <= con_tol && break

		# dual ascent
		augmented_lagrangian_update!(prob.m_data.obj,
			s = ρ_scale, max_penalty = ρ_max)
	end
end


text:
abstract type DeviceEMIS end

abstract type Project{T <: BuildPhase} <: DeviceEMIS end    #Projects, parameterized by buildphase
abstract type GeneratorEMIS{T <: BuildPhase} <: Project{T} end  #Generators - subtype of Projects
abstract type StorageEMIS{T <: BuildPhase} <: Project{T} end  #Storage - subtype of Projects

get_name(project::P) where P <: Project{<: BuildPhase} = project.name
get_tech(project::P) where P <: Project{<: BuildPhase}  = project.tech
get_decision_year(project::P) where P <: Project{<: BuildPhase}  = project.decision_year
get_construction_year(project::P) where P <: Project{<: BuildPhase}  = project.construction_year
get_retirement_year(project::P) where P <: Project{<: BuildPhase}  = project.retirement_year
get_end_life_year(project::P) where P <: Project{<: BuildPhase}  = project.end_life_year
get_finance_data(project::P) where P <: Project{<: BuildPhase}  = project.finance_data
get_products(project::P) where P <: Project{<: BuildPhase} = project.products
get_mincap(project::P) where P <: GeneratorEMIS{<: BuildPhase}  = get_active_power_limits(get_tech(project))[:min]
get_maxcap(project::P) where P <: GeneratorEMIS{<: BuildPhase}  = get_active_power_limits(get_tech(project))[:max]
get_mincap(project::P) where P <: StorageEMIS{<: BuildPhase}  = get_output_active_power_limits(get_tech(project))[:min]
get_maxcap(project::P) where P <: StorageEMIS{<: BuildPhase}  = get_output_active_power_limits(get_tech(project))[:max]

function set_name!(project::P,
                   name::String) where P <: Project{<: BuildPhase}
        project.name = name
    return
end

function set_decision_year!(project::P, year::Int64) where P <: Project{<: BuildPhase}
    project.decision_year = year
    return
end

function set_construction_year!(project::P, year::Int64) where P <: Project{<: BuildPhase}
    project.construction_year = year
    return
end

function set_retirement_year!(project::P, year) where P <: Project{<: BuildPhase}
    project.retirement_year = year
end

function set_end_life_year!(project::P, year::Int64) where P <: Project{<: BuildPhase}
    project.end_life_year = year
    return
end

function set_investment_cost!(project::P,
                            investment_cost::Vector{Float64}) where P <: Project{<: BuildPhase}
    return
end

function set_investment_cost!(project::P,
                             investment_cost::Vector{Float64}) where P <: Project{Option}
   project.finance_data.investment_cost = investment_cost
end

function set_effective_investment_cost!(project::P,
                                      investment_cost::Float64) where P <: Project{<: BuildPhase}
    project.finance_data.effective_investment_cost = investment_cost
end


text:
using LinearAlgebra, PDMats
using PDMats: chol_lower, chol_upper

@testset "chol_lower" begin
    A = rand(100, 100)
    C = A'A
    size_of_one_copy = sizeof(C)
    @assert size_of_one_copy > 100  # ensure the matrix is large enough that few-byte allocations don't matter

    chol_lower(C)
    @test (@allocated chol_lower(C)) < 1.05 * size_of_one_copy  # allow 5% overhead

    for uplo in (:L, :U)
        ch = cholesky(Symmetric(C, uplo))
        chol_lower(ch)
        @test (@allocated chol_lower(ch)) < 50  # allow small overhead for wrapper types
        chol_upper(ch)
        @test (@allocated chol_upper(ch)) < 50  # allow small overhead for wrapper types
    end
end


text:
# Artifact generator for the assets directory.
#
# Requires:
#
#   - A github token stored in `ENV["GITHUB_TOKEN"]` with access to `repo` read/write.
#   - Extra packages: `ghr_jll` for release uploading.
#
# Usage:
#
#   - Update files in `/assets` directory.
#   - Bump the `version` variable below.
#   - Run `rebuild_artifacts()` to create new tarball and update Artifacts.toml.
#   - Commit & push the changes.
#   - Run `release_artifacts()` to create a new release and upload the new tarball.
# Note:
# https://discourse.julialang.org/t/git-tree-sha1-for-artifacts-toml/43879
# can make this process on windows a bit annoying
# TODO figure out if we can use ArtifactUtils or so...

using Pkg.Artifacts
using ghr_jll
using LibGit2

version = v"0.1.2"
user = "JuliaPlots"
repo = "AbstractPlotting.jl"
host = "https://github.com/$user/$repo/releases/download"

build_path = joinpath(@__DIR__, "build")
assets_path = joinpath(@__DIR__, "assets")
artifact_toml = joinpath(@__DIR__, "Artifacts.toml")

function rebuild_artifacts()
    ispath(build_path) && rm(build_path, force=true, recursive=true)
    mkpath(build_path)

    product_hash = create_artifact() do artifact_dir
        cp(assets_path, artifact_dir; force = true)
    end

    archive_filename = "assets-$version.tar.gz"
    path = joinpath(build_path, archive_filename)
    download_hash = archive_artifact(product_hash, path)
    bind_artifact!(
        artifact_toml,
        "assets",
        product_hash,
        force = true,
        download_info = Tuple[
            (
                "$host/assets-$version/$archive_filename",
                download_hash,
            ),
        ],
    )
    return path
end

function release_artifacts(tar_path)
    name = "Release assets $(version)"
    commit = string(LibGit2.GitHash(LibGit2.GitCommit(LibGit2.GitRepo(@__DIR__), "HEAD")))
    token = ENV["GITHUB_TOKEN"]
    tag = "assets-$version"
    ghr() do bin
        run(`$bin -u $user -r $repo -n $name -c $commit -t $token $tag $tar_path`)
    end
end

tar_path = rebuild_artifacts()

release_artifacts(tar_path)


text:
# This file was generated, do not modify it. # hide
lm = LinearRegressor()
mlm = machine(lm, select(X, :Horsepower), y)
fit!(mlm, rows=train)
rms(predict(mlm, rows=test), y[test])^2


text:
using BenchmarkTools
using Pkg
tempdir = mktempdir()
Pkg.activate(tempdir)
Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))
Pkg.add(["BenchmarkTools", "PkgBenchmark","Random"])
Pkg.resolve()

using Random
const SUITE = BenchmarkGroup()
Random.seed!(3)

SUITE = BenchmarkGroup()
SUITE["nodc"] = BenchmarkGroup(
    ["nodc"])
SUITE["dc"] = BenchmarkGroup(["dc"])
# designmatrix generation
SUITE["dc"]["X_gen_lin"] = @benchmarkable 1==1
SUITE["nodc"]["fit_lin"] = @benchmarkable 1==2


text:
#push!(LOAD_PATH,"/home/heltonmc/git-hub/LightPropagation/src")
#push!(LOAD_PATH,"/Users/michaelhelton/LightPropagation/src")

using LightPropagation

import Pkg; Pkg.add("Documenter")
using Documenter

makedocs(
    modules = [LightPropagation],
    sitename = "LightPropagation.jl",
    format = Documenter.HTML(),
    pages = [
                "Home" => "index.md",
                "Getting Started" => "getting-started.md",
                "Tutorials" => Any[
                    "When is semi-infinite semi-infinite?" => "whensemiinf.md"
                ],
                "Forward Models" => Any[
                    "Diffusion Approximation" => Any[
                    "Slab" => "DA_slab_semiinfgeom.md"
                     ],
                ],
    ]
)

deploydocs(
    repo = "github.com/heltonmc/LightPropagation.git",
    devbranch = "main"
)

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using Test

using Bijectors
using Bijectors: LeakyReLU

using LinearAlgebra
using ForwardDiff

true_logabsdetjac(b::Bijector{0}, x::Real) = (log ∘ abs)(ForwardDiff.derivative(b, x))
true_logabsdetjac(b::Bijector{0}, x::AbstractVector) = (log ∘ abs).(ForwardDiff.derivative.(b, x))
true_logabsdetjac(b::Bijector{1}, x::AbstractVector) = logabsdet(ForwardDiff.jacobian(b, x))[1]
true_logabsdetjac(b::Bijector{1}, xs::AbstractMatrix) = mapreduce(z -> true_logabsdetjac(b, z), vcat, eachcol(xs))

@testset "0-dim parameter, 0-dim input" begin
    b = LeakyReLU(0.1; dim=Val(0))
    x = 1.
    @test inv(b)(b(x)) == x
    @test inv(b)(b(-x)) == -x

    # Mixing of types
    # 1. Changes in input-type
    @assert eltype(b(Float32(1.))) == Float64
    @assert eltype(b(Float64(1.))) == Float64

    # 2. Changes in parameter-type
    b = LeakyReLU(Float32(0.1); dim=Val(0))
    @assert eltype(b(Float32(1.))) == Float32
    @assert eltype(b(Float64(1.))) == Float64

    # logabsdetjac
    @test logabsdetjac(b, x) == true_logabsdetjac(b, x)
    @test logabsdetjac(b, Float32(x)) == true_logabsdetjac(b, x)

    # Batch
    xs = randn(10)
    @test logabsdetjac(b, xs) == true_logabsdetjac(b, xs)
    @test logabsdetjac(b, Float32.(x)) == true_logabsdetjac(b, Float32.(x))

    @test logabsdetjac(b, -xs) == true_logabsdetjac(b, -xs)
    @test logabsdetjac(b, -Float32.(xs)) == true_logabsdetjac(b, -Float32.(xs))

    # Forward
    f = forward(b, xs)
    @test f.logabsdetjac ≈ logabsdetjac(b, xs)
    @test f.rv ≈ b(xs)

    f = forward(b, Float32.(xs))
    @test f.logabsdetjac == logabsdetjac(b, Float32.(xs))
    @test f.rv ≈ b(Float32.(xs))
end

@testset "0-dim parameter, 1-dim input" begin
    d = 2

    b = LeakyReLU(0.1; dim=Val(1))
    x = ones(d)
    @test inv(b)(b(x)) == x
    @test inv(b)(b(-x)) == -x

    # Batch
    xs = randn(d, 10)
    @test logabsdetjac(b, xs) == true_logabsdetjac(b, xs)
    @test logabsdetjac(b, Float32.(x)) == true_logabsdetjac(b, Float32.(x))

    @test logabsdetjac(b, -xs) == true_logabsdetjac(b, -xs)
    @test logabsdetjac(b, -Float32.(xs)) == true_logabsdetjac(b, -Float32.(xs))

    # Forward
    f = forward(b, xs)
    @test f.logabsdetjac ≈ logabsdetjac(b, xs)
    @test f.rv ≈ b(xs)

    f = forward(b, Float32.(xs))
    @test f.logabsdetjac == logabsdetjac(b, Float32.(xs))
    @test f.rv ≈ b(Float32.(xs))

    # Mixing of types
    # 1. Changes in input-type
    @assert eltype(b(ones(Float32, 2))) == Float64
    @assert eltype(b(ones(Float64, 2))) == Float64

    # 2. Changes in parameter-type
    b = LeakyReLU(Float32(0.1); dim=Val(1))
    @assert eltype(b(ones(Float32, 2))) == Float32
    @assert eltype(b(ones(Float64, 2))) == Float64
end


text:
using ITensors, Test
using Combinatorics: permutations

@testset "Combiner" begin

i = Index(2,"i")
j = Index(3,"j")
k = Index(4,"k")
l = Index(5,"l")

A = randomITensor(i, j, k, l)

@testset "Two index combiner" begin
    for inds_ij ∈ permutations([i,j])
        C,c = combiner(inds_ij...)
        B = A*C
        @test hasinds(B, l, k, c)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_il ∈ permutations([i,l])
        C,c = combiner(inds_il...)
        B = A*C
        @test hasinds(B, j, k)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_ik ∈ permutations([i,k])
        C,c = combiner(inds_ik...)
        B = A*C
        @test hasinds(B, j, l)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_jk ∈ permutations([j,k])
        C,c = combiner(inds_jk...)
        B = A*C
        @test hasinds(B, i, l)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        B = C*A
        @test hasinds(B, i, l)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_jl ∈ permutations([j,l])
        C,c = combiner(inds_jl...)
        B = A*C
        @test hasinds(B, i, k)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        B = C*A
        @test hasinds(B, i, k)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_kl ∈ permutations([k,l])
        C,c = combiner(inds_kl...)
        B = A*C
        @test hasinds(B, i, j)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        B = C*A
        @test hasinds(B, i, j)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
end

@testset "Three index combiner" begin
    for inds_ijl ∈ permutations([i,j,l])
        C,c = combiner(inds_ijl...)
        B = A*C
        @test hasindex(B, k)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        B = C*A
        @test hasindex(B, k)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_ijk ∈ permutations([i,j,k])
        C,c = combiner(inds_ijk...)
        B = A*C
        @test hasindex(B, l)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        B = C*A
        @test hasindex(B, l)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
    for inds_jkl ∈ permutations([j,k,l])
        C,c = combiner(inds_jkl...)
        B = A*C
        @test hasindex(B, i)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        B = C*A
        @test hasindex(B, i)
        @test c == commonindex(B, C)
        D = B*C
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
        D = C*B
        @test hasinds(D, i, j, k, l)
        @test D ≈ A
    end
end

@testset "SVD/Combiner should play nice" begin
    cmb, ci = combiner(i, j, k)
    Ac = A*cmb
    U,S,V,spec,u,v = svd(Ac, ci)
    Uc = cmb*U
    Ua,Sa,Va,spec,ua,va = svd(A, i, j, k)
    replaceindex!(Ua, ua, u)
    @test A ≈ cmb*Ac
    @test A ≈ Ac*cmb
    @test Ua*cmb ≈ U
    @test cmb*Ua ≈ U
    @test Ua ≈ Uc
    @test Uc*S*V ≈ A
    @test (cmb*Ua)*S*V ≈ Ac
    cmb, ci = combiner(i, j)
    Ac = A*cmb
    U,S,V,spec,u,v = svd(Ac, ci)
    Uc = U*cmb
    Ua,Sa,Va,spec,ua,va = svd(A, i, j)
    replaceindex!(Ua, ua, u)
    @test Ua ≈ Uc
    @test Ua*cmb ≈ U
    @test cmb*Ua ≈ U
    @test Uc*S*V ≈ A
    @test (cmb*Ua)*S*V ≈ Ac
end

@testset "mult/Combiner should play nice" begin
    cmb, ci = combiner(i, j, k)
    Ac = A*cmb
    B = randomITensor(l)
    C = Ac*B
    @test C*cmb ≈ A*B
end

@testset "Replace index combiner" begin
    C,nl = combiner(l, tags="nl")
    B = A*C
    replaceindex!(B, nl, l)
    @test B == A
end

end


text:
@testset "Params" begin
    @testset "Constructor" begin
        params = WebAPI.Params(:a => "1", :b => "2", :c => "3")
        dict = Dict(:a => "1", :b => "2", :c => "3")

        @test params.a == "1"
        @test params.b == "2"
        @test params.c == "3"

        @test params[:a] == "1"
        @test params[:b] == "2"
        @test params[:c] == "3"

        @test params["a"] == "1"
        @test params["b"] == "2"
        @test params["c"] == "3"

        @test Dict(params) == dict
        @test collect(keys(params)) == collect(keys(dict))
        @test collect(values(params)) == collect(values(dict))

        @test get(params, :a, "test") == "1"
        @test get(params, :b, "test") == "2"
        @test get(params, "c", "test") == "3"
        @test get(params, "d", "test") == "test"

        @test get(() -> "test", params, :a) == "1"
        @test get(() -> "test", params, :b) == "2"
        @test get(() -> "test", params, "c") == "3"
        @test get(() -> "test", params, "d") == "test"

        @test haskey(params, :a) == true
        @test haskey(params, "b") == true
        @test haskey(params, :d) == false

        @test propertynames(params) == [:a, :b, :c]
    end

    @testset "Crete Params" begin
        req = HTTP.Request("GET", "/params/1/2/3")
        path = "params/:a/:b/:c"
        params = WebAPI.createparams(req, path)
        dict = Dict(:a => "1", :b => "2", :c => "3")

        @test params.a == "1"
        @test params.b == "2"
        @test params.c == "3"

        @test params[:a] == "1"
        @test params[:b] == "2"
        @test params[:c] == "3"

        @test params["a"] == "1"
        @test params["b"] == "2"
        @test params["c"] == "3"

        @test Dict(params) == dict
        @test collect(keys(params)) == collect(keys(dict))
        @test collect(values(params)) == collect(values(dict))

        @test get(params, :a, "test") == "1"
        @test get(params, :b, "test") == "2"
        @test get(params, "c", "test") == "3"
        @test get(params, "d", "test") == "test"

        @test get(() -> "test", params, :a) == "1"
        @test get(() -> "test", params, :b) == "2"
        @test get(() -> "test", params, "c") == "3"
        @test get(() -> "test", params, "d") == "test"

        @test haskey(params, :a) == true
        @test haskey(params, "b") == true
        @test haskey(params, :d) == false

        @test propertynames(params) == [:a, :b, :c]
    end
end

@testset "Query" begin
    @testset "Constructor" begin
        query = WebAPI.Query(:a => "1", :b => "2", :c => "3")
        dict = Dict(:a => "1", :b => "2", :c => "3")

        @test query.a == "1"
        @test query.b == "2"
        @test query.c == "3"

        @test query[:a] == "1"
        @test query[:b] == "2"
        @test query[:c] == "3"

        @test query["a"] == "1"
        @test query["b"] == "2"
        @test query["c"] == "3"

        @test Dict(query) == dict
        @test collect(keys(query)) == collect(keys(dict))
        @test collect(values(query)) == collect(values(dict))

        @test get(query, :a, "test") == "1"
        @test get(query, :b, "test") == "2"
        @test get(query, "c", "test") == "3"
        @test get(query, "d", "test") == "test"

        @test get(() -> "test", query, :a) == "1"
        @test get(() -> "test", query, :b) == "2"
        @test get(() -> "test", query, "c") == "3"
        @test get(() -> "test", query, "d") == "test"

        @test haskey(query, :a) == true
        @test haskey(query, "b") == true
        @test haskey(query, :d) == false

        @test propertynames(query) == [:a, :b, :c]
    end

    @testset "Crete Query" begin
        req = HTTP.Request("GET", "/query?a=1&b=2&c=3")
        query = WebAPI.createquery(req)
        dict = Dict(:a => "1", :b => "2", :c => "3")

        @test query.a == "1"
        @test query.b == "2"
        @test query.c == "3"

        @test query[:a] == "1"
        @test query[:b] == "2"
        @test query[:c] == "3"

        @test query["a"] == "1"
        @test query["b"] == "2"
        @test query["c"] == "3"

        @test Dict(query) == dict
        @test collect(keys(query)) == collect(keys(dict))
        @test collect(values(query)) == collect(values(dict))

        @test get(query, :a, "test") == "1"
        @test get(query, :b, "test") == "2"
        @test get(query, "c", "test") == "3"
        @test get(query, "d", "test") == "test"

        @test get(() -> "test", query, :a) == "1"
        @test get(() -> "test", query, :b) == "2"
        @test get(() -> "test", query, "c") == "3"
        @test get(() -> "test", query, "d") == "test"

        @test haskey(query, :a) == true
        @test haskey(query, "b") == true
        @test haskey(query, :d) == false

        @test propertynames(query) == [:a, :b, :c]
    end
end


text:
export SingleId, has_index

## ----------  SingleId

# Object without index
SingleId(name :: Symbol) = SingleId(name, Array{Int,1}(), "");

SingleId(name :: Symbol, idx :: T1, descr = "") where T1 <: Integer =
    SingleId(name, [idx], descr);

SingleId(name :: Symbol, idxM :: Array{I1}) where I1 <: Integer =
    SingleId(name, idxM, "");

SingleId(name :: Symbol, descr :: String) = SingleId(name, Array{Int,1}(), descr);

# Make a string of the form "x[2, 1]"
function show_string(s :: SingleId)
    outStr = string(s.name);
    if !isempty(s.index)
        outStr = outStr * "$(s.index)";
    end
    return outStr
end

show(io :: IO,  s :: SingleId) =
    print(io,  "SingleId:  $(show_string(s))");

description(s :: SingleId) = s.description;
name(s :: SingleId) = s.name;
index(s :: SingleId) = s.index;

function has_index(this :: SingleId)
    return !Base.isempty(this.index)
end

function isequal(id1 :: SingleId, id2 :: SingleId)
    return (id1.name == id2.name)  &&  (id1.index == id2.index)
end

isequal(id1V :: Vector{SingleId},  id2V :: Vector{SingleId}) =
    all(isequal.(id1V, id2V));

#     outVal = length(id1V) == length(id2V);
#     if outVal
#         for i1 = 1 : length(id1V)
#             outVal = outVal && isequal(id1V[i1], id2V[i1]);
#         end
#     end
#     return outVal
# end

"""
	$(SIGNATURES)

Make a string from o `SingleId`. Such as "x[2, 1]".
"""
function make_string(id :: SingleId)
    if !has_index(id)
        outStr = "$(id.name)"
    elseif length(id.index) == 1
        outStr = "$(id.name)$(id.index)"
    else
        outStr = "$(id.name)$(id.index)"
    end
    return outStr
end

"""
	$(SIGNATURES)

The inverse of [`make_string`](@ref).
"""
function make_single_id(s :: T1) where T1 <: AbstractString
    if occursin('[', s)
        # Pattern "id1[4, 3]"
        m = match(r"(.+)\[([0-9, ]+)+\]", s);
        idxV = parse.(Int, split(m[2], ","));
        sId = SingleId(Symbol(m[1]),  idxV);
    else
        sId = SingleId(Symbol(s));
    end
    return sId
end

# ----------------


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger

function compute()::Int
    result = 0
    for i ∈ 1:9
        j, n = 1, BigInteger(i)
        while length(n) == j
            result += 1
            j += 1
            n *= i
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
include(raw"utils.jl")
include(raw"invariant_sets.jl")
include(raw"robust_invariant_sets.jl")

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
using Test

@time begin
@testset "PotentialUQ.jl" begin
    include("snap_test_MAP.jl")
    include("snap_test_sample.jl")
end
end


text:
include("d:\\source\\WindFarmCode\\Julia\\FlowFarm.jl\\src\\turbines.jl")
include("d:\\source\\WindFarmCode\\Julia\\FlowFarm.jl\\src\\wake_models.jl")
using Plots

model = Jensen(0.1)

coord = Coord(0.0,0.0,0.0)
rotor_diameter = 150.0
hub_height = 90.0
aI = 1.0/3.0
gamma = 0.0
yaw = 0.0
ct = 0.7 #TODO handle ct and axial induction appropriately

turbine = Turbine(coord, rotor_diameter, hub_height, aI, gamma, yaw, ct)
println(turbine)

num = 10
xlocs = ones(num).*900.0
ylocs = range(-300,stop=300,length=num)
zlocs = ones(num).*turbine.hub_height
locs = zeros(length(xlocs),3)
locs[:,1] = xlocs[:]
locs[:,2] = ylocs[:]
locs[:,3] = zlocs[:]
loss_val = zeros(num)

deflection = [100.0,0.0]
for i = 1:num
    loss_val[i] = wake_model(locs[i,:], deflection, model, turbine)
end

println(loss_val)
# plot(ylocs,loss_val)

model = Multizone([-0.5,0.3,1.0],0.05,[0.5,1.0,5.5],12.0,1.3)
for i = 1:num
    loss_val[i] = wake_model(locs[i,:], deflection, model, turbine)
end
println(loss_val)
plot(ylocs,loss_val)

# test Gauss model
# set params
k_star = 0.075
turbulence_intensity = 0.07
horizontal_spread_rate = k_star
vertical_spread_rate = k_star
alpha_star = 2.32
beta_star = 0.154

# test Gauss model 2014
version = 2014

model = Gauss(version, k_star, turbulence_intensity, horizontal_spread_rate, vertical_spread_rate, alpha_star, beta_star)
for i = 1:num
    loss_val[i] = wake_model(locs[i,:], deflection, model, turbine)
end
println(loss_val)
plot(ylocs,loss_val, show=true)

# test Gauss model 2016
version = 2016

model = Gauss(version, k_star, turbulence_intensity, horizontal_spread_rate, vertical_spread_rate, alpha_star, beta_star)
for i = 1:num
    loss_val[i] = wake_model(locs[i,:], deflection, model, turbine)
end
println(loss_val)
plot(ylocs,loss_val, show=true)


text:
using .BaseAgent
using Params
using Traffic.Simple
using Channel.Gilbert
using Util

export OptHighestSNR

type OptHighestSNR <: Agent
    s :: AgentState
    status :: Status
    OptHighestSNR(i, P, pos) = new(AgentState(i, P, pos), Initialized)
end

#=
I written my own minumum that returns Inf instead of
throwing an error when there is no suitable channel to
prevent unnecessary exception checking and performance
drop.
=#
function minimum_noise(env :: Environment)
    min_noise = Inf
    for i=1:endof(env.channels)
        if env.traffics[i].traffic
            continue
        end
        min_noise = min(min_noise, env.channels[i].noise)
    end
    min_noise
end

function BaseAgent.act(a :: OptHighestSNR, env, t)
    if a.status == Initialized
        if a.s.B_max - a.s.B_empty == 0
            return idle(a)
        end

        # we've done packet checking
        min_noise = minimum_noise(env)
        if min_noise == Inf
            return idle(a)
        end

        min_noise_channels = filter(i -> ! env.traffics[i].traffic && env.channels[i].noise <= min_noise, 1:Params.n_channel)
        chan = int8(min_noise_channels[rand(1:endof(min_noise_channels))])

        return switch(a, chan)
    elseif a.status == Switched
        return sense(a, env, detect_traffic)
    elseif a.status == Sensed
        if detect_traffic(a.s, env.traffics[a.s.chan], a.s.t_remaining)
            return idle(a)
        else
            P_tx = Params.P_levels[end-1]
            pkgs_to_send = min(a.s.B_max - a.s.B_empty, floor((a.s.t_remaining * capacity(env.channels[a.s.chan], P_tx, a.s.pos)) / Params.pkt_size))

            if pkgs_to_send == 0
                return idle(a)
            end
            return transmit!(a, P_tx, env, pkgs_to_send)
        end
    elseif a.status == Transmitted
        return idle(a)
    end
end


text:
module dcraw
"""
Simply provide dcraw 9.26 executable file to julia as `dcraw.DCRAW_EXE` for windows, osx and linux.

The [dcraw website](http://www.cybercom.net/~dcoffin/dcraw/) notes: "Unless otherwise noted in the source code, these programs
are free for all uses, although I would like to receive credit for them."
The program dcraw is written by Dave Coffin.
"""

const RAW_EXT = String[".3fr", ".ari", ".arw", ".bay", ".crw", ".cr2",
    ".cap", ".dcs", ".dcr", ".dng",
    ".drf", ".eip", ".erf", ".fff", ".iiq", ".k25", ".kdc", ".mdc", ".mef", ".mos", ".mrw",
    ".nef", ".nrw", ".obm", ".orf", ".pef", ".ptx", ".pxn", ".r3d", ".raf", ".raw", ".rwl",
    ".rw2", ".rwz", ".sr2", ".srf", ".srw", ".tif", ".x3f"]

const DCRAW_DIR = Pkg.dir("dcraw", "lib")
@static if is_windows()
    const DCRAW_EXE = joinpath(DCRAW_DIR, "dcraw-9.26-ms-64-bit.exe")
end

# On mac, compile dcraw.c with
# 	`llvm-gcc -o dcraw dcraw.c -lm -DNO_JPEG -DNO_LCMS -DNO_JASPER`
# source [http://vkphotoblog.blogspot.be/2014/05/dcraw-921-for-os-x-mavericks-users.html]

@static if is_unix()
    const DCRAW_EXE = joinpath(DCRAW_DIR, "dcraw")
    isexecutable(DCRAW_EXE) || chmod(DCRAW_EXE, 0o755)
end

# From LeafAreaIndex.jl package:
# We use [dcraw](http://www.cybercom.net/~dcoffin/dcraw/) to extract from RAW images
# the blue channel pixels and convert to 16bit pgm without exposure manipulation.

# ### dcraw options
# * -W: avoids stretching ("No matter how dark an image is, dcraw's auto-exposure stretches it so that one percent of its pixels appear white. The "-W" option avoids this behavior.")
# * -4: is for linear 16-bit, same as -6 -W -g 1 1 (with -g for gamma correction)
# * -j: don't stretch or rotate raw pixels
# * -t [0-7]: flip image (0=none)
# * -D: document mode without scaling (totally raw), while -d scales to 16bit (eg from 14bit). Either use -D and then reinterpret in julia or use -d. We now use -d otherwise we need to extract image property for bit depth (12bit, 14bit, ...).
# * -r 0 0 0 1: select only the blue channel (see Brusa and Bunker, 2014). This option selects from a bayer RGGB layout.
# * -v: for verbose output

end # module


text:
using AWSAuth
using Documenter
using HTTP

makedocs(modules=[AWSAuth],
         sitename="AWSAuth.jl",
         authors="The JuliaCloud Developers",
         pages=["index.md"])

deploydocs(repo="github.com/JuliaCloud/AWSAuth.jl.git", target="build")


text:
@testset "Activation Functions" begin

xs = rand(5,5)

@test all(sum(softmax(xs), 1) .≈ 1)

@test sum(softmax(vec(xs))) ≈ 1

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Patch is provided to give a concrete name and type to the Kubernetes PATCH request body.

    IoK8sApimachineryPkgApisMetaV1Patch(;
    )

"""
mutable struct IoK8sApimachineryPkgApisMetaV1Patch <: SwaggerModel

    function IoK8sApimachineryPkgApisMetaV1Patch(;)
        o = new()
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1Patch

const _property_map_IoK8sApimachineryPkgApisMetaV1Patch = Dict{Symbol,Symbol}()
const _property_types_IoK8sApimachineryPkgApisMetaV1Patch = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1Patch }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1Patch))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1Patch }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1Patch[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1Patch }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1Patch[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1Patch)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1Patch }, name::Symbol, val)
end


text:
module nparticles

export NParticleBasis, BosonicNParticleBasis, FermionicNParticleBasis, nparticleoperator_1, nparticleoperator_2

import Base.==

using ..bases, ..operators, ..operators_sparse

function distribute_bosons(particlenumber::Int, singleparticledimension::Int, index::Int=1, occupations::Vector{Int}=zeros(Int,singleparticledimension), results::Vector{Vector{Int}}=Vector{Int}[])
    if index==singleparticledimension
        occupations[index] = particlenumber
        push!(results, copy(occupations))
    else
        for n=particlenumber:-1:0
            occupations[index] = n
            distribute_bosons(particlenumber-n, singleparticledimension, index+1, occupations, results)
        end
    end
    return results
end

function distribute_fermions(particlenumber::Int, singleparticledimension::Int, index::Int=1, occupations::Vector{Int}=zeros(Int,singleparticledimension), results::Vector{Vector{Int}}=Vector{Int}[])
    if (singleparticledimension-index)+1<particlenumber
        return results
    end
    if index==singleparticledimension
        occupations[index] = particlenumber
        push!(results, copy(occupations))
    else
        for n=min(1,particlenumber):-1:0
            occupations[index] = n
            distribute_fermions(particlenumber-n, singleparticledimension, index+1, occupations, results)
        end
    end
    return results
end

"""
Abstract basis for n-particle systems. (Second quantization)

The inheriting bosonic and fermionic bases all have the same fields:

particlenumber
    (Fixed) number of particles in the system.
modenumber
    Number of modes in which the particles can be distributed.
occupations
    Basis states that specify how many particles are in which modes.

The only difference between Fermionic and Bosonic basis is which basis states
are included. These states are created automatically if the occupation
states are ommited when creating a FermionicNParticleBasis or a
BosonicNParticleBasis.
"""
abstract NParticleBasis <: Basis

function check_NParticleBasis_arguments(particlenumber::Int, modenumber::Int, occupations::Vector{Vector{Int}})
    particlenumber < 0 && throw(ArgumentError("Can't have less than zero particles."))
    for occupation=occupations
        modenumber != length(occupation) && throw(ArgumentError("Dimension of single particle basis has to be equal to the dimension of the N-particle basis vector."))
        any(occupation.<0) && throw(ArgumentError("Occupation numbers smaller than zero not possible."))
        sum(occupation) != particlenumber && throw(ArgumentError("Total occupation has to be equal to the particle number."))
    end
end

type BosonicNParticleBasis <: NParticleBasis
    shape::Vector{Int}
    particlenumber::Int
    modenumber::Int
    occupations::Vector{Vector{Int}}
    occupations_hash::UInt

    function BosonicNParticleBasis(particlenumber::Int, modenumber::Int, occupations::Vector{Vector{Int}})
        check_NParticleBasis_arguments(particlenumber, modenumber, occupations)
        new([length(occupations)], particlenumber, modenumber, occupations, hash(occupations))
    end
end

type FermionicNParticleBasis <: NParticleBasis
    shape::Vector{Int}
    particlenumber::Int
    modenumber::Int
    occupations::Vector{Vector{Int}}
    occupations_hash::UInt

    function FermionicNParticleBasis(particlenumber::Int, modenumber::Int, occupations::Vector{Vector{Int}})
        check_NParticleBasis_arguments(particlenumber, modenumber, occupations)
        for occupation=occupations
            any(occupation.>1) && throw(ArgumentError("Occupation numbers greater than zero not possible for Fermions."))
        end
        new([length(occupations)], particlenumber, modenumber, occupations, hash(occupations))
    end
end

BosonicNParticleBasis(particlenumber::Int, modenumber::Int) = BosonicNParticleBasis(particlenumber, modenumber, distribute_bosons(particlenumber, modenumber))
FermionicNParticleBasis(particlenumber::Int, modenumber::Int) = FermionicNParticleBasis(particlenumber, modenumber, distribute_fermions(particlenumber, modenumber))

=={T<:NParticleBasis}(b1::T, b2::T) = (b1.particlenumber==b2.particlenumber && b1.modenumber==b2.modenumber && b1.occupations_hash==b2.occupations_hash)

function nparticleoperator(nparticlebasis::NParticleBasis, op::DenseOperator; rank::Int=1)
    rank<1 && throw(ArgumentError("Rank has to be greater than zero."))
    rank>=2 && throw(ArgumentError("Not yet implemented for ranks greater than one."))
    result = Operator(nparticlebasis)
    N = length(nparticlebasis)
    S = nparticlebasis.modenumber
    for m=1:N, n=1:N
        if m==n
            for i=1:nparticlebasis.modenumber
                result.data[m,m] += op.data[i,i]*basis.occupations[m][i]
            end
        end
        indices = indices_adaggeri_aj(basis.occupations[m], basis.occupations[n])
        if !isnull(indices)
            i, j = get(indices)
            Ni = basis.occupations[n][i]
            Nj = basis.occupations[n][j]+1
            op.data[m,n] += result.data[i,j]*sqrt(Ni*Nj)
        end
    end
    return op
end

"""
Create a N-particle operator from a single particle operator.

The mathematical formalism is described by

.. math::

    X = \\sum_{ij} a_i^\\dagger a_j
                    \\left\\langle u_i \\right|
                    x
                    \\left| u_j \\right\\rangle

where :math:`X` is the N-particle operator, :math:`x` is the single particle operator and
:math:`\\left| u \\right\\rangle` are the single particle states associated to the
different modes of the N-particle basis.

Arguments
---------

nparticlebasis
    NParticleBasis
op
    An operator represented in the single-particle functions associated to the
    modes of the N-particle basis. This means the dimension of this operator has to
    be equal to the modenumber of the N-particle basis.
"""
function nparticleoperator_1(nparticlebasis::NParticleBasis, op::DenseOperator)
    N = length(nparticlebasis)
    S = nparticlebasis.modenumber
    @assert S == length(op.basis_l)
    @assert S == length(op.basis_r)
    result = DenseOperator(nparticlebasis)
    occupations = nparticlebasis.occupations
    for m=1:N, n=1:N
        for i=1:S, j=1:S
            C = coeff(occupations[m], occupations[n], [i], [j])
            if C!=0.
                result.data[m,n] += C*op.data[i,j]
            end
        end
    end
    return result
end

function nparticleoperator_1(nparticlebasis::NParticleBasis, op::SparseOperator)
    N = length(nparticlebasis)
    S = nparticlebasis.modenumber
    @assert S == length(op.basis_l)
    @assert S == length(op.basis_r)
    result = SparseOperator(nparticlebasis)
    occupations = nparticlebasis.occupations
    rows = rowvals(op.data)
    values = nonzeros(op.data)
    for column=1:S, j in nzrange(op.data, column)
        row = rows[j]
        value = values[j]
        for m=1:N, n=1:N
            C = coeff(occupations[m], occupations[n], [row], [column])
            if C!=0.
                result.data[m,n] += C*value
            end
        end
    end
    return result
end

"""
Create a N-particle operator from a two-particle operator.

The mathematical formalism is described by

.. math::

    X = \\sum_{ijkl} a_i^\\dagger a_j^\\dagger a_k a_l
            \\left\\langle u_i \\right| \\left\\langle u_j \\right|
            x
            \\left| u_k \\right\\rangle \\left| u_l \\right\\rangle

where :math:`X` is the N-particle operator, :math:`x` is the two-particle operator and
:math:`\\left| u \\right\\rangle` are the single particle states associated to the
different modes of the N-particle basis.

Arguments
---------

nparticlebasis
    NParticleBasis
op
    A two particle operator represented in the single particle functions associated to the
    modes of the N-particle basis. This means the dimension of this operator has to
    be equal to the modenumber of the N-particle basis.
"""
function nparticleoperator_2(nparticlebasis::NParticleBasis, op::DenseOperator)
    N = length(nparticlebasis)
    S = nparticlebasis.modenumber
    @assert S^2 == length(op.basis_l)
    @assert S^2 == length(op.basis_r)
    result = DenseOperator(nparticlebasis)
    op_data = reshape(op.data, S, S, S, S)
    occupations = nparticlebasis.occupations
    for m=1:N, n=1:N
        for i=1:S, j=1:S, k=1:S, l=1:S
            C = coeff(occupations[m], occupations[n], [i, j], [k, l])
            result.data[m,n] += C*op_data[i, j, k, l]
        end
    end
    return result
end

function nparticleoperator_2(nparticlebasis::NParticleBasis, op::SparseOperator)
    N = length(nparticlebasis)
    S = nparticlebasis.modenumber
    @assert S^2 == length(op.basis_l)
    @assert S^2 == length(op.basis_r)
    result = SparseOperator(nparticlebasis)
    occupations = nparticlebasis.occupations
    rows = rowvals(op.data)
    values = nonzeros(op.data)
    for column=1:S^2, j in nzrange(op.data, column)
        row = rows[j]
        value = values[j]
        for m=1:N, n=1:N
            # println("row:", row, " column:"column, ind_left)
            index = ind2sub((S, S, S, S), (column-1)*S^2 + row)
            C = coeff(occupations[m], occupations[n], index[1:2], index[3:4])
            if C!=0.
                result.data[m,n] += C*value
            end
        end
    end
    return result
end

"""
Calculate the matrix element <{m}|at_1...at_n a_1...a_n|{n}>.
"""
function coeff(occ_m, occ_n, at_indices, a_indices)
    occ_m = deepcopy(occ_m)
    occ_n = deepcopy(occ_n)
    C = 1.
    for i=at_indices
        if occ_m[i] == 0
            return 0.
        end
        C *= sqrt(occ_m[i])
        occ_m[i] -= 1
    end
    for i=a_indices
        if occ_n[i] == 0
            return 0.
        end
        C *= sqrt(occ_n[i])
        occ_n[i] -= 1
    end
    if occ_m == occ_n
        return C
    else
        return 0.
    end
end

"""
Find indices i and j where <{m}|a^t_i a_j |{n}> is not zero.

Returns null if this condition is never fullfilled for the given states.
"""
function indices_adaggeri_aj(occ_m::Vector{Int}, occ_n::Vector{Int})
    null = Nullable{Tuple{Int,Int}}()
    idx_create = 0
    idx_destroy = 0
    for i=1:length(occ_m)
        if occ_m[i]==occ_n[i]
            continue
        end
        delta = occ_n[i] - occ_m[i]
        if delta==-1
            if idx_destroy != 0
                return null
            end
            idx_destroy = i
            N_destroy = occ_n[i]
        elseif delta==1
            if idx_create != 0
                return null
            end
            idx_create = i
            N_create = occ_n[i]+1
        else
            return null
        end
    end
    return Nullable((idx_create, idx_destroy))
end

"""
Find indices i and j where <{m}|a^t_i a_j |{n}> is not zero.

Returns null if this condition is never fullfilled for the given states.
"""
function indices_adaggeri_adaggerj_ak_al(occ_m::Vector{Int}, occ_n::Vector{Int})
    null = Nullable{Tuple{Int,Int,Int,Int}}()
    i = 0
    j = 0
    k = 0
    l = 0
    for i=1:length(occ_m)
        if occ_m[i]==occ_n[i]
            continue
        end
        delta = occ_n[i] - occ_m[i]
        if delta==-1
            if idx_destroy != 0
                return null
            end
            idx_destroy = i
            N_destroy = occ_n[i]
        elseif delta==1
            if idx_create != 0
                return null
            end
            idx_create = i
            N_create = occ_n[i]+1
        else
            return null
        end
    end
    return Nullable((idx_create, idx_destroy))
end

end


text:
#### solver attributes
#MOI.get(m::Union{MosekSolver,MosekModel},::MOI.SupportsDuals) = true
#MOI.get(m::Union{MosekSolver,MosekModel},::MOI.SupportsAddConstraintAfterSolve) = true
#MOI.get(m::Union{MosekSolver,MosekModel},::MOI.SupportsAddVariableAfterSolve) = true
#MOI.get(m::Union{MosekSolver,MosekModel},::MOI.SupportsDeleteConstraint) = true
#MOI.get(m::Union{MosekSolver,MosekModel},::MOI.SupportsDeleteVariable) = true
#MOI.get(m::Union{MosekSolver,MosekModel},::MOI.SupportsConicThroughQuadratic) = false # though actually the solver does

#### objective
MOI.get(m::MosekModel,attr::MOI.ObjectiveValue) = getprimalobj(m.task,m.solutions[attr.resultindex].whichsol)
#MOI.canget(m::MosekSolver,attr::MOI.ObjectiveValue) = true
MOI.canget(m::MosekModel,attr::MOI.ObjectiveValue) = attr.resultindex > 0 && attr.resultindex <= length(m.solutions)

MOI.canget(m::MosekModel,attr::MOI.ObjectiveBound) = getintinf(m.task,MSK_IINF_MIO_NUM_RELAX) > 0
MOI.get(m::MosekModel,attr::MOI.ObjectiveBound) = getdouinf(m.task,MSK_DINF_MIO_OBJ_BOUND)

MOI.canget(m::MosekModel,attr::MOI.RelativeGap) = getdouinf(m.task,MSK_DINF_MIO_OBJ_REL_GAP) > -1.0
MOI.get(m::MosekModel,attr::MOI.RelativeGap) = getdouinf(m.task,MSK_DINF_MIO_OBJ_REL_GAP)

#MOI.canget(m::MosekSolver,attr::MOI.SolveTime) = true
MOI.canget(m::MosekModel,attr::MOI.SolveTime) = true
MOI.get(m::MosekModel,attr::MOI.SolveTime) = getdouinf(m.task,MSK_DINF_OPTIMIZER_TIME)

# NOTE: The MOSEK interface currently only supports Min and Max
function MOI.get(m::MosekModel,attr::MOI.ObjectiveSense)
    sense = getobjsense(m.task)
    if sense == MSK_OBJECTIVE_SENSE_MINIMIZE
        MOI.MinSense
    else
        MOI.MaxSense
    end
end

function MOI.set!(m::MosekModel,attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    if sense == MOI.MinSense
        putobjsense(m.task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    elseif sense == MOI.MaxSense
        putobjsense(m.task,MSK_OBJECTIVE_SENSE_MAXIMIZE)
    else
        @assert sense == MOI.FeasibilitySense
        putobjsense(m.task,MSK_OBJECTIVE_SENSE_MINIMIZE)
        MOI.set!(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.VariableIndex[], Float64[], 0.))
    end
end

#### Solver/Solution information

#MOI.canget(m::MosekSolver,attr::MOI.SimplexIterations) = true
MOI.canget(m::MosekModel,attr::MOI.SimplexIterations) = true
function MOI.get(m::MosekModel,attr::MOI.SimplexIterations)
    miosimiter = getlintinf(m.task,MSK_LIINF_MIO_SIMPLEX_ITER)
    if miosimiter > 0
        Int(miosimiter)
    else
        Int(getintinf(m.task,MSK_IINF_SIM_PRIMAL_ITER) + getintinf(m.task,MSK_IINF_SIM_DUAL_ITER))
    end
end

#MOI.canget(m::MosekSolver,attr::MOI.BarrierIterations) = true
MOI.canget(m::MosekModel,attr::MOI.BarrierIterations) = true
function MOI.get(m::MosekModel,attr::MOI.BarrierIterations)
    miosimiter = getlintinf(m.task,MSK_LIINF_MIO_INTPNT_ITER)
    if miosimiter > 0
        Int(miosimiter)
    else
        Int(getintinf(m.task,MSK_IINF_INTPNT_ITER))
    end
end

#MOI.canget(m::MosekSolver,attr::MOI.NodeCount) = true
MOI.canget(m::MosekModel,attr::MOI.NodeCount) = true
function MOI.get(m::MosekModel,attr::MOI.NodeCount)
        Int(getintinf(m.task,MSK_IINF_MIO_NUM_BRANCH))
end

#MOI.canget(m::MosekSolver,attr::MOI.RawSolver) = true
MOI.canget(m::MosekModel,attr::MOI.RawSolver) = true
MOI.get(m::MosekModel,attr::MOI.RawSolver) = m.task

#MOI.canget(m::MosekSolver,attr::MOI.ResultCount) = true
MOI.canget(m::MosekModel,attr::MOI.ResultCount) = true
MOI.get(m::MosekModel,attr::MOI.ResultCount) = length(m.solutions)

#### Problem information

#MOI.canget(m::MosekSolver,attr::MOI.NumberOfVariables) = true
MOI.canget(m::MosekModel,attr::MOI.NumberOfVariables) = true
MOI.get(m::MosekModel,attr::MOI.NumberOfVariables) = m.publicnumvar

#MOI.canget(m::MosekSolver,attr::MOI.NumberOfConstraints) = true
MOI.canget(m::MosekModel,attr::MOI.NumberOfConstraints{F,D}) where {F,D} = true
MOI.get(m::MosekModel,attr::MOI.NumberOfConstraints{F,D}) where {F,D} = length(select(m.constrmap,F,D))

#MOI.canget(m::MosekSolver,attr::MOI.ListOfVariableIndices) = false
#MOI.canget(m::MosekSolver,attr::MOI.ListOfConstraints) = false

#MOI.canget{F,D}(m::MosekSolver,attr::MOI.ListOfConstraintIndices{F,D}) = true
#MOI.get{F,D}(m::MosekSolver,attr::MOI.ListOfConstraintIndices{F,D}) = keys(select(m.constrmap,F,D))

#### Warm start values

function MOI.set!(m::MosekModel,attr::MOI.VariablePrimalStart, v :: MOI.VariableIndex, val::Float64)
    subj = getindexes(m.x_block,ref2id(v))

    xx = Float64[val]
    for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
        putxxslice(m.task,sol,Int(subj),Int(subj+1),xx)
    end
end

function MOI.set!(m::MosekModel,attr::MOI.VariablePrimalStart, vs::Vector{MOI.VariableIndex}, vals::Vector{Float64})
    subj = Array{Int}(undef,length(vs))
    for i in 1:length(subj)
        getindexes(m.x_block,ref2id(vs[i]),subj,i)
    end

    for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
        if solutiondef(m.task,sol)
            xx = getxx(m.task,sol)
            xx[subj] = vals
            putxx(m.task,sol,xx)
        else
            xx = zeros(Float64,getnumvar(m.task))
            xx[subj] = vals
            putxx(m.task,sol,xx)
        end
    end
end

function MOI.set!(m::MosekModel, attr::MOI.VariableName, index :: MOI.VariableIndex, value :: String)
    subj = getindexes(m.x_block, ref2id(index))
    putvarname(m.task,subj[1],value)
end

# function MOI.set!(m::MosekModel,attr::MOI.ConstraintDualStart, vs::Vector{MOI.ConstraintIndex}, vals::Vector{Float64})
#     subj = Array{Int}(length(vs))
#     for i in 1:length(subj)
#         getindexes(m.x_block,ref2id(vs[i]),subj,i)
#     end

#     for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
#         if solutiondef(m.task,sol)
#             xx = getxx(m.task,sol)
#             xx[subj] = vals
#             putxx(m.task,sol,xx)
#         else
#             xx = zeros(Float64,getnumvar(m.task))
#             xx[subj] = vals
#             putxx(m.task,sol,xx)
#         end
#     end
# end

#### Variable solution values
MOI.canget(m::MosekModel,attr::MOI.VariablePrimal) = attr.N > 0 && attr.N <= length(m.solutions)
#MOI.canget(m::MosekSolver,attr::MOI.VariablePrimal) = true

MOI.canget(m::MosekModel,attr::MOI.VariablePrimal,::Type{MOI.VariableIndex}) = MOI.canget(m,attr)

function MOI.get!(output::Vector{Float64},m::MosekModel,attr::MOI.VariablePrimal, vs::Vector{MOI.VariableIndex})
    subj = Array{Int}(undef,length(vs))
    for i in 1:length(subj)
        getindexes(m.x_block,ref2id(vs[i]),subj,i)
    end

    output[1:length(output)] = m.solutions[attr.N].xx[subj]
end

function MOI.get(m::MosekModel,attr::MOI.VariablePrimal, vs::Vector{MOI.VariableIndex})
    output = Vector{Float64}(undef,length(vs))
    MOI.get!(output,m,attr,vs)
    output
end

function MOI.get(m::MosekModel,attr::MOI.VariablePrimal, vref::MOI.VariableIndex)
    subj = getindexes(m.x_block,ref2id(vref))[1]
    m.solutions[attr.N].xx[subj]
end

#MOI.canget(m::MosekSolver,attr::MOI.VariableBasisStatus) = false # is a solution selector missing?
MOI.canget(m::MosekModel,attr::MOI.VariableBasisStatus)  = false

#### Constraint solution values
function canget(m::MosekModel,attr::MOI.ConstraintPrimal)
    num = 0
    if solutiondef(m.task,MSK_SOL_ITG) num += 1 end
    if solutiondef(m.task,MSK_SOL_BAS) num += 1 end
    if solutiondef(m.task,MSK_SOL_ITR) num += 1 end
    attr.N > 0 && attr.N <= num
end

MOI.canget(m::MosekModel,attr::MOI.ConstraintPrimal, ::Type{<:MOI.ConstraintIndex}) = canget(m,attr)

function MOI.get(
    m     ::MosekModel,
    attr  ::MOI.ConstraintPrimal,
    cref  ::MOI.ConstraintIndex{MOI.SingleVariable,D}) where D

    conid = ref2id(cref)
    idxs  = getindexes(m.xc_block,conid)
    subj  = m.xc_idxs[idxs[1]]

    m.solutions[attr.N].xx[subj]
end

# Semidefinite domain for a variable
function MOI.get!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MOI.ConstraintPrimal,
    cref  ::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.PositiveSemidefiniteConeTriangle})

    whichsol = getsolcode(m,attr.N)
    cid = ref2id(cref)
    @assert(cid < 0)

    # It is in fact a real constraint and cid is the id of an ordinary constraint
    barvaridx = - m.c_block_slack[-cid]
    output[1:length(output)] = sympackedLtoU(getbarxj(m.task,whichsol,barvaridx))
end

# Any other domain for variable vector
function MOI.get!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MOI.ConstraintPrimal,
    cref  ::MOI.ConstraintIndex{MOI.VectorOfVariables,D}) where D

    xcid = ref2id(cref)
    @assert(xcid > 0)

    mask = m.xc_bounds[xcid]
    idxs = getindexes(m.xc_block,xcid)
    subj = m.xc_idxs[idxs]

    output[1:length(output)] = m.solutions[attr.N].xx[subj]
end

function MOI.get(m     ::MosekModel,
                 attr  ::MOI.ConstraintPrimal,
                 cref  ::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D}) where D
    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)[1]
    m.solutions[attr.N].xc[subi] + m.c_constant[subi]
end

function MOI.get!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MOI.ConstraintPrimal,
    cref  ::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},D}) where D

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)

    if     m.c_block_slack[cid] == 0 # no slack
        output[1:length(output)] = m.solutions[attr.N].xc[subi] + m.c_constant[subi]
    elseif m.c_block_slack[cid] >  0 # qcone slack
        xid = m.c_block_slack[cid]
        xsubj = getindexes(m.x_block, xid)
        output[1:length(output)] = m.solutions[attr.N].xx[xsubj]
    else # psd slack
        xid = - m.c_block_slack[cid]
        output[1:length(output)] = sympackedLtoU(getbarxj(m.task,m.solutions[attr.N].whichsol,Int32(xid)))
    end
end

MOI.canget(m::MosekModel,attr::MOI.ConstraintDual) = attr.N > 0 && attr.N <= length(m.solutions)
MOI.canget(m::MosekModel,attr::MOI.ConstraintDual,::Type{<:MOI.ConstraintIndex}) = MOI.canget(m,attr)

function MOI.get(
    m     ::MosekModel,
    attr  ::MOI.ConstraintDual,
    cref  ::MOI.ConstraintIndex{MOI.SingleVariable,D}) where { D <: MOI.AbstractSet }

    xcid = ref2id(cref)
    idxs = getindexes(m.xc_block,xcid) # variable ids

    @assert(blocksize(m.xc_block,xcid) > 0)

    subj  = getindexes(m.x_block, m.xc_idxs[idxs][1])[1]
    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        if     m.xc_bounds[xcid] & boundflag_lower != 0 && m.xc_bounds[xcid] & boundflag_upper != 0
            m.solutions[attr.N].slx[subj] - m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[xcid] & boundflag_lower) != 0
            m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[xcid] & boundflag_upper) != 0
            - m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[xcid] & boundflag_cone) != 0
            m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    else
        if     m.xc_bounds[xcid] & boundflag_lower != 0 && m.xc_bounds[xcid] & boundflag_upper != 0
            m.solutions[attr.N].sux[subj] - m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[xcid] & boundflag_lower) != 0
            - m.solutions[attr.N].slx[subj]
        elseif (m.xc_bounds[xcid] & boundflag_upper) != 0
            m.solutions[attr.N].sux[subj]
        elseif (m.xc_bounds[xcid] & boundflag_cone) != 0
            - m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    end
end

getsolcode(m::MosekModel, N) = m.solutions[N].whichsol

# Semidefinite domain for a variable
function MOI.get!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MOI.ConstraintDual,
    cref  ::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.PositiveSemidefiniteConeTriangle})

    whichsol = getsolcode(m,attr.N)
    cid = ref2id(cref)
    @assert(cid < 0)

    # It is in fact a real constraint and cid is the id of an ordinary constraint
    barvaridx = - m.c_block_slack[-cid]
    dual = sympackedLtoU(getbarsj(m.task,whichsol,barvaridx))
    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        output[1:length(output)] = dual
    else
        output[1:length(output)] = -dual
    end
end

# Any other domain for variable vector
function MOI.get!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MOI.ConstraintDual,
    cref  ::MOI.ConstraintIndex{MOI.VectorOfVariables,D}) where D

    xcid = ref2id(cref)
    @assert(xcid > 0)

    mask = m.xc_bounds[xcid]
    idxs = getindexes(m.xc_block,xcid)
    subj = m.xc_idxs[idxs]

    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        if     mask & boundflag_lower != 0 && mask & boundflag_upper != 0
            output[1:length(output)] = m.solutions[attr.N].slx[subj] - m.solutions[attr.N].sux[subj]
        elseif (mask & boundflag_lower) != 0
            output[1:length(output)] = m.solutions[attr.N].slx[subj]
        elseif (mask & boundflag_upper) != 0
            output[1:length(output)] = - m.solutions[attr.N].sux[subj]
        elseif (mask & boundflag_cone) != 0
            output[1:length(output)] = m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    else
        if     mask & boundflag_lower != 0 && mask & boundflag_upper != 0
            output[1:length(output)] = m.solutions[attr.N].sux[subj] - m.solutions[attr.N].slx[subj]
        elseif (mask & boundflag_lower) != 0
            output[1:length(output)] = - m.solutions[attr.N].slx[subj]
        elseif (mask & boundflag_upper) != 0
            output[1:length(output)] = m.solutions[attr.N].sux[subj]
        elseif (mask & boundflag_cone) != 0
            output[1:length(output)] = - m.solutions[attr.N].snx[subj]
        else
            error("Dual value available for this constraint")
        end
    end
end

function MOI.get(m     ::MosekModel,
                 attr  ::MOI.ConstraintDual,
                 cref  ::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D}) where D

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)[1]

    if getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE
        m.solutions[attr.N].y[subi]
    else
        - m.solutions[attr.N].y[subi]
    end
end

function MOI.get!(
    output::Vector{Float64},
    m     ::MosekModel,
    attr  ::MOI.ConstraintDual,
    cref  ::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},D}) where { D <: MOI.AbstractSet }

    cid = ref2id(cref)
    subi = getindexes(m.c_block,cid)

    if getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE
        if     m.c_block_slack[cid] == 0 # no slack
            output[1:length(output)] = m.solutions[attr.N].y[subi]
        elseif m.c_block_slack[cid] >  0 # qcone slack
            xid = m.c_block_slack[cid]
            xsubj = getindexes(m.x_block, xid)
            output[1:length(output)] = m.solutions[attr.N].snx[xsubj]
        else # psd slack
            whichsol = getsolcode(m,attr.N)
            xid = - m.c_block_slack[cid]
            output[1:length(output)] = sympackedLtoU(getbarsj(m.task,whichsol,Int32(xid)))
        end
    else
        if     m.c_block_slack[cid] == 0 # no slack
            output[1:length(output)] = - m.solutions[attr.N].y[subi]
        elseif m.c_block_slack[cid] >  0 # qcone slack
            xid = m.c_block_slack[cid]
            subj = getindexes(m.x_block, xid)
            output[1:length(output)] = - m.solutions[attr.N].snx[subj]
        else # psd slack
            whichsol = getsolcode(m,attr.N)
            xid = - m.c_block_slack[cid]
            output[1:length(output)] = sympackedLtoU(- getbarsj(m.task,whichsol,Int32(xid)))
        end
    end
end

solsize(m::MosekModel, cref :: MOI.ConstraintIndex{<:MOI.AbstractScalarFunction}) = 1
function solsize(m::MosekModel, cref :: MOI.ConstraintIndex{MOI.VectorOfVariables})
    cid = ref2id(cref)
    if cid < 0
        blocksize(m.c_block,-cid)
    else
        blocksize(m.xc_block,cid)
    end
end

function solsize(m::MosekModel, cref :: MOI.ConstraintIndex{<:MOI.VectorAffineFunction})
    cid = ref2id(cref)
    blocksize(m.c_block,cid)
end

function MOI.get(m::MosekModel, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual}, cref :: MOI.ConstraintIndex{F,D}) where {F <: MOI.AbstractVectorFunction,D}
    cid = ref2id(cref)
    output = Vector{Float64}(undef,solsize(m,cref))
    MOI.get!(output,m,attr,cref)
    output
end

#### Status codes
MOI.canget(m::MosekModel,attr::MOI.TerminationStatus) = true
function MOI.get(m::MosekModel,attr::MOI.TerminationStatus)
    if     m.trm == MSK_RES_OK
        MOI.Success
    elseif m.trm == MSK_RES_TRM_MAX_ITERATIONS
        MOI.IterationLimit
    elseif m.trm == MSK_RES_TRM_MAX_TIME
        MOI.TimeLimit
    elseif m.trm == MSK_RES_TRM_OBJECTIVE_RANGE
        MOI.ObjectiveLimit
    elseif m.trm == MSK_RES_TRM_MIO_NEAR_REL_GAP
        MOI.AlmostSuccess
    elseif m.trm == MSK_RES_TRM_MIO_NEAR_ABS_GAP
        MOI.AlmostSuccess
    elseif m.trm == MSK_RES_TRM_MIO_NUM_RELAXS
        MOI.OtherLimit
    elseif m.trm == MSK_RES_TRM_MIO_NUM_BRANCHES
        MOI.NodeLimit
    elseif m.trm == MSK_RES_TRM_NUM_MAX_NUM_INT_SOLUTIONS
        MOI.SolutionLimit
    elseif m.trm == MSK_RES_TRM_STALL
        MOI.SlowProgress
    elseif m.trm == MSK_RES_TRM_USER_CALLBACK
        MOI.Interrupted
    elseif m.trm == MSK_RES_TRM_MAX_NUM_SETBACKS
        MOI.OtherLimit
    elseif m.trm == MSK_RES_TRM_NUMERICAL_PROBLEM
        MOI.SlowProgress
    elseif m.trm == MSK_RES_TRM_INTERNAL
        MOI.OtherError
    elseif m.trm == MSK_RES_TRM_INTERNAL_STOP
        MOI.OtherError
    else
        MOI.OtherError
    end
end

function MOI.canget(m::MosekModel,attr::MOI.PrimalStatus)
    if attr.N < 0 || attr.N > length(m.solutions)
        false
    else
        solsta = m.solutions[attr.N].solsta
        if     solsta == MSK_SOL_STA_UNKNOWN true
        elseif solsta == MSK_SOL_STA_OPTIMAL true
        elseif solsta == MSK_SOL_STA_PRIM_FEAS true
        elseif solsta == MSK_SOL_STA_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER false
        elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER true
        elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER false
        elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER true
        elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL true
        else false
        end
    end
end
function MOI.get(m::MosekModel,attr::MOI.PrimalStatus)
    solsta = m.solutions[attr.N].solsta
    if     solsta == MSK_SOL_STA_UNKNOWN
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_OPTIMAL
        MOI.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_FEAS
        MOI.FeasiblePoint
    elseif solsta == MSK_SOL_STA_DUAL_FEAS
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS
        MOI.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        MOI.InfeasibilityCertificate
    elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER
        MOI.ReductionCertificate
    elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL
        MOI.FeasiblePoint
    else
        MOI.UnknownResultStatus
    end
end

function MOI.canget(m::MosekModel,attr::MOI.DualStatus)
    if attr.N < 0 || attr.N > length(m.solutions)
        false
    else
        solsta = m.solutions[attr.N].solsta
        if     solsta == MSK_SOL_STA_UNKNOWN true
        elseif solsta == MSK_SOL_STA_OPTIMAL true
        elseif solsta == MSK_SOL_STA_PRIM_FEAS true
        elseif solsta == MSK_SOL_STA_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS true
        elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER true
        elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER false
        elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER true
        elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER false
        elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL false
        else false
        end
    end
end
function MOI.get(m::MosekModel,attr::MOI.DualStatus)
    solsta = m.solutions[attr.N].solsta
    if     solsta == MSK_SOL_STA_UNKNOWN
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_OPTIMAL
        MOI.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_FEAS
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_DUAL_FEAS
        MOI.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS
        MOI.FeasiblePoint
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        MOI.InfeasibilityCertificate
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER
        MOI.ReductionCertificate
    elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER
        MOI.UnknownResultStatus
    elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL
        MOI.FeasiblePoint
    else
        MOI.UnknownResultStatus
    end
end


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

## object-oriented RE2 interface via cre2 ##

export @rr_str, Regex2, Regex2Match
import Base: compile, contains, match, matchall, findfirst, findnext, eachmatch

import Base: getindex, ==, show, print_quoted_literal, ensureroom
import Base: eltype, next, done, start, IteratorSize, SizeUnknown, SubstitutionString

include("cre2.jl")

export        ANCHORED, ENDANCHORED,
              LOG_ERRORS, LATIN1, LONGEST_MATCH, LITERAL, NEVER_NL, DOTALL, NEVER_CAPTURE,
              CASELESS, POSIX_SYNTAX, PERL_CLASSES, WORD_BOUNDARY, MULTILINE

import .CRE2: ANCHORED, ENDANCHORED,
              LOG_ERRORS, LATIN1, LONGEST_MATCH, LITERAL, NEVER_NL, DOTALL, NEVER_CAPTURE,
              CASELESS, POSIX_SYNTAX, PERL_CLASSES, WORD_BOUNDARY, MULTILINE

const DEFAULT_COMPILER_OPTS = PERL_CLASSES | WORD_BOUNDARY
const DEFAULT_MATCH_OPTS = ANCHORED & ~ANCHORED

mutable struct Regex2
    pattern::String
    compile_options::UInt32
    match_options::UInt32
    regex::Ptr{Cvoid}
    match_data::Vector{Ptr{Cvoid}}

    function Regex2(pattern::AbstractString, compile_options::Integer,
                   match_options::Integer)

        pattern = String(pattern)
        compile_options = UInt32(compile_options & ~LOG_ERRORS) # silently ignore
        match_options = UInt32(match_options)
        if (compile_options & ~CRE2.COMPILE_MASK) != 0
            throw(ArgumentError("invalid regex compile options: $compile_options"))
        end
        if (match_options & ~CRE2.EXECUTE_MASK) !=0
            throw(ArgumentError("invalid regex match options: $match_options"))
        end
        if compile_options & POSIX_SYNTAX == 0
            compile_options |= (PERL_CLASSES|WORD_BOUNDARY)
            compile_options &= ~MULTILINE
        end
        re = new(pattern, compile_options, match_options, C_NULL, Ptr{Cvoid}[])
        compile(re)

        finalizer(re) do re
            re.regex == C_NULL || CRE2.free_re(re.regex)
        end
        re
    end
end

using Printf
function Regex2(pattern::AbstractString, flags::AbstractString)
    options = DEFAULT_COMPILER_OPTS
    for f in flags
        options |= f=='i' ? CASELESS  :
                   f=='m' ? MULTILINE :
                   f=='s' ? DOTALL    :
                   f=='l' ? LONGEST_MATCH  :
                   f=='p' ? POSIX_SYNTAX  :
                   throw(ArgumentError("unknown regex flag: $f"))
    end
    Regex2(pattern, options, DEFAULT_MATCH_OPTS)
end
Regex2(pattern::AbstractString) = Regex2(pattern, DEFAULT_COMPILER_OPTS, DEFAULT_MATCH_OPTS)

function compile(re::Regex2)
    if re.regex == C_NULL
        copts = CRE2.new_options(re.compile_options)
        re.regex = CRE2.compile(re.pattern, copts)
        re.match_data = CRE2.create_match_data(re.regex)
        CRE2.free_options(copts)
    end
    re
end

"""
    @rr_str -> Regex2

Construct a regex, such as `rr"^[a-z]*\$"`. The regex also accepts one or more flags,
listed after the ending quote, to change its behaviour:

- `i` enables case-insensitive matching
- `m` treats the `^` and `\$` tokens as matching the start and end of individual lines, as
  opposed to the whole string.
- `s` allows the `.` modifier to match newlines.
- `l` enables "longest match mode".
- `p` enables "POSIX syntax" - like egrep.

For example, this regex has all three flags enabled:

```jldoctest
julia> match(rr"a+.*b+.*?d\$"ism, "Goodbye,\\nOh, angry,\\nBad world\\n")
Regex2Match("angry,\\nBad world")
```
"""
macro rr_str(pattern, flags...) Regex2(pattern, flags...) end

function show(io::IO, re::Regex2)
    opts = re.compile_options
    imsx = CASELESS|MULTILINE|DOTALL|LONGEST_MATCH|POSIX_SYNTAX|LOG_ERRORS
    opts & POSIX_SYNTAX == 0 && ( imsx |= PERL_CLASSES|WORD_BOUNDARY )

    if opts & ~imsx == 0
        print(io, "rr")
        print_quoted_literal(io, re.pattern)
        if opts & CASELESS != 0; print(io, 'i'); end
        if opts & DOTALL != 0; print(io, 's'); end
        if opts & LONGEST_MATCH != 0; print(io, 'l'); end
        if opts & POSIX_SYNTAX != 0; print(io, 'p'); end
        if opts & POSIX_SYNTAX != 0
            if (opts & MULTILINE) != 0; print(io, 'm'); end
        end
    else
        print(io, "Regex2(")
        show(io, re.pattern)
        print(io, ',')
        show(io, opts)
        print(io, ')')
    end
end

# TODO: map offsets into strings in other encodings back to original indices.
# or maybe it's better to just fail since that would be quite slow

struct Regex2Match
    match::SubString{String}
    captures::Vector{Union{Nothing,SubString{String}}}
    offset::Int
    offsets::Vector{Int}
    regex::Regex2
end

function show(io::IO, m::Regex2Match)
    print(io, "Regex2Match(")
    show(io, m.match)
    idx_to_capture_name = CRE2.capture_names(m.regex.regex, m.regex.pattern)
    if !isempty(m.captures)
        print(io, ", ")
        for i = 1:length(m.captures)
            # If the capture group is named, show the name.
            # Otherwise show its index.
            capture_name = get(idx_to_capture_name, i, i)
            print(io, capture_name, "=")
            show(io, m.captures[i])
            if i < length(m.captures)
                print(io, ", ")
            end
        end
    end
    print(io, ")")
end

# Capture group extraction
getindex(m::Regex2Match, idx::Integer) = m.captures[idx]
function getindex(m::Regex2Match, name::AbstractString)
    idx = CRE2.group_number_from_name(m.regex.regex, name)
    idx <= 0 && error("no capture group named $name found in regex")
    m[idx]
end
getindex(m::Regex2Match, name::Symbol) = m[string(name)]

function contains(s::AbstractString, r::Regex2, offset::Integer=0)
    compile(r)
    return CRE2.exec(r.regex, String(s), offset, get_match_options(r), r.match_data)
end

function contains(s::SubString, r::Regex2, offset::Integer=0)
    compile(r)
    return CRE2.exec(r.regex, s, offset, get_match_options(r), r.match_data)
end

(r::Regex2)(s) = contains(s, r)

"""
    match(r::Regex2, s::AbstractString[, idx::Integer[, addopts]])

Search for the first match of the regular expression `r` in `s` and return a `Regex2Match`
object containing the match, or nothing if the match failed. The matching substring can be
retrieved by accessing `m.match` and the captured sequences can be retrieved by accessing
`m.captures` The optional `idx` argument specifies an index at which to start the search.
The optional addopts can be zero, `ANCHORED`, or `ANCHORED | ENDANCHORED`. `ENDANCHORED`
allone is not supported.

# Examples
```jldoctest
julia> rx = rr"a(.)a"
r"a(.)a"

julia> m = match(rx, "cabac")
Regex2Match("aba", 1="b")

julia> m.captures
1-element Array{Union{Nothing, SubString{String}},1}:
 "b"

julia> m.match
"aba"

julia> match(rx, "cabac", 3) == nothing
true
```
"""
function match end

function match(re::Regex2, str::Union{SubString{String}, String},
               idx::Integer, add_opts::Integer=0)

    compile(re)
    opts = get_match_options(re.match_options | UInt32(add_opts))
    CRE2.exec(re.regex, str, idx-1, opts, re.match_data) || return nothing
    md = re.match_data
    pstr = pointer(str)
    n = div(length(md),2) - 1
    ix1, ix2 = substrind(str, md, 0)
    mat = SubString(str, ix1, ix2)
    cap = Union{Nothing,SubString{String}}[md[2i+1] == Ptr{UInt8}(0) ?
                                nothing : substr(str, md, i) for i=1:n]

    off = Int[ substrind(str, md, i)[1] for i=1:n ]
    Regex2Match(mat, cap, ix1, off, re)
end

match(r::Regex2, s::AbstractString) = match(r, s, firstindex(s))
match(r::Regex2, s::AbstractString, i::Integer) = throw(ArgumentError(
    "regex matching is only available for the String type; use String(s) to convert"
))

function get_match_options(optmask::UInt32)
    opt = CRE2.UNANCHORED
    optmask & ANCHORED != 0 && ( opt = max(opt, CRE2.ANCHOR_START) )
    optmask & ENDANCHORED != 0 && ( opt = max(opt, CRE2.ANCHOR_BOTH) )
    opt
end
get_match_options(re::Regex2) = get_match_options(re.match_options)

function substrind(str::Union{String,SubString{String}}, md::Vector{Ptr{Cvoid}}, i::Int)
    ix1 = Int(md[2i+1] + 1 - pointer(str))
    ix1, prevind(str, ix1 + Int(md[2i+2]))
end

function substr(str::Union{String,SubString{String}}, md::Vector{Ptr{Cvoid}}, i::Int)
    SubString(str, substrind(str, md, i)...)
end

function substrr(str::Union{String,SubString{String}}, md::Vector{Ptr{Cvoid}}, i::Int)
    colon(substrind(str, md, i)...)
end

"""
    matchall(r::Regex2, s::AbstractString; overlap::Bool = false]) -> Vector{AbstractString}

Return a vector of the matching substrings from [`eachmatch`](@ref).

# Examples
```jldoctest
julia> rx = r"a.a"
r"a.a"

julia> matchall(rx, "a1a2a3a")
2-element Array{SubString{String},1}:
 "a1a"
 "a3a"

julia> matchall(rx, "a1a2a3a", overlap = true)
3-element Array{SubString{String},1}:
 "a1a"
 "a2a"
 "a3a"
```
"""
function matchall(re::Regex2, str::String; overlap::Bool = false)
    regex = compile(re).regex
    regex2 = regex
    copts = re.compile_options
    longmatch = copts & LONGEST_MATCH != 0
    need2 = !longmatch
    n = nextind(str, lastindex(str))
    offset = Cuint(firstindex(str))
    opts = get_match_options(re)
    opts2 = get_match_options(ANCHORED|re.match_options)
    matches = SubString{String}[]

    while offset <= n && CRE2.exec(regex, str, offset-1, opts, re.match_data, 1)
        ix1, ix2 = substrind(str, re.match_data, 0)
        push!(matches, SubString(str, ix1, ix2))

        if ix1 > ix2 && !longmatch
            if need2
                regex2 = Regex2(re.pattern, copts|LONGEST_MATCH, re.match_options).regex
                need2 = false
            end
            if CRE2.exec(regex2, str, offset-1, opts2, re.match_data, 1)
                ix1, ix2 = substrind(str, re.match_data, 0)
                if ix1 <= ix2
                    push!(matches, SubString(str, ix1, ix2))
                end
            end
        end
        offset == n && break
        offset = Cuint(overlap ? nextind(str, ix1) : nextind(str, max(ix1, ix2)))
    end
    matches
end

matchall(re::Regex2, str::SubString; overlap::Bool = false) =
    matchall(re, String(str), overlap = overlap)

# TODO: return only start index and update deprecation
function findnext(re::Regex2, str::Union{String,SubString}, idx::Integer)
    if idx > nextind(str,lastindex(str))
        throw(BoundsError())
    end
    opts = get_match_options(re.match_options)
    compile(re)
    md = re.match_data
    CRE2.exec(re.regex, str, idx-1, opts, re.match_data) ? substrr(str, md, 0) : (0:-1)
end
findnext(r::Regex2, s::AbstractString, idx::Integer) = throw(ArgumentError(
    "regex search is only available for the String type; use String(s) to convert"
))
findfirst(r::Regex2, s::AbstractString) = findnext(r,s,firstindex(s))

replace_err(repl) = error("Bad replacement string: $repl")

function _write_capture(io, re::Regex2, group)
    len = CRE2.substring_length_bynumber(re.match_data, group)
    ensureroom(io, len+1)
    CRE2.substring_copy_bynumber(re.match_data, group,
        pointer(io.data, io.ptr), len+1)
    io.ptr += len
    io.size = max(io.size, io.ptr - 1)
end

function Base._replace(io, repl_s::SubstitutionString, str, r, re::Regex2)
    SUB_CHAR = '\\'
    GROUP_CHAR = 'g'
    LBRACKET = '<'
    RBRACKET = '>'
    repl = repl_s.string
    i = firstindex(repl)
    e = lastindex(repl)
    while i <= e
        if repl[i] == SUB_CHAR
            next_i = nextind(repl, i)
            next_i > e && replace_err(repl)
            if repl[next_i] == SUB_CHAR
                write(io, SUB_CHAR)
                i = nextind(repl, next_i)
            elseif isdigit(repl[next_i])
                group = parse(Int, repl[next_i])
                i = nextind(repl, next_i)
                while i <= e
                    if isdigit(repl[i])
                        group = 10group + parse(Int, repl[i])
                        i = nextind(repl, i)
                    else
                        break
                    end
                end
                _write_capture(io, re, group)
            elseif repl[next_i] == GROUP_CHAR
                i = nextind(repl, next_i)
                if i > e || repl[i] != LBRACKET
                    replace_err(repl)
                end
                i = nextind(repl, i)
                i > e && replace_err(repl)
                groupstart = i
                while repl[i] != RBRACKET
                    i = nextind(repl, i)
                    i > e && replace_err(repl)
                end
                #  TODO: avoid this allocation
                groupname = SubString(repl, groupstart, prevind(repl, i))
                if all(isdigit, groupname)
                    _write_capture(io, re, parse(Int, groupname))
                else
                    group = CRE2.group_number_from_name(re.regex, groupname)
                    group < 0 && replace_err("Group $groupname not found in regex $re")
                    _write_capture(io, re, group)
                end
                i = nextind(repl, i)
            else
                replace_err(repl)
            end
        else
            write(io, repl[i])
            i = nextind(repl, i)
        end
    end
end

mutable struct Regex2MatchIterator
    regex::Regex2
    string::String
    overlap::Bool
    regex2::Regex2
    longmatch::Bool

    function Regex2MatchIterator(regex::Regex2, string::AbstractString, ovr::Bool=false)
        new(regex, string, ovr, regex, regex.compile_options & LONGEST_MATCH != 0)
    end
end
compile(itr::Regex2MatchIterator) = (compile(itr.regex); itr)
eltype(::Type{Regex2MatchIterator}) = Regex2Match
start(itr::Regex2MatchIterator) = match(itr.regex, itr.string, 1, UInt32(0))
done(itr::Regex2MatchIterator, prev_match) = (prev_match === nothing)
IteratorSize(::Type{Regex2MatchIterator}) = SizeUnknown()

# Assumes prev_match is not nothing
function next(itr::Regex2MatchIterator, prev_match)
    prevempty = isempty(prev_match.match)

    if itr.overlap
        if !prevempty
            offset = nextind(itr.string, prev_match.offset)
        else
            offset = prev_match.offset
        end
    else
        offset = prev_match.offset + lastindex(prev_match.match)
    end

    regex = itr.regex
    if !itr.longmatch && prevempty
        if regex === itr.regex2
            itr.regex2 = Regex2(regex.pattern, regex.compile_options|LONGEST_MATCH, 0)
        end
        mat = match(itr.regex2, itr.string, offset, ANCHORED)
        if mat !== nothing && !isempty(mat.match)
            return (prev_match, mat)
        elseif offset > sizeof(itr.string)
            return (prev_match, nothing)
        end
    end

    offset = nextind(itr.string, min(offset, lastindex(itr.string)))
    mat = match(regex, itr.string, offset, UInt32(0))
    if mat != nothing && isempty(mat.match) &&
        offset == prev_match.offset > lastindex(itr.string)

        mat = nothing
    end
    return (prev_match, mat)
end

"""
    eachmatch(r::Regex2, s::AbstractString; overlap::Bool=false])

Search for all matches of a the regular expression `r` in `s` and return a iterator over the
matches. If overlap is `true`, the matching sequences are allowed to overlap indices in the
original string, otherwise they must be from distinct character ranges.

# Examples
```jldoctest
julia> rx = rr"a.a"
r"a.a"

julia> m = eachmatch(rx, "a1a2a3a")
Base.Regex2MatchIterator(r"a.a", "a1a2a3a", false)

julia> collect(m)
2-element Array{Regex2Match,1}:
 Regex2Match("a1a")
 Regex2Match("a3a")

julia> collect(eachmatch(rx, "a1a2a3a", overlap = true))
3-element Array{Regex2Match,1}:
 Regex2Match("a1a")
 Regex2Match("a2a")
 Regex2Match("a3a")
```
"""
eachmatch(re::Regex2, str::AbstractString; overlap = false) =
    Regex2MatchIterator(re, str, overlap)

## comparison ##

function ==(a::Regex2, b::Regex2)
    a.pattern == b.pattern && a.compile_options == b.compile_options && a.match_options == b.match_options
end

## hash ##
const hashre_seed = UInt === UInt64 ? 0x67e195eb8555e72d : 0xe32373e4
function hash(r::Regex2, h::UInt)
    h += hashre_seed
    h = hash(r.pattern, h)
    h = hash(r.compile_options, h)
    h = hash(r.match_options, h)
end