text:
#=
Given an array of time intervals (start, end) for classroom lectures (possibly overlapping), find the minimum number of rooms required.

For example, given [(30, 75), (0, 50), (60, 150)], you should return 2.
=#

using Test
include("Solutions/problem21_number_of_rooms_for_lectures.jl")

@test number_of_rooms_required([(30, 75), (0, 50), (60, 150)]) == 2

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
## To get this as a string I pasted into vim and then did :%s/\n  //
str = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"

max = 1
for i = 1:(length(str)-13)
    product = 1
    substr = str[i:i+13]
    for j = 1:13
        digit = parse(Int, substr[j])
        product = product * digit
        if product > max
            max = product
        end
    end
end

println(max)

ext: jl

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
# ---
# title: 649. Dota2 Senate
# id: problem649
# author: Indigo
# date: 2020-12-11
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/dota2-senate/description/>
# hidden: true
# ---
#
# In the world of Dota2, there are two parties: the `Radiant` and the `Dire`.
#
# The Dota2 senate consists of senators coming from two parties. Now the senate
# wants to make a decision about a change in the Dota2 game. The voting for this
# change is a round-based procedure. In each round, each senator can exercise
# `one` of the two rights:
#
#   1. `Ban one senator's right`:
# A senator can make another senator lose **all his rights** in this and all the
# following rounds.
#
#   2. `Announce the victory`:
# If this senator found the senators who still have rights to vote are all from
# **the same party** , he can announce the victory and make the decision about
# the change in the game.
#
#
#
# Given a string representing each senator's party belonging. The character 'R'
# and 'D' represent the `Radiant` party and the `Dire` party respectively. Then
# if there are `n` senators, the size of the given string will be `n`.
#
# The round-based procedure starts from the first senator to the last senator in
# the given order. This procedure will last until the end of voting. All the
# senators who have lost their rights will be skipped during the procedure.
#
# Suppose every senator is smart enough and will play the best strategy for his
# own party, you need to predict which party will finally announce the victory
# and make the change in the Dota2 game. The output should be `Radiant` or
# `Dire`.
#
# **Example 1:**
#
#
#
#     Input: "RD"
#     Output: "Radiant"
#     Explanation: The first senator comes from Radiant and he can just ban the next senator's right in the round 1.
#     And the second senator can't exercise any rights any more since his right has been banned.
#     And in the round 2, the first senator can just announce the victory since he is the only guy in the senate who can vote.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: "RDD"
#     Output: "Dire"
#     Explanation:
#     The first senator comes from Radiant and he can just ban the next senator's right in the round 1.
#     And the second senator can't exercise any rights anymore since his right has been banned.
#     And the third senator comes from Dire and he can ban the first senator's right in the round 1.
#     And in the round 2, the third senator can just announce the victory since he is the only guy in the senate who can vote.
#
#
#
#
# **Note:**
#
#   1. The length of the given string will in the range [1, 10,000].
#
#
#
#
## @lc code=start
using LeetCode
function predict_party_victory(senate::String)
    len = length(senate)
    radiant, dire = Queue{Int}(), Queue{Int}()
    for (idx, ch) in enumerate(senate)
        (ch == 'R') ? enqueue!(radiant, idx) : enqueue!(dire, idx)
    end
    while !isempty(radiant) && !isempty(dire)
        r_front, d_fromt = dequeue!(radiant), dequeue!(dire)
        (r_front < d_fromt) ? enqueue!(radiant, r_front + len) : enqueue!(dire, d_fromt + len)

    end
    return isempty(dire) ? "Radiant" : "Dire"
end
## @lc code=end


text:
abstract type AbstractDCARule <: AbstractODRule end

struct DCA{B,R,T} <: AbstractDCARule
    rule::B
    ruleset::R
    states::Int
    radius::T
end

"""
    DCA(rule; states=2, radius=1)

Returns a ```DCA``` object given a specific rule, number of states and radius. The ruleset for the rule is computed and
stored in the struct as well.
"""
function DCA(rule;
    states=2,
    radius=1)

    ruleset = conversion(rule, states, radius)
    DCA(rule, ruleset, states, radius)
end

function (dca::DCA)(starting_array)

    nextgen = evolution(starting_array, dca.ruleset, dca.states, dca.radius)
end

function conversion(rule, states, radius::Int)

    rule_len = states^(2*radius+1)
    rule_bin = parse.(Int, split(string(rule, base=states), ""))
    rule_bin = vcat(zeros(typeof(rule_bin[1]), rule_len-length(rule_bin)), rule_bin)
    reverse!(rule_bin)
end

function conversion(rule, states, radius::Tuple)

    rule_len = states^(sum(radius)+1)
    rule_bin = parse.(Int, split(string(rule, base=states), ""))
    rule_bin = vcat(zeros(typeof(rule_bin[1]), rule_len-length(rule_bin)), rule_bin)
    reverse!(rule_bin)
end

function state_reader(neighborhood, states)

    parse(Int,join(convert(Array{Int},neighborhood)), base=states)+1 #ugly
end

function evolution(cell, ruleset, states, radius::Int)

    neighborhood_size = radius*2+1
    output = zeros(length(cell))
    cell = vcat(cell[end-neighborhood_size÷2+1:end], cell, cell[1:neighborhood_size÷2])

    for i=1:length(cell)-neighborhood_size+1
        output[i] = ruleset[state_reader(cell[i:i+neighborhood_size-1], states)]
    end

    output
end

function evolution(cell, ruleset, states, radius::Tuple)

    neighborhood_size = sum(radius)+1
    output = zeros(length(cell))#da qui in poi da modificare
    cell = vcat(cell[end-radius[1]+1:end], cell, cell[1:radius[2]])

    for i=1:length(cell)-neighborhood_size+1
        output[i] = ruleset[state_reader(cell[i:i+neighborhood_size-1], states)]
    end

    output
end


text:
struct GraphQLError <: Exception
    msg::String
    errs
end
GraphQLError(msg) = GraphQLError(msg, nothing)
GraphQLError(err::HTTP.StatusError) = GraphQLError("Request to server failed.", JSON3.read(err.response.body, GQLResponse{Any}))
function Base.showerror(io::IO, ex::GraphQLError)
    printstyled(io, "GraphQLError: " * ex.msg, color=Base.error_color())
    ex.errs !== nothing && display_errors(io, ex.errs)
end

print_heading(io::IO, heading) = printstyled(io, "\n\n$heading:\n\n", bold=true, color=Base.error_color())

function print_message(io::IO, message)
    print_heading(io, "Message")
    printstyled(io, message, color=Base.error_color())
end

function print_location(io::IO, location::GQLLocation)
    print_heading(io, "Location(s)")
    location = "line: " * string(location.line) * ", column: " * string(location.column)
    printstyled(io, location, color=Base.error_color())
end

function display_errors(io::IO, resp::GQLResponse)
    for err in resp.errors
        print_message(io, err.message)
        !isnothing(err.locations) && foreach(loc->print_location(io, loc), err.locations)
    end
end

struct GraphQLClientException <: Exception
    msg::String
end
function Base.showerror(io::IO, ex::GraphQLClientException)
    printstyled(io, "GraphQLClientException: " * ex.msg, color=Base.error_color())
end


text:
"""
Extract the principal component time series (PCs)
"""
function pcs(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
    scaling::Symbol = :none,
)
    if isnothing(n)
        n = length(eof.eigenvals)
    end

    out = eof.pcs[:, 1:n]

    if scaling == :divide
        # Divide by the square-root of the eigenvalue.
        for i = 1:n
            out[:, i] = out[:, i] ./ sqrt(eof.eigenvals[i])
        end
    elseif scaling == :multiply
        # Multiply by the square root of the eigenvalue.
        for i = 1:n
            out[:, i] = out[:, i] .* sqrt(eof.eigenvals[i])
        end
    elseif scaling != :none
        @error(
            "Could not understand scaling option. Select either :none, :divide, or :multiply"
        )
        return
    end

    return out
end

"""
Extract the empirical orthogonal functions (EOFs)
"""
function eofs(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
    scaling::Symbol = :none,
)
    if isnothing(n)
        n = length(eof.eigenvals)
    end

    out = eof.eofs[:, 1:n]

    if scaling == :divide
        # Divide by the square-root of the eigenvalue.
        for i = 1:n
            out[:, i] = out[:, i] ./ sqrt(eof.eigenvals[i])
        end
    elseif scaling == :multiply
        # Multiply by the square root of the eigenvalue.
        for i = 1:n
            out[:, i] = out[:, i] .* sqrt(eof.eigenvals[i])
        end
    elseif scaling != :none
        @error(
            "Could not understand scaling option. Select either :none, :divide, or :multiply"
        )
        return
    end

    return out
end

"""
Empirical orthogonal functions (EOFs) expressed as the
correlation between the principal component time series (PCs)
and the time series of the `eof` input *dataset* at each grid
point.
"""
function correlationmap(
    eof::EmpiricalOrthogonalFunction;
    outshape::Union{Nothing,Tuple{Int,Int}} = nothing,
    n::Union{Nothing,Int} = nothing,
)

    p = pcs(eof, n = n, scaling = :divide)

    p_center = Statistics.mean(p, dims = 2)

    pcs_cent = p .- p_center
    field_cent = eof.dataset

    pcs_std = Statistics.std(pcs_cent, dims = 1)
    field_std = field_std = Statistics.std(field_cent, dims = 1)

    div = size(pcs_cent, 1)

    cor = transpose(field_cent' * pcs_cent) / div

    cor = cor ./ (pcs_std' * field_std)

    if !isnothing(outshape)
        cor = reshape(cor, (n, outshape...))
    end

    return cor
end

"""
Empirical orthogonal functions (EOFs) expressed as the
covariance between the principal component time series (PCs)
and the time series of the `eof` input *dataset* at each grid
point.
"""
function covariancemap(
    eof::EmpiricalOrthogonalFunction;
    outshape::Union{Nothing,Tuple{Int,Int}} = nothing,
    n::Union{Nothing,Int} = nothing,
    ddof::Int = 1,
    scaling = :none,
)
    p = pcs(eof, n = n, scaling = scaling)

    p_center = Statistics.mean(p, dims = 2)

    pcs_cent = p .- p_center
    field_cent = eof.dataset

    div = size(pcs_cent, 1) - ddof

    cor = transpose(field_cent' * pcs_cent) / div

    if !isnothing(outshape)
        cor = reshape(cor, (n, outshape...))
    end

    return cor
end

"""
Extract the eigenvalues (decreasing variances) associated with each EOF.
"""
function eigenvalues(eof::EmpiricalOrthogonalFunction; n::Union{Nothing,Int} = nothing)
    if isnothing(n)
        n = length(eof.eigenvals)
    end
    return eof.eigenvals[1:n]
end

"""
Fractional EOF mode variances.
The fraction of the total variance explained by each EOF mode,
values between 0 and 1 inclusive.
"""
function variancefraction(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
)
    return eigenvalues(eof; n = n) / totalanomalyvar(eof)
end

"""
Total variance associated with the field of anomalies (the sum
of the eigenvalues).
"""
function totalanomalyvar(eof::EmpiricalOrthogonalFunction)
    return sum(eof.eigenvals)
end

"""
The method of North et al. (1982) is used to compute the typical
error for each eigenvalue. It is assumed that the number of
times in the input data set is the same as the number of
independent realizations. If this assumption is not valid then
the result may be inappropriate.
"""
function northtest(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
    scaled::Bool = false,
)

    nrecords = size(eof.dataset, 1)
    factor = sqrt(2.0 / nrecords)

    if scaled
        factor /= totalanomalyvar(eof)
    end

    return eigenvalues(eof; n = n) * factor
end

"""
Project a field onto the EOFs.
Given a data set, projects it onto the EOFs to generate a
corresponding set of pseudo-PCs.
"""
function projectfield(eof_arr::Array{<:Any,2}, field::Array{<:Any,2})
    # Project the data set onto the EOFs using a matrix multiplication.
    return field * eof_arr
end

function projectfield(
    eof::EmpiricalOrthogonalFunction,
    field::Array{<:Any,2};
    n::Union{Nothing,Int} = nothing,
    scaling::Symbol = :none,
)

    # extract out eofs for projection
    eofs_arr = eofs(eof, n = n, scaling = scaling)

    # check to make sure all data with valid index in EOF is valid in the input field
    @assert all(!ismissing(field[:, eof.valididx]))

    return projectfield(eofs_arr, field)
end

function projectfield(
    eof::EmpiricalOrthogonalFunction,
    field::Array{<:Any,3};
    n::Union{Nothing,Int} = nothing,
    scaling::Symbol = :none,
    timedim::Int = 3,
)

    geodims = [x for x in 1:3 if x != timedim]
    shape = size(dataset)
    spatial_shape = shape[geodims]
    time_shape = shape[timedim]

    if timedim == 2
        field = permutedims(field, (2, 1, 3))
    elseif timedim == 3
        field = permutedims(flattened, (3, 1, 2))
    end

    flat_field = reshape(field, (time_shape, prod(spatial_shape)))

    projected_pcs = projectfield(eof, flat_field; n = n, scaling = scaling)

    # reshape spatial field back to the original
    return reshape(projected_pcs, (spatial_shape..., time_shape))

end

"""
Reconstructed input data field based on a subset of EOFs.
"""
function reconstruct(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
    uncenter::Bool = true,
)

    eofs_arr = eofs(eof, n = n)
    ps = pcs(eof, n = n)

    rval = (ps * eofs_arr')

    if uncenter
        rval = rval .+ eof.center
    end

    return permutedims(rval, (2, 1))
end

function orthorotation(
    components::Array{<:Any,2};
    method::Symbol = :varimax,
    maxiter::Int = 50,
    tol::Float64 = 1e-6,
)
    nrows, ncols = size(components)

    rotation_matrix = I(ncols)
    var = 0
    i = 0

    for _ = 1:maxiter
        comp_rot = components * rotation_matrix
        if method == :varimax
            tmp = comp_rot * transpose(sum((comp_rot .^ 2), dims = 1) / nrows)
        elseif method == :quartimax
            tmp = 0
        end
        x = @. comp_rot^3 - tmp
        u, s, v = svd(transpose(components) * x)
        rotation_matrix = u * v
        var_new = sum(s)
        if var != 0 && var_new < var * (1 + tol)
            break
        end
        var = var_new
        i += 1
    end

    if i == maxiter
        @warn(
            "matrix rotation did not converge, try specifying a larger max iteration value for stable result"
        )
    end

    return components * rotation_matrix
end

"""
Apply orthogonal rotation to EOF data.
After rotation the original dataset will be projected on the rotated EOF
to create new PCs. Additionally new EOFs and PCs are ordered in decreasing
variance
"""
function orthorotation(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
    method::Symbol = :varimax,
    maxiter::Int = 50,
    tol::Float64 = 1e-6,
)
    components = eofs(eof, n = n)

    valididx = eof.valididx

    reofs = Array{Union{Missing,Float64}}(missing, size(components)...)
    reofs[valididx, :] = orthorotation(
        components[valididx, :];
        method = method,
        maxiter = maxiter,
        tol = tol,
    )

    projected_pcs = projectfield(reofs[valididx, :], eof.dataset[:, valididx])

    rotatedvar = Statistics.var(projected_pcs, dims = 1)[1, :]
    rotatedvarnorm =
        (rotatedvar ./ sum(rotatedvar)) .* sum(variancefraction(eof, n = n))

    sortedidx = sortperm(rotatedvarnorm)[end:-1:1]
    rotatedvar = rotatedvar[sortedidx]

    sorted_reofs = reofs[:, sortedidx]
    sorted_pcs = projected_pcs[:, sortedidx]

    return EmpiricalOrthogonalFunction(
        eof.dataset,
        eof.center,
        eof.valididx,
        sorted_reofs,
        rotatedvar,
        sorted_pcs,
    )
end

"""
Inplace orthorotation that will override the data structure of the input EOF object
"""
function orthorotation!(
    eof::EmpiricalOrthogonalFunction;
    n::Union{Nothing,Int} = nothing,
    method::Symbol = :varimax,
    maxiter::Int = 50,
    tol::Float64 = 1e-6,
)
    if isnothing(n)
        n = length(eof.eigenvals)
    else
        truncate!(eof, n)
    end

    valididx = eof.valididx

    eof.eofs[valididx, :] = orthorotation(
        eof.eofs[valididx, :];
        method = method,
        maxiter = maxiter,
        tol = tol,
    )

    eof.pcs[:, :] =
        projectfield(eof.eofs[valididx, :], eof.dataset[:, valididx])

    rotatedvar = Statistics.var(eof.pcs, dims = 1)[1, :]
    rotatedvarnorm =
        rotatedvarnorm =
            (rotatedvar ./ sum(rotatedvar)) .* sum(variancefraction(eof))

    sortedidx = sortperm(rotatedvarnorm)[end:-1:1]
    eof.eigenvals[:] = rotatedvar[sortedidx]

    eof.eofs[:, :] = eof.eofs[:, sortedidx]
    eof.pcs[:, :] = eof.pcs[:, sortedidx]

    return
end

"""
Truncate the EOF data structure to trim unneccesay modes
"""
function truncate(eof::EmpiricalOrthogonalFunction, n::Int)
    pcs_trunc = eof.pcs[:, 1:n]
    eofs_trunc = eof.eofs[:, 1:n]
    ev_trunc = eof.eigenvals[1:n]

    return EmpiricalOrthogonalFunction(
        eof.dataset,
        eof.center,
        eof.valididx,
        eofs_trunc,
        ev_trunc,
        pcs_trunc,
    )
end

"""
Apply truncate inplace on EOF data structure
"""
function truncate!(eof::EmpiricalOrthogonalFunction, n::Int)
    eof.pcs = eof.pcs[:, 1:n]
    eof.eofs = eof.eofs[:, 1:n]
    eof.eigenvals = eof.eigenvals[1:n]
    return
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

module FEMBasis

using LinearAlgebra

abstract type AbstractBasis end

include("subs.jl")
include("vandermonde.jl")

include("create_basis.jl")
export get_reference_element_coordinates, eval_basis!, eval_dbasis!

include("lagrange_segments.jl")
export Seg2, Seg3
include("lagrange_quadrangles.jl")
export Quad4, Quad8, Quad9
include("lagrange_triangles.jl")
export Tri3, Tri6, Tri7
include("lagrange_tetrahedrons.jl")
export Tet4, Tet10
include("lagrange_hexahedrons.jl")
export Hex8, Hex20, Hex27
include("lagrange_wedges.jl")
export Wedge6, Wedge15
include("lagrange_pyramids.jl")
export Pyr5

include("nurbs.jl")
include("nurbs_segment.jl")
export NSeg
include("nurbs_surface.jl")
export NSurf
include("nurbs_solid.jl")
export NSolid

include("math.jl")
export interpolate, interpolate!, jacobian
export grad, grad!
export BasisInfo

end


text:
using ForwardDiff
using DiffRules

import ForwardDiff: Dual, Partials, single_seed, partials, Chunk, Tag, seed!, value,
              gradient, valtype, extract_gradient!, derivative,extract_derivative,
              checktag, vector_mode_gradient, vector_mode_jacobian

const AMBIGUOUS_TYPES = (AbstractFloat, Irrational, Integer, Rational,
                         Real, RoundingMode, ComplexF64)

# The philosophy followed here is a bit different from that used by
# ForwardDiff. Here, the number of partials is limited to 1 or 2, even for
# arrays of vortex elements. For ComplexComplex cases, there are two partials,
# corresponding to differentiation with respect to the real and imaginary parts
# of the differentiation variable. For ComplexReal cases, there is one partial, since the
# differentiation variable is real.
# ComplexComplexDual is designated for derivatives of complex functions
# with respect to complex numbers
# ComplexRealDual is for derivatives of complex functions with respect to
# real numbers
# ComplexDual encompasses all complex dual types
const ComplexDual{T,V,N} = Complex{Dual{T,V,N}}
#const ComplexComplexDual{T,V} = ComplexDual{T,V,2}
#const ComplexRealDual{T,V} = ComplexDual{T,V,1}
#const RealComplexDual{T,V} = Dual{T,V,2}

# Constructor for ComplexDual based on real Partials for dreal and dimag
@inline function ComplexDual{T,V,N}(z::Number,dr::Partials{N,V},di::Partials{N,V}) where {T,N,V<:Real}
    zr, zi = reim(z)
    return Dual{T}(zr,dr) + im*Dual{T}(zi,di)
end

@inline ComplexDual{T}(z::Number,u::Partials{N,V},v::Partials{N,V}) where {T,N,V} =
    ComplexDual{T,V,N}(z,u,v)

# Constructor for ComplexDual based on complex Partials for dz and dz*
@inline function ComplexDual{T}(z::Number,dz::Partials{N,V},dzstar::Partials{N,V}) where {T,N,V<:Complex}
    dr, di = reim_partials(dz,dzstar)
    return ComplexDual{T}(z,dr,di)
end

@inline ComplexDual{T}(z::Number,dz::Complex,dzstar::Complex) where {T} =
    ComplexDual{T}(z,Partials(tuple(dz)),Partials(tuple(dzstar)))

@inline Base.one(::Type{<:ComplexDual{T}},z::Number) where T = ComplexDual{T}(z,one(z),zero(z))

@inline valtype(::ComplexDual{T,V,N}) where {T,V,N} = V
@inline valtype(::Type{ComplexDual{T,V,N}}) where {T,V,N} = V
@inline valtype(::Array{<:ComplexDual{T,V,N}}) where {T,V,N} = V
@inline valtype(::Type{<:Array{<:ComplexDual{T,V,N}}}) where {T,V,N} = V

# ==== extensions of Partials ==== #
@inline Base.conj(v::Partials{N,V}) where {N,V} = Partials{N,V}(conj.(v.values))
@inline Base.real(v::Partials{N,Complex{V}}) where {N,V} = Partials{N,V}(real.(v.values))
@inline Base.imag(v::Partials{N,Complex{V}}) where {N,V} = Partials{N,V}(imag.(v.values))

@inline Base.:*(x::Complex, partials::Partials) = partials*x

@inline function Base.:*(partials::Partials, x::Complex)
        return Partials(ForwardDiff.scale_tuple(partials.values, x))
end

# extract the partial derivatives from a complex dual number
@inline function ForwardDiff.partials(d::ComplexDual)
    dr, di = reim(d)
    return partials(dr), partials(di)
end

# assemble the derivatives d/dz and d/dz* from the partials
@inline function _derivs(dvrx,dvry,dvix,dviy)
    return 0.5*(dvrx+dviy) + 0.5*im*(dvix-dvry), 0.5*(dvrx-dviy) + 0.5*im*(dvix+dvry)
end

# Given a complex dual number, extract Partials of dz and dz*
@inline function dz_partials(d::ComplexDual{T,V,N}) where {T,V,N}
    dr, di = partials(d)
    tmp  = [_derivs(dr[i], dr[i+1], di[i],di[i+1]) for i in 1:2:N]
    return Partials(tuple(first.(tmp)...)), Partials(tuple(last.(tmp)...))
end

# For the dz and dz* partials with respect to the ith element
@inline Base.@propagate_inbounds function dz_partials(d::ComplexDual, i)
    dr, di = partials(d)
    return _derivs(dr[2i-1], dr[2i], di[2i-1],di[2i])
end

@inline Base.@propagate_inbounds dz_partials(::Type{T}, d::ComplexDual{T}, i...) where T = dz_partials(d, i...)

# For dealing with arrays of complex duals, return the dz and dz*
# partials as separate arrays
@inline dz_partials(d::AbstractVector{ComplexDual{T,V,N}}) where {T,V,N} =
  (tmp = dz_partials.(d); return first.(tmp), last.(tmp))

# For the dz and dz* partials with respect to the ith element
@inline Base.@propagate_inbounds dz_partials(d::AbstractVector{ComplexDual{T,V,N}},i) where {T,V,N} =
  (tmp = dz_partials.(d,i); return first.(tmp), last.(tmp))

@inline Base.@propagate_inbounds dz_partials(::Type{T}, d::AbstractVector{ComplexDual{T,V,N}}, i...) where {T,V,N} = dz_partials(d, i...)

# Given a complex dual, return dz and dzstar as Partial types. Ugly way of doing
# this
@inline function extract_gradient!(dz,dzstar,d::ComplexDual{T,V,N}) where {T,V,N}
    tmp, tmpstar = dz_partials(d)
    dz .= tmp
    dzstar .= tmpstar
    return dz,dzstar
end

# This is for making sure that the desired gradient has the same tag as the dual
@inline extract_gradient!(::Type{T},dz,dzstar,d::ComplexDual{T,V,N}) where {T,V,N} =
          extract_gradient!(dz,dzstar,d)

function extract_jacobian!(::Type{T}, dz::AbstractArray, dzstar::AbstractArray,
                ydual::AbstractArray{<:ComplexDual{T,V,N}}, n) where {T,V,N}
    @assert size(dz) == size(dzstar) "Inconsistent result sizes"

    dz_reshaped = reshape(dz, length(ydual), n)
    dzstar_reshaped = reshape(dzstar, length(ydual), n)
    for col in 1:size(dz_reshaped, 2), row in 1:size(dz_reshaped, 1)
        dz_reshaped[row, col], dzstar_reshaped[row, col] = dz_partials(T, ydual[row], col)
    end
    return dz, dzstar
end

"""
    ForwardDiff.extract_derivative(T,d::ComplexDual)

Given a complex dual value `d` and tag `T`, extract the derivatives d/dz
and d/dz* from the partials in `d`.
"""
@inline function extract_derivative(::Type{T},d::ComplexDual{T}) where T
   return dz_partials(d,1)
end

# If the function outputs a tuple of duals, as would be the case for
# nested derivatives, then we deal with each one by one and assemble an overall
# list of derivatives
@inline function extract_derivative(::Type{T},dlist::NTuple{N,ComplexDual{T}}) where {N,T}
   out = ()
   for d in dlist
        out = (out...,extract_derivative(T,d)...)
    end
    out
end

# Given complex Partials for d/dz and d/dz*, return partials for dr/dx,y and di/dx,y
@inline function reim_partials(dz::Partials{N,V},dzstar::Partials{N,V}) where {T,V<:Complex,N}
    drx, dix = reim(dz+dzstar)
    dry, diy = reim(im*dz-im*dzstar)
    return _splice_xy(drx,dry), _splice_xy(dix,diy)
end

# Given partials for the Create separate partials for the real part and the imaginary part.
# These are organized as drx[1],dry[1],drx[2],dry[2],... and dix[1],diy[1],dix[2],diy[2],...
@inline function _splice_xy(dx::Partials{N},dy::Partials{N}) where {N}
    return Partials(tuple([ifelse(mod(i,2) == 0,dy[max(i÷2,1)],dx[(i+1)÷2]) for i in 1:2N]...))
end

# Return the value of a complex dual
@inline function value(d::ComplexDual)
    dr, di = reim(d)
    return value(dr)+im*value(di)
end

@generated function construct_complex_seeds(::Type{Partials{N,V}}) where {N,V}
    return Expr(:tuple, Expr(:tuple,[:(single_seed(Partials{N,V}, Val{2*$i-1}())) for i in 1:N÷2]...),
                        Expr(:tuple,[:(single_seed(Partials{N,V}, Val{2*$i-0}())) for i in 1:N÷2]...))
end

function seed!(duals::AbstractArray{<:ComplexDual{T,V,M}}, x,
               rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}}) where {T,V,N,M}
    for i in 1:N
        #duals[i] = Dual{T,V,M}(real(x[i]), rseeds[i]) + im*Dual{T,V,M}(imag(x[i]), iseeds[i])
        duals[i] = ComplexDual{T,V,M}(x[i],rseeds[i],iseeds[i])
    end
    return duals
end

function seed!(duals::AbstractArray{<:ComplexDual{T,V,M}}, x,
               rseed::Partials{M,V} = zero(Partials{M,V}),iseed::Partials{M,V} = zero(Partials{M,V})) where {T,V,N,M}
    for i in eachindex(duals)
        duals[i] = ComplexDual{T,V,M}(x[i],rseed,iseed)
    end
    return duals
end

function dualize(::Type{T},x::AbstractArray{<:Complex{V}}) where {T,V}
    xdual = similar(x,ComplexDual{T,V,2*length(v)})
    seed!(xdual,x)
end

# =====  Configurations =====  #
"""
    ComplexGradientConfig(f,z)

Create a configuration for calculation of the gradient of a function
`f` with respect to complex array `z` by automatic differentiation.
"""
struct ComplexGradientConfig{T,V,N,D,M} <: ForwardDiff.AbstractConfig{N}
    rseeds::NTuple{N,Partials{M,V}}
    iseeds::NTuple{N,Partials{M,V}}
    duals::D
end

function ComplexGradientConfig(f::F,
                        x::AbstractArray{Complex{V}},
                        ::Chunk{N} = Chunk(x,length(x)),
                        ::T = Tag(f, V)) where {F,V,N,T}
    rseeds, iseeds = construct_complex_seeds(Partials{2N,V})
    duals = similar(x, ComplexDual{T,V,2N})
    return ComplexGradientConfig{T,V,N,typeof(duals),2N}(rseeds, iseeds, duals)
end

# =====  other seed functions ==== ##

"""
    seed(v::Number,cfg::ComplexGradientConfig)

Given a number `v`, create a dual-valued copy of `v` with
a unit seed. The use of complex duals (even if `v` is real-valued)
ensures correct treatment for gradient calculation, since some
calculations will be complex-valued.
"""
function seed(v::Complex,cfg::ComplexGradientConfig)
    vduals = copy(cfg.duals)
    seed!(vduals,[v],cfg.rseeds,cfg.iseeds)
    return vduals[1]
end

function seed(v::Real,cfg::ComplexGradientConfig)
    vduals = copy(cfg.duals)
    seed!(vduals,complex([v]),cfg.rseeds,cfg.iseeds)
    return real(vduals[1])
end

# ===== gradient ===== #

# Given a function f of z and a complex array z, evaluate the gradient wrt z at z
function gradient(f, z::AbstractArray{Complex{V}},
            cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, z)) where {T, V}
    #CHK && checktag(T, f, z)
    checktag(T, f, z)
    vector_mode_gradient(f, z, cfg)
end

function vector_mode_gradient(f::F, z, cfg::ComplexGradientConfig{T},
                evalfcn=ForwardDiff.vector_mode_dual_eval) where {F,T}
    ydual = evalfcn(f, z, cfg)
    dz = similar(z, Complex{valtype(ydual)})
    dzstar = similar(z, Complex{valtype(ydual)})
    return extract_gradient!(T, dz, dzstar,ydual)
end

function ForwardDiff.vector_mode_dual_eval(f::F, z, cfg::ComplexGradientConfig) where {F}
    zduals = cfg.duals
    seed!(zduals,z,cfg.rseeds,cfg.iseeds)
    return f(zduals)
end

# ===== jacobian ===== #

function jacobian(f, z::AbstractArray{Complex{V}},
            cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, z)) where {T, V}
    #CHK && checktag(T, f, z)
    checktag(T, f, z)
    vector_mode_jacobian(f, z, cfg)
end

function vector_mode_jacobian(f::F, z, cfg::ComplexGradientConfig{T,V,N},
            evalfcn=ForwardDiff.vector_mode_dual_eval) where {F,T,V,N}
    ydual = evalfcn(f, z, cfg)
    dz = similar(ydual, Complex{valtype(eltype(ydual))}, length(ydual), N)
    dzstar = similar(ydual, Complex{valtype(eltype(ydual))}, length(ydual), N)
    extract_jacobian!(T, dz, dzstar, ydual, N)
    return dz, dzstar
end

# ===== derivative ===== #

"""
    ForwardDiff.derivative(f,z::Complex)

Compute the derivative of function `f` with respect to `z` and `conj(z)`
"""
@inline function derivative(f::F, z::C) where {F,C<:Complex}
    T = typeof(ForwardDiff.Tag(f, C))
    # making f output complex ensures that it gets dispatched to our extract_derivative
    # for complex Duals rather than the native one in ForwardDiff
    return extract_derivative(T,complex.(f(one(ComplexDual{T},z))))
end

# ======== Extension of diff rules to complex functions ==========#

# preempt other function diffrules. The two entries correspond to d/dz and d/dz*
DiffRules.@define_diffrule Base.abs2(z) = :(conj($z)), :($z)
DiffRules.@define_diffrule Base.sqrt(z) = :(inv(2 * sqrt($z))), :(0)
DiffRules.@define_diffrule Base.log(z) = :(inv($z)), :(0)
DiffRules.@define_diffrule Base.:^(z,p) = :($p * ($z^($p - 1))), :(0)

DiffRules.@define_diffrule Base.abs(z) = :(0.5*conj($z)*inv(abs($z))), :(0.5*$z*inv(abs($z)))

# ======== Train functions to deal with complex duals ==========#

# Extend functions of a single complex argument to accept dual argument
macro extend_unary_dual_to_complex(fcn)
    M = :Base
    f = :($fcn)
    dfz, dfzstar = DiffRules.diffrule(M,f,:v)
    fdef = quote
        function $M.$f(z::ComplexDual{T}) where {T}
            v = value(z)
            dvz, dvzstar = dz_partials(z)

            # chain rule
            dfdz =     $dfz*dvz     + $dfzstar*conj(dvzstar)
            dfdzstar = $dfz*dvzstar + $dfzstar*conj(dvz)

            return ComplexDual{T}($M.$f(v),dfdz,dfdzstar)
        end
    end
    return esc(fdef)
end

# Extend functions of a single complex argument and another (non-dual) argument to accept
# dual argument in place of first argument
macro extend_binary_dual_to_complex(fcn)
    M = :Base
    f = :($fcn)
    dfz, dfzstar = DiffRules.diffrule(M,f,:v,:p)
    defs = quote end
    for R in AMBIGUOUS_TYPES
        expr = quote
            function $M.$f(z::ComplexDual{T},p::$R) where {T}
                v = value(z)
                dvz, dvzstar = dz_partials(z)

                # chain rule
                dfdz =     $dfz*dvz     + $dfzstar*conj(dvzstar)
                dfdzstar = $dfz*dvzstar + $dfzstar*conj(dvz)

                return ComplexDual{T}($M.$f(v,p),dfdz,dfdzstar)

            end
        end
        append!(defs.args, expr.args)
    end
    return esc(defs)
end

@extend_unary_dual_to_complex abs2
@extend_unary_dual_to_complex sqrt
@extend_unary_dual_to_complex log
@extend_unary_dual_to_complex abs

@extend_binary_dual_to_complex ^


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

### truncation
Base.trunc(dt::Date, p::Type{Year}) = Date(UTD(totaldays(year(dt), 1, 1)))
Base.trunc(dt::Date, p::Type{Quarter}) = firstdayofquarter(dt)
Base.trunc(dt::Date, p::Type{Month}) = firstdayofmonth(dt)
Base.trunc(dt::Date, p::Type{Day}) = dt

Base.trunc(dt::DateTime, p::Type{Year}) = DateTime(trunc(Date(dt), Year))
Base.trunc(dt::DateTime, p::Type{Quarter}) = DateTime(trunc(Date(dt), Quarter))
Base.trunc(dt::DateTime, p::Type{Month}) = DateTime(trunc(Date(dt), Month))
Base.trunc(dt::DateTime, p::Type{Day}) = DateTime(Date(dt))
Base.trunc(dt::DateTime, p::Type{Hour}) = dt - Minute(dt) - Second(dt) - Millisecond(dt)
Base.trunc(dt::DateTime, p::Type{Minute}) = dt - Second(dt) - Millisecond(dt)
Base.trunc(dt::DateTime, p::Type{Second}) = dt - Millisecond(dt)
Base.trunc(dt::DateTime, p::Type{Millisecond}) = dt

Base.trunc(t::Time, p::Type{Hour}) = Time(Hour(t))
Base.trunc(t::Time, p::Type{Minute}) = Time(Hour(t), Minute(t))
Base.trunc(t::Time, p::Type{Second}) = Time(Hour(t), Minute(t), Second(t))
Base.trunc(t::Time, p::Type{Millisecond}) = t - Microsecond(t) - Nanosecond(t)
Base.trunc(t::Time, p::Type{Microsecond}) = t - Nanosecond(t)
Base.trunc(t::Time, p::Type{Nanosecond})  = t

"""
    trunc(dt::TimeType, ::Type{Period}) -> TimeType

Truncates the value of `dt` according to the provided `Period` type.

# Examples
```jldoctest
julia> trunc(Dates.DateTime("1996-01-01T12:30:00"), Dates.Day)
1996-01-01T00:00:00
```
"""
Dates.trunc(::Dates.TimeType, ::Type{Dates.Period})

# Adjusters
"""
    firstdayofweek(dt::TimeType) -> TimeType

Adjusts `dt` to the Monday of its week.

# Examples
```jldoctest
julia> Dates.firstdayofweek(DateTime("1996-01-05T12:30:00"))
1996-01-01T00:00:00
```
"""
function firstdayofweek end

firstdayofweek(dt::Date) = Date(UTD(value(dt) - dayofweek(dt) + 1))
firstdayofweek(dt::DateTime) = DateTime(firstdayofweek(Date(dt)))

"""
    lastdayofweek(dt::TimeType) -> TimeType

Adjusts `dt` to the Sunday of its week.

# Examples
```jldoctest
julia> Dates.lastdayofweek(DateTime("1996-01-05T12:30:00"))
1996-01-07T00:00:00
```
"""
function lastdayofweek end

lastdayofweek(dt::Date) = Date(UTD(value(dt) + (7 - dayofweek(dt))))
lastdayofweek(dt::DateTime) = DateTime(lastdayofweek(Date(dt)))

"""
    firstdayofmonth(dt::TimeType) -> TimeType

Adjusts `dt` to the first day of its month.

# Examples
```jldoctest
julia> Dates.firstdayofmonth(DateTime("1996-05-20"))
1996-05-01T00:00:00
```
"""
function firstdayofmonth end

firstdayofmonth(dt::Date) = Date(UTD(value(dt) - day(dt) + 1))
firstdayofmonth(dt::DateTime) = DateTime(firstdayofmonth(Date(dt)))

"""
    lastdayofmonth(dt::TimeType) -> TimeType

Adjusts `dt` to the last day of its month.

# Examples
```jldoctest
julia> Dates.lastdayofmonth(DateTime("1996-05-20"))
1996-05-31T00:00:00
```
"""
function lastdayofmonth end

function lastdayofmonth(dt::Date)
    y, m, d = yearmonthday(dt)
    return Date(UTD(value(dt) + daysinmonth(y, m) - d))
end
lastdayofmonth(dt::DateTime) = DateTime(lastdayofmonth(Date(dt)))

"""
    firstdayofyear(dt::TimeType) -> TimeType

Adjusts `dt` to the first day of its year.

# Examples
```jldoctest
julia> Dates.firstdayofyear(DateTime("1996-05-20"))
1996-01-01T00:00:00
```
"""
function firstdayofyear end

firstdayofyear(dt::Date) = Date(UTD(value(dt) - dayofyear(dt) + 1))
firstdayofyear(dt::DateTime) = DateTime(firstdayofyear(Date(dt)))

"""
    lastdayofyear(dt::TimeType) -> TimeType

Adjusts `dt` to the last day of its year.

# Examples
```jldoctest
julia> Dates.lastdayofyear(DateTime("1996-05-20"))
1996-12-31T00:00:00
```
"""
function lastdayofyear end

function lastdayofyear(dt::Date)
    y, m, d = yearmonthday(dt)
    return Date(UTD(value(dt) + daysinyear(y) - dayofyear(y, m, d)))
end
lastdayofyear(dt::DateTime) = DateTime(lastdayofyear(Date(dt)))

"""
    firstdayofquarter(dt::TimeType) -> TimeType

Adjusts `dt` to the first day of its quarter.

# Examples
```jldoctest
julia> Dates.firstdayofquarter(DateTime("1996-05-20"))
1996-04-01T00:00:00

julia> Dates.firstdayofquarter(DateTime("1996-08-20"))
1996-07-01T00:00:00
```
"""
function firstdayofquarter end

function firstdayofquarter(dt::Date)
    y,m = yearmonth(dt)
    mm = m < 4 ? 1 : m < 7 ? 4 : m < 10 ? 7 : 10
    return Date(y, mm, 1)
end
firstdayofquarter(dt::DateTime) = DateTime(firstdayofquarter(Date(dt)))

"""
    lastdayofquarter(dt::TimeType) -> TimeType

Adjusts `dt` to the last day of its quarter.

# Examples
```jldoctest
julia> Dates.lastdayofquarter(DateTime("1996-05-20"))
1996-06-30T00:00:00

julia> Dates.lastdayofquarter(DateTime("1996-08-20"))
1996-09-30T00:00:00
```
"""
function lastdayofquarter end

function lastdayofquarter(dt::Date)
    y,m = yearmonth(dt)
    mm, d = m < 4 ? (3, 31) : m < 7 ? (6, 30) : m < 10 ? (9, 30) : (12, 31)
    return Date(y, mm, d)
end
lastdayofquarter(dt::DateTime) = DateTime(lastdayofquarter(Date(dt)))

# Temporal Adjusters
struct DateFunction
    f::Function
    # validate boolean, single-arg inner constructor
    function DateFunction(@nospecialize(f), dt::TimeType)
        isa(f(dt), Bool) || throw(ArgumentError("Provided function must take a single TimeType argument and return true or false"))
        return new(f)
    end
end
Base.show(io::IO, df::DateFunction) = println(io, df.f)

# Core adjuster
function adjust(df::DateFunction, start, step, limit)
    for i = 1:limit
        df.f(start) && return start
        start += step
    end
    throw(ArgumentError("Adjustment limit reached: $limit iterations"))
end

function adjust(func::Function, start; step::Period=Day(1), limit::Int=10000)
    return adjust(DateFunction(func, start), start, step, limit)
end

# Constructors using DateFunctions

"""
    Date(f::Function, y[, m, d]; step=Day(1), limit=10000) -> Date

Create a `Date` through the adjuster API. The starting point will be constructed from the
provided `y, m, d` arguments, and will be adjusted until `f::Function` returns `true`.
The step size in adjusting can be provided manually through the `step` keyword.
`limit` provides a limit to the max number of iterations the adjustment API will
pursue before throwing an error (given that `f::Function` is never satisfied).

# Examples
```jldoctest
julia> Date(date -> Dates.week(date) == 20, 2010, 01, 01)
2010-05-17

julia> Date(date -> Dates.year(date) == 2010, 2000, 01, 01)
2010-01-01

julia> Date(date -> Dates.month(date) == 10, 2000, 01, 01; limit = 5)
ERROR: ArgumentError: Adjustment limit reached: 5 iterations
Stacktrace:
[...]
```
"""
function Date(func::Function, y, m=1, d=1; step::Period=Day(1), limit::Int=10000)
    return adjust(DateFunction(func, Date(y, m, d)), Date(y, m, d), step, limit)
end

"""
    DateTime(f::Function, y[, m, d, h, mi, s]; step=Day(1), limit=10000) -> DateTime

Create a `DateTime` through the adjuster API. The starting point will be constructed from
the provided `y, m, d...` arguments, and will be adjusted until `f::Function` returns
`true`. The step size in adjusting can be provided manually through the `step` keyword.
`limit` provides a limit to the max number of iterations the adjustment API will
pursue before throwing an error (in the case that `f::Function` is never satisfied).

# Examples
```jldoctest
julia> DateTime(dt -> Dates.second(dt) == 40, 2010, 10, 20, 10; step = Dates.Second(1))
2010-10-20T10:00:40

julia> DateTime(dt -> Dates.hour(dt) == 20, 2010, 10, 20, 10; step = Dates.Hour(1), limit = 5)
ERROR: ArgumentError: Adjustment limit reached: 5 iterations
Stacktrace:
[...]
```
"""
DateTime(::Function, args...)

function DateTime(func::Function, y, m=1; step::Period=Day(1), limit::Int=10000)
    return adjust(DateFunction(func, DateTime(y, m)), DateTime(y, m), step, limit)
end
function DateTime(func::Function, y, m, d; step::Period=Hour(1), limit::Int=10000)
    return adjust(DateFunction(func, DateTime(y)), DateTime(y, m, d), step, limit)
end
function DateTime(func::Function, y, m, d, h; step::Period=Minute(1), limit::Int=10000)
    return adjust(DateFunction(func, DateTime(y)), DateTime(y, m, d, h), step, limit)
end
function DateTime(func::Function, y, m, d, h, mi; step::Period=Second(1), limit::Int=10000)
    return adjust(DateFunction(func, DateTime(y)), DateTime(y, m, d, h, mi), step, limit)
end
function DateTime(func::Function, y, m, d, h, mi, s; step::Period=Millisecond(1), limit::Int=10000)
    return adjust(DateFunction(func, DateTime(y)), DateTime(y, m, d, h, mi, s), step, limit)
end

"""
    Time(f::Function, h, mi=0; step::Period=Second(1), limit::Int=10000)
    Time(f::Function, h, mi, s; step::Period=Millisecond(1), limit::Int=10000)
    Time(f::Function, h, mi, s, ms; step::Period=Microsecond(1), limit::Int=10000)
    Time(f::Function, h, mi, s, ms, us; step::Period=Nanosecond(1), limit::Int=10000)

Create a `Time` through the adjuster API. The starting point will be constructed from the
provided `h, mi, s, ms, us` arguments, and will be adjusted until `f::Function` returns `true`.
The step size in adjusting can be provided manually through the `step` keyword. `limit`
provides a limit to the max number of iterations the adjustment API will pursue before
throwing an error (in the case that `f::Function` is never satisfied). Note that the default step
will adjust to allow for greater precision for the given arguments; i.e. if hour, minute, and second
arguments are provided, the default step will be `Millisecond(1)` instead of `Second(1)`.

# Examples
```jldoctest
julia> Dates.Time(t -> Dates.minute(t) == 30, 20)
20:30:00

julia> Dates.Time(t -> Dates.minute(t) == 0, 20)
20:00:00

julia> Dates.Time(t -> Dates.hour(t) == 10, 3; limit = 5)
ERROR: ArgumentError: Adjustment limit reached: 5 iterations
Stacktrace:
[...]
```
"""
Time(::Function, args...)

function Time(func::Function, h, mi=0; step::Period=Second(1), limit::Int=10000)
    return adjust(DateFunction(func, Time(h, mi)), Time(h, mi), step, limit)
end
function Time(func::Function, h, mi, s; step::Period=Millisecond(1), limit::Int=10000)
    return adjust(DateFunction(func, Time(h, mi, s)), Time(h, mi, s), step, limit)
end
function Time(func::Function, h, mi, s, ms; step::Period=Microsecond(1), limit::Int=10000)
    return adjust(DateFunction(func, Time(h, mi, s, ms)), Time(h, mi, s, ms), step, limit)
end
function Time(func::Function, h, mi, s, ms, us; step::Period=Nanosecond(1), limit::Int=10000)
    return adjust(DateFunction(func, Time(h, mi, s, ms, us)), Time(h, mi, s, ms, us), step, limit)
end

# Return the next TimeType that falls on dow
ISDAYOFWEEK = Dict(Mon => DateFunction(ismonday, Date(0)),
                   Tue => DateFunction(istuesday, Date(0)),
                   Wed => DateFunction(iswednesday, Date(0)),
                   Thu => DateFunction(isthursday, Date(0)),
                   Fri => DateFunction(isfriday, Date(0)),
                   Sat => DateFunction(issaturday, Date(0)),
                   Sun => DateFunction(issunday, Date(0)))

# "same" indicates whether the current date can be considered or not
"""
    tonext(dt::TimeType, dow::Int; same::Bool=false) -> TimeType

Adjusts `dt` to the next day of week corresponding to `dow` with `1 = Monday, 2 = Tuesday,
etc`. Setting `same=true` allows the current `dt` to be considered as the next `dow`,
allowing for no adjustment to occur.
"""
tonext(dt::TimeType, dow::Int; same::Bool=false) = adjust(ISDAYOFWEEK[dow], same ? dt : dt + Day(1), Day(1), 7)

# Return the next TimeType where func evals true using step in incrementing
"""
    tonext(func::Function, dt::TimeType; step=Day(1), limit=10000, same=false) -> TimeType

Adjusts `dt` by iterating at most `limit` iterations by `step` increments until `func`
returns `true`. `func` must take a single `TimeType` argument and return a [`Bool`](@ref).
`same` allows `dt` to be considered in satisfying `func`.
"""
function tonext(func::Function, dt::TimeType; step::Period=Day(1), limit::Int=10000, same::Bool=false)
    return adjust(DateFunction(func, dt), same ? dt : dt + step, step, limit)
end

"""
    toprev(dt::TimeType, dow::Int; same::Bool=false) -> TimeType

Adjusts `dt` to the previous day of week corresponding to `dow` with `1 = Monday, 2 =
Tuesday, etc`. Setting `same=true` allows the current `dt` to be considered as the previous
`dow`, allowing for no adjustment to occur.
"""
toprev(dt::TimeType, dow::Int; same::Bool=false) = adjust(ISDAYOFWEEK[dow], same ? dt : dt + Day(-1), Day(-1), 7)

"""
    toprev(func::Function, dt::TimeType; step=Day(-1), limit=10000, same=false) -> TimeType

Adjusts `dt` by iterating at most `limit` iterations by `step` increments until `func`
returns `true`. `func` must take a single `TimeType` argument and return a [`Bool`](@ref).
`same` allows `dt` to be considered in satisfying `func`.
"""
function toprev(func::Function, dt::TimeType; step::Period=Day(-1), limit::Int=10000, same::Bool=false)
    return adjust(DateFunction(func, dt), same ? dt : dt + step, step, limit)
end

# Return the first TimeType that falls on dow in the Month or Year
"""
    tofirst(dt::TimeType, dow::Int; of=Month) -> TimeType

Adjusts `dt` to the first `dow` of its month. Alternatively, `of=Year` will adjust to the
first `dow` of the year.
"""
function tofirst(dt::TimeType, dow::Int; of::Union{Type{Year}, Type{Month}}=Month)
    dt = of <: Month ? firstdayofmonth(dt) : firstdayofyear(dt)
    return adjust(ISDAYOFWEEK[dow], dt, Day(1), 366)
end

# Return the last TimeType that falls on dow in the Month or Year
"""
    tolast(dt::TimeType, dow::Int; of=Month) -> TimeType

Adjusts `dt` to the last `dow` of its month. Alternatively, `of=Year` will adjust to the
last `dow` of the year.
"""
function tolast(dt::TimeType, dow::Int; of::Union{Type{Year}, Type{Month}}=Month)
    dt = of <: Month ? lastdayofmonth(dt) : lastdayofyear(dt)
    return adjust(ISDAYOFWEEK[dow], dt, Day(-1), 366)
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
function make_conjugate_symmetric(array)
# conjSymmetric = makeConjugateSymmetric (array)
# array is an array such that imag(array(1)) == 0
# conjSymmetric is a conjugate symmetric array constructed from array

if length(size(array)) > 1
   error("array must be a vector")
end

conjSymmetric = vcat(array, flipdim(conj(array[2:end]),1)); # compatible with FFT
# if mod(length(conjSymmetric),2) ~= 1
#    error('array does not have an odd number of indices');
# end
#
# N = (length(conjSymmetric)-1)/2; # length(array) -   1
#
# first_half = conjSymmetric(1:N);
# second_half = conjSymmetric(N+2:end);
# if any(first_half - conj(flipud(second_half)))
#    error('array is not conjugate symmetric.')
# end

return conjSymmetric

end


text:
using HTTP
using Gumbo
using CSV
using DataFrames

sciNot = r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?"
numNot = r"[-+]?[0-9]*\.?[0-9]+"
shellNot = r"[K-Q] (XII|XI|X|IX|VIII|VII|VI|V|IV|III|II|I)?"

shellmap = (
    "K ", "L I", "L II", "L III",
    "M I", "M II", "M III", "M IV", "M V",
    "N I", "N II", "N III", "N IV", "N V",
    "N VI","N VII", "O I", "O II", "O III",
    "O IV", "O V", "P I", "P II", "P III",
)

function parseEdges(dd)
    lines = split(dd,"\n")
    res = Vector{Any}(missing, length(shellmap))
    for line in lines
        if !isnothing(findfirst("24 edges", line))
            continue
        end
        shells, energies = [], []
        for (sh, num) in zip( eachmatch(shellNot, line), eachmatch(sciNot,line))
            # println("Shell = $(sh.match), Number = $(num.match)")
            res[findfirst(n->n==sh.match,shellmap)] = 1000.0*parse(Float64, num.match)
        end
        @assert length(shells) == length(energies)
    end
    return res
end

"""
    scrape(z::Int)
Scrapes FFAST data from the NIST website by downloading by element the
associated HTML page and parsing the content to extract the various different
values.  Not exported as part of FFAST but provided to assist with downloading
updated values if they should become available.
"""
function scrape(z::Int)
    parseSN(dd) = parse(Float64, match(sciNot, dd).match)
    parseN(dd) = parse(Float64, match(numNot, dd).match)
    multiSN(dd) = [ parse(Float64, m.match) for m in eachmatch(sciNot, dd) ]
    html = HTTP.get("https://physics.nist.gov/cgi-bin/ffast/ffast.pl?Z=$(z)&Formula=&gtype=4&range=S&lower=0.001&upper=1000&density=&frames=no&htmltable=1");
    doc = parsehtml(String(html.body))
    headerNode = doc.root[2][1]
    if z <= 2
        shells = Vector{Any}(missing, length(shellmap))
        shells[1]= 1000.0*multiSN(headerNode[4][1].text)[2]
    else
        shells = parseEdges(headerNode[4][2][1].text)
    end
    data0 = Array{Any}(missing, 1, 7 + length(shells))
    data0[8:end] = shells
    data0[1,1] = parseN(headerNode[1][1][3].text) # A '=\n  54.93800 g mol'
    if z == 32
        data0[1,2] = parseSN(headerNode[2][7].text)
    else
        data0[1,2] = parseSN(headerNode[2][11].text) #  xsec ')  ×  9.12268E+01'
    end
    data0[1,3] = parseSN(headerNode[1][6].text) # density ')  = 7.4300'
    if z == 32
        data0[1,4] = parseSN(headerNode[3][13].text) # EV ')  ×  7.65960E+05'
    else
        data0[1,4] = parseSN(headerNode[3][17].text) # EV ')  ×  7.65960E+05'
    end
    if z <= 2
        rc = multiSN(headerNode[5][4].text)
        data0[1, 5], data0[1,6] = rc[length(rc) - 1:end] # # Reativistic correction '(H82,3/5CL) = -9.8594E-02, -6.1200E-02'
        data0[1, 7] = parseSN(headerNode[6][length(headerNode[6].children) - 3].text)
    else
        rc = multiSN(headerNode[4][6].text)
        data0[1, 5], data0[1, 6] = rc[length(rc) - 1:5] # # Relativistic correction '(H82,3/5CL) = -9.8594E-02, -6.1200E-02'
        data0[1, 7] = parseSN(headerNode[5][length(headerNode[5].children) - 3].text)
    end
    shSym = [ :K, :L1, :L2, :L3, :M1, :M2, :M3, :M4, :M5,
        :N1, :N2, :N3, :N4, :N5, :N6, :N7,
        :O1, :O2, :O3, :O4, :O5, :P1, :P2, :P3   ]
    df0 = DataFrame(data0, [ :A, :xsec, :density, :ev, :rc1, :rc2, :nt, shSym...])
    CSV.write("c:\\Users\\nritchie\\.julia\\dev\\FFAST\\data\\data[$(z)].csv", df0)

    tableNode = doc.root[2][8][1][1]
    data = zeros(Float64, length(tableNode.children) - 1, 8);
    for row in 2:length(tableNode.children)
        tr = tableNode[row]
        idx = 1
        for col in 1:length(tr.children)
            nums = multiSN(tr[col][1].text)
            for num in nums
                data[row - 1,idx] = num
                idx += 1
            end
        end
    end
    df1 = DataFrame(data)
    names!(df1, [:E,:f1,:f2,:μρpe,:μρci,:μρtot,:μρK,:λ])
    CSV.write("c:\\Users\\nritchie\\.julia\\dev\\FFAST\\data\\mac[$(z)].csv", df1)
end

for z in 1:92 scrape(z) end

df0=DataFrame( A=[],xsec=[],density=[],ev=[],rc1=[],rc2=[],nt=[],K=[],
    L1=[], L2=[], L3=[], M1=[], M2=[], M3=[], M4=[], M5=[],
    N1=[], N2=[], N3=[], N4=[], N5=[], N6=[], N7=[],
    O1=[], O2=[], O3=[], O4=[], O5=[], P1=[], P2=[], P3=[]  )
all = mapreduce(z->z == 0 ? df0 : CSV.read("c:\\Users\\nritchie\\.julia\\dev\\FFAST\\data\\data[$(z)].csv"),append!,0:92)
CSV.write("c:\\Users\\nritchie\\.julia\\dev\\FFAST\\data\\shelldata.csv",all)


text:
let
    tol = 1.0e-7

    # Test known position conversions
    geoc1 = [0.0, 0.0, 0.0]
    ecef1 = sGEOCtoECEF(geoc1)

    @test isapprox(ecef1[1], WGS84_a, atol=tol)
    @test isapprox(ecef1[2], 0, atol=tol)
    @test isapprox(ecef1[3], 0, atol=tol)

    geoc2 = [90.0, 0.0, 0.0]
    ecef2 = sGEOCtoECEF(geoc2, use_degrees=true)

    @test isapprox(ecef2[1], 0, atol=tol)
    @test isapprox(ecef2[2], WGS84_a, atol=tol)
    @test isapprox(ecef2[3], 0, atol=tol)

    geoc3 = [0.0, 90.0, 0.0]
    ecef3 = sGEOCtoECEF(geoc3, use_degrees=true)

    @test isapprox(ecef3[1], 0, atol=tol)
    @test isapprox(ecef3[2], 0, atol=tol)
    @test isapprox(ecef3[3], WGS84_a, atol=tol)

    # Test two-input format
    geoc = [0.0, 0.0]
    ecef = sGEOCtoECEF(geoc)

    @test isapprox(ecef[1], WGS84_a, atol=tol)
    @test isapprox(ecef[2], 0, atol=tol)
    @test isapprox(ecef[3], 0, atol=tol)

    geoc = [90.0, 0.0]
    ecef = sGEOCtoECEF(geoc, use_degrees=true)

    @test isapprox(ecef[1], 0, atol=tol)
    @test isapprox(ecef[2], WGS84_a, atol=tol)
    @test isapprox(ecef[3], 0, atol=tol)

    geoc = [0.0, 90.0]
    ecef = sGEOCtoECEF(geoc, use_degrees=true)

    @test isapprox(ecef[1], 0, atol=tol)
    @test isapprox(ecef[2], 0, atol=tol)
    @test isapprox(ecef[3], WGS84_a, atol=tol)

    # Test circularity
    geoc4 = sECEFtoGEOC(ecef1, use_degrees=true)
    geoc5 = sECEFtoGEOC(ecef2, use_degrees=true)
    geoc6 = sECEFtoGEOC(ecef3, use_degrees=true)

    @test isapprox(geoc4[1], geoc1[1], atol=tol)
    @test isapprox(geoc4[2], geoc1[2], atol=tol)
    @test isapprox(geoc4[3], geoc1[3], atol=tol)

    @test isapprox(geoc5[1], geoc2[1], atol=tol)
    @test isapprox(geoc5[2], geoc2[2], atol=tol)
    @test isapprox(geoc5[3], geoc2[3], atol=tol)

    @test isapprox(geoc6[1], geoc3[1], atol=tol)
    @test isapprox(geoc6[2], geoc3[2], atol=tol)
    @test isapprox(geoc6[3], geoc3[3], atol=tol)

    # Random point circularity
    geoc  = [77.875000,    20.975200,     0.000000]
    ecef  = sGEOCtoECEF(geoc, use_degrees=true)
    geocc = sECEFtoGEOC(ecef, use_degrees=true)
    @test isapprox(geoc[1], geocc[1], atol=tol)
    @test isapprox(geoc[2], geocc[2], atol=tol)
    @test isapprox(geoc[3], geocc[3], atol=tol)

    # Test Error Condition
    @test_throws ArgumentError sGEOCtoECEF([0.0,  90.1], use_degrees=true)
    @test_throws ArgumentError sGEOCtoECEF([0.0, -90.1], use_degrees=true)
end

let
    tol = 1.0e-7

    # Test known position conversions
    geod1 = [0, 0, 0]
    ecef1 = sGEODtoECEF(geod1)

    @test isapprox(ecef1[1], WGS84_a, atol=tol)
    @test isapprox(ecef1[2], 0, atol=tol)
    @test isapprox(ecef1[3], 0, atol=tol)

    geod2 = [90.0, 0.0, 0.0]
    ecef2 = sGEODtoECEF(geod2, use_degrees=true)

    @test isapprox(ecef2[1], 0, atol=tol)
    @test isapprox(ecef2[2], WGS84_a, atol=tol)
    @test isapprox(ecef2[3], 0, atol=tol)

    geod3 = [0, 90.0, 0]
    ecef3 = sGEODtoECEF(geod3, use_degrees=true)

    @test isapprox(ecef3[1], 0, atol=tol)
    @test isapprox(ecef3[2], 0, atol=tol)
    @test isapprox(ecef3[3], WGS84_a*(1.0-WGS84_f), atol=tol)

    # Test two input format
    geod = [0.0, 0.0]
    ecef = sGEODtoECEF(geod)

    @test isapprox(ecef[1], WGS84_a, atol=tol)
    @test isapprox(ecef[2], 0, atol=tol)
    @test isapprox(ecef[3], 0, atol=tol)

    geod = [90.0, 0.0]
    ecef = sGEODtoECEF(geod, use_degrees=true)

    @test isapprox(ecef[1], 0, atol=tol)
    @test isapprox(ecef[2], WGS84_a, atol=tol)
    @test isapprox(ecef[3], 0, atol=tol)

    geod = [0.0, 90.0]
    ecef = sGEODtoECEF(geod, use_degrees=true)

    @test isapprox(ecef[1], 0, atol=tol)
    @test isapprox(ecef[2], 0, atol=tol)
    @test isapprox(ecef[3], WGS84_a*(1.0-WGS84_f), atol=tol)

    # Test circularity
    geod4 = sECEFtoGEOD(ecef1, use_degrees=true)
    geod5 = sECEFtoGEOD(ecef2, use_degrees=true)
    geod6 = sECEFtoGEOD(ecef3, use_degrees=true)

    @test isapprox(geod4[1], geod1[1], atol=tol)
    @test isapprox(geod4[2], geod1[2], atol=tol)
    @test isapprox(geod4[3], geod1[3], atol=tol)

    @test isapprox(geod5[1], geod2[1], atol=tol)
    @test isapprox(geod5[2], geod2[2], atol=tol)
    @test isapprox(geod5[3], geod2[3], atol=tol)

    @test isapprox(geod6[1], geod3[1], atol=tol)
    @test isapprox(geod6[2], geod3[2], atol=tol)
    @test isapprox(geod6[3], geod3[3], atol=tol)

    geod  = [77.875000,    20.975200,     0.000000]
    ecef  = sGEODtoECEF(geod, use_degrees=true)
    geodc = sECEFtoGEOD(ecef, use_degrees=true)
    @test isapprox(geod[1], geodc[1], atol=tol)
    @test isapprox(geod[2], geodc[2], atol=tol)
    @test isapprox(geod[3], geodc[3], atol=tol)

    # Test Error Condition
    @test_throws ArgumentError sGEODtoECEF([0.0,  90.1], use_degrees=true)
    @test_throws ArgumentError sGEODtoECEF([0.0, -90.1], use_degrees=true)
end

let
    tol = 1.0e-8

    station_ecef = [0, R_EARTH, 0]

    R_ecef_enz = rECEFtoENZ(station_ecef, conversion="geocentric")
    R_enz_ecef = rENZtoECEF(station_ecef, conversion="geocentric")

    @test R_ecef_enz == R_enz_ecef'

    # State conversion
    epc  = Epoch(2018,1,1,12,0,0)
    oe   = [R_EARTH + 500e3, 1e-3, 97.8, 75, 25, 45]
    ecef = sECItoECEF(epc, sOSCtoCART(oe, use_degrees=true))

    station_ecef = sGEODtoECEF([-122.4056, 37.7716, 0.0], use_degrees=true)

    enz   = sECEFtoENZ(station_ecef, ecef)
    ecef2 = sENZtoECEF(station_ecef, enz)

    @test isapprox(ecef[1], ecef2[1], atol=tol)
    @test isapprox(ecef[2], ecef2[2], atol=tol)
    @test isapprox(ecef[3], ecef2[3], atol=tol)
    @test isapprox(ecef[4], ecef2[4], atol=tol)
    @test isapprox(ecef[5], ecef2[5], atol=tol)
    @test isapprox(ecef[6], ecef2[6], atol=tol)

    ecef         = sGEODtoECEF([-122.4, 37.78, 200.0],    use_degrees=true)
    station_ecef = sGEODtoECEF([-122.4056, 37.7716, 0.0], use_degrees=true)

    enz   = sECEFtoENZ(station_ecef, ecef, conversion="geocentric")
    ecef2 = sENZtoECEF(station_ecef, enz, conversion="geocentric")

    @test isapprox(ecef[1], ecef2[1], atol=tol)
    @test isapprox(ecef[2], ecef2[2], atol=tol)
    @test isapprox(ecef[3], ecef2[3], atol=tol)

    # Test ENZ Error Conditions
    @test_throws ArgumentError rECEFtoENZ([R_EARTH, 0.0])
    @test_throws ArgumentError rECEFtoENZ([R_EARTH, 0.0, 0.0], conversion="unknown")
    @test_throws ArgumentError rENZtoECEF([R_EARTH, 0.0])
    @test_throws ArgumentError sECEFtoENZ([R_EARTH, 0.0], [R_EARTH + 100.0, 0.0, 0.0])
    @test_throws ArgumentError sECEFtoENZ([R_EARTH, 0.0, 0.0], [R_EARTH + 100.0, 0.0])
    @test_throws ArgumentError sENZtoECEF([R_EARTH, 0.0], [0.0, 0.0, 0.0])
    @test_throws ArgumentError sENZtoECEF([R_EARTH, 0.0, 0.0], [0.0, 0.0])

    # Test length of return is 3
    enz = sECEFtoENZ(station_ecef, ecef[1:3], conversion="geocentric")
    @test length(enz) == 3
end

let
    tol = 1.0e-8
    station_ecef = [0, R_EARTH, 0]

    R_ecef_sez = rECEFtoSEZ(station_ecef, conversion="geocentric")
    R_sez_ecef = rSEZtoECEF(station_ecef, conversion="geocentric")

    @test isapprox(R_ecef_sez[1, 1], R_sez_ecef[1, 1], atol=tol)
    @test isapprox(R_ecef_sez[2, 1], R_sez_ecef[1, 2], atol=tol)
    @test isapprox(R_ecef_sez[3, 1], R_sez_ecef[1, 3], atol=tol)

    @test isapprox(R_ecef_sez[1, 2], R_sez_ecef[2, 1], atol=tol)
    @test isapprox(R_ecef_sez[2, 2], R_sez_ecef[2, 2], atol=tol)
    @test isapprox(R_ecef_sez[3, 2], R_sez_ecef[2, 3], atol=tol)

    @test isapprox(R_ecef_sez[1, 3], R_sez_ecef[3, 1], atol=tol)
    @test isapprox(R_ecef_sez[2, 3], R_sez_ecef[3, 2], atol=tol)
    @test isapprox(R_ecef_sez[3, 3], R_sez_ecef[3, 3], atol=tol)

    # State conversion
    epc  = Epoch(2018,1,1,12,0,0)
    oe   = [R_EARTH + 500e3, 1e-3, 97.8, 75, 25, 45]
    ecef = sECItoECEF(epc, sOSCtoCART(oe, use_degrees=true))

    station_ecef = sGEODtoECEF([-122.4056, 37.7716, 0.0], use_degrees=true)

    sez   = sECEFtoSEZ(station_ecef, ecef)
    ecef2 = sSEZtoECEF(station_ecef, sez)

    @test isapprox(ecef[1], ecef2[1], atol=tol)
    @test isapprox(ecef[2], ecef2[2], atol=tol)
    @test isapprox(ecef[3], ecef2[3], atol=tol)
    @test isapprox(ecef[4], ecef2[4], atol=tol)
    @test isapprox(ecef[5], ecef2[5], atol=tol)
    @test isapprox(ecef[6], ecef2[6], atol=tol)

    ecef         = sGEODtoECEF([-122.4, 37.78, 200.0],    use_degrees=true)
    station_ecef = sGEODtoECEF([-122.4056, 37.7716, 0.0], use_degrees=true)

    sez   = sECEFtoSEZ(station_ecef, ecef, conversion="geocentric")
    ecef2 = sSEZtoECEF(station_ecef, sez, conversion="geocentric")

    @test isapprox(ecef[1], ecef2[1], atol=tol)
    @test isapprox(ecef[2], ecef2[2], atol=tol)
    @test isapprox(ecef[3], ecef2[3], atol=tol)

    # Test SEZ Error Conditions
    @test_throws ArgumentError rECEFtoSEZ([R_EARTH, 0.0])
    @test_throws ArgumentError rECEFtoSEZ([R_EARTH, 0.0, 0.0], conversion="unknown")
    @test_throws ArgumentError rSEZtoECEF([R_EARTH, 0.0])
    @test_throws ArgumentError sECEFtoSEZ([R_EARTH, 0.0], [R_EARTH + 100.0, 0.0, 0.0])
    @test_throws ArgumentError sECEFtoSEZ([R_EARTH, 0.0, 0.0], [R_EARTH + 100.0, 0.0])
    @test_throws ArgumentError sSEZtoECEF([R_EARTH, 0.0], [0.0, 0.0, 0.0])
    @test_throws ArgumentError sSEZtoECEF([R_EARTH, 0.0, 0.0], [0.0, 0.0])

    # Test length of return is 3
    sez = sECEFtoSEZ(station_ecef, ecef[1:3], conversion="geocentric")
    @test length(sez) == 3
end

let
    # Test taken from Montenbruck and Gill Exercise 2.4
    # It mixes geodetic and geocentric coordinations in a strange way, but the
    # mixing is retained here for consistenty with the source material test
    epc = Epoch(1997, 1, 1, 0, 0, 0, tsys="UTC")
    oe  = [6378.137e3 + 960e3, 0, 97, 130.7, 0, 0]
    dt  = 15*60

    # Get Satellite position at 15 minutes
    n = mean_motion(oe[1], use_degrees=true)
    oe[6] += n*dt

    sat_eci = sOSCtoCART(oe, use_degrees=true)

    # Low precision ECEF transform
    d = (dt/86400.0 + mjd(epc) - 51544.5)
    # O = 280.4606 + 360.9856473*d
    O = 1.82289510683
    sat_ecef = Rz(0, use_degrees=false) * sat_eci[1:3]

    # Station coordinates
    station_ecef = sGEODtoECEF([48.0, 11.0, 0.0], use_degrees=true)

    # Compute enz and sez state
    enz   = sECEFtoENZ(station_ecef, sat_ecef, conversion="geocentric")
    sez   = sECEFtoSEZ(station_ecef, sat_ecef, conversion="geocentric")

    # Compute azimuth and elevation from topocentric coordinates
    azel_enz = sENZtoAZEL(enz, use_degrees=true)
    azel_sez = sSEZtoAZEL(sez, use_degrees=true)

    @test azel_enz[1] == azel_sez[1]
    @test azel_enz[2] == azel_sez[2]
    @test azel_enz[3] == azel_sez[3]
end

let
    # State conversion
    epc  = Epoch(2018,1,1,12,0,0)
    oe   = [R_EARTH + 500e3, 1e-3, 97.8, 75, 25, 45]
    ecef = sECItoECEF(epc, sOSCtoCART(oe, use_degrees=true))

    station_ecef = sGEODtoECEF([-122.4056, 37.7716, 0.0], use_degrees=true)

    enz = sECEFtoENZ(station_ecef, ecef)
    sez = sECEFtoSEZ(station_ecef, ecef)

    # Compute azimuth and elevation from topocentric coordinates
    azel_enz = sENZtoAZEL(enz, use_degrees=true)
    azel_sez = sSEZtoAZEL(sez, use_degrees=true)

    tol = 1e-7
    @test isapprox(azel_enz[1], azel_sez[1], atol=tol)
    @test isapprox(azel_enz[2], azel_sez[2], atol=tol)
    @test isapprox(azel_enz[3], azel_sez[3], atol=tol)
    @test isapprox(azel_enz[4], azel_sez[4], atol=tol)
    @test isapprox(azel_enz[5], azel_sez[5], atol=tol)
    @test isapprox(azel_enz[6], azel_sez[6], atol=tol)

end

let
    # Test Error Conditions
    enz = [0.0, 0.0, 100, 90.0, 0.0, 0.0]

    # Non-standard input length
    @test_throws ArgumentError sENZtoAZEL(enz[1:2])

    # Cant resolve azimuth without range information
    azel = sENZtoAZEL(enz[1:3])
    @test azel[1] == 0.0

    # Test ability to resolve azimuth ambiguity
    azel = sENZtoAZEL(enz)
    @test azel[1] != 0
    @test azel[2] != 0
    @test azel[3] != 0
end

let
    # Test Error Conditions
    sez = [0.0, 0.0, 100, 90.0, 0.0, 0.0]

    # Non-standard input length
    @test_throws ArgumentError sSEZtoAZEL(sez[1:2])

    # Cant resolve azimuth without range information
    azel = sSEZtoAZEL(sez[1:3])
    @test azel[1] == 0.0

    # Test ability to resolve azimuth ambiguity
    azel = sSEZtoAZEL(sez)
    @test azel[1] != 0
    @test azel[2] != 0
    @test azel[3] != 0
end


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

#Pkg.add("ForwardDiff")
#Pkg.add("SpecialFunctions")
#Pkg.update()
#Pkg.checkout("ForwardDiff")
#Pkg.status()
using Printf
using SpecialFunctions
using ForwardDiff

include("common.jl")

# input should be 1 dimensional
function logsumexp(x)
  mx = maximum(x)
  log.(sum(exp.(x .- mx))) .+ mx
end

@inbounds function gmm_objective(alphas,means,icf,x,wishart::Wishart)
  d = size(x,1)
  n = size(x,2)
  CONSTANT = -n*d*0.5*log(2 * pi)

  sum_qs = sum(icf[1:d,:],dims=1)
  slse = sum(sum_qs)
  Qs = [get_Q(d,icf[:,ik]) for ik in 1:k]
  main_term = zeros(eltype(alphas),1,k)

  slse = 0.
  for ix=1:n
    for ik=1:k
      main_term[ik] = -0.5*sum(abs2, Qs[ik] * (x[:,ix] - means[:,ik]))
    end
    slse += logsumexp(alphas + sum_qs + main_term)
  end

  CONSTANT + slse - n*logsumexp(alphas) + log_wishart_prior(wishart, sum_qs, Qs, icf)
end

# Read instance
if length(ARGS) < 5
  throw("Too few args")
end

dir_in = ARGS[1]
dir_out = ARGS[2]
fn = ARGS[3]
nruns_f = parse(Int,ARGS[4])
nruns_J = parse(Int,ARGS[5])
replicate_point = size(ARGS,1) >= 6 && ARGS[6] == "-rep"

fn_in = string(dir_in, fn)
fn_out = string(dir_out, fn)

alphas,means,icf,x,wishart = read_gmm_instance(string(fn_in,".txt"),replicate_point)
d = size(means,1)
k = size(means,2)
n = size(x,2)

# Objective
# Call once in case of precompilation etc
err = gmm_objective(alphas,means,icf,x,wishart)

tf = @elapsed for i in 1:nruns_f
  gmm_objective(alphas,means,icf,x,wishart)
end
tf = tf/nruns_f;
@printf "tf: %g\n" tf
#@printf "err: %f\n" err

# Gradient helper
function wrapper_gmm_objective(packed)
  alphas,means,icf = unpack(d,k,packed)
  gmm_objective(alphas,means,icf,x,wishart)
end

# Gradient
g = x-> ForwardDiff.gradient(wrapper_gmm_objective, x)

J = g(pack(alphas,means,icf))

tJ = @elapsed for i in 1:nruns_J
  g(pack(alphas,means,icf))
end
tJ = tJ/nruns_J;
@printf "tJ: %g\n" tJ
#println("J:")
#println(J)

name = "Julia"

write_J(string(fn_out,"_J_",name,".txt"),J)
write_times(string(fn_out,"_times_",name,".txt"),tf,tJ)

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
using LinearAlgebra: svd
function beyn(T::Function, A::AbstractMatrix, X::AbstractMatrix, nodes::Integer;
                c=complex(0.0,0.0), r=1.0)
    N, m₀ = size(X)
    if size(A, 1) != size(A, 2)
        error("Incorrect dimensions of A, must be square")
    elseif size(A,1) != N
        error("Incorrect dimensions of X₀, must match A")
    end
    res = zeros(m₀)
    temp = zeros(ComplexF64, N, m₀)

    Q₀, Q₁ = zeros(ComplexF64, N, m₀), zeros(ComplexF64, N, m₀)
    ZmA = similar(A, ComplexF64)
    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    for i=1:nodes
        ZmA .= T(r*exp(θ[i]*im)+c)
        ldiv!(temp, lu!(ZmA), X)
        Q₀ .+= temp .* exp(θ[i]*im)./nodes
        Q₁ .+= (r*exp(θ[i]*im)+c) .* temp .* exp(θ[i]*im)./nodes
    end
    S = svd!(Q₀)
    Aq = S.U' * Q₁ * S.V * inv(Diagonal(S.S))
    Λ, Xq = eigen!(Aq)
    X = S.U * Xq
    for i=1:m₀
            # res[i] = norm(T(Λ[i])*X[:,i])/norm(T(c)*X[:,i])
        res[i] = norm(T(Λ[i])*X[:,i])
    end
    p = sortperm(res)
    Λ[p], X[:, p], res[p]
end

function block_SS!(T, X::AbstractMatrix{ComplexF64}, nodes=2^4, moments=2;
    c=complex(0.0,0.0), r=1.0, debug=false)

    N, m₀ = size(X)
    K = moments*m₀
    X = Matrix(qr(X).Q)

    Temp = similar(X, ComplexF64)
    S = zeros(ComplexF64, N, (2*moments+1)*m₀)
    # Y = copy(X)
    Y = rand(ComplexF64,N, m₀)
    # Y = Matrix(I, N,N)
    l = size(Y,2)

    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    for i=1:nodes
        z = (r*exp(θ[i]*im)+c)
        Temp .= (T(z) \ X) .* (r*exp(θ[i]*im)/nodes)
        for j=1:2*moments+1
            S[:, (j-1)*m₀+1:j*m₀] .+= z^(j-1) .* Temp
        end
    end

    Q₀, Q₁ = zeros(ComplexF64, l*moments, K), zeros(ComplexF64, l*moments, K)

    for i=1:moments, j=1:moments
        Q₀[(i-1)*l+1:i*l, (j-1)*m₀+1:j*m₀] .= Y' * S[:, (i+j-1)*m₀+1:(i+j)*m₀]
        Q₁[(i-1)*l+1:i*l, (j-1)*m₀+1:j*m₀] .= Y' * S[:, (i+j)*m₀+1:(i+j+1)*m₀]
    end

    # D = svd!(Q₀)
    # Aq = D.U' * Q₁ * D.V * inv(Diagonal(D.S))
    # Λ, Xq = eigen!(Aq)
    # # X = D.U[1:N,:] * Xq
    # X = inv(Y') \ (D.U * Xq)[1:l, :]
    # # X = S[:,1:N] * D.U[1:N,:] * Xq
    # # X = S[:,1:K] * D.U * Xq

    V = svd(Q₀)
    n = min(count(V.S/V.S[1] .> 1e-13), K)

    H1 = V.U[:,1:n]' * Q₁ * V.V[:,1:n]
    H0 = V.U[:,1:n]' * Q₀ * V.V[:,1:n]

    Λ, Xq = eigen!(H1, H0)
    X = S[:,1:K] * V.V[:,1:n] *  Xq

    # H1 = V.U[:,1:n]' * Q₁ * V.V[:,1:n] * Diagonal(1 ./ V.S[1:n])
    # Λ, Xq = eigen!(H1)
    # X = S[:,1:K] * V.V[:,1:n] * Diagonal(1 ./ V.S[1:n]) * Xq
    # Λ = Λ .+ c

    res = zeros(n)

    for i=1:size(X, 2)
        X[:,i] ./= norm(X[:,i])
        res[i] = norm(T(Λ[i])*X[:,i])/norm(T(Λ[i]))
    end
    # p = sortperm(res)
    Λ, X, res
end


text:
using QuasiArrays, LinearAlgebra, LazyArrays, Test

@testset "ldiv" begin
    A = QuasiArray(randn(5,5), (0:0.5:2, 1:0.5:3))
    b = QuasiArray(randn(5), (0:0.5:2,))
    L = Ldiv(A, b)
    @test axes(L) == (axes(A,2),)
    @test similar(L) isa QuasiVector
    @test A\b == copy(L) == QuasiArray(parent(A)\parent(b), (1:0.5:3,))

    B = QuasiArray(randn(5,3), (0:0.5:2,Base.OneTo(3)))
    L = Ldiv(A, B)
    @test axes(L) == (axes(A,2),axes(B,2))
    @test similar(L) isa QuasiMatrix
    @test A\B == copy(L) == QuasiArray(parent(A)\parent(B), (1:0.5:3,Base.OneTo(3)))

    A = QuasiArray(randn(5,5), (Base.OneTo(5), 0:0.5:2))
    v = randn(5)
    @test A\v == QuasiArray(parent(A)\v, (0:0.5:2,))
    @test A' \ b == parent(A)' \ parent(b)
end


text:
__precompile__()

module ArtificialBeeColony

export ABC
export search!

abstract type TargetSystem end
struct TimeVariant <: TargetSystem end
struct TimeInvariant <: TargetSystem end

mutable struct Bee
    data :: Vector
    fitness :: AbstractFloat
    count :: Integer
end

function Bee(data::Vector)
    Bee(data, zero(eltype(data)), zero(Integer))
end

const Bees = Vector{Bee}

function Base.copyto!(t::Bee, s::Bee)
    copyto!(t.data, s.data)
    t.fitness = s.fitness
    t.count = s.count
    return t
end

function Base.copy(s::Bee)
    return Bee(copy(s.data), s.fitness, s.count)
end

struct ABC
    bees :: Bees
    best :: Bee
    init :: Function
    function ABC(N::Integer, init::Function)
        dim = length(init())
        bees = Bee[Bee(init()) for i=1:N]
        new(bees, Bee(zeros(dim)), init)
    end
end

@inline function update_fitness!(bee::Bee, g::Function)
    bee.fitness = fitness(g(bee.data))
end

@inline function update_fitness!(bees::Bees, g::Function)
    for bee in bees
        update_fitness!(bee, g)
    end
end

function update_bee!(bee::Bee, beem::Bee, g::Function)
    dim = length(bee.data)
    phi = 2.0*(rand()-0.5)
    j = rand(1:dim)

    val = phi*(bee.data[j] - beem.data[j])

    bee.data[j] += val
    fitnew = fitness(g(bee.data))

    if fitnew > bee.fitness
        bee.fitness = fitnew
        bee.count = 0
    else
        bee.data[j] -= val
        bee.count += 1
    end
end

function update_employed!(bees::Bees, g::Function, mode::TimeInvariant)
    N = length(bees)
    for bee in bees
        update_bee!(bee, bees[rand(1:N)], g)
    end
end

function update_employed!(bees::Bees, g::Function, mode::TimeVariant)
    N = length(bees)
    for bee in bees
        update_fitness!(bee, g)
        update_bee!(bee, bees[rand(1:N)], g)
    end
end

function roulette_select(bees::Bees)
    sf = sum(map(x->x.fitness, bees))
    r = rand()
    rs = zero(r)

    for bee in bees
        rs += bee.fitness/sf
        if r<=rs
            return bee
        end
    end
    return bees[end]
end

function update_outlook!(bees::Bees, g::Function, No::Integer, mode::TimeInvariant)
    N = length(bees)
    for i=1:No
        bee = roulette_select(bees)
        update_bee!(bee, bees[rand(1:N)], g)
    end
end

function update_outlook!(bees::Bees, g::Function, No::Integer, mode::TimeVariant)
    N = length(bees)
    for i=1:No
        bee = roulette_select(bees)
        update_fitness!(bee, g)
        update_bee!(bee, bees[rand(1:N)], g)
    end
end

function update_scout!(bees::Bees, g::Function, init::Function, limit::Integer, mode::TimeInvariant)
    bees_scout = filter(x->x.count >= limit, bees)
    for bee in bees_scout
        bee.data = init()
        bee.count = 0
    end
    update_fitness!(bees_scout, g)
end

function update_scout!(bees::Bees, g::Function, init::Function, limit::Integer, mode::TimeVariant)
    bees_scout = filter(x->x.count >= limit, bees)
    for bee in bees_scout
        bee.data = init()
        bee.count = 0
    end
end

function find_best(bees::Bees)
    best = bees[1]
    for bee in bees
        if bee.fitness > best.fitness
            best = bee
        end
    end
    return best
end

function update_best!(bees::Bees, bee_best::Bee, mode::TimeInvariant)
    bee_cand = find_best(bees)
    if bee_best.fitness < bee_cand.fitness
        copyto!(bee_best, bee_cand)
    end
end

function update_best!(bees::Bees, bee_best::Bee, mode::TimeVariant)
    copyto!(bee_best, find_best(bees))
end

function search!(abc::ABC, g::Function; epoch::Integer=1000, time_invariant::Bool=false)
    mode = time_invariant ? TimeInvariant() : TimeVariant()

    Ne = length(abc.bees)
    No = length(abc.bees)
    dim = length(abc.bees[1].data)
    limit = round(Integer, 0.1*dim*(Ne+No))

    if time_invariant
        update_fitness!(abc.bees, g)
    end
    for i=1:epoch
        update_employed!(abc.bees, g, mode)
        update_outlook!(abc.bees, g, No, mode)
        update_best!(abc.bees, abc.best, mode)
        update_scout!(abc.bees, g, abc.init, limit, mode)
    end
    return abc.best.data
end

function fitness(val::T) where T<:AbstractFloat
    if val >= zero(T)
        return one(T)/(one(T)+val)
    else
        return one(T)+abs(val)
    end
end

end


text:
using Test, FinancialMonteCarlo;
@show "LogNormalMixture Model"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = [0.2, 0.2];
lam = Float64[0.999999999];
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = LogNormalMixture(sigma, lam, Underlying(S0, d));
param_ = FinancialMonteCarlo.get_parameters(Model)

FinancialMonteCarlo.set_parameters!(Model, param_)
display(Model)

@test_throws(ErrorException, FinancialMonteCarlo.set_parameters!(Model, [1, 4]))

@show FwdPrice = pricer(Model, rfCurve, mc, FwdData);
@show EuPrice = pricer(Model, rfCurve, mc, EUData);
@show AmPrice = pricer(Model, rfCurve, mc, AMData);
@show BarrierPrice = pricer(Model, rfCurve, mc, BarrierData);
@show AsianPrice1 = pricer(Model, rfCurve, mc, AsianFloatingStrikeData);
@show AsianPrice2 = pricer(Model, rfCurve, mc, AsianFixedStrikeData);

@test abs(FwdPrice - 99.1078451563562) < toll
@test abs(EuPrice - 8.43005524824866) < toll
@test abs(AmPrice - 8.450489415187354) < toll
@test abs(BarrierPrice - 7.5008664470880735) < toll
@test abs(AsianPrice1 - 4.774451704549382) < toll

@show "Test LogNormalMixture Parameters"
eta = [0.2, 0.2]
lam11 = [0.9999]
@test_throws(ErrorException, LogNormalMixture(-0.2 * ones(3), 0.1 * ones(2), Underlying(S0, d)))
@test_throws(ErrorException, LogNormalMixture(0.2 * ones(3), -0.1 * ones(2), Underlying(S0, d)))
@test_throws(ErrorException, LogNormalMixture(0.2 * ones(3), ones(2), Underlying(S0, d)))
@test_throws(ErrorException, LogNormalMixture(0.2 * ones(3), 0.2 * ones(3), Underlying(S0, d)))


text:
# Activating Julia environment
using Pkg
Pkg.activate(".")

# Packages
using JuMP, Gurobi, Random, Distributions, LinearAlgebra, Plots

T = 25
S = 5
I = 6
Random.seed!(MersenneTwister(314)) # Do not change the seed.
function generate_data(I::Int64, S::Int64, T::Int64)
    mvdist = Multinomial(I, S) # Multinomial w.r.t. server allocation
    raws = Array{Any}(zeros(S)) # Generating tridiagonal multivariate normals for the servers
    for s = 1:S
        randos = 1/5*randperm(S)
        mat = diagm(0 => 1*ones(S), 1 => randos, -1 => randos)
        symmat = mat * mat'
        raws[s] = MvNormal(symmat)
    end
    D = zeros(I, S, T); # Generating demand data
    for t = 1:T
        D[:,:, t] += 2*rand(mvdist,I)'
        for i = 1:I
            for s = 1:S
                D[i, s, t] += abs(rand(raws[s])[1])
            end
        end
    end
    return D
end
D = generate_data(I,S,T);

heatmap(1:I, 1:S, var(D, dims=3)[:,:,1], xlabel = "Resources", ylabel = "Servers", title = "Time-variance of demand")

heatmap(1:I, 1:S, mean(D, dims=3)[:,:,1], xlabel = "Resources", ylabel = "Servers", title = "Time-mean of demand")

# Plotting different kinds of resource costs, i ∈ [1, ... , I]
V = 2*[1.5, 1.3, 0.8, 1.3, 1.2, 0.5]; # expansion costs
F = [1.5, 1.3, 0.8, 1.3, 1.2, 0.9]; # fixed costs
C = [0.4, 0.5, 0.6, 0.7, 0.8, 0.7]; # reallocation costs
# bardata = hcat(V, F, C)
# sx = repeat(["expansion", "fixed", "reallocation"], inner = 6)
# nam = repeat("Resource " .* string.(1:I), outer = 3)
# groupedbar(nam, bardata, group = sx, ylabel = "Costs",
#         title = "Costs for each resource")

# OPTIMIZATION MODEL HERE.
m = JuMP.Model(Gurobi.Optimizer)
@variable(m, r[i=1:I, s=1:S] >= 0)
@variable(m, e[i=1:I, s=1:S, t=1:T] >= 0)
@variable(m, h[s=1:S, t =1:T] >= 0)
@variable(m, u[i=1:I, s1=1:S, s2=1:S, t=1:T] >= 0)

# Demand uncertainty is described by:
# @uncertain(m, d[i=1:I, s=1:S, t =1:T])
# for s=1:S
#     @constraint(m, norm(d[:, s, :], 1) <= Γ)
# end
# Which is a budget uncertainty with ρ = 1 for the ∞ norm, and γ bounding the 1-norm.
# ρ = 0
# Γ = 0
ρ = 1
Γ = sqrt(2log(1/0.05))*sqrt(I*T)

@objective(m, Min, T*sum(sum(F .* r[:, s]) for s=1:S) +
                sum(sum(sum(V .* e[:, s, t]) for s=1:S) for t=1:T) +
                sum(sum(sum(sum(C .* u[:,s1,s2,t]) for s1=1:S) for s2=1:S) for t=1:T));

# Using the robust counterpart here
for i = 1:I
    for s = 1:S
        for t = 1:T
            y = @variable(m, [j=1:I, k=1:T])
            abs_y = @variable(m, [j=1:I, k=1:T])
            max_diff = @variable(m)
            @constraint(m, abs_y .≥ y)
            @constraint(m, abs_y .≥ -y)
            for j = 1:I # The right hand side uncertainty can be tricky!
                for k = 1:T
                    if j == i && k == t
                        @constraint(m, max_diff ≥ (-1 + y[j,k]))
                        @constraint(m, max_diff ≥ -(-1 + y[j,k]))
                    else
                        @constraint(m, max_diff ≥ y[j,k])
                        @constraint(m, max_diff ≥ -y[j,k])
                    end
                end
            end
            @constraint(m, D[i,s,t] + ρ*sum(abs_y) + Γ*max_diff + sum(u[i, s2, s, t] for s2 = 1:S) <= r[i,s] + e[i,s,t] + sum(u[i, s, s2, t] for s2 = 1:S))
        end
    end
end

for s = 1:S
    @constraint(m, h[s, 1] == sum(e[:,s,1] ./3))
    @constraint(m, h[s, 2] == sum(e[:,s,2] ./3) + sum(e[:,s,1] ./3))
    @constraint(m, h[s, 1] <= 1)
    @constraint(m, h[s, 2] <= 1)
    for t = 3:T
        @constraint(m, h[s ,t] == sum(e[:,s,t-2] ./3) + sum(e[:,s,t-1] ./3) + sum(e[:,s,t] ./3))
        @constraint(m, h[s, t] <= 1)
    end
end
optimize!(m)

# Fixed capacities plot
plt1 = heatmap(value.(r)', xlabel = "Resources", ylabel = "Servers", title = "Fixed capacities")

# Time-mean of expansions plot
plt2 = heatmap(mean(value.(e), dims=3)[:,:,1]', xlabel = "Resources",
        ylabel = "Servers", title = "Time-mean of expansions")

# Time-variance of expansions
plt3 = heatmap(var(value.(e), dims=3)[:,:,1]', xlabel = "Resources",
ylabel = "Servers", title = "Time-variance of expansions")

# Time-mean of job transfers out
transfers_out = zeros(I,S,T); # Computing the transfers out of each server
[transfers_out[i,s,t] = sum(value.(u)[i, s, :, t]) for i=1:I, s = 1:S, t = 1:T];
plt4 = heatmap(mean(transfers_out, dims=3)[:,:,1]', xlabel = "Resources",
        ylabel = "Servers", title = "Time-mean of job transfers out")

# Time-variance of job transfers out
plt5 = heatmap(var(transfers_out, dims=3)[:,:,1]', xlabel = "Resources",
ylabel = "Servers", title = "Time-variance of job transfers out")

# Plots of temperature
temps = value.(h)
plt6 = plot(1:T, temps[1,:], label=1)
for s=2:S
    plot!(1:T, temps[s,:], label=s, title = "Server temperatures", xlabel = "Time period (t)", ylabel = "Temperature",
        legend = :bottomright)
end

# Displaying
display(plt6)


text:
using Test

m = test_page_model()
include("../src/components/ClimateTemperature.jl")

climatetemperature = add_comp!(m, ClimateTemperature)

climatetemperature[:y_year_0] = 2008.
climatetemperature[:y_year] = Mimi.dim_keys(m.md, :time)

climatetemperature[:ft_totalforcing] = readpagedata(m, "test/validationdata/ft_totalforcing.csv")
climatetemperature[:fs_sulfateforcing] = readpagedata(m, "test/validationdata/fs_sulfateforcing.csv")

p = load_parameters(m)
set_leftover_params!(m, p)

##running Model
run(m)

rt_g = m[:ClimateTemperature, :rt_g_globaltemperature]
rt_g_compare = readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv")

@test rt_g ≈ rt_g_compare rtol=1e-5


text:
# Test routines for KV_2010.jl

## Mortality probability:
# Visual check
scatter(t,ξ)

## Check κf
plot([0:39],κ./maximum(κ), ylims=(0,1))

## Net labour income function Y()
function linreg(x,y)
    x = [x ones(length(y))]
    x'*x \ x'*y
end

Random.seed!(1234)
plot(labor_income(1,T_ret, 0., 0., 0., κ; y_0 = 0.0)')
any(labor_income(1,T_ret, 0., 0., 0., κ; y_0 = 0.0) .== 1.) != false

Random.seed!(1234)
testy = labor_income(10, T_ret, 0., σ_η, 0., κ; y_0 = 0.0)'
plot(testy)
i = 7
linreg(testy[i,1:end-1],testy[i,2:end])[2]

Random.seed!(1234)
testy = labor_income(10,T_ret, σ_ε, 0., 0., κ; y_0 = 0.0)'
plot!(testy)
i = 6
linreg(testy[i,1:end-1],testy[i,2:end])[2]

## Gross earnings Y_tilde_fn()
Random.seed!(1234)
N = 1000
Y_l = Y(N, T_ret, σ_ε, σ_η, σ_z0, κ)
i = rand(1:N)
plot(Y_l[i,:])
τ_s = 1.45
Y_tilde = Y_tilde_fn.(Y_l,τ_s)
plot!(Y_tilde[i,:])
scatter(Y_tilde[i,:], 1. .- Y_l[i,:] ./ Y_tilde[i,:], ylabel="tax rate", xlabel="gross income", legend=false)
scatter(Y_tilde, 1. .- Y_l ./ Y_tilde, ylabel="tax rate", xlabel="gross income", legend=false)

## Calibration of τ, function G()
Random.seed!(1234)
N = 100
Y_l = Y(N, T_ret, σ_ε, σ_η, σ_z0, κ)
sol = G(Y_l, τ_s_0 = 0.1, Kp = 10, Kd = 1., tol = 1E-3, max_iter = 100, nonconvergence_message = true)
sol.err
sol.iters
sol.τ_s

## Test P(Y_tilde)
## using Y_l instead of Y_tilde for simplicity
Random.seed!(1234)
N = 10000
Y_l = Y(N, T_ret, σ_ε, σ_η, σ_z0, κ)

## Bendpoints
Y_ave = mean(Y_l,dims=2) #lifetime average income
y_tilde_l_ave = mean(Y_ave)
bp1 = 0.18*y_tilde_l_ave
bp2 = 1.1*y_tilde_l_ave
Y_l_SS = P(Y_l)
scatter(Y_ave,Y_l_SS)
vline!([bp1 bp2])
xlims!(0.,1.)
ylims!(0.,0.4)

## Test average income gets 45% of earnings
y_tilde_l_ave = mean(sum(Y_l,dims=2)) #average lifetime net labour earnings
ind_mean = findmin(abs.( sum(Y_l,dims=2) .- y_tilde_l_ave ))
mean_ind = ind_mean[2][1]
sum(Y_l[mean_ind,:])
mean(Y_l[mean_ind,:])
Y_l_SS[mean_ind] - 0.45*mean(Y_l[mean_ind,:]) < 1E-4
plot(Y_l[mean_ind,:])
Y = hcat(Y_l, repeat(Y_l_SS, 1, T-(T_ret-1)))
plot(Y[mean_ind,:], ylims=(0,:))

## Test income generation
Random.seed!(1234)
Y_l = Y(N, T_ret, σ_ε, σ_η, σ_z0, κ)
(Y_tilde, τ_s) = G(Y_l, τ_s_0 = 1.45)
Y_tilde_SS = P(Y_tilde)
Y_SS = Y_tilde_SS - τ.(0.85.*Y_tilde_SS, τ_s)
Y = hcat(Y_l, repeat(Y_SS, 1, T-(T_ret-1)) )
for num = 1:10
    num == 1 ? plot() : nothing
    i = rand(1:N)
    display(plot!(Y[i,:], ylims=(0,8), label=i, title="Net income"))
end

# Discretization
# Test: Simulate the markov chain process and compare to simulated income process
Y_l = Y(N, T_ret, σ_ε, σ_η, σ_z0, κ)
scatter(1:T_ret-1, var(Y_l', dims = 2), label="Y_l", legend=:bottomright)
# Does this match graphs from literature??

# Generate z0's
z_0 = σ_z0 .* randn(N)
var(exp.(z_0))
# Find index of Markov chain corresponding to each z_0
ind_0 = [findmin( abs.(mcs_per[1].state_values .- z_0) )[2] for z_0 in z_0]
z = zeros(length(mcs_per)+1,N)
[z[:,i] = simulate_ns(mcs_per, ind_0[i]) for i in 1:N]
var(exp.(z[1,:]))
z = z[2:end,:]
ε = zeros(T_ret+1,N)
[ε[:,i] = simulate(mc_tran, T_ret+1, init = round(Int,19/2)) for i in 1:N]
ε = ε[2:end,:]
Κ = repeat(κ, 1, N)
Y_mc = exp.(Κ + z + ε)
scatter!(1:T_ret, var(Y_mc, dims = 2), label="Y_mc")
(var(Y_l', dims=2) - var(Y_mc, dims=2)[1:end-1]) ./ var(Y_l', dims=2)*100 |> plot
# Errors are small and don't change systematically over time

scatter(1:T_ret-1, mean(Y_l', dims = 2), label="Y_l", legend=:bottomright)
scatter!(1:T_ret, mean(Y_mc, dims = 2), label="Y_mc")
# Average are very close

# Check initial values are working
scatter(z_0,z_0 - z[1,:])
vline!(mcs_per[1].state_values)

# Boneyard

# Begin from last period
for t in reverse(T_ret:T-1) #reverse(t)
    if t < T_ret
        Ec = zeros(Float64,a_gridpoints,z_gridpoints)
        for (a_t1_ind, a_t1) in enumerate(A_vals), (z_ind, z) in enumerate(mcs_z[t].state_values)
            A_t2 = extrapolate(interpolate(A[t+2], BSpline(Linear())), 0. )(a,z)
            for (ε_ind, ε) in enumerate(mc_ε.state_values), (z_t1_ind, z_t1) in enumerate(mcs_z[t].state_values)
                Ec[a_ind, z_ind] += mcs_z[t].p[z_ind, z_t1_ind] * mc_ε.p[1,ε_ind] * ( (1+r).*a_t1 + exp(κ[t+1] + z_t1 + ε) - A_t2 )^(-γ)
            end
        end
        C[t] = ( β*(1+r) .* Ec).^(-1/γ)
        At = 1/(1+r) .* (C[t] - repeat(PY_tilde_vals', a_gridpoints, 1) + A_grid_l) # ******how to account for all possible combos of z,ε??
        A[t+1] = grid_shift(invert_rule(At, A_vals), A_vals, A_vals, Z_vals_old, Z_vals_new) # ********

    elseif t < (T-1) && t >= T_ret # Be careful with crossover year!
        C_t1 = (1+r).*A_grid_ret + PY_tilde_grid - A[t+2]
        C_t = ( β*ξ[t+1]/ξ[t]*(1+r) )^(-1/γ) .* C_t1
        At = 1/(1+r) .* (C_t - PY_tilde_grid + A_grid_ret)
        A[t+1] = invert_rule(At, A_vals)

    elseif t == T-1
        A[t+1] = zeros(Float64, a_gridpoints, PY_tilde_gridpoints)
        # no need to shift because all zeros
    end
end

# Check 3d interpolation
A = zeros(3,4,2)
A[:,:,1] = [1 2 3 4; 5 6 7 8; 9 10 11 12]
A[:,:,2] = 2*[1 2 3 4; 5 6 7 8; 9 10 11 12]
li = interpolate(A, BSpline(Linear()))
li(1,2.5,1)
mean(A,dims=3)

A = [1 2 3 4; 5 6 7 8]
fill(A, 1,1,2)

A
function testfn()
    A = Float64[]
    for a in 1:3, b in 2:5, c in 1:2
        push!(A, a + b + c)
    end
    return A
end

testfn()

A = [a + b + c for a in 1:3, b in 2:5, c in 1:2 ]
A = [println(a) for a in 1:3, b in 2:5, c in 1:2 ]
A = [(a,b,c) for a in 1:3, b in 4:5, c in 6:7]

# Ec = zeros(Float64,a_gridpoints,z_gridpoints)   # Ec =  E_t[c_{ŧ+1}^{-γ}]
# for (a_t1_ind, a_t1) in enumerate(A_vals), (z_ind, z) in enumerate(mcs_z[t].state_values), ε in enumerate(mc_ε.state_values)
# C_t = ( β*(1+r) .* fill(Ec,1,1,ε_gridpoints).^(-1/γ)
# repeat(mcs_z[1].state_values, 1, ε_gridpoints)
# At = 1/(1+r) .* (C[t] - exp.(κ[t] + z + ε) + A_grid_l) # ******how to account for all possible combos of z,ε??

a_min = 0. # ZBC
a_gridpoints = 10
a_max = 10.
A_vals = exp10.( range(log10(a_min+1), stop = log10(a_max+1), length = a_gridpoints) ) .- 1

for t in reverse(1:T_ret-2) #reverse(t)
    At = zeros(a_gridpoints,z_gridpoints,ε_gridpoints)
    for a_t1 in A_vals, (z_ind, z) in enumerate(mcs_z[t].state_values), ε in mc_ε.state_values
        Ec = 0.
        for (ε_t1_ind, ε_t1) in enumerate(mc_ε.state_values), (z_t1_ind, z_t1) in enumerate(mcs_z[t].state_values)
            A_t2 = extrapolate(interpolate(A[t+2], BSpline(Linear())), 0. )(a_t1, z_t1, ε_t1) # NEED INDICES, NOT VALUES
            Ec += mcs_z[t].p[z_ind, z_t1_ind] * mc_ε.p[1,ε_t1_ind] * ( (1+r).*a_t1 + exp(κ[t+1] + z_t1 + ε_t1) - A_t2 )^(-γ)
        end
        At = 1/(1+r) .* ( ( β*(1+r)*Ec )^(-1/γ) - exp(κ[t] + z + ε) + a_t1 )
    end
    A[t+1] = grid_shift(invert_rule(At, A_vals), A_vals, A_vals, Z_vals_old, Z_vals_new) # ********
end

using ProgressMeter
function testfn()
    for t in 1 #reverse(t)
        i=0
        At = zeros(a_gridpoints,z_gridpoints,ε_gridpoints)
        li = extrapolate(interpolate(A[t+2], BSpline(Linear())), 0. )
        for (ε_ind, ε) in enumerate(mc_ε.state_values), (z_ind, z) in enumerate(mcs_z[t].state_values), (a_t1_ind, a_t1) in enumerate(A_vals)
            Ec = 0.
            for (ε_t1_ind, ε_t1) in enumerate(mc_ε.state_values), (z_t1_ind, z_t1) in enumerate(mcs_z[t].state_values)
                A_t2 = li(a_t1_ind, z_t1_ind, ε_t1_ind) # CHECK THESE INDICES ARE RIGHT
                Ec += mcs_z[t].p[z_ind, z_t1_ind] * mc_ε.p[1,ε_t1_ind] * ( (1+r).*a_t1 + exp(κ[t+1] + z_t1 + ε_t1) - A_t2 )^(-γ)
            end
            At[a_t1_ind,z_ind,ε_ind] = 1/(1+r) .* ( ( β*(1+r)*Ec )^(-1/γ) - exp(κ[t] + z + ε) + a_t1 )
            @show i+= 1
        end
        return At

        # At = Float32[]  #zeros(a_gridpoints,z_gridpoints,ε_gridpoints)
        # i = 0
        # li = extrapolate(interpolate(A[t+2], BSpline(Linear())), 0. )
        # for ε in mc_ε.state_values, (z_ind, z) in enumerate(mcs_z[t].state_values), (a_t1_ind, a_t1) in enumerate(A_vals)
        #     @show i += 1
        #     Ec = 0.
        #     for (ε_t1_ind, ε_t1) in enumerate(mc_ε.state_values), (z_t1_ind, z_t1) in enumerate(mcs_z[t].state_values)
        #         A_t2 = li(a_t1_ind, z_t1_ind, ε_t1_ind) # CHECK THESE INDICES ARE RIGHT
        #         Ec += mcs_z[t].p[z_ind, z_t1_ind] * mc_ε.p[1,ε_t1_ind] * ( (1+r).*a_t1 + exp(κ[t+1] + z_t1 + ε_t1) - A_t2 )^(-γ)
        #     end
        #     push!(At, 1/(1+r) .* ( ( β*(1+r)*Ec )^(-1/γ) - exp(κ[t] + z + ε) + a_t1 ) )
        # end
        # return At
        # # A[t+1] = grid_shift(invert_rule(At, A_vals), A_vals, A_vals, Z_vals_old, Z_vals_new) # ********
    end
end

At = testfn()
A = reshape(At,a_gridpoints,z_gridpoints,ε_gridpoints)
heatmap(A)

function myfunc()
    @showprogress for i in 1:1E9, b in 1:15
        i += i
    end
end

myfunc()

function invert_rule(A,A_vals)
    A_interp = similar(A)
    for col in eachindex(A[1,:])
        # li = LinearInterpolation(A[:,col], A_vals, extrapolation_bc = (0.,Line()))
        # use extrapolate and do different for each side fixed value/line
        A_interp[:,col] = extrapolate(interpolate(A[:,col], BSpline(Linear())), 0.)(A_vals)
    end
    return A_interp
end

A = zeros(11,2)
A_vals = 0:0.1:1
A[:,1] = [(x-5)^3 + 100 for x in A_vals]
A[:,1] = [x^2 for x in A_vals]

plot!(1:25,1:25)
plot()
scatter!(A_vals,A[:,1])
A_inv = invert_rule(A, A[:,1])
scatter!(A[:,1],A_inv[:,1])

plot(x -> (x-5)^3 + 100, 0, 10, ylims=(0,200))

function interpolate_cols(X, Y, X_grid::Array; extrapolation_bc = Throw())
    #linear interpolation of a column of Y on that column of X, evaluated at point x
    Y_interp = similar(X_grid)
    for col in eachindex(X_grid[1,:])
        li = LinearInterpolation(X[:,col], Y[:,col],extrapolation_bc=extrapolation_bc)
        Y_interp[:,col] = li.(X_grid[:,col])
    end
    return Y_interp
end

A_prime_grid = repeat(A_vals, 1, 2)
A_prime = interpolate_cols(A,A_prime_grid,A_prime_grid, extrapolation_bc = Flat())

A_vals = 0:0.1:1
A = [log(x) for x in A_vals]
scatter(A_vals,A)
li = LinearInterpolation(A,A_vals, extrapolation_bc=Flat())
A_inv = li.(A_vals)
scatter!(A_vals,A_inv)

# Check invert_rule
using Interpolations, Plots
A_vals = 0.1:0.01:1.5
Y_vals = 0.1:0.1:0.3
A = [y + x^3 for x in A_vals, y in Y_vals]
minimum(A[:,1])
plot([0; 1],[0; 1], linestyle=:dash, legend=:bottomright)
plot!(A_vals,A[:,2])
# plot!(A_vals,A[:,2])
# plot!(A_vals,A[:,3])
A_inv = invert_rule(A, A_vals, extrapolation_bc = 0.)
plot!(A_vals,A_inv[:,2])
# plot!(A_vals,A_inv[:,2])
# plot!(A_vals,A_inv[:,3])
# ylims!(0,1)

function invert_rule(A, A_vals; extrapolation_bc = Throw())
    #linear interpolation of a column of Y on that column of X, evaluated at point x
    A_vals_interp = similar(A)
    for col in eachindex(A[1,:])
        li = LinearInterpolation(A[:,col], A_vals,extrapolation_bc=extrapolation_bc)
        A_vals_interp[:,col] = li(A_vals)
    end
    return A_vals_interp
end

A
A3d = zeros(size(A)[1],size(A)[2],3)
function myfunc(A)
    for i in 1:size(A3d)[3]
        A3d[:,:,i] = A
    end
    return A3d
end
A3d = myfunc(A)
A3d[:,:,1]==A3d[:,:,3]

function invert_rule_3d(A3d, A_vals; extrapolation_bc = Throw())
    wrapper(A) = invert_rule(A,A_vals,extrapolation_bc = extrapolation_bc)
    A3d_inv = mapslices(wrapper, A3d, dims=1)
end

A3d_inv = invert_rule_3d(A3d, A_vals, extrapolation_bc = Flat())
plot([0; 1],[0; 1], linestyle=:dash, legend=:bottomright)
plot!(A_vals,A3d[:,2,3])
plot!(A_vals,A3d_inv[:,2,3])

# Try to reconstruct construction
A_inv2 = extrapolate(scale(interpolate(A, BSpline(Linear())), A), Flat())(A_vals)
scatter!(A_vals,A_inv2)
A_inv2 = extrapolate()
scale(interpolate(A_vals, BSpline(Linear())), A)
#fail

# Multidimensional LinearInterpolation test (NOT INVERTED)
A_vals = 0:0.1:1
Y_vals = 1:3
A = [y + x^2 for x in A_vals, y in Y_vals]
li = LinearInterpolation((A_vals,Y_vals),A , extrapolation_bc=Flat())
li(0,2)
plotly()
plot(A,st=:surface)
plot!([2 2.5],[1 1],[li(0,2) li(0,2.5)],st=:scatter3d)

function mytestfunc(A_vals, Y_vals)
    for i = 1
        li1 = LinearInterpolation(A_vals,Y_vals)
    end
    for i = 2:3
        @show s = li1(2.)
        li2 = LinearInterpolation(A_vals.+i,Y_vals)
    end
end

mytestfunc([1;2; 3],[1; 4; 9])
# interpolation object doesn't get passed between for loops,
# have to define at the beginning of each

        #             (A_vals, mcs_z[t+1].state_values, mc_ε.state_values)
        # grid_shift(A[t+1], Z_vals_old, Z_vals_new, dims = 2)
        # li = LinearInterpolation((A_vals,mcs_z[t].state_values,),A[t+2] , extrapolation_bc=Flat())

        # function grid_shift(A, A_vals_old, A_vals_new, Y_vals_old, Y_vals_new)
        #     A_grid = zeros( eltype(A[1]), length(A_vals_new), length(Y_vals_new) )
        #
        #     # Find relative coordinates of new basis
        #     li = LinearInterpolation(Y_vals_old, 1:length(Y_vals_old), extrapolation_bc=Line())
        #     Y_ind = li(Y_vals_new)
        #     li = LinearInterpolation(A_vals_old, 1:length(A_vals_old), extrapolation_bc=Line())
        #     A_ind = li(A_vals_new)
        #
        #     # Bi-linear interpolation on new coordinates
        #     li = extrapolate(interpolate(A, BSpline(Linear()) ), Line())  #extrapolate with boundary rules
        #     [ A_grid[i,j] = li(A_ind, Y_ind) for  (i, A_ind) in enumerate(A_ind),
        #                                                 (j, Y_ind) in enumerate(Y_ind)]
        #     return A_grid
        # end

## Income process troubleshoot/calibration

# Net labour income
function Y(N, T_ret, σ_ε, σ_η, σ_z0, κ)
# Takes standard deviations as arguments (not varicances)!
    Y = zeros(Float64, N, T_ret-1)
    z = similar(Y)
    ε = similar(Y)
    z[:,1] = σ_z0 .* randn(N) + σ_η .* randn(N)

    for t in 1:(T_ret-1)
        ε[:,t] = σ_ε .* randn(N)
        Y[:,t] = exp.( κ[t] .+ z[:,t] + ε[:,t] )
        t < (T_ret-1) ? z[:,t+1] = z[:,t] + σ_η .* randn(N) : nothing
    end
    return (Y, z, ε)
end

(Y_l, z_l, ε_l) = Y(10000000, 2, σ_ε, σ_η, σ_z0, κ)
(Y_l, z_l, ε_l) = Y(N, T_ret, σ_ε, 0., 0., κ)
Y_tot = hcat(Y_l, repeat(Y_SS, 1, (T-1)-(T_ret-1)) )
mean(Y_l, dims = 1)
exp(κ[1])
mean(z_l, dims = 1)
exp(κ[1] + mean(z_l[:,1]) + mean(ε_l[:,1]))

li = LinearInterpolation(zeros(10), 1:10)


text:
module Recalls

export @note, @recall, recall

using .Meta: isexpr
using Dates: Second, now
using ExprTools: combinedef, splitdef
using Requires: @require
using UUIDs: UUID, uuid4

include("core.jl")
include("metadata.jl")
include("calls.jl")
include("notes.jl")

# `no_juliainterpreter.jl` includes `juliainterpreter.jl` with a fake
# `Juliainterpreter` module for defining `_call` without a breakpoint.
# `_call` will be re-defined at run-time after `JuliaInterpreter` is
# loaded
include("no_juliainterpreter.jl")

function __init__()
    init!(NOTES)
    @require JuliaInterpreter = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" include("juliainterpreter.jl")
end

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end Recalls

end