text:
# This file contains functions to get both teams that participated in a game

"""
    PublicWyscoutTeams
"""
Base.@kwdef mutable struct PublicWyscoutTeams
    data::DataFrame
end

"""
function teams()

Return a dataframe with both teams that participated in a game.
"""
function teams(event_data::PublicWyscoutLoader, game_id::Int)::PublicWyscoutTeams

    teams_data = PublicWyscoutTeams(data = get_df(joinpath(event_data.path, "teams.json")))
    rename!(teams_data.data, :name => :team_name_short,
                        :officialName => :team_name,
                        :wyId => :team_id)

    teams_data.data = teams_data.data[:,[:team_id, :team_name, :team_name_short]]
    teams_id = collect(keys(lineup(event_data, game_id)))
    teams_id = DataFrame(:team_id => parse.(Int, teams_id))

    teams_data.data = leftjoin(teams_id, teams_data.data, on = :team_id)

    return teams_data
end


text:
using Random
using InvertedIndices
using TensorToolbox

function rapid_tensor_from_theta(theta)
    tensor_shape = size(theta)
    d = ndims(theta)

    rs = []
    for k = 1:d
        theta_k = theta[ ntuple( j->j == k ? (:) : 1, d)... ]
        r = ones( tensor_shape[k] )
        r[1] = exp( theta_k[1] )
        for i = 2:tensor_shape[k]
            r[i] = r[i-1] * exp( theta_k[i] )
        end
        push!(rs, r)
    end

    tensor = ttt( vec(rs[1]), vec(rs[2]) )
    for k = 3:d
        tensor = ttt(tensor, vec(rs[k]) )
    end

    return tensor
end

function rapid_eta_from_tensor(T)
    tensor_shape = size(T)
    eta = zeros( tensor_shape )
    d = ndims(T)
    one_to_N = [1:d;]
    for i = 1:d
        tensor_shape_i = ones(Int, d)
        tensor_shape_i[i] = tensor_shape[i]
        tensor_shape_i_tp = tuple( tensor_shape_i... )
        sums = vec(sum(T, dims=one_to_N[Not(i)]))
        for idxes in reverse(Tuple(CartesianIndices( tensor_shape_i_tp )))
            eta[idxes] = sums[ idxes[i] ]
            if idxes[i] != tensor_shape[i]
                tmp = zeros(Int, d)
                tmp[i] = 1
                eta[idxes] += eta[ idxes + CartesianIndex( tuple(tmp...) ) ]
            end
        end
    end

    return eta
end

function rapid_theta_from_tensor(T)
    tensor_shape = size(T)
    theta = zeros( tensor_shape )
    d = ndims(T)
    for i = 1:d
        tensor_shape_i = ones(Int, d)
        tensor_shape_i[i] = tensor_shape[i]
        tensor_shape_i_tp = tuple( tensor_shape_i... )
        for idxes in CartesianIndices( tensor_shape_i_tp )
            theta[ idxes ] = log( T[ idxes ] )
            if idxes[i] > 1
                tmp = zeros(Int, d)
                tmp[i] = 1
                theta[ idxes ] -= log( T[ idxes - CartesianIndex( tuple(tmp...))] )
            end
        end
    end
    return theta
end

function get_uniform_eta( tensor_shape )
    eta = ones(tensor_shape)
    for idxes in CartesianIndices( tensor_shape )
        tmp = 1.0
        for i = 1:length(idxes)
            tmp *= (tensor_shape[i]-idxes[i]+1)
        end
        eta[ idxes ] = tmp
    end
    eta ./= prod(tensor_shape)
    return eta
end

function get_uniform_theta( tensor_shape )
    theta = zeros(tensor_shape)
    theta[1] = -1.0 * log(prod( tensor_shape ))
    return theta
end

ext: jl

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
function UHFEnergy(H, Dα, Dβ, Fα, Fβ, Vnuc)
    # Calculate energy
    @tensoropt A[:] := 0.5 * (H[i,j] * (Dα[i,j] + Dβ[i,j]) + Fα[i,j] * Dα[i,j] + Fβ[i,j] * Dβ[i,j])
    Ee = A[1]
    return(Ee+Vnuc)
end

function build_fock!(Fα, Fβ, Jα, Jβ, Kα, Kβ, Dα, Dβ, ints::IntegralHelper{Float64,Chonky,AtomicOrbitals})
    # Calculate Fock matrix
    H = ints["T"] + ints["V"]
    Fα .= H
    Fβ .= H
    calcJ!(Jα, Dα, ints)
    calcJ!(Jβ, Dβ, ints)
    calcK!(Kα, Dα, ints)
    calcK!(Kβ, Dβ, ints)
    Fα .+= Jα - Kα + Jβ
    Fβ .+= Jβ - Kβ + Jα
end

function calcJ!(J, D, ints::IntegralHelper{Float64,Chonky,AtomicOrbitals})
    # Calculate Coloumb integrals contracted with D
    ERI = ints["ERI"]
    @tensoropt J[i,j] = ERI[i,j,k,l] * D[l,k]
end

function calcK!(K, D, ints::IntegralHelper{Float64,Chonky,AtomicOrbitals})
    # Calculate Exchange integrals contracted with D
    ERI = ints["ERI"]
    @tensoropt K[i,j] = ERI[i,l,k,j] * D[l,k]
end

function build_fock!(Fα, Fβ, Jα, Jβ, Kα, Kβ, Dα, Dβ, ints::IntegralHelper{Float64,<:AbstractDFERI,AtomicOrbitals})
    H = ints["T"] + ints["V"]
    Fα .= H
    Fβ .= H
    calcJ!(Jα, Dα, ints)
    calcJ!(Jβ, Dβ, ints)
    calcK!(Kα, Dα, ints)
    calcK!(Kβ, Dβ, ints)
    Fα .+= Jα - Kα + Jβ
    Fβ .+= Jβ - Kβ + Jα
end

function calcJ!(J, D, ints::IntegralHelper{Float64,<:AbstractDFERI,AtomicOrbitals})
    # Calculate Coloumb integrals contracted with D
    b = ints["ERI"]
    @tensoropt J[i,j] = b[Q,i,j] * b[Q,k,l] * D[k,l]
end

function calcK!(K, D, ints::IntegralHelper{Float64,<:AbstractDFERI,AtomicOrbitals})
    # Calculate Exchange integrals contracted with D
    b = ints["ERI"]
    @tensoropt K[i,j] = b[Q,i,k] * b[Q,j,l] * D[k,l]
end

function buildD!(D, C, N)
    # Build density matrix
    Co = C[:,1:N]
    @tensoropt D[μ, ν] = Co[μ, i] * Co[ν, i]
end

function build_fock!(Fα, Fβ, Jα, Jβ, Kα, Kβ, Dα, Dβ, ints::IntegralHelper{Float64,<:SparseERI,AtomicOrbitals})
    H = ints["T"] + ints["V"]
    Fα .= H
    Fβ .= H
    calcJK!(Jα, Kα, Dα, ints)
    calcJK!(Jβ, Kβ, Dβ, ints)
    Fα .+= Jα - Kα + Jβ
    Fβ .+= Jβ - Kβ + Jα
end

function calcJK!(J, K, D, ints::IntegralHelper{Float64,<:SparseERI,AtomicOrbitals})
    D = D.data
    J .= 0.0
    K .= 0.0
    eri_vals = ints["ERI"].data
    idxs = ints["ERI"].indexes

    nbas = ints.orbitals.basisset.nbas
    Jarrays = [zeros(Int((nbas^2 + nbas)/2)) for i = 1:Threads.nthreads()]
    Karrays = [zeros(Int((nbas^2 + nbas)/2)) for i = 1:Threads.nthreads()]

    Threads.@threads for z = eachindex(eri_vals)
        i,j,k,l = idxs[z] .+ 1
        ν = eri_vals[z]
        Jt = Jarrays[Threads.threadid()]
        Kt = Karrays[Threads.threadid()]
        @inbounds begin
            @fastmath begin
            ij = Fermi.index2(i-1,j-1)
            ik = Fermi.index2(i-1,k-1) + 1
            il = Fermi.index2(i-1,l-1) + 1

            jk = Fermi.index2(j-1,k-1) + 1
            jl = Fermi.index2(j-1,l-1) + 1

            kl = Fermi.index2(k-1,l-1)

            γij = i !== j
            Xik = i === k ? 2 : 1
            Xjk = j === k ? 2 : 1

            γkl = k !== l
            γab = ij !== kl

            Xil = i === l ? 2 : 1
            Xjl = j === l ? 2 : 1
            if γij && γkl && γab
                # J
                Jt[ij+1] += 2*D[k,l]*ν
                Jt[kl+1] += 2*D[i,j]*ν

                # K
                Kt[ik] += Xik*D[j,l]*ν
                Kt[jk] += Xjk*D[i,l]*ν
                Kt[il] += Xil*D[j,k]*ν
                Kt[jl] += Xjl*D[i,k]*ν

            elseif γkl && γab
                # J
                Jt[ij+1] += 2*D[k,l]*ν
                Jt[kl+1] += 1*D[i,j]*ν

                # K
                Kt[ik] += Xik*D[j,l]*ν
                Kt[il] += Xil*D[j,k]*ν
            elseif γij && γab
                # J
                Jt[ij+1] += 1*D[k,l]*ν
                Jt[kl+1] += 2*D[i,j]*ν

                # K
                Kt[ik] += Xik*D[j,l]*ν
                Kt[jk] += Xjk*D[i,l]*ν

            elseif γij && γkl

                # Only possible if i = k and j = l
                # and i < j ⇒ i < l

                # J
                Jt[ij+1] += 2*D[k,l]*ν

                # K
                Kt[ik] += D[j,l]*ν
                Kt[il] += D[j,k]*ν
                Kt[jl] += D[i,k]*ν
            elseif γab
                # J
                Jt[ij+1] += 1*D[k,l]*ν
                Jt[kl+1] += 1*D[i,j]*ν
                # K
                Kt[ik] += Xik*D[j,l]*ν
            else
                Jt[ij+1] += 1*D[k,l]*ν
                Kt[ik] += D[j,l]*ν
            end
            end
        end
    end

    # Reduce values produces by each thread
    for i::Int16 = 1:nbas
        for j::Int16 = i:nbas
            @inbounds begin
                ij = Fermi.index2(i-1,j-1)
                for k = eachindex(Jarrays)
                    J[i,j] += Jarrays[k][ij+1]
                end

                for k = eachindex(Karrays)
                    K[i,j] += Karrays[k][ij+1]
                end

                J[j,i] = J[i,j]
                K[j,i] = K[i,j]
            end
        end
    end
end

function odadamping(D, Ds, F, Fs)
    dD = D - Ds
    s = tr(Fs * dD)
    c = tr((F - Fs) * (dD))
    if c <= -s/(2*c)
        λ = 1.0
    else
        λ = -s/(2*c)
    end
    Fs .= (1-λ)*Fs + λ*F
    Ds .= (1-λ)*Ds + λ*D
    damp = 1-λ
    F .= Fs
    return damp
end


text:
# This file was generated, do not modify it. # hide
RFR = @load RandomForestRegressor pkg=ScikitLearn

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
struct TerminateCommand <: Exception end

include("git.jl")
include("auth.jl")
include("registry.jl")
include("package.jl")

# Copied from IsURL.jl
# Source: https://github.com/sindresorhus/is-absolute-url (MIT license)
const windowsregex = r"^[a-zA-Z]:[\\]"
const urlregex = r"^[a-zA-Z][a-zA-Z\d+\-.]*:"

"""
    isurl(str)
Checks if the given string is an absolute URL.
# Examples
```julia-repl
julia> isurl("https://julialang.org")
true
julia> isurl("mailto:someone@example.com")
true
julia> isurl("/foo/bar")
false
```
"""
function isurl(str::AbstractString)
    return !occursin(windowsregex, str) && occursin(urlregex, str)
end

function isgithub(url::AbstractString)
    if startswith(url, "http") && startswith(split(url, "//")[2], "github.com")
        return true
    else
        return false
    end
end

function write_project_toml(filepath::String, d::Dict{String, Any})
    open(filepath, "w+") do f # following whatever Pkg does
        TOML.print(f, d; sorted=true, by=key -> (Pkg.Types.project_key_order(key), key))
    end
end


text:
function rrule(::typeof(sparse), I::AbstractVector, J::AbstractVector, V::AbstractVector, m, n, combine::typeof(+))
    project_V = ProjectTo(V)

    function sparse_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        ΔV = project_V(ΔΩ[I .+ m .* (J .- 1)])
        return NoTangent(), NoTangent(), NoTangent(), ΔV, NoTangent(), NoTangent(), NoTangent()
    end

    return sparse(I, J, V, m, n, combine), sparse_pullback
end

function rrule(::Type{T}, A::AbstractMatrix) where T <: AbstractSparseMatrix
    function sparse_pullback(Ω̄)
        return NoTangent(), Ω̄
    end
    return T(A), sparse_pullback
end

function rrule(::Type{T}, v::AbstractVector) where T <: AbstractSparseVector
    function sparse_pullback(Ω̄)
        return NoTangent(), Ω̄
    end
    return T(v), sparse_pullback
end

function rrule(::typeof(findnz), A::AbstractSparseMatrix)
    I, J, V = findnz(A)
    m, n = size(A)

    function findnz_pullback(Δ)
        _, _, V̄ = unthunk(Δ)
        V̄ isa AbstractZero && return (NoTangent(), V̄)
        return NoTangent(), sparse(I, J, V̄, m, n)
    end

    return (I, J, V), findnz_pullback
end

function rrule(::typeof(findnz), v::AbstractSparseVector)
    I, V = findnz(v)
    n = length(v)

    function findnz_pullback(Δ)
        _, V̄ = unthunk(Δ)
        V̄ isa AbstractZero && return (NoTangent(), V̄)
        return NoTangent(), sparsevec(I, V̄, n)
    end

    return (I, V), findnz_pullback
end


text:
using Toxcore, Reactive, GLWindow
println("ima load shit")
window = createwindow("toxclient", 10, 10)

const BOOTSTRAP_ADDRESS = "78.225.128.39"
const BOOTSTRAP_PORT    = 33445
const BOOTSTRAP_KEY     = "7A2306BFBA665E5480AE59B31E116BE9C04DCEFE04D9FE25082316FA34B4DA0C"

typealias FriendID Uint32

immutable FriendRequest
    message::UTF8String
    public_key::ToxPublicKey
end

immutable ToxData
    friends::Vector{FriendID}
    friendrequest::Input{Vector{FriendRequest}}
    messages::Dict{FriendID, Input{UTF8String}}
end
ToxData() = ToxData(Uint8[], Input(FriendRequest[]), Dict{FriendID, Signal{UTF8String}}())

const CLIENT_DICT = Dict{Ptr{Tox}, ToxData}()

function init()
    println("ima init shit")
    my_tox = tox_new()

    # register the callbacks
    tox_callback_friend_request(my_tox, OnFriendRequest, C_NULL)
    tox_callback_friend_message(my_tox, OnFriendMessage, C_NULL)

    # print own address
    address = convert(ASCIIString, tox_self_get_address(my_tox))
    println(address)
    clipboard(address)
    println("key copied to clipboard")

    # define user details
    tox_self_set_name(my_tox, "Julia")
    tox_self_set_status_message(my_tox, "I am Julia, a high-level, high-performance dynamic programming language for technical computing.")

    # bootstrap from the node defined above
    if !tox_bootstrap(my_tox, BOOTSTRAP_ADDRESS, BOOTSTRAP_PORT, ToxPublicKey(BOOTSTRAP_KEY))
        println("Failed to bootstrap.")
        exit()
    end
    println("init done!")

    my_tox
end

function OnFriendRequest(tox::Ptr{Tox}, public_key::Ptr{Uint8}, message::Ptr{Uint8}, len::Csize_t, user_data::Ptr{Void})
    push!(CLIENT_DICT[tox].friendrequest, [FriendRequest(utf8(message, len), ToxPublicKey(public_key))])
    nothing
end
function OnFriendMessage(tox::Ptr{Tox}, friend_number::FriendID, typ::TOX_MESSAGE_TYPE, message::Ptr{Uint8}, len::Csize_t, user_data::Ptr{Void})
    push!(CLIENT_DICT[tox].messages[friend_number], utf8(message, len))
    nothing
end

function main()
    #Create a default Tox
    my_tox = init()
    CLIENT_DICT[my_tox] = ToxData()

    iterate_signal = fpswhen(window.inputs[:open], 60)

    lift(CLIENT_DICT[my_tox].friendrequest) do request_list
        isempty(request_list) && return nothing
        current = pop!(request_list) # get current request
        friend = tox_friend_add_norequest(my_tox, current.public_key) #accept as friend
        push!(CLIENT_DICT[my_tox].friends, friend) # add to friends list

        friendmessages = Input(utf8(""))
        CLIENT_DICT[my_tox].messages[friend] = friendmessages# add to friends list
        println(current.message)
        #register message printout
        lift(println, friendmessages)
        nothing
    end

    while window.inputs[:open].value
        yield()
        tox_iterate(my_tox)
        GLFW.PollEvents()
        sleep(tox_iteration_interval(my_tox)/1000)
    end
    tox_kill(my_tox)
    GLFW.Terminate()
end

main()


text:
using Test, CSV, DelimitedFiles, ModelSelection.Preprocessing

# TODO: Add Seasonal Adjustament to the tests
# TODO: Add panel_data to the tests
# TODO: Add time_data to the tests

function replace_string_to_nothing(data)
    for n=1:size(data[1],1)
        for m=1:size(data[1], 2)
            if isa(data[1][n,m], String) || isa(data[1][n,m], SubString)
                data[1][n,m] = missing
            end
        end
    end
    return data
end

function convert_array_with_missing(data, datatype)
    return convert(Array{Union{Missing, datatype}}, data)
end

data_paneltime_without_missings_dataframe = CSV.read(DATABASE_PANELTIME_WITHOUT_MISSINGS_FILENAME)
data_paneltime_without_missings_tuple = replace_string_to_nothing(readdlm(DATABASE_PANELTIME_WITHOUT_MISSINGS_FILENAME, ',', header=true))
data_paneltime_without_missings_tuple = (convert_array_with_missing(data_paneltime_without_missings_tuple[1], Float64), data_paneltime_without_missings_tuple[2])

data_time_without_missings_tuple = replace_string_to_nothing(readdlm(DATABASE_TIME_WITHOUT_MISSINGS_FILENAME, ',', header=true))
data_time_without_missings_tuple = (convert_array_with_missing(data_time_without_missings_tuple[1], Float64), data_time_without_missings_tuple[2])

data_paneltime_with_missings_dataframe = CSV.read(DATABASE_PANELTIME_WITH_MISSINGS_FILENAME)

data_paneltime_with_missings_tuple = replace_string_to_nothing(readdlm(DATABASE_PANELTIME_WITH_MISSINGS_FILENAME, ',', header=true))
data_paneltime_with_missings_tuple = (convert_array_with_missing(data_paneltime_with_missings_tuple[1], Float64), data_paneltime_with_missings_tuple[2])

data_paneltime_as_missings = replace_string_to_nothing(readdlm(DATABASE_PANELTIME_AS_MISSINGS_FILENAME, ',', header=true))
data_paneltime_as_missings = (convert_array_with_missing(data_paneltime_as_missings[1], Float64), data_paneltime_as_missings[2])

test_time_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_TIME_FILENAME, ','), nothing) )[1], Float32)
test_panel_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_PANEL_FILENAME, ','), nothing) )[1], Float32)
test_paneltime_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_PANELTIME_FILENAME, ','), nothing) )[1], Float32)
test_missing_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_MISSING_FILENAME, ','), nothing) )[1], Float32)
test_removedmissing_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_REMOVEDMISSING_FILENAME, ','), nothing) )[1], Float32)
test_nonmissing_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_NONMISSING_FILENAME, ','), nothing) )[1], Float32)
test_removedoutliers_data = convert_array_with_missing(replace_string_to_nothing( (readdlm(TEST_REMOVEDOUTLIERS_FILENAME, ','), nothing) )[1], Float32)

@testset "Equations" begin
    @testset "Stata like" begin
        data = Preprocessing.input("y x2 x1", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
    end

    @testset "R like" begin
        data = Preprocessing.input("y ~ x2 + x1", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
        data = Preprocessing.input("y ~ x.", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :x10]
    end

    @testset "A String separated by comma with spaces" begin
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
        data = Preprocessing.input("y, x*", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :x10]
    end

    @testset "A String separated by comma without spaces" begin
        data = Preprocessing.input("y,x2,x1", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
        data = Preprocessing.input("y,x*", data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :x10]
    end

    @testset "Array of String" begin
        data = Preprocessing.input(["y", "x2", "x1"], data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
        data = Preprocessing.input(["y" "x2" "x1"], data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
        data = Preprocessing.input(["y" "x*"], data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9, :x10]
    end

    @testset "Array of Symbol" begin
        data = Preprocessing.input([:y, :x2, :x1], data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
        data = Preprocessing.input([:y :x2 :x1], data_paneltime_without_missings_tuple)
        @test data.equation == [:y, :x2, :x1]
    end

    @testset "Variables errors" begin
        data = nothing
        try
            data = Preprocessing.input("a b c", data_paneltime_without_missings_tuple[1])
        catch e
        end
        @test data == nothing

        data = nothing
        try
            data = Preprocessing.input("",  data_paneltime_without_missings_tuple[1])
        catch e
        end
        @test data == nothing

        data = nothing
        try
            data = Preprocessing.input([:a, :b, :c],  data_paneltime_without_missings_tuple)
        catch e
        end
        @test data == nothing
    end
end

@testset "Options" begin
    @testset "Intercept" begin
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple)
        @test data.intercept == true
        @test findfirst(isequal(:_cons), data.expvars) != nothing
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, intercept=true)
        @test data.intercept == true
        @test findfirst(isequal(:_cons), data.expvars) != nothing
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, intercept=false)
        @test data.intercept == false
        @test findfirst(isequal(:_cons), data.expvars) == nothing
    end

    @testset "Method" begin
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple)
        @test data.datatype == Float32
        @test isa(data.depvar_data, Array{Float32}) || isa(data.depvar_data, Array{Union{Missing, Float32}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method="fast")
        @test data.datatype == Float32
        @test isa(data.depvar_data, Array{Float32}) || isa(data.depvar_data, Array{Union{Missing, Float32}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method="FAST")
        @test data.datatype == Float32
        @test isa(data.depvar_data, Array{Float32}) || isa(data.depvar_data, Array{Union{Missing, Float32}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method=:fast)
        @test data.datatype == Float32
        @test isa(data.depvar_data, Array{Float32}) || isa(data.depvar_data, Array{Union{Missing, Float32}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method=:FAST)
        @test data.datatype == Float32
        @test isa(data.depvar_data, Array{Float32}) || isa(data.depvar_data, Array{Union{Missing, Float32}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method="precise")
        @test data.datatype == Float64
        @test isa(data.depvar_data, Array{Float64}) || isa(data.depvar_data, Array{Union{Missing, Float64}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method="PRECISE")
        @test data.datatype == Float64
        @test isa(data.depvar_data, Array{Float64}) || isa(data.depvar_data, Array{Union{Missing, Float64}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method=:precise)
        @test data.datatype == Float64
        @test isa(data.depvar_data, Array{Float64}) || isa(data.depvar_data, Array{Union{Missing, Float64}})
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method=:PRECISE)
        @test data.datatype == Float64
        @test isa(data.depvar_data, Array{Float64}) || isa(data.depvar_data, Array{Union{Missing, Float64}})

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, method="invalid")
        catch e
        end
        @test data == nothing
    end

    @testset "Time" begin
        data = Preprocessing.input("y, x2, x1", data_time_without_missings_tuple, time=:time)
        @test data.time == :time
        @test data.depvar_data == test_time_data[:,1]
        @test data.expvars_data == test_time_data[:,2:end]

        data = Preprocessing.input("y, x2, x1", data_time_without_missings_tuple, time="time")
        @test data.time == :time
        @test data.depvar_data == test_time_data[:,1]
        @test data.expvars_data == test_time_data[:,2:end]

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_time_without_missings_tuple, time="invalid")
        catch e
        end
        @test data == nothing

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_paneltime_as_missings, time=:time)
        catch e
        end
        @test data == nothing
    end

    @testset "Panel" begin
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, panel=:panel)
        @test data.panel == :panel
        @test data.depvar_data == test_panel_data[:,1]
        @test data.expvars_data == test_panel_data[:,2:end]

        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, panel="panel")
        @test data.panel == :panel
        @test data.depvar_data == test_panel_data[:,1]
        @test data.expvars_data == test_panel_data[:,2:end]

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, panel="invalid")
        catch e
        end
        @test data == nothing

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_paneltime_as_missings, panel=:panel)
        catch e
        end
        @test data == nothing
    end

    @testset "Panel Time" begin
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, time=:time, panel=:panel)
        @test data.time == :time
        @test data.panel == :panel

        @test data.depvar_data == test_paneltime_data[:,1]
        @test data.expvars_data == test_paneltime_data[:,2:end]

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_paneltime_with_missings_tuple, time="invalid", panel="invalid")
        catch e
        end
        @test data == nothing

        data = nothing
        try
            data = Preprocessing.input("y, x2, x1", data_paneltime_as_missings, time=:time, panel=:panel)
        catch e
        end
        @test data == nothing
    end

    @testset "Remove outliers" begin
        data = Preprocessing.input("y, x2, x1", data_paneltime_without_missings_tuple, removeoutliers=true)
        t = data.depvar_data == test_removedoutliers_data[:,1]
        @test isa(t, Missing) || t == true
        t = data.expvars_data == test_removedoutliers_data[:,2:end]
        @test isa(t, Missing) || t == true

    end
end

@testset "Data" begin
    @testset "With missings" begin
        @testset "DataFrame" begin
            data = Preprocessing.input("y x2 x1", data_paneltime_without_missings_dataframe)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            @test data.expvars_data == test_nonmissing_data[:,2:end]
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data=data_paneltime_without_missings_dataframe)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            @test data.expvars_data == test_nonmissing_data[:,2:end]
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_dataframe)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data=data_paneltime_with_missings_dataframe)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_dataframe, removemissings=false)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_dataframe, removemissings=true)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_removedmissing_data[:,1]
            t = data.expvars_data == test_removedmissing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 18
        end

        @testset "Tuple" begin
            data = Preprocessing.input("y x2 x1", data_paneltime_without_missings_tuple)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            @test data.expvars_data == test_nonmissing_data[:,2:end]
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data=data_paneltime_without_missings_tuple)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            @test data.expvars_data == test_nonmissing_data[:,2:end]
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_tuple)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data=data_paneltime_with_missings_tuple)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_tuple, removemissings=false)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_tuple, removemissings=true)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_removedmissing_data[:,1]
            @test data.expvars_data == test_removedmissing_data[:,2:end]
            @test data.nobs == 18
        end

        @testset "Data" begin
            data = Preprocessing.input("y x2 x1", data_paneltime_without_missings_tuple[1], datanames=names(data_paneltime_without_missings_dataframe))
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            @test data.expvars_data == test_nonmissing_data[:,2:end]
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data=data_paneltime_without_missings_tuple[1], datanames=names(data_paneltime_without_missings_dataframe))
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            @test data.expvars_data == test_nonmissing_data[:,2:end]
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data=data_paneltime_without_missings_tuple[1], datanames=data_paneltime_without_missings_tuple[2])
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_nonmissing_data[:,1]
            t = data.expvars_data == test_nonmissing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_tuple[1], datanames=data_paneltime_without_missings_tuple[2], removemissings=false)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_missing_data[:,1]
            t = data.expvars_data == test_missing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 20

            data = Preprocessing.input("y x2 x1", data_paneltime_with_missings_tuple[1], datanames=data_paneltime_without_missings_tuple[2], removemissings=true)
            @test data.depvar == :y
            @test data.expvars == [:x2, :x1, :_cons]
            @test data.depvar_data == test_removedmissing_data[:,1]
            t = data.expvars_data == test_removedmissing_data[:,2:end]
            @test isa(t, Missing) || t == true
            @test data.nobs == 18

        end
    end
end


text:
smooth(g::AbstractMatrix, k::IdentityKernel) = g
smooth(g::AbstractMatrix, k::SmoothingKernel) = k.smoother(g, k.ξ)

nobs(mf::MomentFunction) = mf.nobs
nmom(mf::MomentFunction) = mf.nmom
npar(mf::MomentFunction) = mf.npar

function get_mom_deriv(g::Function, dtype::Symbol, k::SmoothingKernel, nobs, nmom, npar)

  s(θ::Vector)             = smooth(g(θ), k)
  gn(θ::Vector, p::Vector) = s(θ)'*p
  g1(θ::Vector, λ::Vector) = s(θ)*λ

  g2(θ::Vector, p::Vector, λ::Vector) = (p'*s(θ)*λ)[1]

  gn_closure(θ::Vector) = gn(θ, __p)
  g1_closure(θ::Vector) = g1(θ, __λ)
  g2_closure(θ::Vector) = g2(θ, __p, __λ)

  gn_closure!(θ::Vector, gg) = gg[:] = gn(θ, __p)
  g1_closure!(θ::Vector, gg) = gg[:] = g1(θ, __λ)

  if dtype==:typed
  	∂gn  = ForwardDiff.forwarddiff_jacobian(gn_closure, Float64, fadtype=:typed)
  	∂g1  = ForwardDiff.forwarddiff_jacobian(g1_closure, Float64, fadtype=:typed)
  	∂²g2 = ForwardDiff.forwarddiff_hessian(g2_closure, Float64, fadtype=:typed)
  elseif dtype==:dual
  	∂gn  = ForwardDiff.forwarddiff_jacobian(gn_closure!, Float64, fadtype=:dual, n = npar, m = nmom)
  	∂g1  = ForwardDiff.forwarddiff_jacobian(g1_closure!, Float64, fadtype=:dual, n = npar, m = nobs)
  	∂²g2 = ForwardDiff.forwarddiff_hessian(g2_closure, Float64, fadtype=:typed)
  elseif dtype==:diff
  	∂gn(θ::Vector)  = Calculus.jacobian(gn_closure, θ, :central)
  	∂g1(θ::Vector)  = Calculus.jacobian(g1_closure, θ, :central)
  	∂²g2(θ::Vector) = Calculus.hessian(g2_closure, θ, :central)
  end
  return (g, s, gn, ∂gn, ∂g1, ∂²g2)
end

function MomentFunction(g::Function, dtype::Symbol; nobs = Nothing, npar = Nothing, nmom = Nothing)
  ## Default is no smoothing
  MomentFunction(get_mom_deriv(g, dtype, IdentityKernel(), nobs, nmom, npar)...,
                                IdentityKernel(), nobs, nmom, npar)
end

function MomentFunction(g::Function, dtype::Symbol, k::SmoothingKernel;
                         nobs = Nothing, npar = Nothing, nmom = Nothing)
  ## Default is no smoothing
  MomentFunction(get_mom_deriv(g, dtype, k, nobs, nmom, npar)..., k, nobs, nmom, npar)
end


text:
#= ercot_stoch_reliability.jl
solves the a stochastic two stage unit-commitment model in which slow generators
are committed in the first stage and fast generators can be committed in real time.
Uses formulation based on Papavasiliou and Oren (2013).
Demand response can function as a slow generator with an advance commmitment
and an advane production schedule if desired.
Uncertainty can come from total demand and from the reliabiltiy of DR. To do
the latter, DR must receive a first-stage production schedule (i.e. be type
3.).

cmd line format should be:
julia ercot_stoch_reliability.jl <date> <inputs_file_name.csv> <multi-runTF> <period_name>

multi-runTF is usually true: it indicates that the current simulation is one time
period of a group of periods that make up a single simulation. This just affects
the file structure.

Written under Julia 0.6.4
Patricia Levi
pjlevi@stanford.edu
NOV 2018
=#

### Packages ###
using JuMP
using Gurobi
using DataFrames
using CSV
include("convert3dto2d.jl")
include("make_scenarios_rejection_sample.jl")
include("writecsvmulti.jl")
include("combine_scenarios.jl")

test = Gurobi.Env() # test that gurobi is working

# -------------------------------------------
# USER CONTROLS & CMD LINE ARGS ##########
# -------------------------------------------

# USER PARAMS #
no_vars = false #stops execution before making variables
debug = false # stops execution before solving model

sherlock_fol = "~/dr_stoch_uc/"
sherlock_input_file = "inputs/"
sherlock_output_file = "outputs/"
sherlock_params_fol = "code/"
laptop_fol = "~/Documents/dr_stoch_uc/"
laptop_input_file = "inputs/"
laptop_output_file = "outputs/"
laptop_params_fol = "code/"

### DR is added to slow/fast categories depending on value of
## DRtype (in inputs csv)
dr_gens = ["DR"]
notdr_gens = ["BIOMASS","COAL","GAS","GAS_CC","GAS_CT","GAS_ICE","GAS_ST","HYDRO",
            "LANDFILL_GAS","NUCLEAR","OIL","SOLAR","WIND"    ]

####### END USER CONTROLS ##########

# --------------------------------------------------------------------------------------
# FILEPATH SETUP
# -------------------------------------------
if split(pwd(),"/")[2] == "Users" #returns true if we're on Patricia's laptop
    Sherlock = false
else
    Sherlock = true
end
@show Sherlock

if Sherlock
    base_fol = sherlock_fol
    input_fol = string(sherlock_fol,sherlock_input_file)
    params_fol = string(sherlock_fol, sherlock_params_fol)
else
    base_fol = laptop_fol
    input_fol = string(laptop_fol,laptop_input_file)
    params_fol = string(laptop_fol, laptop_params_fol)
end
default_data_fol = string(input_fol,"ercot_default/")

# PARSE CMD LINE ARGS #
ARGNAMES = ["date" ,"inputs_file_name","input_verion" ,"multi-runTF", "period_name" ]
defaultARGS = [Dates.format(Dates.now(),"Y-m-d"),"inputs_ercot.csv","rand_o25_00pp_keyDays2","true","periods_1_468_588.csv"]
localARGS = length(ARGS) > 0 ? ARGS : defaultARGS #if ARGS supplied, use those. otherwise, use default
nargs = length(localARGS)
@show localARGS

if nargs == 5
    submitdate = localARGS[1]
    input_file_name = localARGS[2]
    input_version = localARGS[3]
    multiTF = parse(Bool,lowercase(localARGS[4]))
    periodID = localARGS[5]
elseif nargs > 5
    error(string("Too many arguments supplied. Need ",join(ARGNAMES[1,:]," ")))
elseif nargs <5
    warn("not enough arguments supplied. Need ", join(ARGNAMES[1,:]," "))
end

read_inputs = CSV.read(string(params_fol,input_file_name))
read_inputs = read_inputs[:,[:input_name, Symbol(input_version)]]
# add cmd line args to read_inputs file
newdf = DataFrame(input_name = ARGNAMES[1:length(localARGS)], args = localARGS)
names!(newdf.colindex,map(parse,["input_name",input_version]))
read_inputs = vcat(read_inputs, newdf)
@show read_inputs
# ~ transpose read_inputs so that values can be referenced by name
read_inputs[:,:rowkey] = 1
inputs = unstack(read_inputs,:rowkey, :input_name, Symbol(input_version))

# parse out non-string inputs
startlim = parse(Float64,inputs[1,:startlim])
hourlim = parse(Float64,inputs[1,:hourlim])
energylim = parse(Float64,inputs[1,:energylim])
ramplims = parse(Float64,inputs[1,:ramplims])
durationlim = parse(Float64, inputs[1,:durationlim])
dr_override = parse(Bool,lowercase(inputs[1,:dr_override]))
dr_varcost = parse(Float64,inputs[1,:dr_varcost])
randScenarioSel = parse(Bool,lowercase(inputs[1,:randScenarioSel]))
trueBinaryStartup = parse(Bool,lowercase(inputs[1,:trueBinaryStartup]))
DRtype = parse(Int64,inputs[1,:DRtype])
DRrand = parse(Bool,lowercase(inputs[1,:DRrand]))
nd_nrand_o = parse(Int64,inputs[1,:nd_nrand_o])
int_length = parse(Int64,inputs[1,:intlength])
debug = !parse(Bool, lowercase(inputs[1,:solve_model]))
MIPFocusParam = parse(Int64,inputs[1,:MIPFocusParam])
MIPGapParam = parse(Float64,inputs[1,:MIPGapParam])
NodefileStartParam = parse(Float64,inputs[1,:NodefileStartParam])
MethodParam = parse(Int64,inputs[1,:MethodParam])
ThreadsParam = parse(Int64,inputs[1,:ThreadsParam])

if Sherlock
    output_fol = string(sherlock_fol, sherlock_output_file, input_version,"_",submitdate,"/")
else
    output_fol = string(laptop_fol, laptop_output_file, input_version,"_",submitdate,"/")
end
subsel_data_fol = string(input_fol,inputs[1,:timeseriesID],"/")

# setup unique ID for period identification in output
if multiTF
    pbits = split(periodID,r"_|\.|-")
    periodnum = pbits[2]
    periodfirst = pbits[3]
    periodlast = pbits[4]
    periodsave = string("p",periodnum,"_",periodfirst,"_",periodlast)
end

# --------------------------------------------------------------------------------------
### TIME PERIOD DATA
# -------------------------------------------
if multiTF
    hours = CSV.read(string(subsel_data_fol,periodID),datarow=1,types=[Int])[1]
else
    hours = CSV.read(string(subsel_data_fol,"periods_",timeseriesID,".csv"),datarow=1,types=[Int])[1]
end
n_days = convert(Int32,ceil(length(hours)/24))
n_periods = length(hours)

# convert first and notfirst periods to indices
t_firsts = 1
t_notfirst = collect(2:length(hours))

dem2 = CSV.read(string(default_data_fol, inputs[1,:demandFile]),datarow=2,missingstring="NA")
# subselect for just the rows corresponding to 'hours'
dem = dem2[hours,2]

# --------------------------------------------------------------------------------------
# STOCHASTIC VARIABLE DATA
# -------------------------------------------
# TESTING
# vdr = [0.9,1,1.1]
# pro = [0.25,0.5,0.25]
# probs = CSV.read(string(default_data_fol , "dist_input_",stochID,".csv"))
# vdr_in = convert(Array,probs[1,:]) # converts the first row of probs to an Array
# pro_in = rationalize.(convert(Array,probs[2,:])) #to avoid rounding issues later
    # if this becomes a problem, look into https://github.com/JuliaMath/DecFP.jl

## Net Demand uncertainty ##
# if _vdem files do not exist already, assume that the distribution data
# is given by a file describing the probability distribution of deviations
if !isfile(string(subsel_data_fol,"demandScenarios_vdem","_",inputs[1,:stochID],"_",periodsave,".csv"))
    ndprobs = CSV.read(string(default_data_fol , "dist_input_",inputs[1,:stochID],"_nd.csv"))
    ndv_in = convert(Array,ndprobs[1,:]) # converts the first row of probs to an Array
    ndpro_in = rationalize.(convert(Array,ndprobs[2,:])) #to avoid rounding issues later

    demandScenarios = make_scenarios(n_periods, ndv_in, ndpro_in, int_length; randsel = randScenarioSel, nrand = nd_nrand_o)
    vdem_in = demandScenarios[1]
    vdem_p = demandScenarios[2]
    writecsvmulti(DataFrame(vdem_in),subsel_data_fol,string("demandScenarios_vdem","_",inputs[1,:stochID]),multiTF,periodsave)
    writecsvmulti(DataFrame(vdem_p),subsel_data_fol,string("demandScenarios_prob","_",inputs[1,:stochID]),multiTF,periodsave)
else
    vdem_in = convert(Array,CSV.read(string(subsel_data_fol,"demandScenarios_vdem","_",inputs[1,:stochID],"_",periodsave,".csv")))
    vdem_p = convert(Array,CSV.read(string(subsel_data_fol,"demandScenarios_prob","_",inputs[1,:stochID],"_",periodsave,".csv")))
end

## DR Response uncertainty ##
if !isfile(string(subsel_data_fol,"drResponseScenarios_vdr","_",inputs[1,:DRrand_ID],"_",periodsave,".csv"))
    dr_int_length = parse(Int64,inputs[1,:dr_int_length])
    dr_nrand_o = parse(Int64,inputs[1,:dr_nrand_o])

    drprobs = CSV.read(string(default_data_fol , "dist_input_",inputs[1,:DRrand_ID],"_dr.csv"))
    drv_in = convert(Array,drprobs[2,:]) # converts the first row of probs to an Array
    drpro_in = rationalize.(convert(Array,drprobs[1,:])) #to avoid rounding issues later

    demandScenarios = make_scenarios(n_periods, drv_in, drpro_in, dr_int_length; randsel = randScenarioSel, nrand = dr_nrand_o)
    vdr_in = demandScenarios[1]
    vdr_p = demandScenarios[2]
    writecsvmulti(DataFrame(vdr_in),subsel_data_fol,string("drResponseScenarios_vdr","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
    writecsvmulti(DataFrame(vdr_p),subsel_data_fol,string("drResponseScenarios_vdr_p","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
else
    vdr_in = convert(Array,CSV.read(string(subsel_data_fol,"drResponseScenarios_vdr_",inputs[1,:DRrand_ID],"_",periodsave,".csv")))
    vdr_p = convert(Array,CSV.read(string(subsel_data_fol,"drResponseScenarios_vdr_p_",inputs[1,:DRrand_ID],"_",periodsave,".csv")))
end

## Create new vdr, vdem, pro representing combined scenarios from demand
## and DR reliability
println(vdr_p)
println(vdem_p)
scenario_info = combine_scenarios(vdr_in,vdr_p,vdem_in,vdem_p)
vdr = scenario_info[1];
vnd = scenario_info[2];
pro = scenario_info[3];
writecsvmulti(DataFrame(vdr),subsel_data_fol,string("vdr","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
writecsvmulti(DataFrame(vnd),subsel_data_fol,string("vnd","_",inputs[1,:DRrand_ID]),multiTF,periodsave)
writecsvmulti(DataFrame(pro),subsel_data_fol,string("pro","_",inputs[1,:DRrand_ID]),multiTF,periodsave) ## WARNING: passing columns argument with non-AbstractVector entries is deprecated

println("sum of probabilities is ", sum(pro))
n_omega = length(pro) #redefine for new number of scenarios

# Set up demand realizations matrix dem_real[t,o] #
dem_real = dem .* vnd

# --------------------------------------------------------------------------------------
# GENERATOR DATA
# -------------------------------------------

genset = CSV.read(string(default_data_fol,inputs[1,:genFile]), missingstring ="NA")
# genset[Symbol("Plant Name")] # this is how to access by column name if there are spaces
# names(genset) # this is how to get the column names
# anscombe[:,[:X3, :Y1]]  #how to grab several columns by colname

dr_ind = findin(genset[:Fuel],dr_gens)
slow_ind = findin(genset[:Speed],["SLOW"])
fast_ind = findin(genset[:Speed],["FAST"])

# adjust slow/fast defn based on DRtype
if DRtype == 1
    fast_ind = vcat(fast_ind_, dr_ind)
elseif (DRtype == 2) | (DRtype == 3)
    slow_ind = vcat(slow_ind, dr_ind)
else
    error("DRtype must be 1, 2 or 3")
end

#generator min and max
pmin = genset[:PMin]
pmax = genset[:Capacity]
startup = genset[:StartCost]
varcost = genset[:VCost]
rampmax = genset[:ramprate]

# for manual override of DR variable cost
if dr_override
    varcost[dr_ind] = dr_varcost
end

# Generator PARAMS ###
n_gsl= length(slow_ind)# number of slow generators
n_g =nrow(genset)# number of generators
n_gf = length(fast_ind)
n_gdr = length(dr_ind) #number of DR generators
n_t = n_periods # number of timesteps

## GENERATOR AVAILABILITY ##
pf = repeat([1.0], inner = [n_g, n_t])

### Wind and solar
# load wind, solar info
solar_avail = CSV.read(string(default_data_fol,"solar_availability_factors_2016.csv"))
wind_avail = CSV.read(string(default_data_fol,"wind_availability_factors_2016.csv"))
# remove last col of each
solar_avail = solar_avail[:,1:(ncol(solar_avail)-1)]
wind_avail = wind_avail[:,1:(ncol(wind_avail)-1)]
# loop through all colnames, use findin(genset[:plantUnique],XX) to get row
# sub in new info
for i in 1:length(names(solar_avail))
    col = names(solar_avail)[i]
    ind = findin(genset[:plantUnique],[convert(String, col)])
    pf[ind,:] = solar_avail[hours,i]
end

for i in 1:length(names(wind_avail))
    col = names(wind_avail)[i]
    ind = findin(genset[:plantUnique],[convert(String, col)])
    pf[ind,:] = wind_avail[hours,i]
end

### Demand Response
if parse(inputs[1,:availID])!=0
    dr_avail = CSV.read(string(default_data_fol,inputs[1,:availID],".csv"))
    # remove first col
    dr_avail = dr_avail[:,2:ncol(dr_avail)]
    # sub in new info
    for i in 1:length(names(dr_avail))
        col = names(dr_avail)[i]
        ind = findin(genset[:plantUnique],[convert(String, col)])
        if length(ind) > 0
            pf[ind,:] = dr_avail[hours,i]
        else
            error("No corresponding DR unit found for ", names(dr_avail)[i])
        end
    end
end

# --------------------------------------------------------------------------------------
# MODEL SETUP
# -------------------------------------------
### SETS ###
TIME = 1:n_t
SCENARIOS = 1:n_omega
GENERATORS = 1:n_g #all generators
GEN_NODR = findin(genset[:Fuel],notdr_gens)
GF = fast_ind
GSL = slow_ind #slow generators
GDR = dr_ind #DR generators
# need an index for where the DR is in the slow generators
GDR_SL_ind = findin(GSL,GDR)

# -------------------------------------------
### MODEL ###
# m = Model(solver = ClpSolver())
m = Model(solver=GurobiSolver(Threads = ThreadsParam,ImproveStartGap = 0.0002, Method=MethodParam,#1,
                              MIPFocus=MIPFocusParam,#3,
                              MIPGap=MIPGapParam,#0.0003,
                              NodefileStart = NodefileStartParam,#0.05,
                              NodefileDir = "/scratch/users/pjlevi/gurobi_solving_outputs/",
                              Seed = convert(Int64,abs(floor(rand(Float64)*2000000000)))))
                # can also try reducing threadcount, reducing nodefilestart
if no_vars
    error("just testing model so we are stopping here")
else
    println("begin building model")
end

# -----------------------------------------------------------------------------------------------------
### VARIABLES ###
# @variable(m, 0 <= x <= 2 )
@variable(m, z[1:n_gsl,1:n_t]) # slow generator startup
if trueBinaryStartup
    @variable(m, w[1:n_gsl,1:n_t], Bin) # slow generator commitment TRUE BINARY
else
    @variable(m, 0 <= w[1:n_gsl,1:n_t] <= 1) # slow generator commitment RELAXED BINARY
end

#real-time commitment, startup
@variable(m, v[1:n_g,1:n_t,1:n_omega]) # generator startup
if trueBinaryStartup
    @variable(m, u[1:n_g,1:n_t,1:n_omega], Bin) # generator commitment TRUE BINARY
else
    @variable(m, 0 <= u[1:n_g,1:n_t,1:n_omega] <= 1) # generator commitment RELAXED BINARY
end

#production variables
@variable(m, p[1:n_g,1:n_t,1:n_omega] >= 0) #generator production

# if DR requires advance production schedule
if DRtype == 3
    @variable(m, p_dr[1:n_gdr, 1:n_t] >= 0) # DR day-ahead production commitment
end

# track number of startups to determine costs
@variable(m, start_num[1:n_g, 1:n_t, 1:n_omega] >= 0)

# -----------------------------------------------------------------------------------------------------
### CONSTRAINTS ###
# @constraint(m, 1x + 5y <= 3.0 )

#SUPPLY-DEMAND
@constraint(m,supplydemand[t=1:n_t, o=1:n_omega],
    sum( p[g,t,o] for g=1:n_g) >= dem_real[t,o])

#GENMIN
@constraint(m, mingen[g= GEN_NODR, t= 1:n_t, o=1:n_omega ],
    p[g,t,o] >= pmin[g] * u[g,t,o])
#GENMAX - only for GEN_NODR if dr is random
@constraint(m,[g = GEN_NODR, t = 1:n_t, o=1:n_omega],
    p[g,t,o] <= pmax[g] * u[g,t,o] * pf[g,t] )

if DRtype == 3
#GENMAXDR
    @constraint(m,[g = 1:n_gdr, t = 1:n_t],
        p_dr[g,t] <= pmax[GDR[g]] * w[GDR_SL_ind[g],t] * pf[GDR_SL_ind[g],t]) #needed for p_dr
#GENMINDR
    @constraint(m,[g = 1:n_gdr, t = 1:n_t],
        p_dr[g,t] >= pmin[GDR[g]] * w[GDR_SL_ind[g],t])
end

#START_S
@constraint(m,[g = 1:n_gsl, t=t_notfirst],
    z[g,t] == w[g,t] - w[g,t-1])
#START_F
@constraint(m,[g=GENERATORS ,t=t_notfirst, o = SCENARIOS],
    v[g,t,o] == u[g,t,o] - u[g,t-1,o])
#INIT_S
@constraint(m,[g=1:n_gsl,t=t_firsts],
    z[g,t] == w[g,t])
#INIT_F
@constraint(m,[g=GENERATORS,t=t_firsts, o = SCENARIOS],
    v[g,t,o] == u[g,t,o])
#NAN_ST
@constraint(m,[g = 1:n_gsl, t = TIME, o = SCENARIOS],
    v[GSL[g],t,o] == z[g,t])
#NAN_CM
@constraint(m,[g = 1:n_gsl, t = TIME, o = SCENARIOS],
    u[GSL[g],t,o] == w[g,t])

#DR_RAND
if DRtype == 3 && DRrand
    @constraint(m, dr_rand[g=1:n_gdr,t = TIME, o = SCENARIOS],
         p[GDR[g],t,o] == p_dr[g,t] * vdr[t,o])
elseif DRtype != 3 && DRrand
    error("DRrand can only be used with a DRtype of 3")
end
#TODO: make an elseif for DRype ==3 and !DRrand

#STARTUP COUNT
@constraint(m, [g=GENERATORS, t = TIME, o = SCENARIOS],
    start_num[g,t,o] >= v[g,t,o])

#RAMP RATE -
if ramplims !=0
    @constraint(m,ramplim[g=GENERATORS,t=t_notfirst, o = SCENARIOS],
        p[g,t,o] - p[g,t-1,o] <= (rampmax[g] * pmax[g]))
end

# ------ DR USAGE LIMITS
# Require that each period starts at the beginning of a new day
# to simplify computation

#STARTLIM
# number of startups per period
if startlim !=0
    # # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days, o = SCENARIOS],
    #     sum(v[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= startlim)
    # per period
    @constraint(m,limstart[g = 1:n_gdr, o = SCENARIOS],
        sum(start_num[GDR[g],t,o] for t = TIME) <= startlim)
end

#HOURLIM
#number of hours used per period
if hourlim != 0
    # # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days, o = SCENARIOS],
    #     sum(u[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= hourlim)
    # per period
    @constraint(m,limhrs[g = 1:n_gdr, o = SCENARIOS],
        sum(u[GDR[g],t,o] for t = TIME) <= hourlim)
end

#ENERGYLIM
# amount of energy used per period
if energylim != 0
    # per day
    # @constraint(m,[g = 1:n_gdr,d = 1:n_days ,o=SCENARIOS],
    #     sum(p[GDR[g],t,o] for t = (24*(d-1)+1):(24*d)) <= energylim)
    # per period
    @constraint(m,limenergy[g = 1:n_gdr, o=SCENARIOS],
        sum(p[GDR[g],t,o] for t = TIME) <= energylim)
end

#DURATIONLIM
# number of hours in a row that DR can be committed
if durationlim != 0
    l = Int(durationlim + 1)
    s = 1
    WINDOWS = [TIME[i:i+l-1] for i in range(1,s,Int(length(TIME)-l+1))]
    for i in range(1,s,Int(length(TIME)-l+1))
        print(TIME[i:i+l-1])
    end
    for W in WINDOWS
        @constraint(m,[g = 1:n_gdr, o = SCENARIOS],
            sum(u[GDR[g],t,o] for t = W) <= durationlim)
    end
end

### OBJECTIVE ###
@objective(m, Min, sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))

if debug
    error("just testing model so we are stopping here")
end

# -------------------------------------------
# SOLVE MODEL
# -------------------------------------------
@printf("\nSolving:\n")
status = solve(m)
@printf("Status: %s\n", status)

if !isdir(output_fol)
    mkdir(output_fol)
end

## save a copy of inputs & demand scenarios to the output file
writecsvmulti(read_inputs,output_fol,"inputfile",multiTF,periodsave)

# --------------------------------------------------------------------------------------
# SAVE OUTPUT
# -------------------------------------------

# production of DR
y = getvalue(p[GDR,:,:])

y_out = convert3dto2d(y,1, 3,  2,
    vcat([String("o$i") for i in 1:n_omega],"DR_IND","t"),
     genset[dr_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"DR_production",multiTF,periodsave)

# production of slow generators
zs = getvalue(p[GSL,:,:])

y_out = convert3dto2d(zs,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[slow_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"slow_production",multiTF,periodsave)

# production of fast generators
zf = getvalue(p[GF,:,:])
y_out = convert3dto2d(zf,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[fast_ind,:plantUnique])
writecsvmulti(y_out,output_fol,"fast_production",multiTF,periodsave)

# commitment of slow generators
w_out = getvalue(w)
wdf = DataFrame(transpose(w_out))
names!(wdf,[Symbol("$input") for input in genset[slow_ind,:plantUnique]])
writecsvmulti(wdf,output_fol,"slow_commitment",multiTF,periodsave)

# commitment of all generators
u_out = getvalue(u)
y_out = convert3dto2d(u_out,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[:plantUnique])
writecsvmulti(y_out,output_fol,"u_commitment",multiTF,periodsave)

# startup of all generators
v_out = getvalue(v)
y_out = convert3dto2d(v_out,1, 3, 2,
    vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
     genset[:plantUnique])
writecsvmulti(y_out,output_fol,"v_startup",multiTF,periodsave)

# calculate costs
totcost = getvalue(sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] + p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
totstartupcost = getvalue(sum(pro[o] *
    sum(start_num[g,t,o]*startup[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))
totvarcost = getvalue(sum(pro[o] *
    sum(p[g,t,o]*varcost[g] for g = GENERATORS, t = TIME)
    for o = SCENARIOS))

output_summary = DataFrame(TotalCost = totcost, TotStartupCst = totstartupcost,
                            TotVarCst = totvarcost,
                            FracStartupCost = totstartupcost/totcost,
                            FracVarCost = totvarcost/totcost,
                            timeseries = inputs[1,:timeseriesID],
                            stoch_scenario = inputs[1,:stochID])

writecsvmulti(output_summary,output_fol,"summary_stats",multiTF,periodsave)

# Get dual variables and save
# only works if not using true binary startup
if !trueBinaryStartup
    #supplydemand #2D
    sd_shadow = DataFrame(getdual(supplydemand))
    writecsvmulti(sd_shadow,output_fol,"supplydemand_shadow",multiTF,periodsave)

    ### mingen ###
    mingen_shadow = getdual(mingen)
    mingen_sf = convert3dto2d(mingen_shadow,1, 3,  2,
        vcat([String("o$i") for i in 1:n_omega],"GEN_IND","t"),
         genset[:,:plantUnique])
    # only save rows where there is a nonzero shadow price
    shadowsum = sum(convert(Array{Float64},mingen_sf[:,1:n_omega]),2)
    saveind = find(shadowsum .!= 0)
    writecsvmulti(mingen_sf[saveind,:],output_fol,"mingen_shadow",multiTF,periodsave)

    ### type3dr ###
    ## TODO: This does not work.
    if DRtype == 3
        type3dr_shadow = getdual(type3dr)
        type3dr_sf = convert3dto2d(type3dr_shadow,1, 3,  2,
            vcat([String("o$i") for i in 1:n_omega],"DR_IND","t"),
             genset[dr_ind,:plantUnique])
        # only save rows where there is a nonzero shadow price
        shadowsum = sum(convert(Array{Float64},type3dr_sf[:,1:n_omega]),2)
        saveind = find(shadowsum .!= 0)
        writecsvmulti(type3dr_sf[saveind,:],output_fol,"type3dr_shadow",multiTF,periodsave)
    end
end

max_stars_repo_name: pajalevi/Stoch-UC-model-for-DR

max_issues_repo_name: pajalevi/Stoch-UC-model-for-DR

max_forks_repo_name: pajalevi/Stoch-UC-model-for-DR


text:
function chebyshev_weights(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order[i],normalize_node(nodes[i],domain[:,i]))
  end

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      product = one(T)
      @inbounds for j = 1:N
        product *= poly[j][s[j],i[j]]
      end

      numerator   += f[s]*product
      denominator += product^2

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order[i],normalize_node(nodes[i],domain[:,i]))
  end

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        if s[j] === 1 || s[j] === n[j]
          scale = 0.5
        else
          scale = 1.0
        end
        temp = poly[j][s[j],i[j]]
        num *= temp*scale
        den *= (temp^2)*scale
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly       = Array{Array{T,2},1}(undef,N)
  complement = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i]       = chebyshev_polynomial(order[i],normalize_node(nodes[i],domain[:,i]))
    complement[i] = pinv(poly[i])'
  end

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        temp = complement[j][s[j],i[j]]
        num *= temp
        den *= temp*poly[j][s[j],i[j]]
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      product = one(T)
      @inbounds for j = 1:N
        product *= poly[j][s[j],i[j]]
      end

      numerator   += f[s]*product
      denominator += product^2

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        if s[j] == 1 || s[j] == n[j]
          scale = 0.5
        else
          scale = 1.0
        end
        temp = poly[j][s[j],i[j]]
        num *= temp*scale
        den *= (temp^2)*scale
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  complement = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    complement[i] = pinv(poly[i])'
  end

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        temp = complement[j][s[j],i[j]]
        num *= temp
        den *= temp*poly[j][s[j],i[j]]
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order,normalize_node(nodes[i],domain[:,i]))
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        product = one(T)
        @inbounds for j = 1:N
          product *= poly[j][s[j],i[j]]
        end

        numerator   += f[s]*product
        denominator += product^2

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order,normalize_node(nodes[i],domain[:,i]))
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          if s[j] == 1 || s[j] == n[j]
            scale = 0.5
          else
            scale = 1.0
          end
          temp = poly[j][s[j],i[j]]
          num *= temp*scale
          den *= (temp^2)*scale
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly       = Array{Array{T,2},1}(undef,N)
  complement = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i]       = chebyshev_polynomial(order,normalize_node(nodes[i],domain[:,i]))
    complement[i] = pinv(poly[i])'
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          temp = complement[j][s[j],i[j]]
          num *= temp
          den *= temp*poly[j][s[j],i[j]]
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T<:AbstractFloat,N,S<:Integer}

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        product = one(T)
        @inbounds for j = 1:N
          product *= poly[j][s[j],i[j]]
        end

        numerator   += f[s]*product
        denominator += product^2

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          if s[j] == 1 || s[j] == n[j]
            scale = 0.5
          else
            scale = 1.0
          end
          temp = poly[j][s[j],i[j]]
          num *= temp*scale
          den *= (temp^2)*scale
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T<:AbstractFloat,N,S<:Integer}

  complement = Array{Array{T,2},1}(undef,N)
  @inbounds for i = 1:N
    complement[i] = pinv(poly[i])'
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          temp = complement[j][s[j],i[j]]
          num *= temp
          den *= temp*poly[j][s[j],i[j]]
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order[i],normalize_node(nodes[i],domain[:,i]))
  end

  weights = zeros(Tuple(order.+1))

  @inbounds @sync @qthreads for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      product = one(T)
      @inbounds for j = 1:N
        product *= poly[j][s[j],i[j]]
      end

      numerator   += f[s]*product
      denominator += product^2

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order[i],normalize_node(nodes[i],domain[:,i]))
  end

  weights = zeros(Tuple(order.+1))

  @inbounds @sync @qthreads for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        if s[j] == 1 || s[j] == n[j]
          scale = 0.5
        else
          scale = 1.0
        end
        temp = poly[j][s[j],i[j]]
        num *= temp*scale
        den *= (temp^2)*scale
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly       = Array{Array{T,2},1}(undef,N)
  complement = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i]       = chebyshev_polynomial(order[i],normalize_node(nodes[i],domain[:,i]))
    complement[i] = pinv(poly[i])'
  end

  weights = zeros(Tuple(order.+1))

  @inbounds @sync @qthreads for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        temp = complement[j][s[j],i[j]]
        num *= temp
        den *= temp*poly[j][s[j],i[j]]
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds @sync @qthreads for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      product = one(T)
      @inbounds for j = 1:N
        product *= poly[j][s[j],i[j]]
      end

      numerator   += f[s]*product
      denominator += product^2

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds @sync @qthreads for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        if s[j] == 1 || s[j] == n[j]
          scale = 0.5
        else
          scale = 1.0
        end
        temp = poly[j][s[j],i[j]]
        num *= temp*scale
        den *= (temp^2)*scale
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  complement = Array{Array{T,2},1}(undef,N)
  @inbounds for i = 1:N
    complement[i] = pinv(poly[i])'
  end

  weights = Array{T,N}(undef,Tuple(order.+1))

  @inbounds @sync @qthreads for i in CartesianIndices(weights)

    numerator   = zero(T)
    denominator = zero(T)

    @inbounds for s in CartesianIndices(f)

      num = f[s]
      den = one(T)
      @inbounds for j = 1:N
        temp = complement[j][s[j],i[j]]
        num *= temp
        den *= temp*poly[j][s[j],i[j]]
      end

      numerator   += num
      denominator += den

    end

    weights[i] = numerator/denominator

  end

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order,normalize_node(nodes[i],domain[:,i]))
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds @sync @qthreads for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        product = one(T)
        @inbounds for j = 1:N
          product *= poly[j][s[j],i[j]]
        end

        numerator   += f[s]*product
        denominator += product^2

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  poly = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i] = chebyshev_polynomial(order,normalize_node(nodes[i],domain[:,i]))
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds @sync @qthreads for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          if s[j] == 1 || s[j] == n[j]
            scale = 0.5
          else
            scale = 1.0
          end
          temp = poly[j][s[j],i[j]]
          num *= temp*scale
          den *= (temp^2)*scale
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},nodes::NTuple{N,Array{T,1}},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  poly       = Array{Array{T,2},1}(undef,N)
  complement = Array{Array{T,2},1}(undef,N)

  @inbounds for i = 1:N
    poly[i]       = chebyshev_polynomial(order,normalize_node(nodes[i],domain[:,i]))
    complement[i] = pinv(poly[i])'
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds @sync @qthreads for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          temp = complement[j][s[j],i[j]]
          num *= temp
          den *= temp*poly[j][s[j],i[j]]
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T<:AbstractFloat,N,S<:Integer}

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds @sync @qthreads for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        product = one(T)
        @inbounds for j = 1:N
          product *= poly[j][s[j],i[j]]
        end

        numerator   += f[s]*product
        denominator += product^2

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T<:AbstractFloat,N,S<:Integer}

  n = size(f)

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds @sync @qthreads for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          if s[j] == 1 || s[j] == n[j]
            scale = 0.5
          else
            scale = 1.0
          end
          temp = poly[j][s[j],i[j]]
          num *= temp*scale
          den *= (temp^2)*scale
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},poly::NTuple{N,Array{T,2}},order::S) where {T<:AbstractFloat,N,S<:Integer}

  complement = Array{Array{T,2},1}(undef,N)
  @inbounds for i = 1:N
    complement[i] = pinv(poly[i])'
  end

  ord = (order,)
  for i = 2:N
    ord = (ord...,order)
  end

  weights = Array{T,N}(undef,ord.+1)

  @inbounds @sync @qthreads for i in CartesianIndices(weights)
    if sum(Tuple(i)) <= order+N

      numerator   = zero(T)
      denominator = zero(T)

      @inbounds for s in CartesianIndices(f)

        num = f[s]
        den = one(T)
        @inbounds for j = 1:N
          temp = complement[j][s[j],i[j]]
          num *= temp
          den *= temp*poly[j][s[j],i[j]]
        end

        numerator   += num
        denominator += den

      end

      weights[i] = numerator/denominator

    else
      weights[i] = zero(T)
    end

  end

  return weights

end

function chebyshev_weights(cheb::ChebInterpRoots)

  weights = chebyshev_weights_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  return weights

end

function chebyshev_weights(cheb::ChebInterpExtrema)

  weights = chebyshev_weights_extrema_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  return weights

end

function chebyshev_weights(cheb::ChebInterpExtended)

  weights = chebyshev_weights_extended_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  return weights

end

function chebyshev_weights(cheb::ChebInterpVertesi)

  weights = chebyshev_weights_vertesi_threaded(cheb.data,cheb.nodes,cheb.order,cheb.domain)

  return weights

end

# Functions for the one-variable case where the nodes are a vector

function chebyshev_weights(f::AbstractArray{T,1},nodes::Array{T,1},order::Union{S,Array{S,1}},domain=[one(T);-one(T)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(nodes,),order,domain)

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,1},nodes::Array{T,1},order::Union{S,Array{S,1}},domain=[one(T);-one(T)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights_extrema(f,(nodes,),order,domain)

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,1},nodes::Array{T,1},order::Union{S,Array{S,1}},domain=[one(T);-one(T)]) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights_extended(f,(nodes,),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,1},poly::Array{T,2},order::Union{S,Array{S,1}}) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights(f,(poly,),order)

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,1},poly::Array{T,2},order::Union{S,Array{S,1}}) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights_extrema(f,(poly,),order)

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,1},poly::Array{T,2},order::Union{S,Array{S,1}}) where {T<:AbstractFloat,S<:Integer}

  weights = chebyshev_weights_extended(f,(poly,),order)

  return weights

end

# Functions that allow the nodes to be in an array of arrays

function chebyshev_weights(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extrema(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extended(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended(f,tuple(poly...),order)

  return weights

end

# Threaded functions

function chebyshev_weights_threaded(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_threaded(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema_threaded(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::Array{S,1},domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended_threaded(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_threaded(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema_threaded(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::Array{S,1}) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended_threaded(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_threaded(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema_threaded(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},nodes::Array{Array{T,1},1},order::S,domain=[ones(T,1,N);-ones(T,1,N)]) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended_threaded(f,tuple(nodes...),order,domain)

  return weights

end

function chebyshev_weights_threaded(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_threaded(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extrema_threaded(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extrema_threaded(f,tuple(poly...),order)

  return weights

end

function chebyshev_weights_extended_threaded(f::AbstractArray{T,N},poly::Array{Array{T,2},1},order::S) where {T<:AbstractFloat,N,S<:Integer}

  weights = chebyshev_weights_extended_threaded(f,tuple(poly...),order)

  return weights

end

const chebyshev_weights_vertesi          = chebyshev_weights_extended
const chebyshev_weights_vertesi_threaded = chebyshev_weights_extended_threaded

max_stars_repo_name: RJDennis/ChebyshevApprox

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox