text:
export CuSubArray, CuSubVector, CuSubMatrix, CuSubVecOrMat

struct CuSubArray{T,N} <: AbstractCuArray{T,N}
    parent::CuArray{T}
    indices::Tuple
    offset::Int
    contigious::Bool
end

const CuSubVector{T} = CuSubArray{T,1}
const CuSubMatrix{T} = CuSubArray{T,2}
const CuSubVecOrMat{T} = Union{CuSubVector{T},CuSubMatrix{T}}

CuArray(x::CuSubArray) = copyto!(similar(x), x)
Base.Array(x::CuSubArray) = Array(CuArray(x))
Base.size(x::CuSubArray) = map(length, x.indices)
Base.size(x::CuSubArray, i::Int) = length(x.indices[i])
Base.strides(x::CuSubArray) = strides(x.parent)
Base.stride(x::CuSubArray, dim::Int) = stride(x.parent, dim)
Base.convert(::Type{Ptr{T}}, x::CuSubArray) where T = Ptr{T}(pointer(x))
Base.unsafe_convert(::Type{Ptr{T}}, x::CuSubArray) where T = Ptr{T}(pointer(x))
Base.pointer(x::CuSubArray, index::Int=1) = pointer(x.parent, x.offset+index)
iscontigious(x::CuSubArray) = x.contigious

function Base.view(x::CuArray{T,N}, I...) where {T,N}
    @assert length(I) == N
    I = Base.to_indices(x, I)
    offset = sub2ind(size(x), map(first,I)...) - 1
    function iscontigious(x::CuArray, I::Tuple)
        for i = 1:length(I)-1
            length(I[i]) == size(x,i) && continue
            length(I[i+1]) == 1 && continue
            return false
        end
        true
    end
    CuSubArray{T,N}(x, I, offset, iscontigious(x,I))
end

Base.show(io::IO, ::Type{CuSubArray{T,N}}) where {T,N} = print(io, "CuSubArray{$T,$N}")
Base.print_array(io::IO, X::CuSubArray) = Base.showarray(io, Array(X))


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# modules external to our own eco-system, rigorously separate using and import

using Parameters: @with_kw

using Random: shuffle

import Base: ==, length

using LinearAlgebra: norm, dot, mul!, I

using StaticArrays

# -----------------------------------------------------------------------------
# JuLIP, ACE, etc : just use import throughout, this avoids bugs

import JuLIP

import JuLIP: alloc_temp, alloc_temp_d,
              cutoff,
              evaluate, evaluate_d,
              evaluate!, evaluate_d!,
              SitePotential,
              z2i, i2z, numz,
              read_dict, write_dict,
              AbstractCalculator,
              Atoms,
              chemical_symbol,
              fltype, rfltype

import JuLIP.MLIPs: IPBasis, alloc_B, alloc_dB, combine

import JuLIP.Potentials: ZList, SZList, zlist
import JuLIP: JVec, AtomicNumber

ext: jl


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

## dummy stub for https://github.com/JuliaBinaryWrappers/OpenLibm_jll.jl
module OpenLibm_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libopenlibm

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libopenlibm_handle = C_NULL
libopenlibm_path = ""

if Sys.iswindows()
    const libopenlibm = "libopenlibm.dll"
elseif Sys.isapple()
    const libopenlibm = "@rpath/libopenlibm.3.dylib"
else
    const libopenlibm = "libopenlibm.so.3"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libopenlibm_handle = dlopen(libopenlibm)
    global libopenlibm_path = dlpath(libopenlibm_handle)
end

# JLLWrappers API
is_available() = true
get_libopenlibm_path() = libopenlibm_path

end  # module OpenLibm_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
const SINGLE = SingleComponent
const MULTI = MaterialCompounds{MOLAR,FRACTION}

abstract type AbstractVolumeSolver end
struct VolumeBisection <: AbstractVolumeSolver
    pts::Int64
end
VolumeBisection() = VolumeBisection(21)

struct CubicRoots <: AbstractVolumeSolver end
struct SUVA <: AbstractVolumeSolver end

abstract type HelmholtzModel <: ThermoState.ThermoModel end
abstract type CubicModel <: HelmholtzModel end

abstract type SaturationModel <: ThermoState.ThermoModel end
abstract type GasModel <: ThermoState.ThermoModel end
abstract type SatLiquidModel <: ThermoState.ThermoModel end

#catch all to build models that doesnt need automatically
#to stop erroring


text:
#=
Geometry in 2d euclidean space
=#

"""
deltaangle(a::Real, b::Real)

Return the minimum δ such that
    a + δ = mod(b, 2π)
"""
deltaangle(a::Real, b::Real) = atan2(sin(b-a), cos(b-a))

"""
Distance between two angles
"""
angledist( a::Real, b::Real ) = abs(deltaangle(a,b))

"""
Linear interpolation between angles
"""
lerp_angle(a::Real, b::Real, t::AbstractFloat) = a + deltaangle(a, b)*t

"""
True if the values are collinear within a tolerance
"""
function are_collinear(a::AbstractVec, b::AbstractVec, c::AbstractVec, tol::Float64=1e-8)
    # http://mathworld.wolfram.com/Collinear.html
    # if val = 0 then they are collinear
    val = a.x*(b.y-c.y) + b.x*(c.y-a.y)+c.x*(a.y-b.y)
    abs(val) < tol
end

"""
To find orientation of ordered triplet (p, q, r).
The function returns following values
0 --> p, q and r are colinear
1 --> Clockwise
2 --> Counterclockwise
"""
function orientation(P::VecE2, Q::VecE2, R::VecE2)
    val = (Q.y - P.y)*(R.x - Q.x) - (Q.x - P.x)*(R.y - Q.y)
    if val ≈ 0
        return 0  # colinear
    end
    return (val > 0) ? 1 : 2
end

function inertial2body(point::VecE2, reference::VecSE2)

    #=
    Convert a point in an inertial cartesian coordinate frame
    to be relative to a body's coordinate frame

    The body's position is given relative to the same inertial coordinate frame
    =#

    s, c = sin(reference.θ), cos(reference.θ)
    Δx = point.x - reference.x
    Δy = point.y - reference.y
    VecE2(c*Δx + s*Δy, c*Δy - s*Δx)
end
function inertial2body(point::VecSE2, reference::VecSE2)

    #=
    Convert a point in an inertial cartesian coordinate frame
    to be relative to a body's coordinate frame

    The body's position is given relative to the same inertial coordinate frame
    =#

    s, c = sin(reference.θ), cos(reference.θ)
    Δx = point.x - reference.x
    Δy = point.y - reference.y
    VecSE2(c*Δx + s*Δy, c*Δy - s*Δx, point.θ - reference.θ)
end
function body2inertial(point::VecE2, reference::VecSE2)

    #=
    Convert a point in a body-relative cartesian coordinate frame
    to be relative to a the inertial coordinate frame the body is described by
    =#

    c, s = cos(reference.θ), sin(reference.θ)
    VecE2(c*point.x -s*point.y + reference.x, s*point.x +c*point.y + reference.y)
end
function body2inertial(point::VecSE2, reference::VecSE2)

    #=
    Convert a point in a body-relative cartesian coordinate frame
    to be relative to a the inertial coordinate frame the body is described by
    =#

    c, s = cos(reference.θ), sin(reference.θ)
    VecSE2(c*point.x -s*point.y + reference.x, s*point.x +c*point.y + reference.y, reference.θ + point.θ)
end

include("1d.jl")
include("lines.jl")
include("line_segments.jl")
include("rays.jl")
include("projectiles.jl")
include("solids.jl")
include("hyperplanes.jl")


text:
function test_rotation_map(T)
    ϕ = T(pi)/4
    m = rotation_map(ϕ)
    x = [one(T), zero(T)]
    y = m(x)
    @test y[1] ≈ sqrt(T(2))/2
    @test y[2] ≈ sqrt(T(2))/2

    ϕ = T(pi)/4
    m = rotation_map(ϕ, 0, 0)
    x = [zero(T), one(T), zero(T)]
    y = m(x)
    @test y[1] ≈ 0
    @test y[2] ≈ sqrt(T(2))/2
    @test y[3] ≈ sqrt(T(2))/2

    # TODO: add more tests for a 3D rotation

    theta = T(rand())
    phi = T(rand())
    psi = T(rand())
    m2 = rotation_map(theta)
    test_generic_map(m2)
    m3 = rotation_map(phi, theta, psi)
    test_generic_map(m3)

    r = suitable_point_to_map(m2)
    @test norm(m2(r))≈norm(r)

    r = suitable_point_to_map(m3)
    @test norm(m3(r))≈norm(r)
    @test islinear(m3)
end

function test_cart_polar_map(T)
    m1 = CartToPolarMap{T}()
    test_generic_map(m1)
    @test !islinear(m1)
    @test isreal(m1)

    m2 = PolarToCartMap{T}()
    test_generic_map(m2)
    @test !islinear(m2)
    @test isreal(m2)

    @test inverse(m1) == m2
    @test inverse(m2) == m1
end

function test_applications(T)
    test_rotation_map(T)
    test_cart_polar_map(T)
end

@testset "applications" begin
    test_applications(Float64)
    test_applications(BigFloat)
end


text:
export debye_focus

"""
    debye_focus(L, N, NA, z; λ=550e-9, n=1)

Evaluates the Debye integral to get the electrical field at position `z` with
respect to the focus of the lens.
`z>0` means farther away from the lens than the focus
and vice versa.
Lens has a numerical aperture `NA`.
`λ` is the wavelength and `n` the refractive index.

Returns a quadratic array with size `(N, N)`.
The side length is `L`, specifically `fftpos(L, N)` are the coordinates.
"""
function debye_focus(L, N, NA, z; λ=550e-9, n=1)
    # half opening angle
    α = asin(NA / n)
    # wavevector
    k = calc_k(λ, n)

    # positions we evaluate the integral at
    # the sqrt(2) because we evaluate the integral over the whole
    # diagonal of the quadratic output array
    # later we rotate the values to get the full solution over a
    # quadratic
    # we only evaluate at positive values, because of symmetry

    # we choose more data points  (we actually would only need sqrt(2)* N/2
    # but since we do an interpolation later, we go for more data points
    xs = range(0, sqrt(2) * L / 2, length=10 * N)

    # 1D array storing the diagonal integral values
    out = zeros(ComplexF64, length(xs))

    # integrand
    fi(θ, r) = besselj0(k*r*sin(θ))*exp(-1im*k*z*cos(θ))*sin(θ)
    # loop over the xs positions we want to evaluate at
    Threads.@threads for j in eachindex(xs)
        integral, err = quadgk(θ -> fi(θ, xs[j]), zero(eltype(α)), α)
        out[j] = integral
    end
    # some prefactor
    pre_c = 2*π*im / λ
    out .*= pre_c

    # output coordinates where we finally want to have the values
    xs_out = fftpos(L, N)
    ys_out = copy(xs_out)

    out_2D = apply_rot_symmetry(out, xs, xs_out, ys_out)

    # normalize to 1
    return out_2D
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Add literature references to AdaptiveMHTuning docstring.

"""
    MCMCNoOpTuning <: MCMCTuningAlgorithm

No-op tuning, marks MCMC chains as tuned without performing any other changes
on them. Useful if chains are pre-tuned or tuning is an internal part of the
MCMC sampler implementation.
"""
struct MCMCNoOpTuning <: MCMCTuningAlgorithm end
export MCMCNoOpTuning

struct MCMCNoOpTuner <: AbstractMCMCTunerInstance end

(tuning::MCMCNoOpTuning)(chain::MCMCIterator) = MCMCNoOpTuner()

function MCMCNoOpTuning(tuning::MCMCNoOpTuning, chain::MCMCIterator)
    MCMCNoOpTuner()
end

function tuning_init!(tuner::MCMCNoOpTuner, chain::MCMCIterator)
    chain.info = MCMCIteratorInfo(chain.info, cycle = chain.info.cycle + 1)
    nothing
end

function tuning_postinit!(tuner::MCMCNoOpTuner, chain::MCMCIterator, samples::DensitySampleVector)
end

function tuning_update!(tuner::MCMCNoOpTuner, chain::MCMCIterator, samples::DensitySampleVector)
    chain.info = MCMCIteratorInfo(chain.info, tuned = true)
    nothing
end


text:
struct LinearPolicyWithBasis{TP, TB} <: AbstractPolicy where {TP, TB}
    π::TP
    ϕ::TB
end

function params(π::LinearPolicyWithBasis)
    return params(π.π)
end

function paramsvec(π::LinearPolicyWithBasis)
    return paramsvec(π.π)
end

function paramsfromvec!(π::LinearPolicyWithBasis, θ)
    paramsfromvec!(π.π, θ)
    return π
end

function rrule(::typeof(paramsfromvec!), π::LinearPolicyWithBasis, θ)
    _, pback = rrule(paramsfromvec!, π.π, θ)
    function paramsfromvec_linearwithbasis!(ȳ)
        dθ = pback(ȳ.π)[3]
        return NoTangent(), ȳ, dθ
    end
    return π, paramsfromvec_linearwithbasis!
end

function (π::LinearPolicyWithBasis)(s)
    feats = π.ϕ(s)
    return π.π(feats)
end

function logpdf(π::LinearPolicyWithBasis, s, a)
    feats = π.ϕ(s)
    return logpdf(π.π, feats, a)
end

function grad_logpdf!(ψ, π::LinearPolicyWithBasis, s, a)
    feats = π.ϕ(s)
    return grad_logpdf!(ψ, π.π, feats, a)
end

