text:
## Orthogonal polynomials

abstract PolynomialSpace{D} <: RealUnivariateSpace{D}

Multiplication{U<:PolynomialSpace}(f::Fun{U},sp::PolynomialSpace)=ConcreteMultiplication(f,sp)
bandinds{U<:PolynomialSpace,V<:PolynomialSpace}(M::ConcreteMultiplication{U,V})=(1-ncoefficients(M.f),ncoefficients(M.f)-1)
rangespace{U<:PolynomialSpace,V<:PolynomialSpace}(M::ConcreteMultiplication{U,V})=domainspace(M)

# All polynomials contain constant
union_rule(A::ConstantSpace,B::PolynomialSpace)=B
Base.promote_rule{T<:Number,S<:PolynomialSpace,V}(::Type{Fun{S,V}},::Type{T})=Fun{S,promote_type(V,T)}
Base.promote_rule{T<:Number,S<:PolynomialSpace}(::Type{Fun{S}},::Type{T})=Fun{S,T}

## Evaluation

function evaluate(f::AbstractVector,S::PolynomialSpace,x)
    if x in domain(S)
        clenshaw(S,f,tocanonical(S,x))
    else
        zero(eltype(f))
    end
end

evaluate(f::AbstractVector,S::PolynomialSpace,x::AbstractArray)=map(y->evaluate(f,S,y),x)

# we need the ... for multi-dimensional
evaluate(f::AbstractVector,S::PolynomialSpace,x,y,z...)=
    evaluate(f,S,Vec(x,y,z...))

function evaluate(f::AbstractVector,S::PolynomialSpace,x::Fun)
    if issubset(Interval(minimum(x),maximum(x)),domain(S))
        clenshaw(S,f,tocanonical(S,x))
    else
        error("Implement splitatpoints for evaluate ")
    end
end

## Extrapolation

extrapolate(f::AbstractVector,S::PolynomialSpace,x) = clenshaw(S,f,tocanonical(S,x))
extrapolate(f::AbstractVector,S::PolynomialSpace,x::AbstractArray)=map(y->extrapolate(f,S,y),x)

######
# Recurrence encodes the recurrence coefficients
# or equivalently multiplication by x
######
immutable Recurrence{S,T} <: TridiagonalOperator{T}
    space::S
end

Recurrence(sp)=Recurrence{typeof(sp),promote_type(eltype(sp),eltype(domain(sp)))}(sp)

Base.convert{T,S}(::Type{Operator{T}},J::Recurrence{S})=Recurrence{S,T}(J.space)

#####
# recα/β/γ are given by
#       x p_{n-1} =γ_n p_{n-2} + α_n p_{n-1} +  p_n β_n
#####

function getindex{S,T}(R::Recurrence{S,T},k::Integer,j::Integer)
    if j==k-1
        recβ(T,R.space,k-1)
    elseif j==k
        recα(T,R.space,k)
    elseif j==k+1
        recγ(T,R.space,k+1)
    else
        zero(T)
    end
end

######
# JacobiZ encodes [BasisFunctional(1);(J-z*I)[2:end,:]]
# where J is the Jacobi operator
######
immutable JacobiZ{S,T} <: TridiagonalOperator{T}
    space::S
    z::T
end

JacobiZ(sp,z) =
    (T = promote_type(eltype(sp),eltype(domain(sp)),typeof(z)); JacobiZ{typeof(sp),T}(sp,T(z)))

Base.convert{T,S}(::Type{Operator{T}},J::JacobiZ{S}) = JacobiZ{S,T}(J.space,J.z)

domainspace(::JacobiZ) = ℓ⁰
rangespace(::JacobiZ) = ℓ⁰

#####
# recα/β/γ are given by
#       x p_{n-1} =γ_n p_{n-2} + α_n p_{n-1} +  p_n β_n
#####

function getindex{S,T}(J::JacobiZ{S,T},k::Integer,j::Integer)
    if j==k-1
        recγ(T,J.space,k)
    elseif j==k
        k == 1 ? one(T) : recα(T,J.space,k)-J.z
    elseif j==k+1 && k > 1
        recβ(T,J.space,k)
    else
        zero(T)
    end
end

#####
# Multiplication can be built from recurrence coefficients
#  multiplication by J=M[x] is Recurrence operator
#  and assume p_0 = 1, then we have
#
#   M[p_0] = 1
#   M[p_1] = (J/β_1 - α_1/β_1)*M[p_0]
#   M[p_k] = (J/β_k - α_k/β_k)*M[p_{k-1}] - γ_k/β_k*M[p_{k-2}]
#####

getindex{PS<:PolynomialSpace,T,C<:PolynomialSpace}(M::ConcreteMultiplication{C,PS,T},k::Integer,j::Integer) = M[k:k,j:j][1,1]

function Base.convert{PS<:PolynomialSpace,V,T,C<:PolynomialSpace}(::Type{BandedMatrix},S::SubOperator{T,ConcreteMultiplication{C,PS,V,T},
                                                                            Tuple{UnitRange{Int},UnitRange{Int}}})
    M=parent(S)
    kr,jr=parentindexes(S)

    A=bzeros(S)

    a=coefficients(M.f)

    shft=bandshift(A)

    for k=kr ∩ jr
        A[k-kr[1]+1,k-jr[1]+1]=a[1]
    end

    if length(a) > 1
        sp=M.space
        fsp=space(M.f)
        jkr=max(1,min(kr[1],jr[1])-length(a)+1):max(kr[end],jr[end])+length(a)-1

        #Multiplication is transpose
        J=Recurrence(sp)[jkr,jkr]

        # the sub ranges of jkr that correspond to kr, jr
        kr2,jr2=kr-jkr[1]+1,jr-jkr[1]+1

        C0=beye(size(J,1),size(J,2),0,0)
        C1=(1/recβ(T,fsp,1))*J-(recα(T,fsp,1)/recβ(T,fsp,1))*C0

        BLAS.axpy!(a[2],view(C1,kr2,jr2),A)

        for k=1:length(a)-2
            rβ=recβ(T,fsp,k+1)
            C1,C0=(1/rβ)*J*C1-(recα(T,fsp,k+1)/rβ)*C1-(recγ(T,fsp,k+1)/rβ)*C0,C1
            BLAS.axpy!(a[k+2],view(C1,kr2,jr2),A)
        end
    end
    A
end

## All polynomial spaces can be converted provided spaces match

isconvertible(a::PolynomialSpace,b::PolynomialSpace)=domain(a)==domain(b)
union_rule{D}(a::PolynomialSpace{D},b::PolynomialSpace{D})=domainscompatible(a,b)?(a<b?a:b):NoSpace()   # the union of two polys is always a poly

## General clenshaw
clenshaw(sp::PolynomialSpace,c::AbstractVector,x::AbstractArray) = clenshaw(c,x,
                                            ClenshawPlan(promote_type(eltype(c),eltype(x)),sp,length(c),length(x)))
clenshaw(sp::PolynomialSpace,c::AbstractMatrix,x::AbstractArray) = clenshaw(c,x,ClenshawPlan(promote_type(eltype(c),eltype(x)),sp,size(c,1),length(x)))
clenshaw(sp::PolynomialSpace,c::AbstractMatrix,x) = clenshaw(c,x,ClenshawPlan(promote_type(eltype(c),typeof(x)),sp,size(c,1),size(c,2)))

clenshaw!(sp::PolynomialSpace,c::AbstractVector,x::AbstractVector)=clenshaw!(c,x,ClenshawPlan(promote_type(eltype(c),eltype(x)),sp,length(x)))

function clenshaw(sp::PolynomialSpace,c::AbstractVector,x)
    N,T = length(c),promote_type(eltype(c),typeof(x))
    TT = eltype(T)
    if isempty(c)
        return zero(T)
    end

    bk1,bk2 = zero(T),zero(T)
    A,B,C = recA(TT,sp,N-1),recB(TT,sp,N-1),recC(TT,sp,N)
    for k = N:-1:2
        bk2, bk1 = bk1, muladd(muladd(A,x,B),bk1,muladd(-C,bk2,c[k])) # muladd(-C,bk2,muladd(muladd(A,x,B),bk1,c[k])) # (A*x+B)*bk1+c[k]-C*bk2
        A,B,C = recA(TT,sp,k-2),recB(TT,sp,k-2),recC(TT,sp,k-1)
    end
    muladd(muladd(A,x,B),bk1,muladd(-C,bk2,c[1])) # muladd(-C,bk2,muladd(muladd(A,x,B),bk1,c[1])) # (A*x+B)*bk1+c[1]-C*bk2
end

# Supports constants in operators
promoterangespace{CS<:ConstantSpace}(M::ConcreteMultiplication{CS,UnsetSpace},
                                                ps::PolynomialSpace) =
                        promoterangespace(Multiplication(M.f,space(M.f)),ps)


text:
using ModelingToolkit

@variables x1 x2
@parameters dt

eq = [
    -1 + 9 *x1 - 2*x1^3 + 9 *x2 - 2*x2^3,
    1 - 11*x1 + 2*x1^3 + 11*x2 - 2*x2^3
]

ModelingToolkit.gradient.(eq, [[x1,x2]])

#=
os = OptimizationSystem(eq, [], [x1,x2])

calculate_gradient(os)
=#

ext: jl


text:
using MuZero
include("game.jl")

g = Connect4()
s1 = step!(g, 4)
s2 = step!(g, 4)
s3 = step!(g, 3)
s2 = get_state(g)
s2s = s2 |> svg
g |> show
g |> svg

get_state(g) |> svg

g = Connect4()
step!(g, 4)
step!(g, 7)
step!(g, 1)


text:
# News Boy Problem solved by L-Shaped Benders Decomposition
# Author: Edward J. Xu, edxu96@outlook.com
# Date: April 5th, 2019

function get_data()
    demand = [12, 14, 16, 18, 20, 22, 24, 26, 28, 30] # Demand of newspapers in each scenario
    lengthS = length(demand)
    vec_prob = [0.05, 0.10, 0.10, 0.10, 0.15, 0.15, 0.10, 0.10, 0.10, 0.05 ] # probability of scenario
    c = 20 # purchase price
    p = 70 # selling price
    h = 10 # scrap value
    y = 20
    # 2. Data Transfer -----------------------------------------------------------------------------------------------------
    n_x = 10  # Number of all the x
    n_y = 1
    vec_min_y = hcat([0])
    vec_max_y = hcat([30])
    vec_f = hcat(c - 10)
    vec_pi = hcat(vec_prob)
    mat_c = zeros(10, 1, 1)
    for i = 1:10
        mat_c[i, :, :] = hcat([- 70 + 10])
    end
    mat_h = zeros(10, 2, 1)
    for i = 1: 10
        mat_h[i, 1, :] = hcat(- demand[i])
        mat_h[i, 2, :] = hcat(0)
    end
    mat3_t = zeros(10, 2, 1)
    for i = 1:10
        mat3_t[i, 1, :] = hcat(0)
        mat3_t[i, 2, :] = hcat(1)
    end
    mat3_w = - ones(10, 2, 1)
    epsilon = 0.001
    timesIterationMax = 500
    return
end

function main()

    BendersLshaped(; n_x = n_x, n_y = n_y, vec_min_y = vec_min_y, vec_max_y = vec_max_y, vec_f = vec_f,
        vec_pi = vec_pi, mat_c = mat_c, mat_h = mat_h, mat3_t = mat3_t, mat3_w = mat3_w,
        epsilon = epsilon, timesIterationMax = timesIterationMax)
    # answer: obj = - 976
end


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

using Documenter
using OpticSim

makedocs(
    sitename = "OpticSim.jl",
    format = Documenter.HTML(
        # prettyurls = get(ENV, "CI", nothing) == "true",
        assets = [asset("assets/logo.svg", class = :ico, islocal = true)],
    ),
    modules = [OpticSim],
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Geometry" => [
            "Basic Types" => "basic_types.md",
            "Primitives" => "primitives.md",
            "CSG" => "csg.md"
        ],
        "Optical" => [
            "Systems" => "systems.md",
            "Emitters" => "emitters.md",
            "Interfaces" => "interfaces.md",
            "Lenses" => "lenses.md"
        ],
        "Visualization" => "vis.md",
        "Glass Functions" => "glasscat.md",
        "Optimization" => "optimization.md",
        "Cloud Execution" => "cloud.md",
        "Notebook utilities" => "notebooksutils.md",
        "Reference" => "ref.md",
        "Roadmap" => "roadmap.md"
    ],
    expandfirst = ["glasscat.md", "systems.md", "vis.md"])

deploydocs(
    repo = "github.com/microsoft/OpticSim.jl.git",
    devbranch = "main",
    push_preview = true,
)

# function children(m::Module)
#     ns = names(m, imported = false, all = true)
#     ms = []
#     for n in ns
#         try
#             x = Core.eval(m, n)
#             if x isa Module
#                 if(x != OpticSim.GlassCat)
#                 println("x $x")
#                 push!(ms, x)
#             end
#         end
#         catch
#         end
#     end
#     return ms
# end

#WARNING: I think this code creates a type for each glass name, which overwrites the definition in the src files where each glass name corresponds to an integer indexing into a glass table. Obviously nothing works after this. Needs major surgery to generate glass documentation at the same time as the documentation for everything else.

# # write a code file for the catalog with docstrings
# io = open(joinpath(@__DIR__, "../src/AGFGlassCatDocs.jl"), "w")
# catalogs = children(OpticSim.GlassCat)
# cat_pages = []
# for catname in catalogs
#     println(catalogs)
#     eval_string = ["module $(nameof(catname))"]
#     escape_catalog_name = replace(string(nameof(catname)), "_" => "\\_")
#     push!(cat_pages, escape_catalog_name => "$(nameof(catname)).md")
#     iomd = open(joinpath(@__DIR__, "src/$(nameof(catname)).md"), "w")
#     write(iomd, "# $escape_catalog_name\n\n")
#     write(iomd, "```@raw html\n")
#     write(iomd, "<style>article p {display:flex;justify-content:space-between;}</style>\n")
#     write(iomd, "```\n")
#     write(iomd, "```@docs\n")
#     catalog_module = eval(catname)
#     glass_names = names(catalog_module, all = true, imported = false)
#     for glass_name in glass_names
#         glass_name_str = string(glass_name)
#         if !occursin("#", glass_name_str) && glass_name_str != "eval" && glass_name_str != "include" && glass_name != nameof(catname) && !in(glass_name,(:MODEL, :MIL, :AGF, :OTHER, :AIR))
#             glass = Core.eval(catalog_module, glass_name)

#             temp = typeof(glass)
# println("type of glass $temp")

#             let io = IOBuffer()
#                 OpticSim.GlassCat.docstring(io, glass)
#                 infostr = String(take!(io))
#                 push!(eval_string, "\"\"\"\n$infostr\n\"\"\"")
#             end
#             push!(eval_string, "function $glass_name_str()\nend")
#             write(iomd, "$catname.$glass_name_str\n")
#         end
#     end
#     write(iomd, "```\n")
#     push!(eval_string, "end") # module
#     eval_string = join(eval_string, "\n")
#     write(io, eval_string * "\n")
#     close(iomd)
# end
# close(io)

# # include source with docstrings
# OpticSim.GlassCat.include(joinpath(@__DIR__, "../src/AGFGlassCatDocs.jl"))

# clean up
# rm(joinpath(@__DIR__, "../src/AGFGlassCatDocs.jl"))
# for p in cat_pages
#     rm(joinpath(@__DIR__, "src/" * p[2]))
# end


text:
module TimeScales

export TimeScale

"""
All timescales are subtypes of the abstract type `TimeScale`.
The following timescales are defined:

* `UTC` — Coordinated Universal Time
* `UT1` — Universal Time
* `TAI` — International Atomic Time
* `TT` — Terrestrial Time
* `TCG` — Geocentric Coordinate Time
* `TCB` — Barycentric Coordinate Time
* `TDB` — Barycentric Dynamical Time
"""
abstract type TimeScale end

const scales = (
    :CoordinatedUniversalTime,
    :UniversalTime,
    :InternationalAtomicTime,
    :TerrestrialTime,
    :GeocentricCoordinateTime,
    :BarycentricCoordinateTime,
    :BarycentricDynamicalTime,
)

const acronyms = (
    :UTC,
    :UT1,
    :TAI,
    :TT,
    :TCG,
    :TCB,
    :TDB,
)

for (acronym, scale) in zip(acronyms, scales)
    name = String(acronym)
    @eval begin
        struct $scale <: TimeScale end
        const $acronym = $scale()
        Base.show(io::IO, ::$scale) = print(io, $name)
        export $scale, $acronym
    end
end

end


text:
module CUDA

	import Base.length, Base.size
	import Base.copy!

	export

	# errors
	CuDriverError, description,

	# base
	@cucall, driver_version,

	# devices
	CuDevice, CuCapability, dev_count, name, totalmem, attribute, capability,
	list_devices,

	# context
	CuContext, create_context, destroy, push, pop,

	# module
	CuModule, CuFunction, unload,

	# stream
	CuStream, synchronize,

	# execution
	launch,

	# arrays
	CuPtr, CuArray, free, to_host

	include("errors.jl")

	include("base.jl")
	include("devices.jl")
	include("context.jl")
	include("module.jl")
	include("stream.jl")
	include("execution.jl")

	include("arrays.jl")
end