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.
# *********************************************************************************
"""
    Financial

Financial data struct with inner constructor:
```julia
function Financial(;
    om_cost_escalation_pct::Float64 = 0.025,
    elec_cost_escalation_pct::Float64 = 0.023,
    offtaker_tax_pct::Float64 = 0.26,
    offtaker_discount_pct = 0.083,
    third_party_ownership::Bool = false,
    owner_tax_pct::Float64 = 0.26,
    owner_discount_pct::Float64 = 0.083,
    analysis_years::Int = 25,
    VoLL::Union{Array{R,1}, R} where R<:Real = 1.00,
    microgrid_premium_pct::Float64 = 0.3,
    macrs_five_year::Array{Float64,1} = [0.2, 0.32, 0.192, 0.1152, 0.1152, 0.0576],  # IRS pub 946
    macrs_seven_year::Array{Float64,1} = [0.1429, 0.2449, 0.1749, 0.1249, 0.0893, 0.0892, 0.0893, 0.0446]
)
```

!!! note
    When `third_party_ownership` is `false` the offtaker's discount and tax percentages are used throughout the model:
    ```julia
        if !third_party_ownership
            owner_tax_pct = offtaker_tax_pct
            owner_discount_pct = offtaker_discount_pct
        end
    ```
"""
struct Financial
    om_cost_escalation_pct::Float64
    elec_cost_escalation_pct::Float64
    offtaker_tax_pct::Float64
    offtaker_discount_pct
    third_party_ownership::Bool
    owner_tax_pct::Float64
    owner_discount_pct::Float64
    analysis_years::Int
    VoLL::Union{Array{R,1}, R} where R<:Real
    microgrid_premium_pct::Float64
    macrs_five_year::Array{Float64,1}
    macrs_seven_year::Array{Float64,1}

    function Financial(;
        om_cost_escalation_pct::Float64 = 0.025,
        elec_cost_escalation_pct::Float64 = 0.023,
        offtaker_tax_pct::Float64 = 0.26,
        offtaker_discount_pct = 0.083,
        third_party_ownership::Bool = false,
        owner_tax_pct::Float64 = 0.26,
        owner_discount_pct::Float64 = 0.083,
        analysis_years::Int = 25,
        VoLL::Union{Array{R,1}, R} where R<:Real = 1.00,
        microgrid_premium_pct::Float64 = 0.3,
        macrs_five_year::Array{Float64,1} = [0.2, 0.32, 0.192, 0.1152, 0.1152, 0.0576],  # IRS pub 946
        macrs_seven_year::Array{Float64,1} = [0.1429, 0.2449, 0.1749, 0.1249, 0.0893, 0.0892, 0.0893, 0.0446]
    )
        if !third_party_ownership
            owner_tax_pct = offtaker_tax_pct
            owner_discount_pct = offtaker_discount_pct
        end

        return new(
            om_cost_escalation_pct,
            elec_cost_escalation_pct,
            offtaker_tax_pct,
            offtaker_discount_pct,
            third_party_ownership,
            owner_tax_pct,
            owner_discount_pct,
            analysis_years,
            VoLL,
            microgrid_premium_pct,
            macrs_five_year,
            macrs_seven_year
        )
    end
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


text:
module TwoPointFunctions

using ForwardDiff
using StaticArrays

export cosχ
export ∇ϕ₁cosχ
export ∇ϕ₂cosχ
export ∂²θ₁cosχ
export ∂²θ₂cosχ
export ∂²ϕ₁cosχ
export ∂²ϕ₂cosχ
export ∇₁cosχ
export ∇₂cosχ

_cosχ2(x) = sin(x[1])sin(x[3])cos(x[2] - x[4])
_cosχ2((θ₁, ϕ₁), (θ₂, ϕ₂)) = sin(θ₁)sin(θ₂)cos(ϕ₁ - ϕ₂)
cosχ(x::AbstractVector{<:Real}) = cos(x[1])cos(x[3]) + _cosχ2(x)

∂cosχ(x) = ForwardDiff.gradient(cosχ, x)
∂²cosχ(x) = ForwardDiff.hessian(cosχ, x)

"""
	cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Cosine of the angle between the unit vectors denoted by the spherical coordinates
``(\\theta_1,\\phi_1)`` and ``(\\theta_2,\\phi_2)``. This is equivalent to the cosine of the
angular distance along a great circle between these two points on a sphere.

Numerically this is given by
```math
\\cos\\chi((\\theta_1,\\phi_1), (\\theta_2,\\phi_2)) = \\cos(\\theta_1)\\cos(\\theta_2) +
\\sin(\\theta_1)\\sin(\\theta_2)\\cos(\\phi_1 - \\phi_2)
```

The unicode character at the end of the function's name is `\\chi`, and not `x`.

# Examples
```jldoctest
julia> cosχ((0,0), (0,0))
1.0

julia> cosχ((0,0), (pi,0))
-1.0
```
"""
function cosχ end

"""
	∂θ₁cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Partial derivative of [`cosχ`](@ref) with respect to `θ₁`.

# Examples
```jldoctest
julia> ∂θ₁cosχ((0,0), (0,0))
0.0

julia> ∂θ₁cosχ((0,0), (pi/2,0))
1.0
```

See also: [`∂θ₂cosχ`](@ref)
"""
function ∂θ₁cosχ end

"""
	∂θ₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Partial derivative of [`cosχ`](@ref) with respect to `θ₂`.

# Examples
```jldoctest
julia> ∂θ₂cosχ((0,0), (0,0))
0.0

julia> ∂θ₂cosχ((0,0), (pi/2,pi/2))
-1.0
```

See also: [`∂θ₁cosχ`](@ref)
"""
function ∂θ₂cosχ end

"""
	∂ϕ₁cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Partial derivative of [`cosχ`](@ref) with respect to `ϕ₁`.

# Examples
```jldoctest
julia> ∂ϕ₁cosχ((pi/2, 0), (pi/2, pi/3)) ≈ sin(pi/3)
true

julia> ∂ϕ₁cosχ((pi/2, 0), (pi/2, pi/2))
1.0
```

See also: [`∂ϕ₂cosχ`](@ref)
"""
function ∂ϕ₁cosχ end

"""
	∂ϕ₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Partial derivative of [`cosχ`](@ref) with respect to `ϕ₂`.

# Examples
```jldoctest
julia> ∂ϕ₂cosχ((pi/2, 0), (pi/2, pi/3)) ≈ -sin(pi/3)
true

julia> ∂ϕ₂cosχ((pi/2, 0), (pi/2, pi/2))
-1.0
```

See also: [`∂ϕ₁cosχ`](@ref)
"""
function ∂ϕ₂cosχ end

for fn in ["cosχ","∂cosχ","∂²cosχ"]
	@eval $(Symbol(fn))((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = $(Symbol(fn))(@SVector [θ₁,ϕ₁,θ₂,ϕ₂])
end

# use subscripted variables
subscripts = Dict(1=>'\u2081',2=>'\u2082') # \_1 and \_2
powers = Dict(2=>'\u00B2',3=>'\u00B3')
for n in 4:9
	powers[n] = '\u2070' + n
end

# First derivatives
for pt in 1:2, (coord_ind,coord) in enumerate(("θ","ϕ"))
	fname = Symbol("∂$(coord)$(subscripts[pt])cosχ")
	flat_ind = (pt-1)*2 + coord_ind
	@eval ($fname)(n1, n2) = ∂cosχ(n1, n2)[$flat_ind]
	@eval export $fname
end

# Second derivatives
for pt1 in 1:2, (coord1_ind,coord1) in enumerate(("θ","ϕ"))
	for pt2 in 1:2, (coord2_ind,coord2) in enumerate(("θ","ϕ"))
		fname = Symbol("∂$(coord1)$(subscripts[pt1])∂$(coord2)$(subscripts[pt2])cosχ")
		coord1_flat_ind = (pt1-1)*2 + coord1_ind
		coord2_flat_ind = (pt2-1)*2 + coord2_ind
		@eval ($fname)(n1, n2) = ∂²cosχ(n1, n2)[$coord1_flat_ind,$coord2_flat_ind]
	end
end

"""
	∂²θ₁cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Second derivative of [`cosχ`](@ref) with respect to `θ₁`.

# Examples
```jldoctest
julia> ∂²θ₁cosχ((0,0), (0,0))
-1.0

julia> ∂²θ₁cosχ((0,0), (pi,0))
1.0
```
"""
∂²θ₁cosχ(n1, n2) = -cosχ(n1, n2)

"""
	∂²θ₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Second derivative of [`cosχ`](@ref) with respect to `θ₂`.

# Examples
```jldoctest
julia> ∂²θ₂cosχ((0,0), (0,0))
-1.0

julia> ∂²θ₂cosχ((0,0), (pi,0))
1.0
```
"""
∂²θ₂cosχ(n1, n2) = -cosχ(n1, n2)

"""
	TwoPointFunctions.∂θ₁∂θ₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Mixed derivative of [`cosχ`](@ref) with respect to `θ₁` and `θ₂`.

# Examples
```jldoctest
julia> TwoPointFunctions.∂θ₁∂θ₂cosχ((0,0), (0,0))
1.0

julia> TwoPointFunctions.∂θ₁∂θ₂cosχ((0,0), (pi,0))
-1.0
```
"""
function ∂θ₁∂θ₂cosχ end

# The following functions are not defined if the gradients are evaluated at the poles

"""
	∇ϕ₁cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

The ``{\\hat{\\phi}}_1`` component of ``\\nabla_1 \\cos \\chi`` in the spherical polar basis.
Numerically this is equivalent to

```math
\\nabla \\phi_1 \\cos\\chi ((\\theta_1, \\phi_1), (\\theta_2, \\phi_2)) =
\\frac{1}{\\sin(\\theta_1)} \\partial \\phi_1 \\cos\\chi ((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))
```

except `∇ϕ₁cosχ` is defined to avoid coordinate singularities at the poles.

# Examples
```jldoctest
julia> ∇ϕ₁cosχ((0, 0), (0, 0))
-0.0

julia> ∇ϕ₁cosχ((pi/2, 0), (pi/2, pi/2))
1.0
```

See also: [`∂ϕ₁cosχ`](@ref), [`∇ϕ₂cosχ`](@ref)
"""
∇ϕ₁cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = -sin(θ₂)sin(ϕ₁-ϕ₂)

"""
	∇ϕ₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

The ``{\\hat{\\phi}}_2`` component of ``\\nabla_2 \\cos \\chi`` in the spherical polar basis.
Numerically this is equivalent to

```math
\\nabla \\phi_2 \\cos\\chi ((\\theta_1, \\phi_1), (\\theta_2, \\phi_2)) =
\\frac{1}{\\sin(\\theta_2)} \\partial \\phi_2 \\cos\\chi ((\\theta_1, \\phi_1), (\\theta_2, \\phi_2))
```

except `∇ϕ₂cosχ` is defined to avoid coordinate singularities at the poles.

# Examples
```jldoctest
julia> ∇ϕ₂cosχ((0, 0), (0, 0))
0.0

julia> ∇ϕ₂cosχ((pi/2, 0), (pi/2, pi/2))
-1.0
```

See also: [`∂ϕ₂cosχ`](@ref), [`∇ϕ₁cosχ`](@ref)
"""
∇ϕ₂cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = sin(θ₁)sin(ϕ₁-ϕ₂)

"""
	∂²ϕ₁cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Second derivative of [`cosχ`](@ref) with respect to `ϕ₁`.

# Examples
```jldoctest
julia> ∂²ϕ₁cosχ((0,0), (0,0))
-0.0

julia> ∂²ϕ₁cosχ((pi/2,0), (pi/2, 0))
-1.0
```
"""
function ∂²ϕ₁cosχ end

"""
	∂²ϕ₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Second derivative of [`cosχ`](@ref) with respect to `ϕ₂`.

# Examples
```jldoctest
julia> ∂²ϕ₂cosχ((0,0), (0,0))
-0.0

julia> ∂²ϕ₂cosχ((pi/2,0), (pi/2, 0))
-1.0
```
"""
function ∂²ϕ₂cosχ end

for pt in 1:2
	∂²ϕ = "∂²ϕ$(subscripts[pt])cosχ"
	∂ϕ∂ϕ = "∂ϕ$(subscripts[pt])∂ϕ$(subscripts[pt])cosχ"
	@eval $(Symbol(∂²ϕ))(n1, n2) = -_cosχ2(n1, n2)
end

∇ϕ₁∂ϕ₁cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = -sin(θ₂)cos(ϕ₁-ϕ₂)
∇ϕ₁∂ϕ₂cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = sin(θ₂)cos(ϕ₁-ϕ₂)
∂ϕ₁∇ϕ₂cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = sin(θ₁)cos(ϕ₁-ϕ₂)
∇ϕ₂∂ϕ₂cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = -sin(θ₂)cos(ϕ₁-ϕ₂)
∂ϕ₁∇ϕ₁cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = ∇ϕ₁∂ϕ₁cosχ((θ₁,ϕ₁),(θ₂,ϕ₂))
∂ϕ₂∇ϕ₂cosχ((θ₁,ϕ₁)::Tuple{<:Real,<:Real}, (θ₂,ϕ₂)::Tuple{<:Real,<:Real}) = ∇ϕ₂∂ϕ₂cosχ((θ₁,ϕ₁),(θ₂,ϕ₂))

# Higher-order derivatives wrt. ϕ
for pt in 1:2
	for order in 3:2:9
		fn = "∂$(powers[order])ϕ$(subscripts[pt])cosχ"
		@eval $(Symbol(fn))(n1, n2) = $(Symbol("∂ϕ$(subscripts[pt])cosχ"))(n1, n2) * real(im^($order-1))
	end
	for order in 4:2:8
		fn = "∂$(powers[order])ϕ$(subscripts[pt])cosχ"
		@eval $(Symbol(fn))(n1, n2) = $(Symbol("∂²ϕ$(subscripts[pt])cosχ"))(n1, n2) * real(im^($order-2))
	end
end

"""
	∇₁cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Return
``(\\hat{\\theta}_1 \\cdot \\nabla_1 \\cos \\chi,\\,
\\hat{\\phi}_1 \\cdot \\nabla_1 \\cos \\chi)``, where ``{\\nabla}_1`` represents
the gradient on the unit sphere with respect to ``(\\theta_1, \\phi_1)``.

# Examples
```jldoctest
julia> ∇₁cosχ((0,0), (pi/2, pi/3)) == (cos(pi/3), sin(pi/3))
true
```
"""
function ∇₁cosχ end

"""
	∇₂cosχ((θ₁,ϕ₁), (θ₂,ϕ₂))

Return
``(\\hat{\\theta}_2 \\cdot \\nabla_2 \\cos \\chi,\\,
\\hat{\\phi}_2 \\cdot \\nabla_2 \\cos \\chi)``, where ``{\\nabla}_2`` represents
the gradient on the unit sphere with respect to ``(\\theta_2, \\phi_2)``.

# Examples
```jldoctest
julia> ∇₂cosχ((0,0), (pi/2, pi/3))
(-1.0, -0.0)
```
"""
function ∇₂cosχ end

for pt in 1:2
	fname = "∇$(subscripts[pt])cosχ"
	∇θ = "∂θ$(subscripts[pt])cosχ"
	∇ϕ = "∇ϕ$(subscripts[pt])cosχ"
	@eval $(Symbol(fname))(n1, n2) = ($(Symbol(∇θ))(n1, n2), $(Symbol(∇ϕ))(n1, n2))
end

end

ext: jl


text:
import Base.sin
import Base.cos
import Base.exp
import Base.log
import Base.sec
import Base.sqrt
import Base.tan
import Base.tanh
import Base.inv

begin

Xarray = Array

function sin!(x::Xarray)
    @. x = sin(x)
end

function Base.:sin(x::Xarray)
    return sin.(x)
end

function cos!(x::Xarray)
    @. x = cos(x)
end

function Base.:cos(x::Xarray)
    return cos.(x)
end

function sec!(x::Xarray)
    @. x = sec(x)
end

function Base.:sec(x::Xarray)
    return sec.(x)
end

function tan!(x::Xarray)
    @. x = tan(x)
end

function Base.:tan(x::Xarray)
    return tan.(x)
end

function tanh!(x::Xarray)
    @. x = tanh(x)
end

function Base.:tanh(x::Xarray)
    return tanh.(x)
end

function exp!(x::Xarray)
    @. x = exp(x)
end

function Base.:exp(x::Xarray)
    return exp.(x)
end

function log!(x::Xarray)
    @. x = log(x)
end

function Base.:log(x::Xarray)
    return log.(x)
end

function sqrt!(x::Xarray)
    @. x = sqrt(x)
end

function Base.:sqrt(x::Xarray)
    return sqrt.(x)
end

function inv!(x::Xarray)
    @. x = inv(x)
end

function inv(x::Xarray)
    return inv.(x)
end

end # begin

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
include("Viability.jl")

dim = 2
num_per_dim = unsigned(50)
num_total = num_per_dim ^ 2
bounds = COORDINATE_TYPE[-1 1; -1 1]

sd_boundary = 0.5

estimation_args = Dict{String,Any}(
    "num-per-dim" => num_per_dim,
    "bounds" => bounds
)

model_params = Dict{String,Any}(
    "sd-boundary" => 0.5
)

model_info = Dict{String,Any}(
    "model" => "simple-test",
    "dimension" => 2,
    "estimation-args" => estimation_args,
    "model-params" => model_params
)

function step_func(x::POINT_TYPE, a::Float64, k::Float64)
    (x-a)/k
end
f1(x::POINT_TYPE) = step_func(x, 0.0, 2.0)
f2(x::POINT_TYPE) = step_func(x, 1.0, 2.0)
step_functions = [f1, f2]

function is_sunny(x::POINT_TYPE)
    squared_norm = sum(x.^2)
    squared_norm < sd_boundary^2 && ( x[1] > -0.01 || x[1] < -0.3 )
end

##################################################
# start the actual stuff for the run
##################################################

ret = 2
points = create_2d_grid(num_per_dim, bounds)
states = map_over_points(STATE_TYPE, is_sunny, points)
for i=1:2
    points = create_2d_grid(num_per_dim, bounds)
    states = map_over_points(STATE_TYPE, is_sunny, points)
    @time ret = sp_viability_kernel(
        points,
        states,
        step_functions,
        delta=0.05
    )
end

if ret == 0
    println("computed viability kernel")
else:
    println("error value during computation $ret")
end

fname = "output.txt"
print("writing output to '$fname' ... ")
write_result_file(fname, model_info, points, states)
println("done")


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

struct WithForwardDiff{F} <: Function
    f::F
end

@inline (wrapped_f::WithForwardDiff)(xs...) = wrapped_f.f(xs...)

# Desireable for consistent behavior?
# Base.broadcasted(wrapped_f::WithForwardDiff, xs...) = broadcast(wrapped_f.f, xs...)

"""
    fwddiff(f::Base.Callable)::Function

Use `ForwardDiff` dual numbers to implement `ChainRulesCore` pullbacks For

* `fwddiff(f)(args...)
* `fwddiff(f).(args...)

Example:

```
using ForwardDiffPullbacks, StaticArrays

f = (xs...) -> (sum(map(x -> sum(map(x -> x^2, x)), xs)))
xs = (2, (3, 4), SVector(5, 6, 7))
f(xs...) == 139

using ChainRulesCore

y, back = rrule(fwddiff(f), xs...)
y == 139
map(unthunk, back(1)) == (NoTangent(), 4, (6, 8), [10, 12, 14])

using Zygote

Zygote.gradient(fwddiff(f), xs...) == Zygote.gradient(f, xs...)

Xs = map(x -> fill(x, 100), xs)
Zygote.gradient((Xs...) -> sum(fwddiff(f).(Xs...)), Xs...) ==
    Zygote.gradient((Xs...) -> sum(f.(Xs...)), Xs...)
```

The gradient is the same with and without `fwddiff`, but `fwddiff` makes the
gradient calculation a lot faster here.
"""
fwddiff(f::Base.Callable) = WithForwardDiff(f)
export fwddiff


text:
# Run in .julia/registries directory

# find each package w/ a direct dependency on Tables.jl
pkgs = dirname.(readlines(`grep -rl ^Tables General/. --include=Deps.toml`))

pkgnames = [splitpath(x)[end] for x in pkgs]

function parseurl(file)
    repo = readlines(file)[end]
    return strip(split(repo, " = ")[end], '"')
end

urls = [parseurl(joinpath("General", string(first(nm)), nm, "Package.toml")) for nm in pkgnames]

open("/Users/jacobquinn/.julia/dev/Tables/INTEGRATIONS.md", "w+") do io
    println(io, "Packages currently integrating with Tables.jl:")
    for (nm, url) in zip(pkgnames, urls)
        println(io, "* [$nm]($url)")
    end
end


text:
# QCQP example
#    maximize x + y
#
#    s.t.  x, y >= 0
#          x^2 + y^2 <= 1
#
#    solution: (0.71, 0.71) objv = 1.414

using Gurobi, Compat.Test

@testset "QCQP1" begin

    env = Gurobi.Env()
    setparam!(env, "OutputFlag", 0)

    model = Gurobi.Model(env, "qcqp_01", :maximize)

    add_cvars!(model, [1., 1.], 0., Inf)
    update_model!(model)

     # add_qpterms!(model, linearindices, linearcoeffs, qrowinds, qcolinds, qcoeffs, sense, rhs)
    add_qconstr!(model, [], [], [1, 2], [1, 2], [1, 1.], '<', 1.0)
    update_model!(model)

    optimize(model)

    @test isapprox(get_solution(model), [0.707107, 0.707107], atol=1e-5)
    @test isapprox(get_objval(model), 1.414214, atol=1e-5)
end


text:
# This file was generated, do not modify it. # hide
hot = machine(OneHotEncoder(), Xs)

W = transform(hot, Xs)
z = log(ys);

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
import Base: first, last, length,
    getindex, lastindex, eachindex, iterate,
    eltype, similar,
    show
abstract type AbstractKnotSet{T} end

order(t::AbstractKnotSet) = t.k
numintervals(t::AbstractKnotSet) = length(t.t)-1

first(t::AbstractKnotSet) = first(t.t)
last(t::AbstractKnotSet) = last(t.t)
length(t::AbstractKnotSet) = length(t.t) + 2t.k - 2

function getindex(t::AbstractKnotSet, i::Integer)
    if i < order(t)
        first(t)
    elseif i < order(t) + numintervals(t)
        t.t[i-order(t)+1]
    else
        last(t)
    end
end
lastindex(t::AbstractKnotSet) = length(t)
eachindex(t::AbstractKnotSet) = 1:length(t)
function iterate(t::AbstractKnotSet, state=(t[1],0))
    element, i = state
    if i >= length(t)
        return nothing
    end
    element, (t[i+2], i+1)
end

eltype(t::AbstractKnotSet{T}) where T = T
similar(t::AbstractKnotSet{T}) where T = Vector{T}(undef, length(t))

function show(io::IO, t::AbstractKnotSet)
    write(io, "$(typeof(t)) of order k=$(order(t)) on [$(first(t)),$(last(t))] ($(numintervals(t)) intervals)")
end

struct LinearKnotSet{T} <: AbstractKnotSet{T}
    k::Integer
    t::AbstractRange{T}
end
LinearKnotSet(k::Integer, a::Integer, b::Integer, N::Integer) =
    LinearKnotSet(k, range(a, stop=b, length=N+1))

struct ExpKnotSet{T} <: AbstractKnotSet{T}
    k::Integer
    exponents::AbstractRange{T}
    base::T
    t::AbstractVector{T}
    include0::Bool
end
function ExpKnotSet(k::Integer, a::T, b::T, N::Integer;
                    base::T=T(10), include0::Bool=true) where T
    exponents = range(a, stop=b, length=include0 ? N : N+1)
    t = base .^ exponents
    ExpKnotSet(k, exponents, eltype(t)(base), include0 ? vcat(0,t) : t, include0)
end

function show(io::IO, t::ExpKnotSet)
    write(io, "$(typeof(t)) of order k=$(order(t)) on [")
    t.include0 && write(io, "0,")
    write(io, "$(t.base^first(t.exponents))..$(t.base^last(t.exponents))] ($(numintervals(t)) intervals)")
end

# function arcsin_knot_set(k::Integer, a::Integer, b::Integer, N::Integer)
#     N2 = N/2
#     arcsin.(range(-N2, stop=N2, lengthN+1)/N2)*(b-a)/π+(a+b)/2
# end

# arcsin_half_knot_set(a, b, N) = arcsin(linspace(0,N,N+1)/N)*(b-a)*2.0/π+a

# function exp_linear_knot_set(a, b, N)
#     ap = a != 0 ? a : 1e-1
#     @assert t[2][1]+t[2][2] == N
#     [logspace(log10(ap),log10(t[1]), t[2][1]); linspace(t[1],b,t[2][2]+1)[2:end]]
# end

@recipe function plot(t::AbstractKnotSet)
    markershape --> :circle
    y = similar(t)
    y′ = 1
    t′ = -Inf
    for i in eachindex(t)
        if t[i] == t′
            y′ += 1
        else
            y′ = 1
        end
        t′ = t[i]
        y[i] = y′
    end
    collect(t),y
end

export LinearKnotSet, ExpKnotSet, order, numintervals


text:
# ---
# title: 1177. Can Make Palindrome from Substring
# id: problem1177
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, String
# link: <https://leetcode.com/problems/can-make-palindrome-from-substring/description/>
# hidden: true
# ---
#
# Given a string `s`, we make queries on substrings of `s`.
#
# For each query `queries[i] = [left, right, k]`, we may **rearrange**  the
# substring `s[left], ..., s[right]`, and then choose **up to** `k` of them to
# replace with any lowercase English letter.
#
# If the substring is possible to be a palindrome string after the operations
# above, the result of the query is `true`. Otherwise, the result is `false`.
#
# Return an array `answer[]`, where `answer[i]` is the result of the `i`-th
# query `queries[i]`.
#
# Note that: Each letter is counted **individually** for replacement so if for
# example `s[left..right] = "aaa"`, and `k = 2`, we can only replace two of the
# letters.  (Also, note that the initial string `s` is never modified by any
# query.)
#
#
#
# **Example :**
#
#
#
#     Input: s = "abcda", queries = [[3,3,0],[1,2,0],[0,3,1],[0,3,2],[0,4,1]]
#     Output: [true,false,false,true,true]
#     Explanation:
#     queries[0] : substring = "d", is palidrome.
#     queries[1] : substring = "bc", is not palidrome.
#     queries[2] : substring = "abcd", is not palidrome after replacing only 1 character.
#     queries[3] : substring = "abcd", could be changed to "abba" which is palidrome. Also this can be changed to "baab" first rearrange it "bacd" then replace "cd" with "ab".
#     queries[4] : substring = "abcda", could be changed to "abcba" which is palidrome.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length, queries.length <= 10^5`
#   * `0 <= queries[i][0] <= queries[i][1] < s.length`
#   * `0 <= queries[i][2] <= s.length`
#   * `s` only contains lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
"""
`solve!` - solve a single instance of a RandomODEProblem
"""
function solve!(problem::RandomODEProblem)
    problem.u = problem.A\problem.rhs
    return problem
end

"""
`solve_ensemble!` - solve an ensemble of RandomODEProblems
"""
function solve_ensemble!(problems::Vector)
    num_problems = length(problems)

    Threads.@threads for i ∈ 1:num_problems
        potential = sample_field(problems[i])
        problems[i].V = potential
        assemble_system!(problems[i])
        solve!(problems[i])
    end

    return problems
end


text:
@testset "LinearACEModel"  begin

##

using ACE, ACEbase
using Printf, Test, LinearAlgebra, ACE.Testing, Random
using ACE: evaluate, evaluate_d, SymmetricBasis, NaiveTotalDegree, PIBasis
using ACEbase.Testing: fdtest

##

@info("Basic test of LinearACEModel construction and evaluation")

# construct the 1p-basis
D = NaiveTotalDegree()
maxdeg = 6
ord = 3

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, D = D)

# generate a configuration
nX = 10
Xs = rand(EuclideanVectorState, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

φ = ACE.Invariant()
pibasis = PIBasis(B1p, ord, maxdeg; property = φ)
basis = SymmetricBasis(pibasis, φ)

BB = evaluate(basis, cfg)
c = rand(length(BB)) .- 0.5
naive = ACE.LinearACEModel(basis, c, evaluator = :naive)
standard = ACE.LinearACEModel(basis, c, evaluator = :standard)

##

# evaluate(naivemodel, cfg)
# evaluate(standard, cfg)
# ACE.grad_config(standard, cfg)

evaluate_ref(basis, cfg, c) = sum(evaluate(basis, cfg) .* c)

grad_config_ref(basis, cfg, c) = permutedims(evaluate_d(basis, cfg)) * c

grad_params_ref(basis, cfg, c) = evaluate(basis, cfg)

grad_params_config_ref(basis, cfg, c) = evaluate_d(basis, cfg)

for (fun, funref, str) in [
         (evaluate, evaluate_ref, "evaluate"),
         (ACE.grad_config, grad_config_ref, "grad_config"),
         (ACE.grad_params, grad_params_ref, "grad_params"),
         (ACE.grad_params_config, grad_params_config_ref, "grad_params_config"),
      ]
   @info("Testing `$str` for different model evaluators")
   for ntest = 1:30
      cgf = rand(EuclideanVectorState, B1p.bases[1], nX) |> ACEConfig
      c = rand(length(basis)) .- 0.5
      ACE.set_params!(naive, c)
      ACE.set_params!(standard, c)
      val = funref(basis, cfg, c)
      val_naive = fun(naive, cfg)
      val_standard = fun(standard, cfg)
      print_tf(@test( val ≈ val_naive ≈ val_standard ))
   end
   println()
end

##

end


text:
# Example of vOptSpecific, for the 2LAP
# when the numerical instance is read on a file, according to the LAP format
# July 2017

using vOptSpecific
m = load2LAP("2ap03.dat") # comment: set up the path to the file if required
z1, z2, S = vSolve( m )

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
"""
    read(
        npd :: NASAPrecipitationDataset,
        geo :: GeoRegion,
        dt  :: TimeType;
        lonlat :: Bool = false
    ) -> NCDataset           (if lonlat = false)
      -> lon, lat, NCDataset (if lonlat = true)

Reads a NASA Precipitation dataset specified by `npd` for a GeoRegion specified by `geo` at a date specified by `dt`.

Arguments
=========
- `npd` : a `NASAPrecipitationDataset` specifying the dataset details and date download range
- `geo` : a `GeoRegion` (see [GeoRegions.jl](https://github.com/JuliaClimate/GeoRegions.jl)) that sets the geographic bounds of the data array in lon-lat
- `dt`  : A specified date. The NCDataset retrieved may will contain data for the date, although it may also contain data for other dates depending on the `NASAPrecipitationDataset` specified by `npd`

Keyword Arguments
=================
- `lonlat` : if `true`, then return the longitude and latitude vectors for the dataset. Otherwise only the NCDataset type will be returned.
"""
function read(
	npd :: NASAPrecipitationDataset,
	geo :: GeoRegion,
    dt  :: TimeType;
    lonlat :: Bool = false
)

    pnc = npdfnc(npd,geo,dt)
    if !isfile(pnc)
        error("$(modulelog()) - The $(npd.lname) Dataset for the $(geo.regID) GeoRegion at Date $dt does not exist at $(pnc).  Check if files exist at $(npd.sroot) or download the files here")
    end
    @info "$(modulelog()) - Opening the $(npd.lname) NCDataset in the $(geo.regID) GeoRegion for $dt"
    pds = NCDataset(pnc)

    if !lonlat
          return pds
    else; return pds, pds["longitude"][:], pds["latitude"][:]
    end

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
using Perceptron
using Base.Test

# write your own tests here
@test Perceptron.accuracy([1,2,2], [1,2,2])==1

@test Perceptron.accuracy([1,2,2,3], [1,2,2,2])==0.75


text:
include("../src/includes.jl")

user = "guest"  # default is usually "guest"
password = "guest"  # default is usually "guest"
auth_params = compose_auth(user,password)
host="127.0.0.1"
conn_def = AmqpConnectionDef(host,auth_params)

amqp_conn_define!(conn_def)

queue!(chan -> source!(chan, "testQueue", String) |>
safe_map_ack!(String, msg -> msg != "bye" ? msg : error("bye"), (e) -> ()) |>
sink!(msg -> ()))

queue!(chan -> source!(chan, "test1Queue", Dict{String,Any}) |>
sink!(msg -> ()))

tasks = execute_queues!()

readline()

println(fetch(tasks[1]))
println(fetch(tasks[2]))


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
n, t=10, 100
X=genDataMatrix(n, t)
Y=genDataMatrix(n, t)
Cx=Symmetric((X*X')/t)
Cy=Symmetric((Y*Y')/t)
Cxy=(X*Y')/t
mC=mca(Cxy, simple=true)
@test Cxy≈mC.F[1]*mC.D*mC.F[2]'
D=mC.F[1]'Cxy*mC.F[2]
@test norm(D-Diagonal(D))+1. ≈ 1.

# Method (1) complex
Xc=genDataMatrix(ComplexF64, n, t)
Yc=genDataMatrix(ComplexF64, n, t)
Cxc=Symmetric((Xc*Xc')/t)
Cyc=Symmetric((Yc*Yc')/t)
Cxyc=(Xc*Yc')/t
mCc=mca(Cxyc, simple=true)
@test Cxyc≈mCc.F[1]*mCc.D*mCc.F[2]'
Dc=mCc.F[1]'Cxyc*mCc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.

# Method (2) real
mXY=mca(X, Y, simple=true)
D=mXY.F[1]'*Cxy*mXY.F[2]
@test norm(D-Diagonal(D))+1≈1.
@test mXY==mC

# Method (2) complex
mXYc=mca(Xc, Yc, simple=true)
Dc=mXYc.F[1]'*Cxyc*mXYc.F[2]
@test norm(Dc-Diagonal(Dc))+1. ≈ 1.
@test mXYc==mCc

# Method (3) real
# maximum covariance analysis of the average covariance and cross-covariance
k=10
Xset=[genDataMatrix(n, t) for i=1:k]
Yset=[genDataMatrix(n, t) for i=1:k]

m=mca(Xset, Yset)

# ... selecting subspace dimension allowing an explained variance = 0.5
m=mca(Xset, Yset; eVar=0.5)

# ... subtracting the mean from the matrices in Xset and Yset
m=mca(Xset, Yset; meanX=nothing, meanY=nothing, eVar=0.5)

# mca on the average of the covariance and cross-covariance matrices
# computed along dims 1
m=mca(Xset, Yset; dims=1, eVar=0.5)

# name of the filter
m.name

using Plots
# plot regularized accumulated eigenvalues
plot(m.arev)

# plot the original cross-covariance matrix and the rotated
# cross-covariance matrix
 Cmax=maximum(abs.(Cxy));
 h1 = heatmap(Cxy, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="Cxy");
 D=mC.F[1]'*Cxy*mC.F[2];
 Dmax=maximum(abs.(D));
 h2 = heatmap(D, clim=(0, Dmax), yflip=true, c=:amp, title="F[1]'*Cxy*F[2]");
 📈=plot(h1, h2, size=(700,300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigMCA.png")

# Method (3) complex
# maximum covariance analysis of the average covariance and cross-covariance

k=10
Xcset=[genDataMatrix(ComplexF64, n, t) for i=1:k]
Ycset=[genDataMatrix(ComplexF64, n, t) for i=1:k]

mc=mca(Xcset, Ycset)

# ... selecting subspace dimension allowing an explained variance = 0.5
mc=mca(Xcset, Ycset; eVar=0.5)

# ... subtracting the mean from the matrices in Xset and Yset
mc=mca(Xcset, Ycset; meanX=nothing, meanY=nothing, eVar=0.5)

# mca on the average of the covariance and cross-covariance matrices
# computed along dims 1
mc=mca(Xcset, Ycset; dims=1, eVar=0.5)

# name of the filter
mc.name


text:
using AxisArrays, IntervalSets

export probe_field, poynting, flux_surface_integral

"     probe_field(field::Field, component::Symbol, pt:Point)"
function probe_field(field::Field, component::Symbol, pt::Point)
    xind = AxisArrays.axisindexes(axes(field)[XX], atvalue(pt.x, atol=dx(field.grid)/2))
    yind = AxisArrays.axisindexes(axes(field)[YY], atvalue(pt.y, atol=dy(field.grid)/2))
    isa(field, FieldTM) && ~in(component, [:Ez, :Hx, :Hy]) && error("$component is invalid for the TM polarization. Valid options are :Ez, :Hx, :Hy")
    isa(field, FieldTE) && ~in(component, [:Hz, :Ex, :Ey]) && error("$component is invalid for the TE polarization. Valid options are :Hz, :Ex, :Ey")
    return field[xind, yind, component]
end

"    poynting(field::Field)"
function poynting(field::Field)
    if isa(field, FieldTM)
        Ez_x = grid_average(field[:, :, :Ez], x̂);
        Ez_y = grid_average(field[:, :, :Ez], ŷ);
        Sx = -0.5*real.(Ez_x.*conj(field[:, :, :Hy]));
        Sy =  0.5*real.(Ez_y.*conj(field[:, :, :Hx]));
        return Flux(field.grid, field.ω, Sx, Sy)
    end
    if isa(field, FieldTE)
        Hz_x = grid_average(field[:, :, :Hz], x̂);
        Hz_y = grid_average(field[:, :, :Hz], ŷ);
        Sx =  0.5*real.(field[:, :, :Ey].*conj(field.Hz_x));
        Sy = -0.5*real.(field[:, :, :Ex].*conj(field.Hz_y));
        return Flux(field.grid, field.ω, Sx, Sy)
    end
end

""""
    flux_surface_integral(field::Field, center::Point, width::Real, nrm::Direction)
Compute integral of flux normal to a surface defined by `center`, `width`, and
 `normal`.
"""
function flux_surface_integral(field::Field, center::Point, width::Real, normal::Direction)
    if normal == x̂
        if isa(field, FieldTM)
            xind = AxisArrays.axisindexes(AxisArrays.axes(field)[XX], atvalue(center.x, atol=dx(field.grid)/2))
            yint = (center.y-width)..(center.y+width)
            yind = AxisArrays.axisindexes(AxisArrays.axes(field)[YY], yint)
            Ez   = view(field.data, xind:xind+1, yind, :Ez)
            Ez_x̂ = grid_average(Ez, x̂)[1,:]
            Hy   = view(field.data, xind, yind, :Hy);
            return sum(-0.5*real.(Ez_x̂.*conj(Hy)))*dy(field.grid)
        end
        if isa(field, FieldTE)
            xind = AxisArrays.axisindexes(AxisArrays.axes(field)[XX], atvalue(pt.x, atol=dx(field.grid)/2))
            yint = (center.y-width)..(center.y+width)
            yind = AxisArrays.axisindexes(AxisArrays.axes(field)[YY], yint)
            Hz   = view(field.data, xind:xind+1, yind, :Hz)
            Hz_x̂ = grid_average(Hz, x̂)[1,:]
            Ey   = view(field.data, xind, yind, :Hy);
            return sum(-0.5*real.(Ey.*conj(Hz_x̂)))*dy(field.grid)
        end
    end
    if normal == ŷ
        if isa(field, FieldTM)
            error("ŷ normal TM calculation not yet implemented")
        end
        if isa(field, FieldTE)
            error("ŷ normal TE calculation not yet implemented")
        end
    end
end

# "    scattering_parameters(fields::Array{<:Field}, d::AbstractDevice)"
# function scattering_parameters(fields::Array{<:Field}, d::AbstractDevice)
#     @assert length(fields) == length(d.ω)
#     (ϵ₀, μ₀, c₀) = normalize_parameters(d)
#     S = zeros(Complex, length(d.ω), length(d.modes))
#     for i = 1:length(d.ω)
#         ω = d.ω[i];
#         for j = 1:length(d.modes)
#             mode = d.modes[j];
#             (β, Ej, indx, indy) = get_modes(d, mode.pol, ω, mode.neff, 1, mode.coor, mode.dir, mode.width)
#             Ec = fields[i].Ez[indx,indy] # computed field
#             Z = β/ω/ϵ₀
#
#             if j == 1
#                 numer = (Ec-Ej).*conj.(Ej)
#             else
#                 numer = Ec.*conj.(Ej)
#             end
#             denom = Ej.*conj.(Ej)
#             S[i, j] = sum(numer)/sum(denom)
#         end
#     end
#     return S
# end