function grad_logpdf(π::LinearPolicyWithBasis, s, a)
    feats = π.ϕ(s)
    logp, ψ = grad_logpdf(π.π, feats, a)
    return logp, ψ
end

function sample_with_trace!(ψ, action, π::LinearPolicyWithBasis, s)
    feats = π.ϕ(s)
    return sample_with_trace!(ψ, action, π.π, feats)
end

function sample_with_trace(π::LinearPolicyWithBasis, s)
    feats = π.ϕ(s)
    return sample_with_trace(π.π, feats)
end


text:
"""
`pbuild_coarse_transition_matrix`: Contruct a transition matrix amongst the
bins in parallel.  It assumed that a pool of workers has already been
contructed.  Returns a sparse matrix.

### Arguments
* `mutation!` - an in place mutation function
* `bin_id` - bin identification function
* `x0_vals` - an array of starting values
* `bin0_vals` - an array of the bins corresponding to `x0_vals`
* `n_bins` - total number of bins
* `n_trials` - number of independent trials for each x0 starting value
"""

function pbuild_coarse_transition_matrix(mutation!, bin_id, x0_vals, n_bins, n_trials)
   n_x0 = length(x0_vals);
   row_vals = SharedArray{Int}(n_x0*n_trials);
   col_vals = SharedArray{Int}(n_x0*n_trials);
   X = similar(x0_vals[1]);

   @sync @distributed for k in 1:n_trials
      for l in 1:n_x0
         X .= deepcopy(x0_vals[l]);
         row_vals[n_x0*(k-1) + l] = bin_id(X);
         mutation!(X);
         col_vals[n_x0*(k-1) + l] = bin_id(X);
      end
   end

   K = sparse(row_vals,col_vals,ones(size(row_vals)),n_bins, n_bins);
   for i in 1:n_bins
      rowZ = sum(K[i,:]);
      @. K[i,:] = K[i,:]/rowZ;
   end
   return K
end


text:
module Day12

export day12

mutable struct Position{T}
    long::T # East-West
    lat::T  # North-South
end

function parseinput(T, filename)
    input = Vector{Pair{Char, T}}()
    for line in readlines(filename)
        push!(input, line[1] => parse(T, line[2:end]))
    end
    input
end

turn(facing, angle) = mod(facing + div(angle, 90), 4)

function forward!(position, facing, distance)
    if facing == 0
        position.lat += distance
    elseif facing == 1
        position.long += distance
    elseif facing == 2
        position.lat -= distance
    else
        position.long -= distance
    end
    nothing
end

function day12(::Val{:a}, filename)
    dir2num = Dict{Char, Int}(['N' => 0, 'E' => 1, 'S' => 2, 'W' => 3])
    input = parseinput(Int, filename)
    position = Position{Int}(0,0)
    facing = dir2num['E']
    for instruction in input
        if first(instruction) == 'N'
            position.lat += last(instruction)
        elseif first(instruction) == 'S'
            position.lat -= last(instruction)
        elseif first(instruction) == 'E'
            position.long += last(instruction)
        elseif first(instruction) == 'W'
            position.long -= last(instruction)
        elseif first(instruction) == 'R'
            facing = turn(facing, last(instruction))
        elseif first(instruction) == 'L'
            facing = turn(facing, -last(instruction))
        elseif first(instruction) == 'F'
            forward!(position, facing, last(instruction))
        end
    end
    abs(position.long) + abs(position.lat)
end

function turn_waypoint!(waypoint, bycos, bysin)
    x = waypoint.long
    y = waypoint.lat
    waypoint.long = x * bycos - y * bysin
    waypoint.lat = y * bycos + x * bysin
    nothing
end

function move!(position, waypoint, times)
    position.long += times * waypoint.long
    position.lat += times * waypoint.lat
    nothing
end

function day12(::Val{:b}, filename)
    input = parseinput(Float64, filename)
    waypoint = Position(10.0,1.0)
    position = Position(0.0,0.0)
    for instruction in input
        if first(instruction) == 'N'
            waypoint.lat += last(instruction)
        elseif first(instruction) == 'S'
            waypoint.lat -= last(instruction)
        elseif first(instruction) == 'E'
            waypoint.long += last(instruction)
        elseif first(instruction) == 'W'
            waypoint.long -= last(instruction)
        elseif first(instruction) == 'R'
            facing = turn_waypoint!(waypoint, cos(last(instruction)*pi/180.0), -sin(last(instruction)*pi/180.0))
        elseif first(instruction) == 'L'
            facing = facing = turn_waypoint!(waypoint, cos(last(instruction)*pi/180.0), sin(last(instruction)*pi/180.0))
        elseif first(instruction) == 'F'
            move!(position, waypoint, last(instruction))
        end
    end
    abs(position.long) + abs(position.lat)
end

end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Rays

using FromFile: @from
using StaticArrays: SVector, SMatrix
using LinearAlgebra: normalize
@from "Transforms.jl" import Transforms: Transform, rotate

# export AbstractRay, Ray, origin, direction, closestpointonray

abstract type AbstractRay{T <: Real, N} end

"""
    Ray{T,N} <: AbstractRay{T,N}

Purely geometric ray, defined as `origin + alpha * direction`.

```julia
Ray(origin::SVector{N,T}, direction::SVector{N,T})
```

Has the following accessor methods:
```julia
direction(ray::Ray{T,N}) -> SVector{N,T}
origin(ray::Ray{T,N}) -> SVector{N,T}
```
"""
struct Ray{T, N} <: AbstractRay{T, N}
    origin::SVector{N, T}
    direction::SVector{N, T}

    function Ray(
        origin::SVector{N, T},
        direction::SVector{N, T},
    ) where {T <: Real, N}
        return new{T, N}(origin, normalize(direction))
    end

    function Ray(
        origin::AbstractArray{T, 1},
        direction::AbstractArray{T, 1},
    ) where {T <: Real}
        @assert length(origin) == length(direction)
        N = length(origin)
        return new{T, N}(
            SVector{N, T}(origin),
            normalize(SVector{N, T}(direction)),
        )
    end
end

direction(ray::Ray{T, N}) where {T <: Real, N} = ray.direction
origin(ray::Ray{T, N}) where {T <: Real, N} = ray.origin

function Base.print(io::IO, a::Ray{T, N}) where {T, N}
    println(io, "$(rpad("Origin:", 20)) $(origin(a))")
    return println(io, "$(rpad("Direction:", 20)) $(direction(a))")
end

"""
    closestpointonray(r::Ray{T,N}, point::SVector{N,T}) -> SVector{T,N

Returns the point on the ray closest to point.
"""
function closestpointonray(
    r::AbstractRay{T, N},
    point::SVector{N, T},
) where {T, N}
    #find t that gives the smallest distance between the ray and the point
    #ray = p0 + tr̂, distance = ||line - point||
    o = origin(r)
    r̂ = direction(r)
    t = dot(point .- o, r̂)

    if t >= 0.0
        return o .+ t .* r̂
    else
        return o
    end
end

"""
    distance(r::Ray{T,N}, point::SVector{N,T}) -> Union{Nothing,T}

Returns distance to the position on the ray closest to point. If t < 0 returns nothing.
"""
function distance(r::AbstractRay{T, N}, point::SVector{N, T}) where {T, N}
    #find t that gives the smallest distance between the ray and the point
    #ray = p0 + tr̂, distance = ||line - point||
    # this has redundant code with closestpointonray but it is faster to do this way than to call closestpointonray and then comput norm(result .- point)

    o = origin(r)
    r̂ = direction(r)
    t = dot(point .- o, r̂)

    if t >= 0.0
        temp = o .+ t .* r̂
        return norm(temp .- point)
    else
        return nothing
    end
end

"""
    α(ray::AbstractRay{T,N}, point::SVector{N,T}) -> T

Computes the alpha corresponding to the closest position on the ray to point
"""
function α(ray::AbstractRay{T, N}, point::SVector{N, T}) where {T <: Real, N}
    return dot(direction(ray), (point .- origin(ray)))
end

"""
Apply a Transform to a Ray object
"""
function Base.:*(a::Transform{T}, r::Ray{T, 3})::Ray{T, 3} where {T}
    return Ray(a * origin(r), rotate(a, direction(r)))
end

end # module


text:
# nuclear Norm (times a constant)

"""
  NuclearNorm(λ::Real=1.0)

Returns the function `λ∑σ_i(X)`, where `σ_i(X)` is i-th singular value of matrix X.
"""

immutable NuclearNorm{R <: Real} <: ProximableConvex
  lambda::R
  function NuclearNorm(lambda::R)
    if lambda < 0
      error("parameter λ must be nonnegative")
    else
      new(lambda)
    end
  end
end

NuclearNorm{R <: Real}(lambda::R=1.0) = NuclearNorm{R}(lambda)

function (f::NuclearNorm){T <: RealOrComplex}(X::AbstractArray{T,2})
  U, S, V = svd(X);
  return f.lambda * sum(S);
end

function prox!{T <: RealOrComplex}(f::NuclearNorm, X::AbstractArray{T,2}, Y::AbstractArray{T,2}, gamma::Real=1.0)
  U, S, V = svd(X)

  for i in eachindex(S)
    S[i] = max(0, S[i] - f.lambda*gamma);
  end

  Y[:] = U * diagm(S) * V'
  return f.lambda * sum(S);
end

fun_name(f::NuclearNorm) = "nuclear norm"
fun_dom(f::NuclearNorm) = "AbstractArray{Real,2}, AbstractArray{Complex,2}"
fun_expr(f::NuclearNorm) = "X ↦ λ∑σ_i(X)"
fun_params(f::NuclearNorm) = "λ = $(f.lambda)"

function prox_naive{T <: RealOrComplex}(f::NuclearNorm, X::AbstractArray{T,2}, gamma::Real=1.0)
  U,S,V = svd(X)
  ftemp = NormL1(1.0)
  S_γ, fS_γ =  prox(ftemp,S,f.lambda*gamma)
  Y = U * diagm(S_γ) * V'
  return Y, f.lambda * sum(S_γ)
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8]);F = [7625*x1^4+3080*x1^2*x2^2+4760*x1^2*x2*x6+200*x1^2*x2*x7-280*x1^2*x5*x8+600*x1^2*x6*x8+5929*x2^4+770*x2^3*x7+9*x2^2*x4^2-1078*x2^2*x5*x8+784*x2^2*x6^2+2310*x2^2*x6*x8+25*x2^2*x7^2+36*x2*x3*x4*x7+360*x2*x4^2*x7+102*x2*x4*x7*x8-70*x2*x5*x7*x8+150*x2*x6*x7*x8+36*x3^2*x7^2+720*x3*x4*x7^2+204*x3*x7^2*x8+3600*x4^2*x7^2+2040*x4*x7^2*x8+49*x5^2*x8^2-210*x5*x6*x8^2+225*x6^2*x8^2+289*x7^2*x8^2-2550*x1^2*x3+9690*x1^2*x5-15810*x1^2*x6-13940*x1^2*x7+2240*x1^2*x8+8624*x2^2*x8-840*x2*x3*x6-306*x2*x4^2+234*x2*x4*x5+3192*x2*x5*x6-5208*x2*x6^2-4592*x2*x6*x7+560*x2*x7*x8-612*x3*x4*x7+468*x3*x5*x7-6120*x4^2*x7+4680*x4*x5*x7-1734*x4*x7*x8+1326*x5*x7*x8-784*x5*x8^2+1680*x6*x8^2+225*x3^2-1710*x3*x5+2790*x3*x6+2460*x3*x7+2601*x4^2-3978*x4*x5+4770*x5^2-10602*x5*x6-9348*x5*x7+8649*x6^2+15252*x6*x7+6724*x7^2+3136*x8^2, 6160*x1^2*x2+4760*x1^2*x6+200*x1^2*x7+23716*x2^3+2310*x2^2*x7+18*x2*x4^2-2156*x2*x5*x8+1568*x2*x6^2+4620*x2*x6*x8+50*x2*x7^2+36*x3*x4*x7+360*x4^2*x7+102*x4*x7*x8-70*x5*x7*x8+150*x6*x7*x8+17248*x2*x8-840*x3*x6-306*x4^2+234*x4*x5+3192*x5*x6-5208*x6^2-4592*x6*x7+560*x7*x8, 36*x2*x4*x7+72*x3*x7^2+720*x4*x7^2+204*x7^2*x8-2550*x1^2-840*x2*x6-612*x4*x7+468*x5*x7+450*x3-1710*x5+2790*x6+2460*x7, 18*x2^2*x4+36*x2*x3*x7+720*x2*x4*x7+102*x2*x7*x8+720*x3*x7^2+7200*x4*x7^2+2040*x7^2*x8-612*x2*x4+234*x2*x5-612*x3*x7-12240*x4*x7+4680*x5*x7-1734*x7*x8+5202*x4-3978*x5, -280*x1^2*x8-1078*x2^2*x8-70*x2*x7*x8+98*x5*x8^2-210*x6*x8^2+9690*x1^2+234*x2*x4+3192*x2*x6+468*x3*x7+4680*x4*x7+1326*x7*x8-784*x8^2-1710*x3-3978*x4+9540*x5-10602*x6-9348*x7, 4760*x1^2*x2+600*x1^2*x8+1568*x2^2*x6+2310*x2^2*x8+150*x2*x7*x8-210*x5*x8^2+450*x6*x8^2-15810*x1^2-840*x2*x3+3192*x2*x5-10416*x2*x6-4592*x2*x7+1680*x8^2+2790*x3-10602*x5+17298*x6+15252*x7, 200*x1^2*x2+770*x2^3+50*x2^2*x7+36*x2*x3*x4+360*x2*x4^2+102*x2*x4*x8-70*x2*x5*x8+150*x2*x6*x8+72*x3^2*x7+1440*x3*x4*x7+408*x3*x7*x8+7200*x4^2*x7+4080*x4*x7*x8+578*x7*x8^2-13940*x1^2-4592*x2*x6+560*x2*x8-612*x3*x4+468*x3*x5-6120*x4^2+4680*x4*x5-1734*x4*x8+1326*x5*x8+2460*x3-9348*x5+15252*x6+13448*x7, -280*x1^2*x5+600*x1^2*x6-1078*x2^2*x5+2310*x2^2*x6+102*x2*x4*x7-70*x2*x5*x7+150*x2*x6*x7+204*x3*x7^2+2040*x4*x7^2+98*x5^2*x8-420*x5*x6*x8+450*x6^2*x8+578*x7^2*x8+2240*x1^2+8624*x2^2+560*x2*x7-1734*x4*x7+1326*x5*x7-1568*x5*x8+3360*x6*x8+6272*x8];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
name(as::AbstractFeature) = as.name
name(as::AbstractFeature) = as.name

const _ranks = (
    domain     = 0,
    kingdom    = 1,
    phylum     = 2,
    class      = 3,
    order      = 4,
    family     = 5,
    genus      = 6,
    species    = 7,
    subspecies = 8,
    strain     = 9
)

const _shortranks = (
    d = :domain,
    k = :kingdom,
    p = :phylum,
    c = :class,
    o = :order,
    f = :family,
    g = :genus,
    s = :species,
    t = :subspecies,
    u = missing
)

"""
    Taxon(name::String, rank::Union{Missing, Symbol, Int}) <: AbstractFeature
    Taxon(name::String)

Microbial taxon with a name and a rank that can be one of

0. `:domain`
1. `:kingom`
2. `:phylum`
3. `:class`
4. `:order`
5. `:faamily`
6. `:genus`
7. `:species`
8. `:subspecies`
9. `:strain`

or `missing`. Contructors can also use numbers 0-9, or pass a string alone
(in which case the `taxon` will be stored as `missing`).

See also [`taxon`](@ref Microbiome.taxon).
"""
struct Taxon <: AbstractFeature
    name::String
    rank::Union{Missing, Symbol}

    Taxon(s::AbstractString, ::Missing) = new(s, missing)
    Taxon(s::AbstractString, rank::Symbol) = in(rank, keys(_ranks)) ?
                                                        new(s, rank)  :
                                                        error("Invalid rank $rank, must be one of $(keys(_ranks))")
end

Taxon(n::AbstractString, rank::Int) = 0 <= rank <= 9 ?
                                            Taxon(n, keys(_ranks)[rank+1]) :
                                            error("Invalid rank $rank, must be one of $_ranks")
Taxon(n::AbstractString) = Taxon(n, missing)

function Base.String(t::Taxon)
    if hasrank(t)
        return taxrank(t) == :strain ? string("t__", name(t)) : string(first(string(taxrank(t))), "__", name(t))
    else
        return string("u__", name(t))
    end
end

"""
    taxon(::AbstractString)

Return a [`Taxon`](@ref) from a string representation.
If the string contains taxonomic rank information in the form
`"x__Thename"` where `x` is the first letter of the rank,
this information will be used.

## Examples

```julia-repl
julia> taxon("Unknown")
Taxon("Unknown", missing)

julia> taxon("s__Prevotella_copri")
Taxon("Prevotella_copri", :species)
```
"""
function taxon(n::AbstractString)
    m = match(r"^([dkpcofgstu])__(.+)", n)
    isnothing(m) && return Taxon(n)
    return Taxon(string(m.captures[2]), _shortranks[Symbol(m.captures[1])])
end

"""
    taxrank(t::Union{Taxon, missing})

Get the `rank` field from a [`Taxon`](@ref) `t`.
Returns `missing` if the rank is not set.
"""
taxrank(t::Taxon) = t.rank
taxrank(::Missing) = missing

"""
    hasrank(t::Taxon)::Bool

Boolean function that returns `true` if the `rank`
field in [`Taxon`](@ref) `t` is not `missing`,
or `false` if it is `missing`
"""
hasrank(t::Taxon) = !ismissing(taxrank(t))

"""
    GeneFunction(name::String, taxon::Union{Taxon, String, Missing}) <: AbstractFeature
    GeneFunction(name::String)

Microbial gene function object with optional stratification (`taxon`).
"""
struct GeneFunction <: AbstractFeature
    name::String
    taxon::Union{Missing, Taxon}
end

GeneFunction(n::AbstractString) = GeneFunction(n, missing)
GeneFunction(n::AbstractString, t::AbstractString) = GeneFunction(n, taxon(t))

"""
    taxon(gf::GeneFunction)

Get the `taxon` field from a [`GeneFunction`](@ref), `gf`.
Returns `missing` if the taxon is not set.
"""
taxon(gf::GeneFunction) = gf.taxon

"""
    hastaxon(gf::GeneFunction)::Bool

Boolean function that returns `true` if the `taxon`
field in a [`GeneFunction`](@ref) `gf` is not `missing`,
or `false` if it is `missing`
"""
hastaxon(gf::GeneFunction) = !ismissing(taxon(gf))

"""
    taxrank(gf::GeneFunction)

Get the `rank` field from the `taxon` field of a [`GeneFunction`](@ref) `gf`
if it has one.
Returns `missing` if the `taxon` or `rank` is not set.
"""
taxrank(gf::GeneFunction) = taxrank(taxon(gf))

"""
    hasrank(gf::GeneFunction)::Bool

Boolean function that returns:

- `true` if `gf` has a [`Taxon`](@ref) with a non-missing `rank` field,
- `false` if there's no `Taxon`, or
- `false` if the `Taxon` has no `rank`
"""
hasrank(gf::GeneFunction) = hastaxon(gf) && !ismissing(taxrank(gf))

Base.String(gf::GeneFunction) = hastaxon(gf) ? string(name(gf), '|', String(taxon(gf))) : name(gf)

Base.:(==)(g1::GeneFunction, g2::GeneFunction) = String(g1) == String(g2)

@testset "GeneFunction Equality" begin
    @test GeneFunction("test") == GeneFunction("test")
    @test GeneFunction("test", Taxon("taxon")) == GeneFunction("test", Taxon("taxon"))
    @test GeneFunction("test") != GeneFunction("test", Taxon("taxon"))
end

"""
    genefunction(n::AbstractString)

Make a [`GeneFunction`](@ref) from a string,
Converting anything after an initial `|` as a [`Taxon`](@ref).
"""
function genefunction(n::AbstractString)
    if contains(n, '|')
        spl = split(n, '|')
        return GeneFunction(string(spl[1]), taxon(join(spl[2:end], '|')))
    else
        return GeneFunction(n)
    end
end

Base.string(f::AbstractFeature) = String(f)

"""
    Metabolite(name::String, commonname::Union{Missing, String}, mz::Union{Missing, Float64}, rt::Union{Missing, Float64}) <: AbstractFeature
    Metabolite(name::String)

Represents a small-molecule metabolite coming from an LCMS.
The fields are

- `name`: required, this should be a unique identifier
- `commonname`: might refer to a chemical name like "proprionate"
- `mz`: The mass/charge ratio
- `rt`: The retention time
"""
struct Metabolite <: AbstractFeature
    name::String
    commonname::Union{Missing, String}
    mz::Union{Missing, Float64}
    rt::Union{Missing, Float64}
end

Metabolite(n::AbstractString) = Metabolite(n, missing, missing, missing)

name(m::Metabolite) = m.name

"""
    commonname(m::Metabolite)

Accessor function for the `commonname` field of a [`Metabolite`](@ref).
"""
commonname(m::Metabolite) = m.commonname

"""
    masscharge(m::Metabolite)

Accessor function for the `mz` field of a [`Metabolite`](@ref).
"""
masscharge(m::Metabolite) = m.mz

"""
    retentiontime(m::Metabolite)

Accessor function for the `rt` field of a [`Metabolite`](@ref).
"""
retentiontime(m::Metabolite) = m.rt

@testset "Metabolites" begin
    m1 = Metabolite("name", "common", 1., 1.)
    @test name(m1) == "name"
    @test commonname(m1) == "common"
    @test masscharge(m1) == 1
    @test retentiontime(m1) == 1
    m2 = Metabolite("name")
    @test name(m2) == "name"
    @test ismissing(commonname(m2))
    @test ismissing(masscharge(m2))
    @test ismissing(masscharge(m2))

end


text:
#=
Authors:
Darius Russell Kish

Created: April 14, 2020
Modified: April 14, 2020

This implementation is based strongly in the Disjoint Set implementation
in JuliaCollections' DataStructures which is under active development. It
represents proper implementation of data structures for Julia and much of the
typing and integration into the language is taken from their implementation,
however the function implementations are my own.
https://github.com/JuliaCollections/DataStructures.jl/blob/master/src/disjoint_set.jl

Implementations for collections of arbitrary inputs are *not* provided here
and it is assumed this low-level data structure will be wrapped for such
functionality. This avoids issues of key-value association time complexity for
conversion between elements and their index value in analyses of this data
structure.
=#

#=
Notation T <: Integer :: T is a subtype of Integer

An example type tree:
abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end
=#
mutable struct QuickUnion{T<:Integer}
    ids::Vector{T}
    nGroups::T
end

# Construction type 1
# Notation: collect(Base.OneTo(n)) == list(range(n)) in Python
QuickUnion(n::T) where {T <: Integer} =
    QuickUnion{T}(collect(Base.OneTo(n)), n)

# Construction type 2
QuickUnion{T}(n::Integer) where {T <: Integer} =
    QuickUnion{T}(collect(Base.OneTo(n)), n)

# The length or "size" of a UnionFind is the length of either of its vectors
length(uf::QuickUnion) = Base.length(uf.ids)

"""
    num_groups(uf::QuickUnion)
Get a number of groups.
"""
num_groups(uf::QuickUnion) = uf.nGroups

"""
    root(uf::QuickUnion{T}, x::T)
Finds the root node of node x
"""
function root(uf::QuickUnion{T}, x::T) where {T <: Integer}
    root = uf.ids[x]
    # extracting one iteration of the loop
    # allows ensuring inbounds for the rest
    @inbounds while root != uf.ids[root]
        root = uf.ids[root]
    end
    return root
end

"""
    find(uf::QuickUnion{T}, x::T, y::T)
Determine if x and y belong to the same group
"""
function find(uf::QuickUnion{T}, x::T, y::T) where {T <: Integer}
    return root(uf, x) == root(uf, y)
end

"""
    unite(uf::QuickUnion{T}, x::T, y::T)
Merges two groups together
"""
function unite!(uf::QuickUnion{T}, x::T, y::T) where {T <: Integer}
    pid = root(uf, x)
    tid = root(uf, y)

    uf.ids[pid] = tid

    uf.nGroups -= 1
end


text:
abstract type DistributionsPreMetric <: Distances.PreMetric end
abstract type DistributionsSemiMetric <: Distances.SemiMetric end
abstract type DistributionsMetric <: Distances.Metric end

const DistributionsDistance = Union{
    DistributionsPreMetric,DistributionsSemiMetric,DistributionsMetric
}


text:
using Quantica: nbands, nvertices, nedges, nsimplices, Subspace
using Arpack

@testset "basic bandstructures" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1))
    b = bandstructure(h, subticks = 13, showprogress = false)
    @test nbands(b)  == 1

    h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(2)
    b = bandstructure(h, cuboid(0.99 .* (-π, π), 0.99 .* (-π, π), subticks = 13), showprogress = false)
    @test nbands(b)  == 3

    h = LatticePresets.honeycomb() |>
        hamiltonian(onsite(0.5, sublats = :A) + onsite(-0.5, sublats = :B) +
                    hopping(-1, range = 1/√3))
    b = bandstructure(h, subticks = (13, 15), showprogress = false)
    @test nbands(b)  == 2

    h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2)
    b = bandstructure(h, subticks = (5, 7, 5), showprogress = false)
    @test nbands(b)  == 1

    b = bandstructure(h, :Γ, :X; subticks = 4, showprogress = false)
    @test nbands(b)  == 1

    b = bandstructure(h, :Γ, :X, (0, π), :Z, :Γ; subticks = 4, showprogress = false)
    @test nbands(b)  == 1
    @test nvertices(b) == 73

    b = bandstructure(h, :Γ, :X, (0, π), :Z, :Γ; subticks = (4,5,6,7), showprogress = false)
    @test nbands(b) == 1
    @test nvertices(b) == 113

    # complex spectra
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(im) + hopping(-1)) |> unitcell(2)
    b = bandstructure(h, cuboid((-π, π), (-π, π), subticks = 7), showprogress = false)
    @test nbands(b)  == 1

    # spectrum sorting
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(10)
    b = bandstructure(h, :K, :K´, subticks = 7, method = ArpackPackage(sigma = 0.1im, nev = 12), showprogress = false)
    @test nbands(b) == 1

    # number of simplices contant across BZ (computing their degeneracy with the projs at each vertex j)
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1)) |> unitcell(3)
    nev = size(h, 1)
    b = bandstructure(h, splitbands = false)
    @test nbands(b) == 1
    band = bands(b, 1)
    @test all(j -> all(rng -> sum(i -> degeneracy(band.sbases[i][j]), rng) == nev, band.sptrs), 1:3)
end

@testset "functional bandstructures" begin
    hc = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, sublats = :A=>:B, plusadjoint = true))
    matrix = similarmatrix(hc, LinearAlgebraPackage())
    hf((x,)) = bloch!(matrix, hc, (x, -x))
    m = cuboid((0, 1))
    b = bandstructure(hf, m, showprogress = false)
    @test nbands(b) == 2
    @test nsimplices(b)  == 24
    @test b[(1,), around = 0] isa Subspace

    hc2 = LatticePresets.honeycomb() |> hamiltonian(hopping(-1))
    hp2 = parametric(hc2, @hopping!((t; s) -> s*t))
    matrix2 = similarmatrix(hc2, LinearAlgebraPackage())
    hf2((s, x)) = bloch!(matrix2, hp2(s = s), (x, x))
    m2 = cuboid((0, 1), (0, 1))
    b = bandstructure(hf2, m2, showprogress = false)
    @test nbands(b) == 1
    @test nsimplices(b)  == 576
    @test degeneracy.(b[(0,0), around = (0, 2)]) == [1, 1]
end

@testset "bandstructures lifts & transforms" begin
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(2) + hopping(-1, range = 1/√3))
    mesh1D = cuboid((0, 2π))
    b = bandstructure(h, mesh1D, mapping = φ -> (φ, -φ), transform = inv, showprogress = false)
    b´ = transform!(inv, bandstructure(h, mesh1D, mapping = φ -> (φ, -φ), showprogress = false))
    @test nbands(b)  == nbands(b´) == 1
    @test vertices(bands(b)[1]) == vertices(bands(b´)[1])
    h´ = unitcell(h)
    s1 = spectrum(h´, transform = inv)
    s2 = transform!(inv, spectrum(h´))
    @test s1.energies == s2.energies
    # no automatic mapping from 2D to 3D
    h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2)
    @test_throws DimensionMismatch bandstructure(h, cuboid((0, 2pi), (0, 2pi)), showprogress = false)
end

@testset "parametric bandstructures" begin
    ph = LatticePresets.linear() |> hamiltonian(onsite(0I) + hopping(-I), orbitals = Val(2)) |> unitcell(2) |>
         parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k = 2, p = [1,2])-> t - k*I .+ p'p))
    mesh2D = cuboid((0, 1), (0, 2π), subticks = 15)
    b = bandstructure(ph, mesh2D, mapping = (x, k) -> (x, (;k = k)), showprogress = false)
    @test nbands(b)  == 4
    b = bandstructure(ph, mesh2D, mapping = (x, k) -> ((x,), (;k = k)), showprogress = false)
    @test nbands(b)  == 4
    b = bandstructure(ph, mesh2D, mapping = (x, k) -> (SA[x], (;k = k)), showprogress = false)
    @test nbands(b)  == 4
    b = bandstructure(ph, mesh2D, mapping = (k, φ) -> (1, (;k = k, p = SA[1, φ])), showprogress = false)
    @test nbands(b)  == 1

    ph = LatticePresets.linear() |> hamiltonian(onsite(0I) + hopping(-I), orbitals = Val(2)) |> unitcell(2) |>
        unitcell |> parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k = 2, p = [1,2])-> t - k*I .+ p'p))
    b = bandstructure(ph, mesh2D, mapping = (k, φ) -> (;k = k, p = SA[1, φ]), showprogress = false)
    @test nbands(b)  == 1
    b = bandstructure(ph, mesh2D, mapping = (k, φ) -> ((;k = k, p = SA[1, φ]),), showprogress = false)
    @test nbands(b)  == 1
    @test_throws UndefKeywordError bandstructure(ph, mesh2D, mapping = (k, φ) -> ((;p = SA[1, φ]),), showprogress = false)
end

@testset "unflatten" begin
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = (Val(2), Val(1))) |> unitcell(2) |> unitcell
    sp = spectrum(h).states[:,1]
    sp´ = Quantica.unflatten_orbitals_or_reinterpret(sp, orbitalstructure(h))
    l = size(h, 1)
    @test length(sp) == 1.5 * l
    @test length(sp´) == l
    @test all(x -> iszero(last(x)), sp´[l+1:end])
    @test sp´ isa Vector
    @test sp´ !== sp

    h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = (Val(2), Val(1))) |> unitcell(2)
    b = bandstructure(h)
    psi = b[(0,0), around = 0]
    psi´ = flatten(psi)
    psi´´ = unflatten(psi´, orbitalstructure(psi))
    @test eltype(psi.basis) == eltype(psi´´.basis) == SVector{2, ComplexF64}
    @test eltype(psi´.basis) == ComplexF64
    @test psi.basis == psi´´.basis

    h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = Val(2)) |> unitcell(2) |> unitcell
    sp = spectrum(h).states[:,1]
    sp´ = Quantica.unflatten_orbitals_or_reinterpret(sp, orbitalstructure(h))
    l = size(h, 1)
    @test length(sp) == 2 * l
    @test length(sp´) == l
    @test sp´ isa Base.ReinterpretArray

    h = LatticePresets.honeycomb() |> hamiltonian(onsite(2I) + hopping(I, range = 1), orbitals = Val(2)) |> unitcell(2) |> unitcell
    sp = spectrum(h).states[:,1]
    sp´ = Quantica.unflatten_orbitals_or_reinterpret(sp, orbitalstructure(h))
    @test sp === sp
end

@testset "spectrum indexing" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(2) |> wrap
    s = spectrum(h)
    ϵv, ψm = s
    @test ϵv == first(s) == s.energies
    @test ψm == last(s) == s.states
    @test s[1] isa Subspace
    @test degeneracy(s[2]) == 3
    @test s[1:3] isa Vector{<:Subspace}
    @test s[[1,3]] isa Vector{<:Subspace}
    ϵ, ψs = s[1]
    @test ϵ isa Number
    @test ψs isa SubArray{<:Complex, 2}
end

@testset "bandstructure indexing" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |> unitcell(2)
    bs = bandstructure(h, subticks = 13)
    @test sum(degeneracy, bs[(1,2)]) == size(h,1)
    @test sum(degeneracy, bs[(0.2,0.3)]) == size(h,1)
    @test size(bs[(1,2), around = 0].basis, 1) == size(h, 2)
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1f0*I), orbitals = (Val(1),Val(2))) |> unitcell(2)
    bs = bandstructure(h, subticks = 13)
    @test sum(degeneracy, bs[(1,2)]) == size(h,1) * 1.5
    @test sum(degeneracy, bs[(0.2,0.3)]) == size(h,1) * 1.5
    @test bs[(1,2), around = 0] |> last |> eltype == SVector{2, ComplexF64}
    @test size(bs[(1,2), around = 0].basis, 1) == size(h, 2)
end

@testset "subspace flatten" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(I), orbitals = (Val(1), Val(2))) |> unitcell(2) |> unitcell
    sub = spectrum(h)[1]
    @test size(sub.basis) == (8,1) && size(flatten(sub).basis) == (12,1)
end

@testset "diagonalizer" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1))
    d = diagonalizer(h)
    @test first(d((0, 0))) ≈ [-3,3]
    h = wrap(h)
    d = diagonalizer(h)
    @test first(d()) ≈ [-3,3]
    @test first(d(())) ≈ [-3,3]
    @test all(d() .≈ Tuple(spectrum(h)))
end

@testset "bandstructure extrema" begin
    h = LP.honeycomb() |>
        hamiltonian(hopping(I) + onsite(0.1, sublats = :A) - onsite(0.1, sublats = :B)) |>
        unitcell(3) |>
        Quantica.wrap(1)
    b = bandstructure(h, subticks = 20)
    @test !isapprox(gap(b, 0; refinesteps = 0), 0.2)
    @test isapprox(gap(b, 0; refinesteps = 1), 0.2)
    @test isapprox(gap(b, 0.3; refinesteps = 1), 0)
    @test isapprox(gap(b, 4; refinesteps = 1), Inf)
    @test all(gapedge(b, 0, +; refinesteps = 1) .≈ (0, 0.1))
    @test all(gapedge(b, 0, -; refinesteps = 1) .≈ (0, -0.1))

    b = bandstructure(h, subticks = 20; method = ArpackPackage(nev = 6, sigma = 0.1*im))
    @test !isapprox(gap(b, 0; refinesteps = 0), 0.2)
    @test isapprox(gap(b, 0; refinesteps = 1), 0.2)
    @test isapprox(gap(b, 0.3; refinesteps = 1), 0)
    @test isapprox(gap(b, 4; refinesteps = 1), Inf)
    @test all(gapedge(b, 0, +; refinesteps = 1) .≈ (0, 0.1))
    @test all(gapedge(b, 0, -; refinesteps = 1) .≈ (0, -0.1))

    h = LP.honeycomb() |> hamiltonian(hopping(I)) |> Quantica.wrap(1)
    b = bandstructure(h, subticks = 20)
    length.(minima(b, refinesteps = 0)) == [2, 0]
    length.(maxima(b, refinesteps = 0)) == [0, 2]
    length.(minima(b, refinesteps = 1)) == [1, 0]
    length.(maxima(b, refinesteps = 1)) == [0, 1]
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
# -------------- Core functions for Supplementary Information --------------

# Packages required
using LinearAlgebra, Parameters, DifferentialEquations, RecursiveArrayTools
using Statistics: mean, std
using PyPlot

function local_max(ts)
    lmaxs = Float64[]
    for i in 2:(length(ts) - 3)
        if ts[i - 1] < ts[i] && ts[i] > ts[i + 1]
            push!(lmaxs, ts[i])
            # @show ts[i]
        end
    end

    if length(lmaxs) == 0
        push!(lmaxs, ts[length(ts) - 3])
    end

    return lmaxs
end

function local_min(ts)
    lmins = Float64[]
    for i in 2:(length(ts) - 3)
        if ts[i - 1] > ts[i] && ts[i] < ts[i + 1]
            push!(lmins, ts[i])
            #@show ts[i]
        end
    end
    if length(lmins) == 0
        push!(lmins, ts[length(ts) - 3])
    end

    return lmins
end

max_stars_repo_name: McCannLab/NDID

max_issues_repo_name: McCannLab/NDID

max_forks_repo_name: McCannLab/NDID


text:
push!(LOAD_PATH, string(pwd(),"/src"))

using ModifiedStokesSolver
import AnalyticDomains
using Base.Test
using PyPlot
import IterativeSolvers

@testset "MultiplyConnected" begin

    # Setup discretization
    R1 = 1.0
    R2 = 0.25
    alpha = 3.0
    panelorder = 16
    numpanels1 = 40
    numpanels2 = 8
    tol = 1e-8

    # Setup problem
    xsrc1 = [0.01, 0.02]
    xsrc2 = [1.0, 1.0]
    srand(1)
    qsrc = 1-2*rand(2)
    nsrc = 1-2*rand(2)
    fsrc = 1-2*rand(2)
    ubc(x) = stresslet(x-xsrc1, qsrc, nsrc, alpha) +
        stokeslet(x-xsrc1, fsrc, alpha) +
        stresslet(x-xsrc2, qsrc, nsrc, alpha) +
        stokeslet(x-xsrc2, fsrc, alpha)

    # Discretize
    curve1 = AnalyticDomains.starfish(;amplitude = 0.1, radius = R1)
    curve2 = AnalyticDomains.starfish(;amplitude = 0.0, radius = R2, interior=true)
    dcurve = discretize([curve1, curve2], [numpanels1, numpanels2], panelorder)
    N = dcurve.numpoints

    # Right hand side
    rhs = zeros(2*N)
    for i=1:dcurve.numpoints
        rhs[2(i-1)+(1:2)] = ubc(dcurve.points[:,i])
    end

    # Dense matrix
    println("* Matrix assembly and solve")
    @time LHS = system_matrix(dcurve, alpha)
    @time sol = LHS\rhs
    density = reshape(sol, 2, N)

    # Interior test point
    zt = [(R1+R2)/2, 0.01]
    ref = ubc(zt)
    dlp = doublelayer(dcurve, density, zt, alpha)
    relerr = maximum(abs.(ref.-dlp)) / maximum(abs.(ref))
    println("* Interior eval")
    @show relerr
    @test relerr < 1e-15

    # Field test
    Ngrid = 50
    xmax = maximum(dcurve.points[1,:])
    xmin = minimum(dcurve.points[1,:])
    ymax = maximum(dcurve.points[2,:])
    ymin = minimum(dcurve.points[2,:])
    x = linspace(xmin, xmax, Ngrid)
    y = linspace(ymin, ymax, Ngrid)
    X, Y = ndgrid(x, y)
    zt = [vec(X) vec(Y)]'
    interior = interior_points(dcurve, zt)
    numeval = sum(interior)
    zt_interior = zt[:, interior]
    @time u = doublelayer_fast(dcurve, density, zt_interior, alpha, specquad=true)
    uref = Array{Float64}(2, numeval)
    for i=1:numeval
        uref[:,i] = ubc(zt_interior[:, i])
    end
    unorm = norm(vec(uref), Inf)
    # Error
    err = u .- uref
    maxerr = Array{Float64}(numeval)
    for i=1:numeval
        maxerr[i] = norm(err[:, i], Inf)
    end

    max_relerr_grid = norm(vec(maxerr), Inf) / unorm
    @show max_relerr_grid
    @test max_relerr_grid < 1e-13

    # Plot
    E = zeros(X)
    E[interior] = maxerr
    clf()
    pcolor(X, Y, log10.(E))
    colorbar()
end


text:
import Base: ==, ≈, copy
import Random.rand

export LazySet,
       ρ, support_function,
       σ, support_vector,
       dim,
       norm,
       radius,
       diameter,
       an_element,
       isbounded, isbounded_unit_dimensions,
       neutral,
       absorbing,
       tosimplehrep,
       isuniversal,
       translate

"""
    LazySet{N}

Abstract type for convex sets, i.e., sets characterized by a (possibly infinite)
intersection of halfspaces, or equivalently, sets ``S`` such that for any two
elements ``x, y ∈ S`` and ``0 ≤ λ ≤ 1`` it holds that ``λ·x + (1-λ)·y ∈ S``.

### Notes

`LazySet` types should be parameterized with a type `N`, typically `N<:Real`,
for using different numeric types.

Every concrete `LazySet` must define the following functions:
- `σ(d::AbstractVector{N}, S::LazySet{N}) where {N<:Real}` -- the support vector
    of `S` in a given direction `d`; note that the numeric type `N` of `d` and
    `S` must be identical; for some set types `N` may be more restrictive than
    `Real`
- `dim(S::LazySet)::Int` -- the ambient dimension of `S`

The subtypes of `LazySet` (including abstract interfaces):

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(LazySet, false)
17-element Array{Any,1}:
 AbstractCentrallySymmetric
 AbstractPolyhedron
 CacheMinkowskiSum
 CartesianProduct
 CartesianProductArray
 ConvexHull
 ConvexHullArray
 EmptySet
 ExponentialMap
 ExponentialProjectionMap
 Intersection
 IntersectionArray
 LinearMap
 MinkowskiSum
 MinkowskiSumArray
 ResetMap
 Translation
```

If we only consider *concrete* subtypes, then:

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(LazySet, true)
37-element Array{Type,1}:
 Ball1
 Ball2
 BallInf
 Ballp
 CacheMinkowskiSum
 CartesianProduct
 CartesianProductArray
 ConvexHull
 ConvexHullArray
 Ellipsoid
 EmptySet
 ExponentialMap
 ExponentialProjectionMap
 HPolygon
 HPolygonOpt
 HPolyhedron
 HPolytope
 HalfSpace
 Hyperplane
 Hyperrectangle
 Intersection
 IntersectionArray
 Interval
 Line
 LineSegment
 LinearMap
 MinkowskiSum
 MinkowskiSumArray
 ResetMap
 Singleton
 SymmetricIntervalHull
 Translation
 Universe
 VPolygon
 VPolytope
 ZeroSet
 Zonotope
 ```
"""
abstract type LazySet{N} end

# --- common LazySet functions ---

"""
    ρ(d::AbstractVector{N}, S::LazySet{N})::N where {N<:Real}

Evaluate the support function of a set in a given direction.

### Input

- `d` -- direction
- `S` -- convex set

### Output

The support function of the set `S` for the direction `d`.

### Notes

The numeric type of the direction and the set must be identical.
"""
function ρ(d::AbstractVector{N}, S::LazySet{N})::N where {N<:Real}
    return dot(d, σ(d, S))
end

"""
    support_function

Alias for the support function ρ.
"""
const support_function = ρ

"""
    σ

Function to compute the support vector σ.
"""
function σ end

"""
    support_vector

Alias for the support vector σ.
"""
const support_vector = σ

"""
    isbounded(S::LazySet)::Bool

Determine whether a set is bounded.

### Input

- `S` -- set

### Output

`true` iff the set is bounded.

### Algorithm

We check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(S::LazySet)::Bool
    return isbounded_unit_dimensions(S)
end

"""
    isbounded_unit_dimensions(S::LazySet{N})::Bool where {N<:Real}

Determine whether a set is bounded in each unit dimension.

### Input

- `S` -- set

### Output

`true` iff the set is bounded in each unit dimension.

### Algorithm

This function performs ``2n`` support function checks, where ``n`` is the
ambient dimension of `S`.
"""
function isbounded_unit_dimensions(S::LazySet{N})::Bool where {N<:Real}
    n = dim(S)
    @inbounds for i in 1:n
        for o in [one(N), -one(N)]
            d = SingleEntryVector(i, n, o)
            if ρ(d, S) == N(Inf)
                return false
            end
        end
    end
    return true
end

"""
    norm(S::LazySet, [p]::Real=Inf)

Return the norm of a convex set.
It is the norm of the enclosing ball (of the given ``p``-norm) of minimal volume
that is centered in the origin.

### Input

- `S` -- convex set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the norm.
"""
function norm(S::LazySet, p::Real=Inf)
    if p == Inf
        return norm(Approximations.ballinf_approximation(S), p)
    else
        error("the norm for this value of p=$p is not implemented")
    end
end

"""
    radius(S::LazySet, [p]::Real=Inf)

Return the radius of a convex set.
It is the radius of the enclosing ball (of the given ``p``-norm) of minimal
volume with the same center.

### Input

- `S` -- convex set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the radius.
"""
function radius(S::LazySet, p::Real=Inf)
    if p == Inf
        return radius(Approximations.ballinf_approximation(S)::BallInf, p)
    else
        error("the radius for this value of p=$p is not implemented")
    end
end

"""
    diameter(S::LazySet, [p]::Real=Inf)

Return the diameter of a convex set.
It is the maximum distance between any two elements of the set, or,
equivalently, the diameter of the enclosing ball (of the given ``p``-norm) of
minimal volume with the same center.

### Input

- `S` -- convex set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the diameter.
"""
function diameter(S::LazySet, p::Real=Inf)
    return radius(S, p) * 2
end

"""
    an_element(S::LazySet{N}) where {N<:Real}

Return some element of a convex set.

### Input

- `S` -- convex set

### Output

An element of a convex set.
"""
function an_element(S::LazySet{N}) where {N<:Real}
    return σ(sparsevec([1], [one(N)], dim(S)), S)
end

"""
    ==(X::LazySet, Y::LazySet)

Return whether two LazySets of the same type are exactly equal.

### Input

- `X` -- any `LazySet`
- `Y` -- another `LazySet` of the same type as `X`

### Output

- `true` iff `X` is equal to `Y`.

### Notes

The check is purely syntactic and the sets need to have the same base type.
For instance, `X::VPolytope == Y::HPolytope` returns `false` even if `X` and `Y`
represent the same polytope.
However `X::HPolytope{Int64} == Y::HPolytope{Float64}` is a valid comparison.

### Algorithm

We recursively compare the fields of `X` and `Y` until a mismatch is found.

### Examples

```jldoctest
julia> HalfSpace([1], 1) == HalfSpace([1], 1)
true

julia> HalfSpace([1], 1) == HalfSpace([1.0], 1.0)
true

julia> Ball1([0.], 1.) == Ball2([0.], 1.)
false
```
"""
function ==(X::LazySet, Y::LazySet)
    # if the common supertype of X and Y is abstract, they cannot be compared
    if isabstracttype(promote_type(typeof(X), typeof(Y)))
        return false
    end

    for f in fieldnames(typeof(X))
        if getfield(X, f) != getfield(Y, f)
            return false
        end
    end

    return true
end
"""
    ≈(X::LazySet, Y::LazySet)

Return whether two LazySets of the same type are approximately equal.

### Input

- `X` -- any `LazySet`
- `Y` -- another `LazySet` of the same type as `X`

### Output

- `true` iff `X` is equal to `Y`.

### Notes

The check is purely syntactic and the sets need to have the same base type.
For instance, `X::VPolytope ≈ Y::HPolytope` returns `false` even if `X` and `Y`
represent the same polytope.
However `X::HPolytope{Int64} ≈ Y::HPolytope{Float64}` is a valid comparison.

### Algorithm

We recursively compare the fields of `X` and `Y` until a mismatch is found.

### Examples

```jldoctest
julia> HalfSpace([1], 1) ≈ HalfSpace([1], 1)
true

julia> HalfSpace([1], 1) ≈ HalfSpace([1.00000001], 0.99999999)
true

julia> HalfSpace([1], 1) ≈ HalfSpace([1.0], 1.0)
true

julia> Ball1([0.], 1.) ≈ Ball2([0.], 1.)
false
```
"""
function ≈(X::LazySet, Y::LazySet)
    # if the common supertype of X and Y is abstract, they cannot be compared
    if isabstracttype(promote_type(typeof(X), typeof(Y)))
        return false
    end

    for f in fieldnames(typeof(X))
        if !(getfield(X, f) ≈ getfield(Y, f))
            return false
        end
    end

    return true
end

# hook into random API
function rand(rng::AbstractRNG, ::SamplerType{T}) where T<:LazySet
    rand(T, rng=rng)
end

"""
    copy(S::LazySet)

Return a deep copy of the given set by copying its values recursively.

### Input

- `S` -- any `LazySet`

### Output

A copy of `S`.

### Notes

This function performs a `deepcopy` of each field in `S`, resulting in a
completely independent object. See the documentation of `?deepcopy` for further
details.
"""
copy(S::LazySet) = deepcopy(S)

"""
    tosimplehrep(S::LazySet)

Return the simple H-representation ``Ax ≤ b`` of a set from its list of linear
constraints.

### Input

- `S` -- set

### Output

The tuple `(A, b)` where `A` is the matrix of normal directions and `b` is the
vector of offsets.

### Notes

This function only works for sets that can be represented exactly by a finite
list of linear constraints.
This fallback implementation relies on `constraints_list(S)`.
"""
tosimplehrep(S::LazySet) = tosimplehrep(constraints_list(S))

"""
    isuniversal(X::LazySet{N}, [witness]::Bool=false
               )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a given convex set is universal, and otherwise optionally compute
a witness.

### Input

- `X`       -- convex set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X`` is universal
* If `witness` option is activated:
  * `(true, [])` iff ``X`` is universal
  * `(false, v)` iff ``X`` is not universal and ``v ∉ X``

### Notes

This is a naive fallback implementation.
"""
function isuniversal(X::LazySet{N}, witness::Bool=false
                    )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    if isbounded(X)
        result = false
    else
        error("cannot determine universality of the set")
    end

    if result
        return witness ? (true, N[]) : true
    elseif witness
        error("witness production is currently not supported")
    else
        return false
    end
end

"""
    plot_recipe(X::LazySet{N}, [ε]::N=N(PLOT_PRECISION)) where {N<:Real}

Convert a convex set to a pair `(x, y)` of points for plotting.

### Input

- `X` -- convex set
- `ε` -- (optional, default: `PLOT_PRECISION`) approximation error bound

### Output

A pair `(x, y)` of points that can be plotted.

### Notes

Plotting of unbounded sets is not implemented yet (see
[#576](https://github.com/JuliaReach/LazySets.jl/issues/576)).

### Algorithm

We first assert that `X` is bounded.

One-dimensional sets are converted to an `Interval`.
We do not support three-dimensional or higher-dimensional sets at the moment.

For two-dimensional sets, we first compute a polygonal overapproximation.
The second argument, `ε`, corresponds to the error in Hausdorff distance between
the overapproximating set and `X`.
The default value `PLOT_PRECISION` is chosen such that the unit ball in the
2-norm is approximated with reasonable accuracy.
On the other hand, if you only want to produce a fast box-overapproximation of
`X`, pass `ε=Inf`.
Finally, we use the plot recipe for polygons.
"""
function plot_recipe(X::LazySet{N}, ε::N=N(PLOT_PRECISION)) where {N<:Real}
    @assert dim(X) <= 2 "cannot plot a $(dim(X))-dimensional $(typeof(X))"
    @assert isbounded(X) "cannot plot an unbounded $(typeof(X))"

    if dim(X) == 1
        Y = convert(Interval, X)
    else
        Y = overapproximate(X, ε)
    end
    return plot_recipe(Y, ε)
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

"""
    REoptInputs(d::Dict)

Return REoptInputs(s) where s in `Scenario` defined in dict `d`.
"""

function REoptInputs(d::Dict)
    REoptInputs(Scenario(d))
end

"""
	run_reopt(m::JuMP.AbstractModel, fp::String)

Solve the model using the `Scenario` defined in JSON file stored at the file path `fp`.
"""
function run_reopt(m::JuMP.AbstractModel, fp::String)
	s = Scenario(JSON.parsefile(fp))
	run_reopt(m, REoptInputs(s))
end

"""
	run_reopt(m::JuMP.AbstractModel, d::Dict)

Solve the model using the `Scenario` defined in dict `d`.
"""
function run_reopt(m::JuMP.AbstractModel, d::Dict)
	s = Scenario(d)
	run_reopt(m, REoptInputs(s))
end

"""
	run_reopt(m::JuMP.AbstractModel, s::AbstractScenario)

Solve the model using a `Scenario` or `BAUScenario`.
"""
function run_reopt(m::JuMP.AbstractModel, s::AbstractScenario)
	run_reopt(m, REoptInputs(s))
end

"""
    run_reopt(t::Tuple{JuMP.AbstractModel, AbstractScenario})

Method for use with Threads when running BAU in parallel with optimal scenario.
"""
function run_reopt(t::Tuple{JuMP.AbstractModel, AbstractInputs})
	run_reopt(t[1], t[2]; organize_pvs=false)
    # must organize_pvs after adding proforma results
end

"""
    run_reopt(ms::AbstractArray{T, 1}, fp::String) where T <: JuMP.AbstractModel

Solve the `Scenario` and `BAUScenario` in parallel using the first two (empty) models in `ms` and inputs defined in the
JSON file at the filepath `fp`.
"""
function run_reopt(ms::AbstractArray{T, 1}, fp::String) where T <: JuMP.AbstractModel
    d = JSON.parsefile(fp)
    run_reopt(ms, d)
end

"""
    run_reopt(ms::AbstractArray{T, 1}, d::Dict) where T <: JuMP.AbstractModel

Solve the `Scenario` and `BAUScenario` in parallel using the first two (empty) models in `ms` and inputs from `d`.
"""
function run_reopt(ms::AbstractArray{T, 1}, d::Dict) where T <: JuMP.AbstractModel
    s = Scenario(d)
    if s.settings.off_grid_flag
        @warn "Only using first Model and not running BAU case because Settings.off_grid_flag == true. The BAU scenario is not applicable for off-grid microgrids."
	    results = run_reopt(ms[1], s)
        return results
    end

    run_reopt(ms, REoptInputs(s))
end

"""
    run_reopt(ms::AbstractArray{T, 1}, p::REoptInputs) where T <: JuMP.AbstractModel

Solve the `Scenario` and `BAUScenario` in parallel using the first two (empty) models in `ms` and inputs from `p`.
"""
function run_reopt(ms::AbstractArray{T, 1}, p::REoptInputs) where T <: JuMP.AbstractModel
    bau_inputs = BAUInputs(p)
    inputs = ((ms[1], bau_inputs), (ms[2], p))
    rs = Any[0, 0]
    Threads.@threads for i = 1:2
        rs[i] = run_reopt(inputs[i])
    end
    # TODO when a model is infeasible the JuMP.Model is returned from run_reopt (and not the results Dict)
    results_dict = combine_results(p, rs[1], rs[2], bau_inputs.s)
    results_dict["Financial"] = merge(results_dict["Financial"], proforma_results(p, results_dict))
    if !isempty(p.techs.pv)
        organize_multiple_pv_results(p, results_dict)
    end
    return results_dict
end

"""
	build_reopt!(m::JuMP.AbstractModel, fp::String)

Add variables and constraints for REopt model.
`fp` is used to load in JSON file to construct REoptInputs.
"""
function build_reopt!(m::JuMP.AbstractModel, fp::String)
	s = Scenario(JSON.parsefile(fp))
	build_reopt!(m, REoptInputs(s))
	nothing
end

"""
	build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)
Add variables and constraints for REopt model.
"""
function build_reopt!(m::JuMP.AbstractModel, p::REoptInputs)

	add_variables!(m, p)

	for ts in p.time_steps_without_grid

		for tier in 1:p.s.electric_tariff.n_energy_tiers
			fix(m[:dvGridPurchase][ts, tier] , 0.0, force=true)
		end

		for t in p.s.storage.types.elec
			fix(m[:dvGridToStorage][t, ts], 0.0, force=true)
		end

        if !isempty(p.s.electric_tariff.export_bins)
            for t in p.techs.elec, u in p.export_bins_by_tech[t]
                fix(m[:dvProductionToGrid][t, u, ts], 0.0, force=true)
            end
        end
	end

	for b in p.s.storage.types.all
		if p.s.storage.attr[b].max_kw == 0 || p.s.storage.attr[b].max_kwh == 0
			@constraint(m, [ts in p.time_steps], m[:dvStoredEnergy][b, ts] == 0)
			@constraint(m, m[:dvStorageEnergy][b] == 0)
			@constraint(m, [ts in p.time_steps], m[:dvDischargeFromStorage][b, ts] == 0)
			if b in p.s.storage.types.elec
				@constraint(m, m[:dvStoragePower][b] == 0)
				@constraint(m, [ts in p.time_steps], m[:dvGridToStorage][b, ts] == 0)
				@constraint(m, [t in p.techs.elec, ts in p.time_steps_with_grid],
						m[:dvProductionToStorage][b, t, ts] == 0)
			end
		else
			add_storage_size_constraints(m, p, b)
			add_general_storage_dispatch_constraints(m, p, b)
			if b in p.s.storage.types.elec
				add_elec_storage_dispatch_constraints(m, p, b)
			elseif b in p.s.storage.types.hot
				add_hot_thermal_storage_dispatch_constraints(m, p, b)
			elseif b in p.s.storage.types.cold
				add_cold_thermal_storage_dispatch_constraints(m, p, b)
			else
				@error("Invalid storage does not fall in a thermal or electrical set")
			end
		end
	end

	if any(max_kw->max_kw > 0, (p.s.storage.attr[b].max_kw for b in p.s.storage.types.elec))
		add_storage_sum_constraints(m, p)
	end

	add_production_constraints(m, p)

    m[:TotalTechCapCosts] = 0.0
    m[:TotalPerUnitProdOMCosts] = 0.0
    m[:TotalPerUnitHourOMCosts] = 0.0
    m[:TotalFuelCosts] = 0.0
    m[:TotalProductionIncentive] = 0
	m[:dvComfortLimitViolationCost] = 0.0
	m[:TotalCHPStandbyCharges] = 0
	m[:OffgridOtherCapexAfterDepr] = 0.0

	if !isempty(p.techs.all)
		add_tech_size_constraints(m, p)

        if !isempty(p.techs.no_curtail)
            add_no_curtail_constraints(m, p)
        end

        if !isempty(p.techs.gen)
            add_gen_constraints(m, p)
            m[:TotalPerUnitProdOMCosts] += m[:TotalGenPerUnitProdOMCosts]
            m[:TotalFuelCosts] += m[:TotalGenFuelCosts]
        end

        if !isempty(p.techs.chp)
            add_chp_constraints(m, p)
            m[:TotalPerUnitProdOMCosts] += m[:TotalCHPPerUnitProdOMCosts]
            m[:TotalFuelCosts] += m[:TotalCHPFuelCosts]
            m[:TotalPerUnitHourOMCosts] += m[:TotalHourlyCHPOMCosts]

			if p.s.chp.standby_rate_us_dollars_per_kw_per_month > 1.0e-7
				m[:TotalCHPStandbyCharges] += sum(p.s.financial.pwf_e * 12 * p.s.chp.standby_rate_us_dollars_per_kw_per_month * m[:dvSize][t] for t in p.techs.chp)
			end

			if !isempty(p.techs.thermal)
				m[:TotalTechCapCosts] += sum(p.s.chp.supplementary_firing_capital_cost_per_kw * m[:dvSupplementaryFiringSize][t] for t in p.techs.chp)
			end
        end

        if !isempty(p.techs.boiler)
            add_boiler_tech_constraints(m, p)
        end

		if !isempty(p.techs.cooling)
            add_cooling_tech_constraints(m, p)
        end

        if !isempty(p.techs.thermal)
            add_thermal_load_constraints(m, p)  # split into heating and cooling constraints?
        end

        if !isempty(p.techs.pbi)
            @warn "adding binary variable(s) to model production based incentives"
            add_prod_incent_vars_and_constraints(m, p)
        end
	end

	add_elec_load_balance_constraints(m, p)

	if p.s.settings.off_grid_flag
		add_operating_reserve_constraints(m, p)
	end

	if !isempty(p.s.electric_tariff.export_bins)
		add_export_constraints(m, p)
	end

	if !isempty(p.s.electric_tariff.monthly_demand_rates)
		add_monthly_peak_constraint(m, p)
	end

	if !isempty(p.s.electric_tariff.tou_demand_ratchet_time_steps)
		add_tou_peak_constraint(m, p)
	end

	if !(p.s.electric_utility.allow_simultaneous_export_import) & !isempty(p.s.electric_tariff.export_bins)
		add_simultaneous_export_import_constraint(m, p)
	end

	if p.s.electric_tariff.n_energy_tiers > 1
		add_energy_tier_constraints(m, p)
	end

    if p.s.electric_tariff.demand_lookback_percent > 0
        add_demand_lookback_constraints(m, p)
    end

    if !isempty(p.s.electric_tariff.coincpeak_periods)
        add_coincident_peak_charge_constraints(m, p)
    end

    if !isempty(setdiff(p.techs.all, p.techs.segmented))
        m[:TotalTechCapCosts] += p.third_party_factor *
            sum( p.cap_cost_slope[t] * m[:dvPurchaseSize][t] for t in setdiff(p.techs.all, p.techs.segmented))
    end

    if !isempty(p.techs.segmented)
        @warn "adding binary variable(s) to model cost curves"
        add_cost_curve_vars_and_constraints(m, p)
        for t in p.techs.segmented  # cannot have this for statement in sum( ... for t in ...) ???
            m[:TotalTechCapCosts] += p.third_party_factor * (
                sum(p.cap_cost_slope[t][s] * m[Symbol("dvSegmentSystemSize"*t)][s] +
                    p.seg_yint[t][s] * m[Symbol("binSegment"*t)][s] for s in 1:p.n_segs_by_tech[t])
            )
        end
    end

	@expression(m, TotalStorageCapCosts, p.third_party_factor * (
		sum( p.s.storage.attr[b].net_present_cost_per_kw * m[:dvStoragePower][b] for b in p.s.storage.types.elec) +
		sum( p.s.storage.attr[b].net_present_cost_per_kwh * m[:dvStorageEnergy][b] for b in p.s.storage.types.all )
	))

	@expression(m, TotalPerUnitSizeOMCosts, p.third_party_factor * p.pwf_om *
		sum( p.om_cost_per_kw[t] * m[:dvSize][t] for t in p.techs.all )
	)

	add_elec_utility_expressions(m, p)

	if !isempty(p.s.electric_utility.outage_durations)
		add_dv_UnservedLoad_constraints(m,p)
		add_outage_cost_constraints(m,p)
		add_MG_production_constraints(m,p)
		if !isempty(p.s.storage.types.elec)
			add_MG_storage_dispatch_constraints(m,p)
		else
			fix_MG_storage_variables(m,p)
		end
		add_cannot_have_MG_with_only_PVwind_constraints(m,p)
		add_MG_size_constraints(m,p)

		if !isempty(p.techs.gen)
			add_MG_fuel_burn_constraints(m,p)
			add_binMGGenIsOnInTS_constraints(m,p)
		else
			m[:ExpectedMGFuelUsed] = 0
			m[:ExpectedMGFuelCost] = 0
			@constraint(m, [s in p.s.electric_utility.scenarios, tz in p.s.electric_utility.outage_start_time_steps, ts in p.s.electric_utility.outage_time_steps],
				m[:binMGGenIsOnInTS][s, tz, ts] == 0
			)
		end

		if p.s.site.min_resil_time_steps > 0
			add_min_hours_crit_ld_met_constraint(m,p)
		end
	end

	if p.s.settings.off_grid_flag
		offgrid_other_capex_depr_savings = get_offgrid_other_capex_depreciation_savings(p.s.financial.offgrid_other_capital_costs,
			p.s.financial.owner_discount_pct, p.s.financial.analysis_years, p.s.financial.owner_tax_pct)
		m[:OffgridOtherCapexAfterDepr] = p.s.financial.offgrid_other_capital_costs - offgrid_other_capex_depr_savings
	end

	#################################  Objective Function   ########################################
	@expression(m, Costs,
		# Capital Costs
		m[:TotalTechCapCosts] + TotalStorageCapCosts +

		# Fixed O&M, tax deductible for owner
		TotalPerUnitSizeOMCosts * (1 - p.s.financial.owner_tax_pct) +

		# Variable O&M, tax deductible for owner
		(m[:TotalPerUnitProdOMCosts] + m[:TotalPerUnitHourOMCosts]) * (1 - p.s.financial.owner_tax_pct) +

		# Total Fuel Costs, tax deductible for offtaker
        m[:TotalFuelCosts] * (1 - p.s.financial.offtaker_tax_pct) +

		# CHP Standby Charges
		m[:TotalCHPStandbyCharges] * (1 - p.s.financial.offtaker_tax_pct) +

		# Utility Bill, tax deductible for offtaker
		m[:TotalElecBill] * (1 - p.s.financial.offtaker_tax_pct) -

        # Subtract Incentives, which are taxable
		m[:TotalProductionIncentive] * (1 - p.s.financial.owner_tax_pct) +

		# Comfort limit violation costs
		m[:dvComfortLimitViolationCost] +

		# Additional annual costs, tax deductible for owner (only applies when off_grid_flag is true)
		p.s.financial.offgrid_other_annual_costs * p.pwf_om * (1 - p.s.financial.owner_tax_pct) +

		# Additional capital costs, depreciable (only applies when off_grid_flag is true)
		m[:OffgridOtherCapexAfterDepr]

	);
	if !isempty(p.s.electric_utility.outage_durations)
		add_to_expression!(Costs, m[:ExpectedOutageCost] + m[:mgTotalTechUpgradeCost] + m[:dvMGStorageUpgradeCost] + m[:ExpectedMGFuelCost])
	end

	@objective(m, Min, m[:Costs])

	if !(isempty(p.s.storage.types.elec)) && p.s.settings.add_soc_incentive # Keep SOC high
		@objective(m, Min, m[:Costs] -
		sum(m[:dvStoredEnergy][b, ts] for b in p.s.storage.types.elec, ts in p.time_steps) /
			(8760. / p.hours_per_time_step)
		)

	end

	for b in p.s.storage.types.elec
		if p.s.storage.attr[b].model_degradation
			add_degradation(m, p; b=b)
			if p.s.settings.add_soc_incentive
				@warn "Settings.add_soc_incentive is set to true but no incentive will be added because it conflicts with the battery degradation model."
			end
		end
	end

	nothing
end

function run_reopt(m::JuMP.AbstractModel, p::REoptInputs; organize_pvs=true)

	build_reopt!(m, p)

	@info "Model built. Optimizing..."
	tstart = time()
	optimize!(m)
	opt_time = round(time() - tstart, digits=3)
	if termination_status(m) == MOI.TIME_LIMIT
		status = "timed-out"
    elseif termination_status(m) == MOI.OPTIMAL
        status = "optimal"
    else
		status = "not optimal"
		@warn "REopt solved with " termination_status(m), ", returning the model."
		return m
	end
	@info "REopt solved with " termination_status(m)
	@info "Solving took $(opt_time) seconds."

	tstart = time()
	results = reopt_results(m, p)
	time_elapsed = time() - tstart
	@info "Results processing took $(round(time_elapsed, digits=3)) seconds."
	results["status"] = status
	results["solver_seconds"] = opt_time
    if organize_pvs && !isempty(p.techs.pv)  # do not want to organize_pvs when running BAU case in parallel b/c then proform code fails
        organize_multiple_pv_results(p, results)
    end
	return results
end

"""
    add_variables!(m::JuMP.AbstractModel, p::REoptInputs)

Add JuMP variables to the model.
"""
function add_variables!(m::JuMP.AbstractModel, p::REoptInputs)
	@variables m begin
		dvSize[p.techs.all] >= 0  # System Size of Technology t [kW]
		dvPurchaseSize[p.techs.all] >= 0  # system kW beyond existing_kw that must be purchased
		dvGridPurchase[p.time_steps, 1:p.s.electric_tariff.n_energy_tiers] >= 0  # Power from grid dispatched to meet electrical load [kW]
		dvRatedProduction[p.techs.all, p.time_steps] >= 0  # Rated production of technology t [kW]
		dvCurtail[p.techs.all, p.time_steps] >= 0  # [kW]
		dvProductionToStorage[p.s.storage.types.all, p.techs.all, p.time_steps] >= 0  # Power from technology t used to charge storage system b [kW]
		dvDischargeFromStorage[p.s.storage.types.all, p.time_steps] >= 0 # Power discharged from storage system b [kW]
		dvGridToStorage[p.s.storage.types.elec, p.time_steps] >= 0 # Electrical power delivered to storage by the grid [kW]
		dvStoredEnergy[p.s.storage.types.all, 0:p.time_steps[end]] >= 0  # State of charge of storage system b
		dvStoragePower[p.s.storage.types.all] >= 0   # Power capacity of storage system b [kW]
		dvStorageEnergy[p.s.storage.types.all] >= 0   # Energy capacity of storage system b [kWh]
		dvPeakDemandTOU[p.ratchets, 1:p.s.electric_tariff.n_tou_demand_tiers] >= 0  # Peak electrical power demand during ratchet r [kW]
		dvPeakDemandMonth[p.months, 1:p.s.electric_tariff.n_monthly_demand_tiers] >= 0  # Peak electrical power demand during month m [kW]
		MinChargeAdder >= 0
	end

	if !isempty(p.techs.gen)  # Problem becomes a MILP
		@warn """Adding binary variable to model gas generator.
				 Some solvers are very slow with integer variables"""
		@variables m begin
			binGenIsOnInTS[p.techs.gen, p.time_steps], Bin  # 1 If technology t is operating in time step h; 0 otherwise
		end
	end

    if !isempty(p.techs.fuel_burning)
		@variable(m, dvFuelUsage[p.techs.fuel_burning, p.time_steps] >= 0) # Fuel burned by technology t in each time step [kWh]
    end

    if !isempty(p.s.electric_tariff.export_bins)
        @variable(m, dvProductionToGrid[p.techs.elec, p.s.electric_tariff.export_bins, p.time_steps] >= 0)
    end

	if !(p.s.electric_utility.allow_simultaneous_export_import) & !isempty(p.s.electric_tariff.export_bins)
		@warn """Adding binary variable to prevent simultaneous grid import/export.
				 Some solvers are very slow with integer variables"""
		@variable(m, binNoGridPurchases[p.time_steps], Bin)
	end

    if !isempty(p.techs.thermal)
        @variables m begin
			dvThermalProduction[p.techs.thermal, p.time_steps] >= 0
			dvSupplementaryThermalProduction[p.techs.chp, p.time_steps] >= 0
			dvSupplementaryFiringSize[p.techs.chp] >= 0  #X^{\sigma db}_{t}: System size of CHP with supplementary firing [kW]
		end
        if !isempty(p.techs.chp)
            @variable(m, dvProductionToWaste[p.techs.chp, p.time_steps] >= 0)
        end
    end

	if !isempty(p.s.electric_utility.outage_durations) # add dvUnserved Load if there is at least one outage
		@warn """Adding binary variable to model outages.
				 Some solvers are very slow with integer variables"""
		max_outage_duration = maximum(p.s.electric_utility.outage_durations)
		outage_time_steps = p.s.electric_utility.outage_time_steps
		tZeros = p.s.electric_utility.outage_start_time_steps
		S = p.s.electric_utility.scenarios
		# TODO: currently defining more decision variables than necessary b/c using rectangular arrays, could use dicts of decision variables instead
		@variables m begin # if there is more than one specified outage, there can be more othan one outage start time
			dvUnservedLoad[S, tZeros, outage_time_steps] >= 0 # unserved load not met by system
			dvMGProductionToStorage[p.techs.elec, S, tZeros, outage_time_steps] >= 0 # Electricity going to the storage system during each time_step
			dvMGDischargeFromStorage[S, tZeros, outage_time_steps] >= 0 # Electricity coming from the storage system during each time_step
			dvMGRatedProduction[p.techs.elec, S, tZeros, outage_time_steps]  # MG Rated Production at every time_step.  Multiply by ProdFactor to get actual energy
			dvMGStoredEnergy[S, tZeros, 0:max_outage_duration] >= 0 # State of charge of the MG storage system
			dvMaxOutageCost[S] >= 0 # maximum outage cost dependent on number of outage durations
			dvMGTechUpgradeCost[p.techs.elec] >= 0
			dvMGStorageUpgradeCost >= 0
			dvMGsize[p.techs.elec] >= 0

			dvMGFuelUsed[p.techs.elec, S, tZeros] >= 0
			dvMGMaxFuelUsage[S] >= 0
			dvMGMaxFuelCost[S] >= 0
			dvMGCurtail[p.techs.elec, S, tZeros, outage_time_steps] >= 0

			binMGStorageUsed, Bin # 1 if MG storage battery used, 0 otherwise
			binMGTechUsed[p.techs.elec], Bin # 1 if MG tech used, 0 otherwise
			binMGGenIsOnInTS[S, tZeros, outage_time_steps], Bin
		end
	end

	if p.s.settings.off_grid_flag
		@variables m begin
			dvOpResFromBatt[p.s.storage.types.elec, p.time_steps_without_grid] >= 0 # Operating reserves provided by the electric storage [kW]
			dvOpResFromTechs[p.techs.providing_oper_res, p.time_steps_without_grid] >= 0 # Operating reserves provided by techs [kW]
			1 >= dvOffgridLoadServedFraction[p.time_steps_without_grid] >= 0 # Critical load served in each time_step. Applied in off-grid scenarios only. [fraction]
		end
	end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
function create_lp_model!(model)
    model.options.lp_optimizer === nothing && return
    com = model.inner
    com.sense == MOI.FEASIBILITY_SENSE && return
    lp_model = Model()

    set_optimizer(lp_model, model.options.lp_optimizer)
    lp_x = Vector{VariableRef}(undef, length(com.search_space))
    for variable in com.search_space
        lp_x[variable.idx] = @variable(
            lp_model,
            lower_bound = variable.lower_bound,
            upper_bound = variable.upper_bound
        )
    end
    lp_backend = backend(lp_model)
    # iterate through all constraints and add all supported constraints
    for constraint in com.constraints
        if MOI.supports_constraint(
            model.options.lp_optimizer.optimizer_constructor(),
            typeof(constraint.fct),
            typeof(constraint.set),
        )
            MOI.add_constraint(lp_backend, constraint.fct, constraint.set)
        elseif constraint.set isa CPE.Strictly # transform a < into a <= by changing the rhs
            if MOI.supports_constraint(
                model.options.lp_optimizer.optimizer_constructor(),
                typeof(constraint.fct),
                typeof(constraint.set.set),
            )
                rhs = get_rhs_from_strictly(model.inner, constraint, constraint.fct, constraint.set)
                computed_set = typeof(constraint.set.set)(rhs)
                MOI.add_constraint(lp_backend, constraint.fct, computed_set)
            end
        end
    end
    # add objective
    !MOI.supports(lp_backend, MOI.ObjectiveSense()) &&
        @error "The given lp solver doesn't allow objective functions"
    typeof_objective = typeof(com.objective.fct)
    if MOI.supports(lp_backend, MOI.ObjectiveFunction{typeof_objective}())
        MOI.set(lp_backend, MOI.ObjectiveFunction{typeof_objective}(), com.objective.fct)
    else
        @error "The given `lp_optimizer` doesn't support the objective function $(typeof_objective)"
    end
    if MOI.supports(lp_backend, MOI.ObjectiveSense())
        MOI.set(lp_backend, MOI.ObjectiveSense(), com.sense)
    else
        @error "The given `lp_optimizer` doesn't support setting `ObjectiveSense`"
    end
    com.lp_x = lp_x
    com.lp_model = lp_model
end

function create_lp_variable!(lp_model, lp_x; lb = typemin(Int64), ub = typemax(Int64))
    v = @variable(lp_model, lower_bound = lb, upper_bound = ub)
    push!(lp_x, v)
    return length(lp_x)
end


text:
function construct_device!(canonical::Canonical,
                           sys::PSY.System,
                           model::DeviceModel{L, D},
                           ::Type{S};
                           kwargs...) where {L<:PSY.ControllableLoad,
                                             D<:AbstractControllablePowerLoadFormulation,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(L, sys)

    if validate_available_devices(devices, L)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    reactivepower_variables!(canonical, devices)

    #Constraints
    activepower_constraints!(canonical, devices, D, S)

    reactivepower_constraints!(canonical, devices, D, S)

    feedforward!(canonical, L, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

function construct_device!(canonical::Canonical,
                           sys::PSY.System,
                           model::DeviceModel{L, D},
                           ::Type{S};
                           kwargs...) where {L<:PSY.ControllableLoad,
                                             D<:AbstractControllablePowerLoadFormulation,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(L, sys)

    if validate_available_devices(devices, L)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    #Constraints
    activepower_constraints!(canonical, devices, D, S)

    feedforward!(canonical, L, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, D, S)

    return

end

function construct_device!(canonical::Canonical,
                           sys::PSY.System,
                           model::DeviceModel{L, InterruptiblePowerLoad},
                           ::Type{S};
                           kwargs...) where {L<:PSY.ControllableLoad,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(L, sys)

    if validate_available_devices(devices, L)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    reactivepower_variables!(canonical, devices)

    commitment_variables!(canonical, devices)

    #Constraints
    activepower_constraints!(canonical, devices, model.formulation, S)

    reactivepower_constraints!(canonical, devices, model.formulation, S)

    feedforward!(canonical, L, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, model.formulation, S)

    return

end

function construct_device!(canonical::Canonical,
                           sys::PSY.System,
                           model::DeviceModel{L, InterruptiblePowerLoad},
                           ::Type{S};
                           kwargs...) where {L<:PSY.ControllableLoad,
                                             S<:PM.AbstractActivePowerModel}

    devices = PSY.get_components(L, sys)

    if validate_available_devices(devices, L)
        return
    end

    #Variables
    activepower_variables!(canonical, devices)

    commitment_variables!(canonical, devices)

    #Constraints
    activepower_constraints!(canonical, devices, model.formulation, S)

    feedforward!(canonical, L, model.feedforward)

    #Cost Function
    cost_function(canonical, devices, model.formulation, S)

    return

end

function construct_device!(canonical::Canonical,
                           sys::PSY.System,
                           model::DeviceModel{L, StaticPowerLoad},
                           ::Type{S};
                           kwargs...) where {L<:PSY.ElectricLoad,
                                             S<:PM.AbstractPowerModel}

    devices = PSY.get_components(L, sys)

    if validate_available_devices(devices, L)
        return
    end

    nodal_expression!(canonical, devices, S)

    return

end

function construct_device!(canonical::Canonical,
                           sys::PSY.System,
                           model::DeviceModel{L, D},
                           ::Type{S};
                           kwargs...) where {L<:PSY.StaticLoad,
                                             D<:AbstractControllablePowerLoadFormulation,
                                             S<:PM.AbstractPowerModel}

    if D != StaticPowerLoad
        @warn("The Formulation $(D) only applies to FormulationControllable Loads, \n Consider Changing the Device Formulation to StaticPowerLoad")
    end

    construct_device!(canonical,
                      sys,
                      DeviceModel(L, StaticPowerLoad),
                      S;
                      kwargs...)

end


text:
@testset "CrystalViscoPlastic" begin
    slipsystems = MaterialModels.slipsystems(MaterialModels.FCC(), rand(RodriguesParam))
    m = MaterialModels.CrystalViscoPlasticRed(E=200e3, ν=0.3, τ_y=400., H_iso=1e3, H_kin=1e3, q=0.0, α_∞=100., t_star=20., σ_c=50., m=10., slipsystems=slipsystems)
    S = MaterialModels.get_n_slipsystems(m)
    state = initial_material_state(m)

    # strain at yield point for uniaxial stress
    ε11_yield_all = [m.τ_y/ms[1,1]/m.E for ms in m.MS]
    ε11_yield = minimum(abs.(ε11_yield_all))
    i = findfirst(v->v==ε11_yield, abs.(ε11_yield_all))

    # elastic response (yield point)
    Δε = SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? ε11_yield : (i == j ? -ε11_yield*m.ν : 0.0))
    σ, ∂σ∂ε, temp_state = material_response(m, Δε, state)
    @test σ ≈ SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? ε11_yield*m.E : 0.0)
    @test abs(σ ⊡ m.MS[i]) ≈ m.τ_y
    @test ∂σ∂ε ≈ m.Eᵉ
    @test temp_state.κ ≈ zeros(S)
    @test temp_state.α ≈ zeros(S)
    @test temp_state.μ ≈ zeros(S)

    # plastic branch
    Δε = SymmetricTensor{2,3,Float64}((i,j) -> i==1 && j==1 ? 2ε11_yield : (i == j ? -2ε11_yield*m.ν : 0.0))
    σ, ∂σ∂ε, temp_state = material_response(m, Δε, state)
    @test abs(σ ⊡ m.MS[i]) > m.τ_y
    # these are rather lousy tests, we should come up with something better
    @test !(temp_state.κ[i] ≈ 0.0)
    @test !(temp_state.α[i] ≈ 0.0)
    @test !(temp_state.μ[i] ≈ 0.0)
end

function get_visco_plastic_loading()
    strain1 = range(0.0,  0.01, length=5)
    strain2 = range(0.005, 0.001, length=5)
    strain3 = range(0.001, 0.007, length=5)

    # _C = [strain1..., strain2..., strain3...]
    _C = [strain1[i]-strain1[i-1] for i=2:length(strain1)]
    ε = [SymmetricTensor{2,3}((x, x/10, 0.0, 0.0, 0.0, 0.0)) for x in _C]

    return ε
end

@testset "CrystalViscoPlasticRed jld2" begin
    slipsystems = MaterialModels.slipsystems(MaterialModels.FCC(),RodriguesParam(-0.665267, 0.0203875, 1.08633))
    m = MaterialModels.CrystalViscoPlasticRed(E=200e3, ν=0.3, τ_y=400., H_iso=1e3, H_kin=1e3, q=0.0, α_∞=100., t_star=20., σ_c=50., m=10., slipsystems=slipsystems)

    loading = get_visco_plastic_loading()
    # want the same results as from the other implementation
    check_jld2(m, loading, "CrystalViscoPlastic1")#, debug_print=true, OVERWRITE_JLD2=true)
end


text:
# ∀ ci ∈ C(ci/c(i+1) ∈ N) , where C=⟨c1,…,cn⟩, ci>ci+1 og 1≤i≤n−1
function can_use_greedy(c)
    n = length(c)
    for i = 1:(n-1) # 1≤i≤n−1
        # check if there is a remainder of modulo, i.e. the division does not produce a natural number
        # Also if ci>ci+1, the modulo returns a remainder..
        if (c[i] % c[i+1] != 0)
            return false
        end
    end
    return true
end

function min_coins_greedy(coins, value)
    i = 1 # current coin
    n_coins = 0
    while (value > 0)
        if coins[i] > value
            i = i+1
        else
            value = value - coins[i]
            n_coins = n_coins+1
        end
    end
    return n_coins
end

function min_coins_dynamic(coins,value)
    results = fill(typemax(Int),value)
    new_coins = [] # holds the coins we'll be using
    for i = 1:length(coins)
        if (coins[i] <= value) # check if the coins are actually usable
            results[coins[i]] = 1
            push!(new_coins,coins[i])
        end
    end
    println(new_coins)
    for val = 1:(value) # element [1] is 0
        if (results[val] < typemax(Int))
            continue
        end
        println(val)
        best = typemax(Int)
        for i = 1:length(new_coins)
            if (coins[i] <= val)
                current = 1 + results[val-coins[i]]

                if (current < best)
                    best = current
                end
            end
        end
        results[val] = best
    end
    return results[value]
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
using NonconvexAugLagLab
using Test, LinearAlgebra

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

m = Model(f)
addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
add_ineq_constraint!(
    m,
    FunctionWrapper(x -> [g(x, 2, 0), g(x, -1, 1)], 2),
)
x0 = [2.0, 2.0]

#alg = AugLag2(primaloptimizer = MMA87())
alg = AugLag2()
options = AugLag2Options(alg)
r = optimize(m, alg, x0, options = options)
@test abs(r.minimum - sqrt(8/27)) < 1e-4
@test norm(r.minimizer - [1/3, 8/27]) < 1e-4


text:
function initialize_pss!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, A, TG, PSY.PSSFixed},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov} end

#Currently not working properly.
#=
function initialize_pss!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, A, TG, PSY.PSSSimple},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov}

    ω0 = PSY.get_ext(dynamic_device)[CONTROL_REFS][ω_ref_index]
    P0 = PSY.get_active_power(static)
    τe = get_inner_vars(dynamic_device)[τe_var]

    #Get parameters
    pss = PSY.get_pss(dynamic_device)
    K_ω = PSY.get_K_ω(pss)
    K_p = PSY.get_K_p(pss)

    #It assumes that the system is balanced so:
    get_inner_vars(dynamic_device)[V_pss_var] = K_ω * (ω0 - ω0) + K_p * (ω0 * τe - P0)

end
=#


text:
## activate package environment
cd(@__DIR__)
using Pkg
Pkg.activate("Project.toml")

## percentage x of Al, wavelengths λ and index of refraction n from publication Aspnes 1986
x = [0.491 0.590]
λ = [0.520]
n = [3.8264 + 0.19880im;
     3.7477 + 0.15232im]
N = real(n)
k = imag(n)

## interpolation of n for wanted x
x_wanted = [0.509]
l = length(λ)
m = length(x_wanted)
N_wanted = zeros(m,l)
k_wanted = zeros(m,l)
for i = 1:l
    for j = 0:m-1
        local coeffs_N = [x[j+1,:] ones(2)]\N[(1+j*2):(1+j*2)+1,i]
        local coeffs_k = [x[j+1,:] ones(2)]\k[(1+j*2):(1+j*2)+1,i]
        N_wanted[i,j+1] = coeffs_N[1]*x_wanted[j+1] + coeffs_N[2]
        k_wanted[i,j+1] = coeffs_k[1]*x_wanted[j+1] + coeffs_k[2]
    end
end
N_wanted = round.(N_wanted,sigdigits=4)
k_wanted = round.(k_wanted,sigdigits=4)

## output solution
println("")
for j = 1:m
    for i = 1:l
        println("n(x = $(x_wanted[j]), λ = $(lambda[i])) = $(N_wanted[i,j]) + $(k_wanted[i,j])i")
    end
end

max_stars_repo_name: MaxGrossmann/TRPL

max_issues_repo_name: MaxGrossmann/TRPL

max_forks_repo_name: MaxGrossmann/TRPL


text:
using ds_DistributedBoosting

# Install the R-Packages opal, dsBaseClient, dsStatsClient, and dsModellingClient in your R-library

# Please change to your settings!
YourRpath = "C:/Users/Username/Documents/R/win-library/R-Version"
YourUrls = ["https://server_url1:port", "https://server_url2:port"]
YourUser = "user"
YourPassword = "password"
YourTable = ["Projectname.Tablename"]
YourServernames = ["Server1","Server2"]

# Basic parameters
YourStepno = 30
YourEndpoint = "Y"
YourStartBuffer = 2
YourRunBuffer = 2
YourCheck = true
YourTolerance=0.00000000000005

YourWarning = true
YourShrink = 0.1
YourMaxvar = 15
YourLabels = ["X1","X2","X3","X4","X5","X6","X7","X8","X9","X10"]

# Perform Boosting for all available variables
result = ds_boosting(YourStepno, YourEndpoint, YourStartBuffer, YourRunBuffer,
            YourRpath, YourUrls,YourUser,YourPassword,YourTable,
            YourServernames,YourCheck)
# Selected variables are:
unique(result.selections)
# For the example data, covariances are called 6 times.
# Of a=2 variables before starting, for additional 5 after 1 boosting step,
# for additional 4 after 7 boosting steps, 4 after 10, 5 after 14 and 3 after 25.
# 9 unique variables are selected:
# X24, X9, X21, X30, X27, X18, X6, X12, X15

# Perform Boosting for selected variables (YourLabels). Check will be performed only for these.
result_sel = ds_boosting(YourStepno, YourEndpoint, YourStartBuffer, YourRunBuffer,
            YourRpath, YourUrls,YourUser,YourPassword,YourTable,
            YourServernames,YourCheck,YourWarning, YourTolerance,
            YourShrink, YourMaxvar, YourLabels)
# For the example data, covariances are called 3 times.
# Of a=2 variables before starting, for additional 4 after 4 boosting steps,
# and 3 after 6 boosting steps.
# 6 unique variables are selected:
# X9, X6, X3, X5, X10, X8


text:
using VortexDistributions
using Test

@testset "Creation and detection " begin include("tests.jl") end


text:
module SubModule2
    import AbstractModule.foo
    foo(x::Float64) = println("Hello there Float")
end


text:
"""
Abstract type for probability inference
"""
abstract type InferenceMethod end

@inline function _ensure_query_nodes_in_pgm_and_not_in_evidence(query::NodeNames, nodenames::NodeNames, evidence::Assignment)
    isempty(query) && return

    q = first(query)
    (q ∈ nodenames) || throw(ArgumentError("Query $q is not in the probabilistic graphical model"))
    haskey(evidence, q) && throw(ArgumentError("Query $q is part of the evidence"))

    return _ensure_query_nodes_in_pgm_and_not_in_evidence(query[2:end], nodenames, evidence)
end

"""
Type for capturing the inference state
"""
struct InferenceState{PGM<:ProbabilisticGraphicalModel}
    pgm::PGM
    query::NodeNames
    evidence::Assignment

    function InferenceState{PGM}(pgm::PGM, query::NodeNames, evidence::Assignment) where PGM<:ProbabilisticGraphicalModel

        _ensure_query_nodes_in_pgm_and_not_in_evidence(query, names(pgm), evidence)
        return new(pgm, query, evidence)
    end
end
function InferenceState(pgm::PGM, query::NodeNameUnion, evidence::Assignment=Assignment()) where {PGM<:ProbabilisticGraphicalModel}
    query = unique(nodeconvert(NodeNames, query))
    return InferenceState{PGM}(pgm, query, evidence)
end

function Base.show(io::IO, inf::InferenceState)
    println(io, "Query: $(inf.query)")
    println(io, "Evidence:")
    for (k, v) in inf.evidence
        println(io, "  $k => $v")
    end
end

"""
    infer(InferenceMethod, InferenceState)
Infer p(query|evidence)
"""
infer(im::InferenceMethod, inf::InferenceState) = error("infer not implemented for $(typeof(im)) and $(typeof(inf))")
infer(im::InferenceMethod, pgm::ProbabilisticGraphicalModel, query::NodeNameUnion; evidence::Assignment=Assignment()) = infer(im, InferenceState(pgm, query, evidence))


text:
using SparseArrays
using LinearAlgebra

mutable struct SRLearner{F<:Number, LU<:LearningUpdate} <: Learner
    ψ::Matrix{F}
    r_w::Matrix{F}
    update::LU

    num_demons::Int
    num_actions::Int

    feature_size::Int

    num_tasks::Int
    feature_projector::FeatureCreator

end

function get_weights(learner::SRLearner)
    return [learner.ψ, learner.r_w]
end

function SRLearner{F}(lu, feature_size, num_demons, num_actions, num_tasks, fp, w_init) where {F<:Number}
    ψ_init = if w_init == 0 0 else 1 end
    SRLearner(ones(F, num_demons-num_tasks, feature_size * num_actions)*ψ_init,
              ones(F, num_tasks, size(fp))*w_init,
              lu,
              num_demons,
              num_actions,
              feature_size,
              num_tasks, fp)
end

SRLearner(lu, feature_size, num_demons, num_actions, num_tasks, fp, w_init) =
    SRLearner{Float64}(lu, feature_size, num_demons, num_actions, num_tasks, fp, w_init)

update(l::SRLearner) = l.update

Base.size(learner::SRLearner) = learner.num_demons

function predict_SF(learner::SRLearner, ϕ::SparseVector, action)
    active_state_action = get_active_action_state_vector(ϕ, action, length(ϕ), learner.num_actions)
    learner.ψ[:, active_state_action.nzind] * active_state_action.nzval
end

function predict(learner::SRLearner, ϕ::SparseVector, action)

    SF = predict_SF(learner, ϕ, action)

    #Column is SF per task
    reshaped_SF = reshape(SF, :, learner.num_tasks)

    Q = learner.r_w * reshaped_SF
    return Q[diagind(Q)]
end

(l::SRLearner)(ϕ, a) = predict(l, ϕ, a)
(l::SRLearner)(ϕ) = [predict(l, ϕ, a) for a ∈ 1:num_actions]

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
const STRING = "String"
const INT = "Int"
const FLOAT = "Float"
const BINARY = "Binary"
const DATE = "Date"
const TIMESTAMP = "Timestamp"
const BINARYCATEGORY = "Binary Category"
const CATEGORY = "Category"
const IMAGE = "Image"
const TEXT = "Text"
const ARRAY = "Array"
const DATATYPES = [STRING, INT, FLOAT, BINARY, DATE, TIMESTAMP,
                   BINARYCATEGORY, CATEGORY, IMAGE, TEXT, ARRAY]

struct Config;
    inputs::Array{Tuple{String, String}, 1};
    outputs::Array{Tuple{String, String}, 1};
end

Config() = Config(Array{Tuple{String, String}, 1}(),
                  Array{Tuple{String, String}, 1}())

function getfeatures(c::Config; ftype="all")
    if ftype == "all"
        features = vcat(c.inputs, c.outputs)
    elseif ftype == "input"
        features = c.inputs
    elseif ftype == "output"
        features = c.outputs
    end
end

function getftypes(c::Config; ftype="all")
    ftypes = [ftype for (fname, ftype) in getfeatures(c; ftype=ftype)]
end

function getfnames(c::Config; ftype="all")
    ftypes = [fname for (fname, ftype) in getfeatures(c; ftype=ftype)]
end

function getfdict(c::Config; ftype="all")
    ftypes = getftypes(c::Config; ftype=ftype)
    Dict(ftype => count(i->i==ftype, ftypes) for ftype in unique(ftypes))
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML