text:
mutable struct SlotInfo
    defs::Vector{Int}
    uses::Vector{Int}
    any_newvar::Bool
end
SlotInfo() = SlotInfo(Int[], Int[], false)

function scan_entry!(result::Vector{SlotInfo}, idx::Int, @nospecialize(stmt))
    # NewVarNodes count as defs for the purpose
    # of liveness analysis (i.e. they kill use chains)
    if isa(stmt, NewvarNode)
        result[slot_id(stmt.slot)].any_newvar = true
        push!(result[slot_id(stmt.slot)].defs, idx)
        return
    elseif isexpr(stmt, :(=))
        if isa(stmt.args[1], SlotNumber)
            push!(result[slot_id(stmt.args[1])].defs, idx)
        end
        stmt = stmt.args[2]
    end
    if isa(stmt, Union{SlotNumber, TypedSlot})
        push!(result[slot_id(stmt)].uses, idx)
        return
    end
    for op in userefs(stmt)
        val = op[]
        if isa(val, Union{SlotNumber, TypedSlot})
            push!(result[slot_id(val)].uses, idx)
        end
    end
end

function lift_defuse(cfg::CFG, defuse)
    map(defuse) do slot
        SlotInfo(
            Int[block_for_inst(cfg, x) for x in slot.defs],
            Int[block_for_inst(cfg, x) for x in slot.uses],
            slot.any_newvar
        )
    end
end

@inline slot_id(s) = isa(s, SlotNumber) ? (s::SlotNumber).id : (s::TypedSlot).id
function scan_slot_def_use(nargs, ci::CodeInfo)
    nslots = length(ci.slotnames)
    result = SlotInfo[SlotInfo() for i = 1:nslots]
    # Set defs for arguments
    for var in result[1:(1+nargs)]
        push!(var.defs, 0)
    end
    for (idx, stmt) in Iterators.enumerate(ci.code)
        scan_entry!(result, idx, stmt)
    end
    result
end

function renumber_ssa(stmt::SSAValue, ssanums::Vector{Any}, new_ssa::Bool=false, used_ssa::Union{Nothing, Vector{Int}}=nothing)
    id = stmt.id + (new_ssa ? 0 : 1)
    if id > length(ssanums)
        return stmt
    end
    val = ssanums[id]
    if isa(val, SSAValue) && used_ssa !== nothing
        used_ssa[val.id] += 1
    end
    return val
end

function renumber_ssa!(@nospecialize(stmt), ssanums::Vector{Any}, new_ssa::Bool=false, used_ssa::Union{Nothing, Vector{Int}}=nothing)
    isa(stmt, SSAValue) && return renumber_ssa(stmt, ssanums, new_ssa, used_ssa)
    return ssamap(val->renumber_ssa(val, ssanums, new_ssa, used_ssa), stmt)
end

function make_ssa!(ci::CodeInfo, idx, slot, @nospecialize(typ))
    (idx == 0) && return Argument(slot)
    stmt = ci.code[idx]
    @assert isexpr(stmt, :(=))
    push!(ci.ssavaluetypes, typ)
    ssa = length(ci.ssavaluetypes)-1
    stmt.args[1] = SSAValue(ssa)
    ssa
end

struct UndefToken
end
const undef_token = UndefToken()

function new_to_regular(@nospecialize(stmt))
    if isa(stmt, NewSSAValue)
        return SSAValue(stmt.id)
    end
    urs = userefs(stmt)
    urs === () && return stmt
    for op in urs
        val = op[]
        if isa(val, NewSSAValue)
            op[] = SSAValue(val.id)
        end
    end
    urs[]
end

function fixup_slot!(ir::IRCode, ci::CodeInfo, idx::Int, slot::Int, @nospecialize(stmt::Union{SlotNumber, TypedSlot}), @nospecialize(ssa))
    # We don't really have the information here to get rid of these.
    # We'll do so later
    if ssa === undef_token
        insert_node!(ir, idx, Any, Expr(:throw_undef_if_not, ci.slotnames[slot], false))
        return undef_token
    end
    if !isa(ssa, Argument) && !(ssa === nothing) && ((ci.slotflags[slot] & SLOT_USEDUNDEF) != 0)
        insert_node!(ir, idx, Any, Expr(:undefcheck, ci.slotnames[slot], ssa))
    end
    if isa(stmt, SlotNumber)
        return ssa
    elseif isa(stmt, TypedSlot)
        return NewSSAValue(insert_node!(ir, idx, stmt.typ, PiNode(ssa, stmt.typ)).id)
    end
end

function fixemup!(cond, rename, ir::IRCode, ci::CodeInfo, idx::Int, @nospecialize(stmt))
    if isa(stmt, Union{SlotNumber, TypedSlot}) && cond(stmt)
        return fixup_slot!(ir, ci, idx, slot_id(stmt), stmt, rename(stmt))
    end
    if isexpr(stmt, :(=))
        stmt.args[2] = fixemup!(cond, rename, ir, ci, idx, stmt.args[2])
        return stmt
    end
    if isa(stmt, PhiNode)
        for i = 1:length(stmt.edges)
            isassigned(stmt.values, i) || continue
            val = stmt.values[i]
            isa(val, Union{SlotNumber, TypedSlot}) || continue
            cond(val) || continue
            bb_idx = block_for_inst(ir.cfg, stmt.edges[i])
            from_bb_terminator = last(ir.cfg.blocks[bb_idx].stmts)
            stmt.values[i] = fixup_slot!(ir, ci, from_bb_terminator, slot_id(val), val, rename(val))
        end
        return stmt
    end
    if isexpr(stmt, :isdefined)
        val = stmt.args[1]
        if isa(val, Union{SlotNumber, TypedSlot})
            slot = slot_id(val)
            if (ci.slotflags[slot] & SLOT_USEDUNDEF) == 0
                return true
            else
                ssa = rename(val)
                if ssa === undef_token
                    return false
                elseif !isa(ssa, SSAValue) && !isa(ssa, NewSSAValue)
                    return true
                end
            end
            stmt.args[1] = ssa
        end
        return stmt
    end
    urs = userefs(stmt)
    urs === () && return stmt
    for op in urs
        val = op[]
        if isa(val, Union{SlotNumber, TypedSlot}) && cond(val)
            x = fixup_slot!(ir, ci, idx, slot_id(val), val, rename(val))
            # We inserted an undef error node. Delete subsequent statement
            # to avoid confusing the optimizer
            if x === undef_token
                return nothing
            end
            op[] = x
        end
    end
    urs[]
end

function fixup_uses!(ir::IRCode, ci::CodeInfo, uses::Vector{Int}, slot, @nospecialize(ssa))
    for use in uses
        ci.code[use] = fixemup!(stmt->slot_id(stmt)==slot, stmt->ssa, ir, ci, use, ci.code[use])
    end
end

function rename_uses!(ir::IRCode, ci::CodeInfo, idx::Int, @nospecialize(stmt), renames::Vector{Any})
    return fixemup!(stmt->true, stmt->renames[slot_id(stmt)], ir, ci, idx, stmt)
end

function strip_trailing_junk!(code::Vector{Any}, lines::Vector{Int})
    # Remove `nothing`s at the end, we don't handle them well
    # (we expect the last instruction to be a terminator)
    for i = length(code):-1:1
        if code[i] !== nothing
            resize!(code, i)
            resize!(lines, i)
            break
        end
    end
    # If the last instruction is not a terminator, add one. This can
    # happen for implicit return on dead branches.
    term = code[end]
    if !isa(term, GotoIfNot) && !isa(term, GotoNode) && !isa(term, ReturnNode)
        push!(code, ReturnNode{Any}())
        push!(lines, 0)
    end
    return code
end

struct DelayedTyp
    phi::NewSSAValue
end

# maybe use expr_type?
function typ_for_val(@nospecialize(val), ci::CodeInfo)
    isa(val, Expr) && return val.typ
    isa(val, GlobalRef) && return abstract_eval_global(val.mod, val.name)
    isa(val, SSAValue) && return ci.ssavaluetypes[val.id+1]
    isa(val, Argument) && return ci.slottypes[val.n]
    isa(val, NewSSAValue) && return DelayedTyp(val)
    isa(val, QuoteNode) && return Const(val.value)
    isa(val, Union{Symbol, PiNode, PhiNode, SlotNumber, TypedSlot}) && error("unexpected val type")
    return Const(val)
end

# Run iterated dominance frontier
function idf(cfg::CFG, defuse, domtree::DomTree, slot::Int)
    # This should be a priority queue, but TODO - sorted array for now
    defs = defuse[slot].defs
    pq = Tuple{Int, Int}[(defs[i], domtree.nodes[defs[i]].level) for i in 1:length(defs)]
    sort!(pq, by=x->x[2])
    phiblocks = Int[]
    processed = IdSet{Int}()
    while !isempty(pq)
        node, level = pop!(pq)
        worklist = Int[]
        visited = IdSet{Int}()
        push!(worklist, node)
        while !isempty(worklist)
            active = pop!(worklist)
            for succ in cfg.blocks[active].succs
                succ_level = domtree.nodes[succ].level
                succ_level > level && continue
                succ in processed && continue
                push!(processed, succ)
                # <- TODO: Use liveness here
                push!(phiblocks, succ)
                if !(succ in defuse[slot].defs)
                    push!(pq, (succ, succ_level))
                    sort!(pq, by=x->x[2])
                end
            end

            for child in domtree.nodes[active].children
                child in visited && continue
                push!(visited, child)
                push!(worklist, child)
            end
        end
    end
    phiblocks
end

function construct_ssa!(ci::CodeInfo, ir::IRCode, domtree::DomTree, defuse, nargs::Int)
    cfg = ir.cfg
    left = Int[]
    defuse_blocks = lift_defuse(ir.cfg, defuse)
    phi_slots = Vector{Int}[Vector{Int}() for _ = 1:length(ir.cfg.blocks)]
    phi_nodes = Vector{Pair{Int,PhiNode}}[Vector{Pair{Int,PhiNode}}() for _ = 1:length(cfg.blocks)]
    phi_ssas = SSAValue[]
    for (idx, slot) in Iterators.enumerate(defuse)
        # No uses => no need for phi nodes
        isempty(slot.uses) && continue
        # TODO: Restore this optimization
        if false # length(slot.defs) == 1 && slot.any_newvar
            if slot.defs[] == 0
                typ = ci.slottypes[idx]
                ssaval = Argument(idx)
                fixup_uses!(ir, ci, slot.uses, idx, ssaval)
            elseif isa(ci.code[slot.defs[]], NewvarNode)
                typ = MaybeUndef(Union{})
                ssaval = nothing
                for use in slot.uses[]
                    insert_node!(ir, use, Union{}, Expr(:throw_undef_if_not, ci.slotnames[idx], false))
                end
                fixup_uses!(ir, ci, slot.uses, idx, nothing)
            else
                val = ci.code[slot.defs[]].args[2]
                typ = typ_for_val(val, ci)
                ssaval = SSAValue(make_ssa!(ci, slot.defs[], idx, typ))
                fixup_uses!(ir, ci, slot.uses, idx, ssaval)
            end
            continue
        end
        # TODO: Perform liveness here to eliminate dead phi nodes
        phiblocks = idf(cfg, defuse_blocks, domtree, idx)
        for block in phiblocks
            push!(phi_slots[block], idx)
            node = PhiNode()
            ssa = insert_node!(ir, first_insert_for_bb(ci.code, cfg, block), Union{}, node)
            push!(phi_nodes[block], ssa.id=>node)
        end
        push!(left, idx)
    end
    # Perform SSA renaming
    initial_incoming_vals = Any[
        if 0 in defuse[x].defs
            Argument(x)
        elseif !defuse[x].any_newvar
            undef_token
        else
            SSAValue(-1)
        end for x in 1:length(ci.slotnames)
    ]
    worklist = Any[(1, 0, initial_incoming_vals)]
    visited = IdSet{Int}()
    type_refine_phi = IdSet{Int}()
    while !isempty(worklist)
        (item, pred, incoming_vals) = pop!(worklist)
        # Insert phi nodes if necessary
        for (idx, slot) in Iterators.enumerate(phi_slots[item])
            ssaval, node = phi_nodes[item][idx]
            incoming_val = incoming_vals[slot]
            if incoming_val == SSAValue(-1)
                # Optimistically omit this path.
                # Liveness analysis would probably have prevented us from inserting this phi node
                continue
            end
            push!(node.edges, pred)
            if incoming_val == undef_token
                resize!(node.values, length(node.values)+1)
            else
                push!(node.values, incoming_val)
            end
            # TODO: Remove the next line, it shouldn't be necessary
            push!(type_refine_phi, ssaval)
            if isa(incoming_val, NewSSAValue)
                push!(type_refine_phi, ssaval)
            end
            typ = incoming_val == undef_token ? MaybeUndef(Union{}) : typ_for_val(incoming_val, ci)
            new_node_id = ssaval - length(ir.stmts)
            old_insert, old_typ, _, old_line = ir.new_nodes[new_node_id]
            if isa(typ, DelayedTyp)
                push!(type_refine_phi, ssaval)
            end
            new_typ = isa(typ, DelayedTyp) ? Union{} : tmerge(old_typ, typ)
            ir.new_nodes[new_node_id] = (old_insert, new_typ, node, old_line)
            incoming_vals[slot] = NewSSAValue(ssaval)
        end
        (item in visited) && continue
        push!(visited, item)
        for idx in cfg.blocks[item].stmts
            stmt = ci.code[idx]
            if isa(stmt, NewvarNode)
                incoming_vals[slot_id(stmt.slot)] = undef_token
                ci.code[idx] = nothing
            else
                stmt = rename_uses!(ir, ci, idx, stmt, incoming_vals)
                if stmt === nothing && idx == last(cfg.blocks[item].stmts)
                    # preserve the CFG
                    stmt = ReturnNode{Any}()
                end
                ci.code[idx] = stmt
                # Record a store
                if isexpr(stmt, :(=)) && isa(stmt.args[1], SlotNumber)
                    id = slot_id(stmt.args[1])
                    val = stmt.args[2]
                    typ = typ_for_val(val, ci)
                    incoming_vals[id] = SSAValue(make_ssa!(ci, idx, id, typ))
                end
            end
        end
        for succ in cfg.blocks[item].succs
            push!(worklist, (succ, item, copy(incoming_vals)))
        end
    end
    # Delete any instruction in unreachable blocks
    for bb in setdiff(IdSet{Int}(1:length(cfg.blocks)), visited)
        for idx in cfg.blocks[bb].stmts
            ci.code[idx] = nothing
        end
    end
    # Convert into IRCode form
    code = ir.stmts
    ssavalmap = Any[SSAValue(-1) for _ in 1:(length(ci.ssavaluetypes)+1)]
    types = Any[Any for _ in 1:length(code)]
    # Detect statement positions for assignments and construct array
    for (idx, stmt) in Iterators.enumerate(ci.code)
        if isexpr(stmt, :(=)) && isa(stmt.args[1], SSAValue)
            ssavalmap[stmt.args[1].id + 1] = SSAValue(idx)
            types[idx] = ci.ssavaluetypes[stmt.args[1].id + 1]
            stmt = stmt.args[2]
            if isa(stmt, PhiNode)
                edges = Any[block_for_inst(cfg, edge) for edge in stmt.edges]
                code[idx] = PhiNode(edges, stmt.values)
            else
                code[idx] = stmt
            end
        # Convert GotoNode/GotoIfNot/PhiNode to BB addressing
        elseif isa(stmt, GotoNode)
            code[idx] = GotoNode(block_for_inst(cfg, stmt.label))
        elseif isa(stmt, GotoIfNot)
            code[idx] = GotoIfNot(stmt.cond, block_for_inst(cfg, stmt.dest))
        else
            code[idx] = stmt
        end
    end
    # This is a bit awkward, because it basically duplicates what type
    # inference does. Ideally, we'd just use this representation earlier
    # to make sure phi nodes have accurate types
    changed = true
    while changed
        changed = false
        for phi in type_refine_phi
            new_idx = phi - length(ir.stmts)
            old_insert, old_typ, node, old_line = ir.new_nodes[new_idx]
            new_typ = Union{}
            for i = 1:length(node.values)
                if !isassigned(node.values, i)
                    if !isa(new_typ, MaybeUndef)
                        new_typ = MaybeUndef(new_typ)
                    end
                    continue
                end
                typ = typ_for_val(node.values[i], ci)
                if isa(typ, DelayedTyp)
                    typ = ir.new_nodes[typ.phi.id - length(ir.stmts)][2]
                end
                new_typ = tmerge(new_typ, typ)
            end
            if !(old_typ ⊑ new_typ) || !(new_typ ⊑ old_typ)
                ir.new_nodes[new_idx] = (old_insert, new_typ, node, old_line)
                changed = true
            end
        end
    end
    types = Any[isa(types[i], DelayedTyp) ? ir.new_nodes[types[i].phi.id - length(ir.stmts)][2] : types[i] for i in 1:length(types)]
    new_nodes = NewNode[let (pos, typ, node, line) = ir.new_nodes[i]
            typ = isa(typ, DelayedTyp) ? ir.new_nodes[typ.phi.id - length(ir.stmts)][2] : typ
            (pos, typ, node, line)
        end for i in 1:length(ir.new_nodes)]
    # Renumber SSA values
    code = Any[new_to_regular(renumber_ssa!(code[i], ssavalmap)) for i in 1:length(code)]
    new_nodes = NewNode[let (pt, typ, stmt, line) = new_nodes[i]
            (pt, typ, new_to_regular(renumber_ssa!(stmt, ssavalmap)), line)
        end for i in 1:length(new_nodes)]
    return IRCode(ir, code, types, ir.lines, ir.cfg, new_nodes)
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
using Test
using LegendrePolynomials
using OffsetArrays
using HyperDualNumbers
using QuadGK
using Symbolics
using LinearAlgebra

tohyper(x) = Hyper(x, one(x), one(x), zero(x))

@testset "LegendrePolynomialIterator" begin
    iter = LegendrePolynomials.LegendrePolynomialIterator(0)
    @test Base.IteratorSize(iter) === Base.IsInfinite()
    @test eltype(iter) == Float64
    iter = LegendrePolynomials.LegendrePolynomialIterator(0, 1)
    @test Base.IteratorSize(iter) === Base.IsInfinite()
    iter = LegendrePolynomials.LegendrePolynomialIterator(big(0))
    @test eltype(iter) == BigFloat
end

@testset "Pl and collectPl" begin
	x = 0.5
    lmax = 5
    lmin = max(0, lmax - 1)

    @testset "norm" begin
        @testset "standard norm" begin
            P = @inferred collectPl(x, lmax = lmax)
            T = Union{OffsetVector{Float64, Vector{Float64}}, OffsetVector{BigFloat, Vector{BigFloat}}}
            P3 = @inferred collectdnPl(x, lmax = lmax, n = 0)
            @test P[0] ≈ P3[0] ≈ 1
            @test P[1] ≈ P3[1] ≈ x
            @test P[2] ≈ P3[2] ≈ (3x^2 - 1)/2
            @test P[3] ≈ P3[3] ≈ (5x^3 - 3x)/2
            @test P[4] ≈ P3[4] ≈ (35x^4 - 30x^2 + 3)/8
            @test P[5] ≈ P3[5] ≈ (63x^5 - 70x^3 + 15x)/8
        end
    end

    @testset "x = 0" begin
        for l = 1:2:101
            P1 = @inferred Pl(0, l)
            P2 = @inferred Union{BigFloat, Float64} Plm(0, l, 0)
            P3 = @inferred dnPl(0, l, 0)
            @test P1 == P2 == P3 == 0
            P1 = @inferred Pl(0, l, norm = Val(:normalized))
            P2 = @inferred Plm(0, l, 0, norm = Val(:normalized))
            @test P1 == P2 == 0
        end
    end

    @testset "x = ±1" begin
        lmax = 10

        P = collectPl(1, lmax = lmax)
        @test all(==(1), P)
        P2 = collectPlm(1, lmax = lmax, m = 0)
        @test all(≈(1), P2)
        P3 = collectdnPl(1, lmax = lmax, n = 0)
        @test all(==(1), P3)

        P = collectPl(-1, lmax = lmax)
        P2 = collectPlm(-1, lmax = lmax, m = 0)
        P3 = collectdnPl(-1, lmax = lmax, n = 0)
        for l in axes(P, 1)
            @test P[l] ≈ P2[l] ≈ P3[l] ≈ (-1)^l
        end
    end

    @testset "large degree" begin
        l = 200
        x = 0.5
        expval = -big"50516977220606216108218169174406199263421124121643988159364881659548
                    89526732539052701378087401221240954262583383510627"/
                    big"3227812347608635737069898965003764842912132241036529391038324195675
                    80952752105149328705669160017228929487896496593436672"

        @test Pl(x, l) ≈ expval rtol=1e-14
        @test Pl(big(x), big(l)) ≈ expval rtol=1e-75
    end

    @testset "parity" begin
        for x in range(0, 1, length=10), l in 0:4
            @test Pl(-x, l) ≈ (-1)^l * Pl(x, l) atol=1e-12 rtol=1e-8
        end
    end

    @test_throws DomainError collectPl(-2, lmax = lmax)
    @test_throws ArgumentError collectPl(0, lmax = -1)

    @test_throws ArgumentError collectPl(0, lmax=1, norm = Val(:undefined))
end

@testset "Plm and collectPlm" begin
    x = 0.5
    lmax = 5
    lmin = max(0, lmax - 1)

    @testset "double factorial" begin
        @test Plm(0.5, 2, 2) ≈ 9/4
        @test Plm(0.5, 3, 3) ≈ -45/8*sqrt(3)
        @test Plm(0.5, 3, 3) ≈ -45/8*sqrt(3)
        @test Plm(0.5, 10, 0) ≈ -49343/262144
        @test Plm(0.5, 10, 6) ≈ -674999325/8192
        @test Plm(0.5, 10, 10) ≈ 159099165225/1024
    end
    @testset "norm" begin
        @testset "standard norm" begin
            P = @inferred collectPl(x, lmax = lmax)
            T = Union{OffsetVector{Float64, Vector{Float64}}, OffsetVector{BigFloat, Vector{BigFloat}}}
            P2 = @inferred T collectPlm(x, lmax = lmax, m = 0)
            @test P ≈ P2

            Q = @inferred collectPl(x, lmax = lmax, lmin = lmin, norm = Val(:standard))
            Q2 = @inferred T collectPlm(x, lmax = lmax, lmin = lmin, m = 0, norm = Val(:standard))
            @test Q ≈ Q2
            @test @views Q[lmin:lmax] == P[lmin:lmax]
        end

        @testset "normalized" begin
            P = @inferred collectPl(x, lmax = lmax, norm = Val(:normalized))
            Q = @inferred collectPlm(x, lmax = lmax, m = 0, norm = Val(:normalized))
            @test P ≈ Q
            @test P[0] ≈ √(1/2)

            P2 = @inferred collectPl(x, lmax = lmax, lmin = lmin, norm = Val(:normalized))
            Q2 = @inferred collectPlm(x, lmax = lmax, lmin = lmin, m = 0, norm = Val(:normalized))
            @test P2 ≈ Q2
            @test @views P2[lmin:lmax] == P[lmin:lmax]

            for m in 1:3
                Q = @inferred collectPlm(x, lmax = lmax, m = m, norm = Val(:normalized))
                @testset for l in m:lmax
                    @test Q[l] ≈ √((2l+1)/2*(factorial(l-m)/factorial(l+m))) * Plm(x, l, m)
                end
            end
        end

        @testset "schmidtquasi" begin
            P = @inferred collectPl(x, lmax = lmax, norm = Val(:schmidtquasi))
            Q = @inferred collectPlm(x, lmax = lmax, m = 0, norm = Val(:schmidtquasi))
            @test P ≈ Q

            for m in 1:3
                Q = @inferred collectPlm(x, lmax = lmax, m = m, norm = Val(:schmidtquasi))
                @testset for l in m:lmax
                    @test Q[l] ≈ √(2*(factorial(l-m)/factorial(l+m))) * Plm(x, l, m)
                end
            end
        end

        @testset "schmidt" begin
            P = @inferred collectPl(x, lmax = lmax, norm = Val(:schmidt))
            Q = @inferred collectPlm(x, lmax = lmax, m = 0, norm = Val(:schmidt))
            @test P ≈ Q

            for m in 1:3
                Q = @inferred collectPlm(x, lmax = lmax, m = m, norm = Val(:schmidt))
                @testset for l in m:lmax
                    @test Q[l] ≈ √(2*(2l+1)*(factorial(l-m)/factorial(l+m))) * Plm(x, l, m)
                end
            end
        end
    end
    @testset "non-standard integers" begin
        l = 2
        m = 1
        x = 0.5
        p_Int = Plm(x, l, m)
        p_Int8 = Plm(x, Int8(l), Int8(m))
        @test p_Int == p_Int8
        p_BigInt = Plm(big(x), big(l), big(m))
        @test p_BigInt ≈ p_Int rtol=1e-12
    end
    @testset "Condon-Shortley phase" begin
        p1 = Plm(0.5, 2, 1)
        p2 = Plm(0.5, 2, 1, csphase = false)
        @test p1 ≈ -p2
        p1 = Plm(0.5, 2, 2)
        p2 = Plm(0.5, 2, 2, csphase = false)
        @test p1 ≈ p2
        for norm in [:standard, :normalized, :schmidtquasi]
            @test all(Plm(0.5, 3, m, csphase = false, norm = Val(norm)) ≈
                    (-1)^m * Plm(0.5, 3, m, norm = Val(norm)) for m in -3:3)
        end
    end
    @testset "avoid overflow" begin
        l = 200
        m = 200
        x = 0.5
        # copied from Mathematica
        expval = big"2604065338240773053120460757061528950680509925586185542679571264679180
           1169571890502639513862292536574563617492830604806347049485953379217258
           5643579361782965615917473831726956215382240884618736130956683834034044
           2039232933380686854433175974549155182434563952079326518988454706912117
           0695345051735139971477695522848840659242654598631485913523486512584654
           3918383779732853064471286369286966147111339555668732856505146906124873
           14193734385778184852110786895451610689633525907993316650390625"/
           big"1606938044258990275541962092341162602522202993782792835301376"
        @test Plm(x, l, m) ≈ expval rtol=1e-12
        @test Plm(big(x), big(l), m) ≈ expval rtol=1e-70
    end
    @testset "negative order" begin
        x = 0.5
        pm1 = collectPlm(x, lmax = 2, m = -1)
        @test Plm(x, 1, -1) ≈ pm1[1] ≈ √3/4
        @test Plm(x, 2, -1) ≈ pm1[2] ≈ √3/8
        pm2 = collectPlm(x, lmax = 3, m = -2)
        @test Plm(x, 2, -2) ≈ pm2[2] ≈ 3/32
        @test Plm(x, 3, -2) ≈ pm2[3] ≈ 3/64
    end
    @testset "parity" begin
        for x in range(0, 1, length=10), l in 0:4, m in -l:l
            @test Plm(-x, l, m) ≈ (-1)^(l+m) * Plm(x, l, m) atol=1e-12 rtol=1e-8
        end
    end
    @test_throws ArgumentError Plm(0.5, 3, 10)
    @test_throws ArgumentError Plm(0.5, -1, 3)
    @test_throws ArgumentError Plm(0.5, 1, -3)
    @test_throws DomainError Plm(2.5, 1, 3)

    @test_throws ArgumentError collectPlm(0, lmax=1, m=1, norm = Val(:undefined))
end

@testset "dnPl" begin
    @test dnPl(0.5, 3, 10) == 0
    @test dnPl(0.5, 3, 10, zeros(100)) == 0
    @testset "double factorial" begin
        @test dnPl(0.5, 2, 2) == 3
        @test dnPl(0.5, 3, 3) == 15
        @test dnPl(0.5, 3, 3, zeros(1)) == 15
    end
    @test_throws ArgumentError dnPl(0.5, 10, 3, zeros(1))
    @test_throws ArgumentError dnPl(0.5, -1, 3)
    @test_throws ArgumentError dnPl(0.5, 1, -3)
    @test_throws DomainError dnPl(2.5, 1, 3)
end

@testset "collectPl!" begin
    v = zeros(10);
    x = 0.6

    collectPl!(v, x)
    for l in 0:length(v) - 1
        @test v[firstindex(v) + l] == Pl(x, l)
    end

    v .= 0
    lmax = 6
    collectPl!(v, x, lmax = lmax)
    for l in 0:lmax
        @test v[firstindex(v) + l] == Pl(x, l)
    end

    @test_throws DomainError collectPl!(v, -2, lmax = lmax)
    @test_throws ArgumentError collectPl!(v, 0, lmax = -1)
end

@testset "collectPlm!" begin
    v = zeros(10);
    @test_throws DomainError collectPlm!(v, -2, lmax = 3, m = 2)
    @test_throws ArgumentError collectPlm!(v, 0, lmax = -1, m = 3)
    @test_throws ArgumentError collectPlm!(v, 0.5, lmax = 1, m = -3)
end

@testset "collectPlm" begin
    @test_throws DomainError collectPlm(-2, lmax = 3, m = 2)
    @test_throws ArgumentError collectPlm(0, lmax = -1, m = 3)
    @test_throws ArgumentError collectPlm(0.5, lmax = 1, m = -3)
end

@testset "collectdnPl!" begin
    v = zeros(10);
    @test_throws DomainError collectdnPl!(v, -2, lmax = 3, n = 2)
    @test_throws ArgumentError collectdnPl!(v, 0, lmax = -1, n = 3)
    @test_throws ArgumentError collectdnPl!(v, 0.5, lmax = 1, n = -3)
    @test_throws ArgumentError collectdnPl!(zeros(1), 0, lmax = 3, n = 3)
end

@testset "collectdnPl" begin
    @test_throws DomainError collectdnPl(-2, lmax = 3, n = 2)
    @test_throws ArgumentError collectdnPl(0, lmax = -1, n = 3)
    @test_throws ArgumentError collectdnPl(0.5, lmax = 1, n = -3)
end

@testset "dPl/dx" begin
    x = 2rand() - 1
    xh = tohyper(x)
    lmax = 5
    P = eps1.(collectPl(xh, lmax = lmax))
    P2 = collectPlm(x, lmax = lmax, m = 1) / (- sqrt(1 - x^2))
    P3 = collectdnPl(x, lmax = lmax, n = 1)
    @test P[0] == P3[0] == 0
    @test P[1] ≈ P2[1] ≈ P3[1] ≈ 1
    @test P[2] ≈ P2[2] ≈ P3[2] ≈ 3x
    @test P[3] ≈ P2[3] ≈ P3[3] ≈ (-3 + 15x^2)/2
    @test P[4] ≈ P2[4] ≈ P3[4] ≈ 1/8*(-60x + 140x^3)
    @test P[5] ≈ P2[5] ≈ P3[5] ≈ 1/8*(15 - 210x^2 + 315x^4)

    @testset "x = ±1" begin
        lmax = 10
        x = 1
        P = collectPl(tohyper(x), lmax = lmax)
        dP1h = eps1.(P)
        dP1 = collectdnPl(x, lmax = lmax, n = 1)

        for l in 1:lmax
            @test dP1h[l] == dP1[l] == l*(l+1)/2
        end

        x = -1
        P = collectPl(tohyper(x), lmax = lmax)
        dPm1h = eps1.(P)
        dPm1 = collectdnPl(x, lmax = lmax, n = 1)

        for l in 1:lmax
            @test dPm1h[l] == dPm1[l] == (-1)^(l+1) * l*(l+1)/2
        end
    end
end

@testset "d^2Pl/dx^2" begin
    x = 2rand() - 1
    xh = tohyper(x)
    lmax = 5
    P = eps1eps2.(collectPl(xh, lmax = lmax))
    P2 = collectPlm(x, lmax = lmax, m = 2) / (1 - x^2)
    P3 = collectdnPl(x, lmax = lmax, n = 2)

    @test P[0] == P3[0] == 0
    @test P[1] == P3[1] == 0
    @test P[2] ≈ P2[2] ≈ P3[2] ≈ 3
    @test P[3] ≈ P2[3] ≈ P3[3] ≈ 15x
    @test P[4] ≈ P2[4] ≈ P3[4] ≈ 1/8*(-60 + 420x^2)
    @test P[5] ≈ P2[5] ≈ P3[5] ≈ 1/8*(-420x + 1260x^3)
end

@testset "Askey-Gasper" begin
    x = 0.5
    for lmax=1:100:1_000
        @test sum(collectPl(x, lmax = lmax)) >= 0
        @test sum(collectPlm(x, lmax = lmax, m = 0)) >= 0
        @test sum(collectdnPl(x, lmax = lmax, n = 0)) >= 0
    end
end

@testset "norm" begin
    @testset "Pl" begin
        @testset "normalized" begin
            for l in 0:4
                f = x -> Pl(x, l, norm = Val(:normalized))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 1 rtol=1e-2 atol=E
            end
        end

        @testset "standard" begin
            for l in 0:4
                f = x -> Pl(x, l, norm = Val(:standard))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 2/(2l+1) rtol=1e-2 atol=E
            end
        end

        @testset "schmidtquasi" begin
            for l in 0:4
                f = x -> Pl(x, l, norm = Val(:schmidtquasi))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 2/(2l+1) rtol=1e-2 atol=E
            end
        end

        @testset "schmidt" begin
            for l in 0:4
                f = x -> Pl(x, l, norm = Val(:schmidt))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 2 rtol=1e-2 atol=E
            end
        end
    end
    @testset "Plm" begin
        @testset "normalized" begin
            @testset for m in -4:4, l in abs(m):4
                f = x -> Plm(x, l, m, norm = Val(:normalized))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 1 rtol=1e-2 atol=E
            end
        end

        @testset "schmidtquasi" begin
            @testset for m in -4:4, l in abs(m):4
                f = x -> Plm(x, l, m, norm = Val(:schmidtquasi))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 2*(2 - (m==0))/(2l + 1) rtol=1e-2 atol=E
            end
        end

        @testset "schmidt" begin
            @testset for m in -4:4, l in abs(m):4
                f = x -> Plm(x, l, m, norm = Val(:schmidt))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ 2*(2 - (m==0)) rtol=1e-2 atol=E
            end
        end

        @testset "standard" begin
            for m in -4:4, l in abs(m):4
                f = x -> Plm(x, l, m, norm = Val(:standard))^2
                I, E = quadgk(f, -1, 1, rtol=1e-2)
                @test I ≈ (2/(2l+1)*factorial(l+m)/factorial(l-m)) rtol=1e-2 atol=E
            end
        end
    end
end

@testset "custom types" begin
    @testset "Symbolics" begin
        @variables x
        f_expr = build_function(Pl(x, 3), [x])
        myf = eval(f_expr)
        @test myf(0.4) ≈ Pl(0.4, 3)
        f_expr = build_function(Plm(x, 3, 2), [x])
        myf2 = eval(f_expr)
        @test myf2(0.4) ≈ Plm(0.4, 3, 2)
    end

    @testset "Matrix" begin
        A = reshape([0.1:0.1:0.4;], 2, 2)
        @test Pl(A, 1) == A
        @test Pl(A, 2) ≈ (3A^2 - one(A))/2
        @test Plm(A, 2, 1) ≈ -3A*√(I - A^2)

        D = Diagonal([1.0, 0.5])
        @test Pl(D, 1) == D
        @test Pl(D, 2) ≈ (3D^2 - one(D))/2
        @test Plm(D, 2, 1) ≈ -3D*√(I - D^2)
    end
end

ext: jl


text:
module FrameFunInterface

using StaticArrays

using ..Platforms, ..ApproximationProblems, BasisFunctions
using ..ParameterPaths: ParametrizedPlatform
import BasisFunctions: discretemeasure, measure, interpolation_grid, elements

export SamplingStrategy
abstract type SamplingStrategy end

export DefaultSamplingStrategy
"""
The default strategy tries to keep the oversamplingfactor correct by distributing equally over all degrees of freedom.

"""
struct DefaultSamplingStrategy <: SamplingStrategy end

default_aznormalization(a...) = false

export smoothingoperator
smoothingoperator(ap::ApproximationProblem; options...) =
    WeightedSmoothingOperator(dictionary(ap); options...)

export normalizationoperator
normalizationoperator(::DictionaryOperatorStyle, ap::ApproximationProblem; sstyle=SamplingStyle(ap), opts...) =
    normalizationoperator(sstyle, ap; opts...)

function normalizationoperator(sstyle::DiscreteStyle, ap::ApproximationProblem; T=coefficienttype(ap), options...)
    grid = sampling_grid(sstyle, ap; options...)
    sampling_normalization(GridBasis{T}(grid),
        measure(sstyle, ap; options...))
end

# Utility functions to create tensor FrameFunInterface functions

function elements(fun, ss::ProductSamplingStyle, ap::ProductPlatformApproximation, args...; options...)
    samplingparameter(ss, ap; options...)
    @assert length(elements(ss))==length(elements(ap))
    map((ssi, api, argi)->fun(ssi, api, argi...; options...), elements(ss), elements(ap), zip(args...))
end

function elements(fun, ss::ProductSamplingStyle, ap::ProductPlatformApproximation; options...)
    samplingparameter(ss, ap; options...)
    @assert length(elements(ss))==length(elements(ap))
    map((ssi, api)->fun(ssi, api; options...), elements(ss), elements(ap))
end

INTERFACEFUNCTIONS = (:samplingoperator, :dualsamplingoperator, :samplingparameter,
    :sampling_grid, :interpolation_grid, :platform_grid, :discretization, :dualdiscretization,
    :azdual_dict, :measure, :discretemeasure, :solver, :AZ_A, :AZ_Z, :AZ_Zt, :plungeoperator,
    :plungematrix, :plungerank, :smoothingoperator)

# parameter                                     # rename to
# platform: the platform to approximate with    #
# dict: the dict to approximate with            #
# M: grid size of problem                       # M
# L: grid size of bounding box                  # sparam
# param: index of platform                      # pparam
# N: (length(dict))

"""
    @aptoplatform `function`

Implement the function at the level of the dictionary.
"""
macro aptoplatform(ex)
    def = Meta.parse("default_"*string(ex))
    intermediate = Meta.parse("_"*string(ex))
    ret = quote
        $(ex)(ap::ApproximationProblem, args...; samplingstyle=SamplingStyle(ap),options...) =
            $(intermediate)(samplingstyle, ap, samplingparameter(samplingstyle, ap; options...), args...; options...)
        $(ex)(ss::SamplingStyle, ap::ApproximationProblem, args...; options...) =
            $(intermediate)(ss, ap, samplingparameter(ss, ap; options...), args...; options...)
        $(intermediate)(ss::SamplingStyle, ap::ApproximationProblem, L, args...; options...) =
            $(ex)(ss, platform(ap), parameter(ap), L, args...; options...)
        $(ex)(ss::SamplingStyle, platform::Platform, param, L, args...; dict = dictionary(platform, param), options...) =
            $(ex)(ss, dict, L, args...)
        $(ex)(ss::SamplingStyle, platform::ParametrizedPlatform, param, L, args...; options...) =
            $(ex)(ss, platform.platform, platform.path[param], L, args...; options...)
        $(ex)(ss::SamplingStyle, dict::Dictionary, L, args...; options...) =
            $(def)(ss, dict, L, args...)
    end
    esc(ret)
end

macro trial(ex)
    def = Meta.parse("default_"*string(ex))
    ret = quote
        $(ex)(dict::Dictionary, args...; options...) = $(ex)(approximationproblem(dict, args...); options...)
        $(ex)(platform::Platform, args...; options...) = $(ex)(approximationproblem(platform, args...); options...)

        function $(ex)(ap::ApproximationProblem, args...; samplingstyle=SamplingStyle(ap),options...)
            samplingparameter(samplingstyle, ap; options...)
            $(ex)(samplingstyle, ap, args...; options...)
        end

        $(ex)(ss::SamplingStyle, ap::ApproximationProblem, args...; options...) =
            $(ex)(ss, platform(ap), parameter(ap), args...; options...)
        $(ex)(ss::SamplingStyle, platform::Platform, param, args...; dict = dictionary(platform, param), options...) =
            $(ex)(ss, dict, args...; options...)
        $(ex)(ss::SamplingStyle, platform::ParametrizedPlatform, param, args...; options...) =
            $(ex)(ss, platform.platform, platform.path[param], args...; options...)
        $(ex)(ss::SamplingStyle, dict::Dictionary, args...; options...) =
            $(def)(ss, dict, args...; options...)
    end
    esc(ret)
end

macro platformtoap(ex)
    ret = quote
        $(ex)(dict::Dictionary, args...; options...) = $(ex)(approximationproblem(dict, args...); options...)
        $(ex)(platform::Platform, args...; options...) = $(ex)(approximationproblem(platform, args...); options...)

        function $(ex)(ap::ApproximationProblem, args...; samplingstyle=SamplingStyle(ap),options...)
            samplingparameter(samplingstyle, ap; options...)
            $(ex)(samplingstyle, ap, args...; options...)
        end
    end
    esc(ret)
end

macro addsamplinparameter(ex)
    intermediate = Meta.parse("_"*string(ex))
    ret = quote
        $(ex)(ss::SamplingStyle, ap::ApproximationProblem, args...; options...) =
            $(intermediate)(ss, ap, samplingparameter(ss, ap; options...), args...; options...)
        $(intermediate)(ss::SamplingStyle, ap::ApproximationProblem, L, args...; options...) =
            $(ex)(ss, platform(ap), parameter(ap), L, args...; options...)
    end
    esc(ret)
end

macro addsamplingstyle(ex)
    ret = quote
        $(ex)(ap::ApproximationProblem; samplingstyle = SamplingStyle(ap), options...) =
            $(ex)(samplingstyle, ap; options...)
    end
    esc(ret)
end

macro addproblemstyle(ex)
    ret = quote
        $(ex)(ap::ApproximationProblem; problemstyle = ProblemStyle(ap), options...) =
            $(ex)(problemstyle, ap; options...)
    end
    esc(ret)
end

# We delegate some functionality from the platform to the dictionary corresponding
# to a value of the platform parameter param.
# Platforms may choose to override these to enable different behaviour.
export interpolation_grid
@trial interpolation_grid
function default_interpolation_grid(ss::SamplingStyle, dict::Dictionary; oversamplingfactor=1, options...)
    oversamplingfactor != 1 && @warn "InterpolationStyle does not support option `oversamplingfactor=$oversamplingfactor`, use `OversamplingStyle()` instead."
    interpolation_grid(dict)
end

export oversampling_grid
@trial oversampling_grid
oversampling_grid(ss::SamplingStyle, ap::ApproximationProblem; options...) =
    oversampling_grid(ss, platform(ap), parameter(ap), samplingparameter(ss, ap; options...); options...)
oversampling_grid(samplingstyle::SamplingStyle, platform::Platform, param, L; dict = dictionary(platform, param), options...) =
    oversampling_grid(dict, L)
include("oversampling_grid.jl")

export samplingparameter
@platformtoap samplingparameter
@addsamplingstyle samplingparameter
include("samplingparameter.jl")

export samplingoperator, dualsamplingoperator
@platformtoap samplingoperator
@platformtoap dualsamplingoperator
include("samplingoperator.jl")

export sampling_grid, platform_grid
@platformtoap sampling_grid
@trial platform_grid
include("grids.jl")
default_sampling_grid(ss::SamplingStyle, dict::Dictionary, L; options...) = oversampling_grid(dict, L)
default_platform_grid(ss::SamplingStyle, dict::Dictionary, args...; grid, options...) = grid

export discretemeasure
@trial discretemeasure
discretemeasure(ss::ProductSamplingStyle, ap::ApproximationProblem; options...) =
    productmeasure(elements(discretemeasure, ss, ap; options...)...)
discretemeasure(ss::SamplingStyle, ap::ApproximationProblem; options...) =
    discretemeasure(ss, platform(ap), parameter(ap), ap; options...)
discretemeasure(ss::SamplingStyle, platform::Platform, param, ap; options...) =
    discretemeasure(sampling_grid(ss, ap; options...))

export measure
@trial measure
default_measure(::SamplingStyle, dict::Dictionary; options...) =
    measure(dict)
measure(ss::DiscreteGramStyle, ap::ApproximationProblem; options...) =
    discrete_gram_measure(ss, ap; options...)
@trial discrete_gram_measure
discrete_gram_measure(ss::SamplingStyleSuper{DiscreteGramStyle}, ap::ApproximationProblem; options...) =
    discretemeasure(ss, ap; options...)
discrete_gram_measure(ss::SamplingStyle, ap::ApproximationProblem; options...) =
    error("Only with `DiscreteGramStyle`")

export azdual_dict
@trial azdual_dict
include("azdual_dict.jl")

export discretization, dualdiscretization
@trial discretization
@trial dualdiscretization
include("discretization.jl")

export normalized_discretization, normalized_dualdiscretization
@trial normalized_discretization
@trial normalized_dualdiscretization
include("normalized_discretization.jl")

export solver
@platformtoap solver
include("solver.jl")

## The AZ algorithm
export AZ_A
@platformtoap AZ_A
@addproblemstyle AZ_A
AZ_A(pstyle::ProblemStyle, ap::ApproximationProblem; options...) =
    normalized_discretization(ap; options...)

AZ_A(::GenericOperatorStyle, ap::ApproximationProblem; options...) =
    SynthesisOperator(dictionary(ap), measure(ap; options...))

export AZ_Z
@platformtoap AZ_Z
@addproblemstyle AZ_Z
AZ_Z(pstyle::ProblemStyle, ap::ApproximationProblem; options...) =
    normalized_dualdiscretization(ap; options...)
AZ_Z(::GenericOperatorStyle, ap::ApproximationProblem; options...) =
    SynthesisOperator(azdual_dict(ap; options...), measure(ap; options...))

export AZ_Zt
@platformtoap AZ_Zt
@addproblemstyle AZ_Zt
AZ_Zt(pstyle::ProblemStyle, ap::ApproximationProblem; options...) = AZ_Z(pstyle, ap; options...)'
samplingoperator(pstyle::GenericOperatorStyle, ap; options...) = AZ_Zt(pstyle, ap; options...)
samplingoperator(pstyle::DictionaryOperatorStyle, ap; options...) = samplingoperator(ap; options...)

export plungeoperator
@platformtoap plungeoperator
@addproblemstyle plungeoperator
function plungeoperator(problemstyle::ProblemStyle, ap::ApproximationProblem; options...)
    A = AZ_A(problemstyle, ap; options...)
    Zt = AZ_Zt(problemstyle, ap; options...)
    I - A*Zt
end

export plungematrix, firstAZstepoperator
@platformtoap plungematrix
@addproblemstyle plungematrix
function plungematrix(pstyle::ProblemStyle, ap::ApproximationProblem; options...)
    A = AZ_A(pstyle, ap; options...)
    P = plungeoperator(pstyle, ap; options...)
    P * A
end
function plungematrix(pstyle::GenericOperatorStyle, ap::ApproximationProblem; options...)
    S = samplingoperator(pstyle, ap; options...)
    A = AZ_A(pstyle, ap; options...)
    Zt = AZ_Zt(pstyle, ap; options...)
    tmp1 = apply(Zt, A; options...)
    tmp2 = apply(S, A; options...)
    tmp2 - tmp2*tmp1
end
const firstAZstepoperator = plungematrix

export plungerank
@platformtoap plungerank
function plungerank(ap::ApproximationProblem;
            rankestimate = 40,
            options...)
    C = plungematrix(ap; options...)
    Q = rSVD_solver(C; rankestimate = rankestimate, options...)
    length(Q.Sinv)
end

import BasisFunctions: approximate
export Fun
include("approximate.jl")

end


text:
#######
# Unscented Kalman filter types and filtering interface

abstract UnscentedKalmanFilter <: KalmanFilter
abstract AbstractUnscentedState <: AbstractState

type UnscentedState{T} <: AbstractUnscentedState
    x::Vector{T}
    p::Matrix
    α::Real
    β::Real
    κ::Real
end

type AdditiveUnscentedObservationModel <: ObservationModel
    h::Function
    r::Matrix
end

type AdditiveUnscentedModel <: Model
    f::Function
    q::Matrix
end

## ap function for AdditiveUnscentedModels and UnscentedStates

function ap(f::AdditiveUnscentedModel,s::UnscentedState)
    σn,wm,wc = sigma(s)
    σs = zeros(σn)
    xn = zeros(s.x)
    for i in 1:size(σn,2)
        σs[:,i] = f.f(σn[:,i])
        xn += wm[i] * σs[:,i]
    end
    pn = zeros(s.p)
    for i in 1:size(σn,2)
        pn += wc[i] * (σs[:,i]-xn)*(σs[:,i]-xn)'
    end
    pn += f.q
    UnscentedState(xn,pn,s.α,s.β,s.κ)
end

type AdditiveUnscentedKalmanFilter <: UnscentedKalmanFilter
    x::UnscentedState
    f::AdditiveUnscentedModel
    z::AdditiveUnscentedObservationModel
end

function AdditiveUnscentedKalmanFilter(x::Vector,p::Matrix,f::AdditiveUnscentedModel,z::AdditiveUnscentedObservationModel,α::Real,β::Real,κ::Real)
    s = UnscentedState(x,p,α,β,κ)
    AdditiveUnscentedKalmanFilter(s,f,z)
end

### covs function for additive unscented Kalman filters

function covs(kf::AdditiveUnscentedKalmanFilter,y::Observation)
    n = size(kf.x.x,1)
    m = size(kf.z.r,1)
    σn,wm,wc = sigma(kf.x.x,kf.x.p,kf.x.α,kf.x.β,kf.x.κ)
    k = size(σn,2)
    yp = zeros(m,k)
    yhat = zeros(m)
    for i in 1:k
        yp[:,i] = kf.z.h(σn[:,i])
        yhat += wm[i] * yp[:,i]
    end

    resx = zeros(σn)
    resy = zeros(yp)

    for i in 1:k
        resx[:,i] = σn[:,i]-kf.x.x
        resy[:,i] = yp[:,i]-yhat
    end

    res = y.y-yhat

    ph = zeros(n,m)
    s = zeros(m,m)
    for i in 1:k
        ph += wc[i]*(resx[:,i]*resy[:,i]')
        s += wc[i]*(resy[:,i]*resy[:,i]')
    end

    s += kf.z.r
    return res,ph,s
end
#######
# Sigma point generating functions

function sigma(x::Vector,p::Matrix,α,β,κ)
    n = length(x)
    σ = zeros(n,2n+1)
    wm = zeros(2n+1)
    wc = zeros(2n+1)
    λ = α^2*(n+κ)-n
    σ[:,1] = x
    wm[1] = λ/(n+λ)
    wc[1] = wm[1] + (1-α^2+β)
    wm[2:end] = 1/(2*(n+λ))
    wc[2:end] = wm[2:end]
    sp = sqrt(n+λ)*full(chol(p))
    for i = 2:n+1
        σ[:,i] = x + sp[:,i-1]
        σ[:,i+n] = x - sp[:,i-1]
    end
    (σ,wm,wc)
end

function sigma(s::UnscentedState)
    n = length(s.x)
    σ = zeros(n,2n+1)
    wm = zeros(2n+1)
    wc = zeros(2n+1)
    λ = s.α^2*(n+s.κ)-n
    σ[:,1] = s.x
    wm[1] = λ/(n+λ)
    wc[1] = wm[1] + (1-s.α^2+s.β)
    wm[2:end] = 1/(2*(n+λ))
    wc[2:end] = wm[2:end]
    sp = sqrt(n+λ)*full(chol(s.p))
    for i = 2:n+1
        σ[:,i] = s.x + sp[:,i-1]
        σ[:,i+n] = s.x - sp[:,i-1]
    end
    (σ,wm,wc)
end


text:
module RSCG
    using LinearAlgebra
    using LinearMaps
    using SparseArrays

    export greensfunctions

    @doc raw"""
        greensfunctions(i::Integer,j::Integer,σ::Array{ComplexF64,1},A)

    Calculate the matrix element of the Green's function:
    ```math
    G_{ij}(\sigma_k) = \left[ (\sigma_k I - A)^{-1} \right]_{ij},
    ```
    with the use of the reduced-shifted conjugate gradient method
    (See, Y. Nagai, Y. Shinohara, Y. Futamura, and T. Sakurai,[arXiv:1607.03992v2 or DOI:10.7566/JPSJ.86.014708]).
    One can obtain ``G_{ij}(\sigma_k)`` with different frequencies ``\sigma_k``, simultaneously.

    Inputs:

    * `i` :index of the Green's function

    * `j` :index of the Green's function

    * `σ` :frequencies

    * `A` :hermitian matrix. We can use Arrays,LinearMaps, SparseArrays

    * `eps` :residual (optional) Default:`1e-12`

    * `maximumsteps` : maximum number of steps (optional) Default:`20000`

    Output:
    * `Gij[1:M]`: the matrix element Green's functions at M frequencies defined by ``\sigma_k``.
    """
    function greensfunctions(i::Integer,j::Integer,σ,A;eps=1e-12,maximumsteps=20000)

        M=length(σ)
        N=size(A,1)
        atype = eltype(A)
        σtype = eltype(σ)
        b = zeros(atype,N)
        b[j] = 1.0

    #--Line 2 in Table III.
        x = zeros(atype,N)
        r = copy(b)
        p = copy(b)
        αm = 1.0
        βm = 0.0

        m = 1
    #--
        Σ =zeros(atype,m)
        Σ[1] = b[i] #Line 3

        Θ = zeros(σtype,M)
        Π = ones(σtype,M,m)
        for mm=1:m
            Π[:,mm] *= Σ[mm]
        end
        ρk = ones(σtype,M)
        ρkm = copy(ρk)
        ρkp = copy(ρk)
        Ap = similar(p)

        for k=0:maximumsteps
            mul!(Ap,A,-p)
            #A_mul_B!(Ap,A,-p)
            rnorm = dot(r,r)
            α = rnorm/dot(p,Ap)
            @. x += α*p #Line 7
            @. r += -α*Ap #Line 8
            β = r'*r/rnorm #Line9
            @. p = r + β*p #Line 10

            Σ[1] = r[i] #Line 11

            for j = 1:M
                if abs(ρk[j]) > eps
                    ρkp[j] = ρk[j]*ρkm[j]*αm/(ρkm[j]*αm*(1.0+α*σ[j])+α*βm*(ρkm[j]-ρk[j]))#Line 13
                    αkj = α*ρkp[j]/ρk[j]#Line 14
                    Θ[j] += αkj*Π[j,1] #Line 15
                    βkj = ((ρkp[j]/ρk[j])^2)*β #Line 16
                    Π[j,:] = ρkp[j]*Σ+ βkj*Π[j,:] #Line 17
                    ρkm[j] = ρk[j]
                    ρk[j] = ρkp[j]
                end

            end
            αm = α
            βm = β
            hi = real(rnorm)*maximum(abs.(ρk))
            if hi < eps
                return Θ
            end
        end

        println("After $maximumsteps steps, this was not converged")
        return Θ

    end

    @doc raw"""
        greensfunctions(vec_left::Array{<:Integer,1},j::Integer,σ::Array{ComplexF64,1},A)

    Calculate the matrix element of the Green's function:
    ```math
    G_{ij}(\sigma_k) = \left[ (\sigma_k I - A)^{-1} \right]_{ij},
    ```
    with the use of the reduced-shifted conjugate gradient method
    (See, Y. Nagai, Y. Shinohara, Y. Futamura, and T. Sakurai,[arXiv:1607.03992v2 or DOI:10.7566/JPSJ.86.014708]).
    One can obtain ``G_{ij}(\sigma_k)`` with different frequencies ``\sigma_k``, simultaneously.

    Inputs:

    * `vec_left` :i indices of the Green's function

    * `j` :index of the Green's function

    * `σ` :frequencies

    * `A` :hermitian matrix. We can use Arrays,LinearMaps, SparseArrays

    * `eps` :residual (optional) Default:`1e-12`

    * `maximumsteps` : maximum number of steps (optional) Default:`20000`

    Output:
    * `Gij[1:M,1:length(vec_left)]`: the matrix element Green's functions at M frequencies defined by ``\sigma_k``.
    """
    function greensfunctions(vec_left::Array{<:Integer,1},j::Integer,σ,A;eps=1e-12,maximumsteps=20000)

        M=length(σ)
        N=size(A,1)
        atype = eltype(A)
        σtype = eltype(σ)
        b = zeros(atype,N)
        b[j] = 1.0
        m = length(vec_left)

    #--Line 2 in Table III.
        x = zeros(atype,N)
        r = copy(b)
        p = copy(b)
        αm = 1.0
        βm = 0.0

    #--
        Σ =zeros(atype,m)
        for mm=1:m
            Σ[mm] = b[vec_left[mm]] #Line 3
        end

        Θ = zeros(σtype,M,m)
        Π = ones(σtype,M,m)
        for mm=1:m
            Π[:,mm] *= Σ[mm]
        end
        ρk = ones(σtype,M)
        ρkm = copy(ρk)
        ρkp = copy(ρk)
        Ap = similar(p)

        for k=0:maximumsteps
            mul!(Ap,A,-p)
            #A_mul_B!(Ap,A,-p)
            rnorm = dot(r,r)
            α = rnorm/dot(p,Ap)
            @. x += α*p #Line 7
            @. r += -α*Ap #Line 8
            β = r'*r/rnorm #Line9
            @. p = r + β*p #Line 10

            for mm=1:m
                Σ[mm] = r[vec_left[mm]] #Line 11
            end

            for j = 1:M
                update = ifelse(abs(ρk[j]) > eps,true,false)
                if update
                    ρkp[j] = ρk[j]*ρkm[j]*αm/(ρkm[j]*αm*(1.0+α*σ[j])+α*βm*(ρkm[j]-ρk[j]))#Line 13
                    αkj = α*ρkp[j]/ρk[j]#Line 14
                    @. Θ[j,:] += αkj*Π[j,:] #Line 15
                    βkj = ((ρkp[j]/ρk[j])^2)*β #Line 16
                    Π[j,:] = ρkp[j]*Σ+ βkj*Π[j,:] #Line 17

                end
                ρkm[j] = ρk[j]
                ρk[j] = ρkp[j]
            end
            αm = α
            βm = β
            hi = real(rnorm)*maximum(abs.(ρk))
            if hi < eps
                return Θ
            end
        end

        println("After $maximumsteps steps, this was not converged")
        return Θ

    end

end # module


text:
function coreVIfull(logl::Function, μarray::Array{Array{Float64,1},1}, Σarray::Array{Array{Float64, 2},1}; gradlogl = gradlogl, seed = 1, S = 100, optimiser=Optim.LBFGS(), iterations=1, numerical_verification = false, Stest=0, show_every=-1, inititerations=0)

    D = length(μarray[1])

    @assert(D == size(Σarray[1], 1) == size(Σarray[1], 2))

    @assert(length(μarray) == length(Σarray))

    @printf("Running VI with full covariance: seed=%d, S=%d, Stest=%d, D=%d for %d iterations\n", seed, S, Stest, D, iterations)

    #----------------------------------------------------
    # generate latent variables
    #----------------------------------------------------

    Ztrain = generatelatentZ(S = S, D = D, seed = seed)

    Ztest  = generatelatentZ(S = Stest, D = D, seed = seed+1)

    #----------------------------------------------------
    function unpack(param)
    #----------------------------------------------------

        @assert(length(param) == D+D*D)

        local μ = param[1:D]

        local C = reshape(param[D+1:D+D*D], D, D)

        return μ, C

    end

    #----------------------------------------------------
    function minauxiliary(param)
    #----------------------------------------------------

        local μ, C = unpack(param)

        return -1.0 * elbo(μ, C, Ztrain)

    end

    #----------------------------------------------------
    function minauxiliary_grad(param)
    #----------------------------------------------------

        local μ, C = unpack(param)

        return -1.0 * elbo_grad(μ, C, Ztrain)

    end

    #----------------------------------------------------
    function getcov(C)
    #----------------------------------------------------

        return C*C'

    end

    #----------------------------------------------------
    function getcovroot(C)
    #----------------------------------------------------

        return C

    end

    #----------------------------------------------------
    function elbo(μ, C, Z)
    #----------------------------------------------------

        mean(map(z -> logl(μ .+ C*z), Z)) + ApproximateVI.entropy(C)

    end

    #----------------------------------------------------
    function elbo_grad(μ, C, Z)
    #----------------------------------------------------

        local Σroot = getcovroot(C)
        local gradC = (Σroot\I)' # entropy contribution
        local gradμ = zeros(eltype(μ), D)
        local S     = length(Z)

        for s=1:S
            g = gradlogl(μ .+ Σroot*Z[s])
            gradC += (1/S)*g*Z[s]'
            gradμ += g/S
        end

        return [vec(gradμ); vec(gradC)]

    end

    gradhelper(storage, param) = copyto!(storage, minauxiliary_grad(param))

    #----------------------------------------------------
    # Numerically verify gradient
    #----------------------------------------------------

    if numerical_verification

        local C = Matrix(cholesky(Σarray[1]).L)
        adgrad = ForwardDiff.gradient(minauxiliary, [μarray[1]; vec(C)])
        angrad = minauxiliary_grad([μarray[1];vec(C)])
        @printf("gradient from AD vs analytical gradient\n")
        display([vec(adgrad) vec(angrad)])
        @printf("maximum absolute difference is %f\n", maximum(abs.(vec(adgrad) - vec(angrad))))

    end

    #----------------------------------------------------
    # Evaluate initial solutions for few iterations
    #----------------------------------------------------

    initoptimise(μ, Σ) = Optim.optimize(minauxiliary, gradhelper, [μ; vec(Matrix(cholesky(Σ).L))], optimiser, Optim.Options(iterations = inititerations))

    results = if inititerations>0
        @showprogress "Initial search with random start " map(initoptimise, μarray, Σarray)
    else
        map(initoptimise, μarray, Σarray)
    end

    bestindex = argmin(map(r -> r.minimum, results))

    bestinitialsolution = results[bestindex].minimizer

    #----------------------------------------------------
    # Call optimiser
    #----------------------------------------------------

    options = Optim.Options(extended_trace = false, store_trace = false, show_trace = show_every > 0, show_every = show_every, iterations = iterations, g_tol = 1e-6)

    result  = Optim.optimize(minauxiliary, gradhelper, bestinitialsolution, optimiser, options)

    μopt, Copt = unpack(result.minimizer)

    #----------------------------------------------------
    # Return results
    #----------------------------------------------------

    Σopt = getcov(Copt)

    return MvNormal(μopt, Σopt), elbo(μopt, Copt, generatelatentZ(S = 10*S, D = D, seed = seed+2))

end


text:
include("FourierWavelet.jl")
include("Collect.jl")
include("Multiplication.jl")
include("NFFT.jl")
include("misc.jl")
include("Freq2Wave.jl")
include("CGNR.jl")


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(sec(c+d*x)^6*(A+C*sec(c+d*x)^2), x) == :(1//7*d^-1*(6C+7A)*tan(c+d*x)+1//21*d^-1*tan(c+d*x)^3*(12C+14A)+1//35*d^-1*tan(c+d*x)^5*(6C+7A)+1//7*C*d^-1*sec(c+d*x)^6*tan(c+d*x))
@test integrate(sec(c+d*x)^5*(A+C*sec(c+d*x)^2), x) == :(1//16*d^-1*(5C+6A)*arctanh(sin(c+d*x))+1//6*C*d^-1*sec(c+d*x)^5*tan(c+d*x)+1//16*d^-1*(5C+6A)*sec(c+d*x)*tan(c+d*x)+1//24*d^-1*sec(c+d*x)^3*(5C+6A)*tan(c+d*x))
@test integrate(sec(c+d*x)^4*(A+C*sec(c+d*x)^2), x) == :(1//5*d^-1*(4C+5A)*tan(c+d*x)+1//15*d^-1*tan(c+d*x)^3*(4C+5A)+1//5*C*d^-1*sec(c+d*x)^4*tan(c+d*x))
@test integrate(sec(c+d*x)^3*(A+C*sec(c+d*x)^2), x) == :(1//8*d^-1*(3C+4A)*arctanh(sin(c+d*x))+1//4*C*d^-1*sec(c+d*x)^3*tan(c+d*x)+1//8*d^-1*(3C+4A)*sec(c+d*x)*tan(c+d*x))
@test integrate(sec(c+d*x)^2*(A+C*sec(c+d*x)^2), x) == :(1//3*d^-1*(2C+3A)*tan(c+d*x)+1//3*C*d^-1*sec(c+d*x)^2*tan(c+d*x))
@test integrate((A+C*sec(c+d*x)^2)*sec(c+d*x), x) == :((1/2)*d^-1*(C+2A)*arctanh(sin(c+d*x))+(1/2)*C*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate(A+C*sec(c+d*x)^2, x) == :(A*x+C*d^-1*tan(c+d*x))
@test integrate((A+C*sec(c+d*x)^2)*cos(c+d*x), x) == :(A*d^-1*sin(c+d*x)+C*d^-1*arctanh(sin(c+d*x)))
@test integrate(cos(c+d*x)^2*(A+C*sec(c+d*x)^2), x) == :(x*(C+(1/2)*A)+(1/2)*A*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^3*(A+C*sec(c+d*x)^2), x) == :(d^-1*(A+C)*sin(c+d*x)+-1//3*A*d^-1*sin(c+d*x)^3)
@test integrate(cos(c+d*x)^4*(A+C*sec(c+d*x)^2), x) == :(x*((1/2)*C+3//8*A)+1//4*A*d^-1*cos(c+d*x)^3*sin(c+d*x)+1//8*d^-1*(3A+4C)*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^5*(A+C*sec(c+d*x)^2), x) == :(d^-1*(A+C)*sin(c+d*x)+-1//3*d^-1*sin(c+d*x)^3*(C+2A)+1//5*A*d^-1*sin(c+d*x)^5)
@test integrate(cos(c+d*x)^6*(A+C*sec(c+d*x)^2), x) == :(x*(3//8*C+5//16*A)+1//6*A*d^-1*cos(c+d*x)^5*sin(c+d*x)+1//16*d^-1*(5A+6C)*cos(c+d*x)*sin(c+d*x)+1//24*d^-1*cos(c+d*x)^3*(5A+6C)*sin(c+d*x))
@test integrate(sec(c+d*x)^m*(C*sec(c+d*x)^2+-1*C*m*(1+m)^-1), x) == :(C*d^-1*(1+m)^-1*sec(c+d*x)^(1+m)*sin(c+d*x))
@test integrate(sec(c+d*x)^m*(A+-1*A*m^-1*sec(c+d*x)^2*(1+m)), x) == :(-1*A*d^-1*m^-1*sec(c+d*x)^(1+m)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^5//2*(A+C*sec(c+d*x)^2), x) == :(2//7*C*d^-1*(b*sec(c+d*x))^5//2*tan(c+d*x)+2//21*b*d^-1*(b*sec(c+d*x))^3//2*(5C+7A)*sin(c+d*x)+2//21*b^2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(5C+7A)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*(A+C*sec(c+d*x)^2), x) == :(2//5*C*d^-1*(b*sec(c+d*x))^3//2*tan(c+d*x)+2//5*b*d^-1*(b*sec(c+d*x))^(1/2)*(3C+5A)*sin(c+d*x)+-2//5*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(3C+5A)*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*(A+C*sec(c+d*x)^2), x) == :(2//3*C*d^-1*(b*sec(c+d*x))^(1/2)*tan(c+d*x)+1//3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(2C+6A)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*(A+C*sec(c+d*x)^2), x) == :(2*C*d^-1*(b*sec(c+d*x))^-1//2*tan(c+d*x)+d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(-2C+2A)*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*(A+C*sec(c+d*x)^2), x) == :(2//3*A*d^-1*(b*sec(c+d*x))^-3//2*tan(c+d*x)+1//3*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(2A+6C)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*(A+C*sec(c+d*x)^2), x) == :(2//5*A*d^-1*(b*sec(c+d*x))^-5//2*tan(c+d*x)+1//5*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(6A+10C)*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-7//2*(A+C*sec(c+d*x)^2), x) == :(2//7*A*d^-1*(b*sec(c+d*x))^-7//2*tan(c+d*x)+1//21*b^-3*d^-1*(b*sec(c+d*x))^-1//2*(10A+14C)*sin(c+d*x)+1//21*b^-4*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(10A+14C)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-9//2*(A+C*sec(c+d*x)^2), x) == :(2//9*A*d^-1*(b*sec(c+d*x))^-9//2*tan(c+d*x)+1//45*b^-3*d^-1*(b*sec(c+d*x))^-3//2*(14A+18C)*sin(c+d*x)+1//15*b^-4*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(14A+18C)*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate(sec(c+d*x)^-1//2*(3+3*sec(c+d*x)^2), x) == :(6*d^-1*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate(sec(e+f*x)^m*(m+-1*sec(e+f*x)^2*(1+m)), x) == :(-1*f^-1*sec(e+f*x)^(1+m)*sin(e+f*x))
@test integrate(sec(e+f*x)^5*(5+-6*sec(e+f*x)^2), x) == :(-1*f^-1*sec(e+f*x)^5*tan(e+f*x))
@test integrate(sec(e+f*x)^4*(4+-5*sec(e+f*x)^2), x) == :(-1*f^-1*sec(e+f*x)^4*tan(e+f*x))
@test integrate(sec(e+f*x)^3*(3+-4*sec(e+f*x)^2), x) == :(-1*f^-1*sec(e+f*x)^3*tan(e+f*x))
@test integrate(sec(e+f*x)^2*(2+-3*sec(e+f*x)^2), x) == :(-1*f^-1*sec(e+f*x)^2*tan(e+f*x))
@test integrate((1+-2*sec(e+f*x)^2)*sec(e+f*x), x) == :(-1*f^-1*sec(e+f*x)*tan(e+f*x))
@test integrate(-1*sec(e+f*x)^2, x) == :(-1*f^-1*tan(e+f*x))
@test integrate(-1*cos(e+f*x), x) == :(-1*f^-1*sin(e+f*x))
@test integrate(cos(e+f*x)^2*(-2+sec(e+f*x)^2), x) == :(-1*f^-1*cos(e+f*x)*sin(e+f*x))
@test integrate(cos(e+f*x)^3*(-3+2*sec(e+f*x)^2), x) == :(-1*f^-1*cos(e+f*x)^2*sin(e+f*x))
@test integrate(cos(e+f*x)^4*(-4+3*sec(e+f*x)^2), x) == :(-1*f^-1*cos(e+f*x)^3*sin(e+f*x))
@test integrate(cos(e+f*x)^5*(-5+4*sec(e+f*x)^2), x) == :(-1*f^-1*cos(e+f*x)^4*sin(e+f*x))
@test integrate(sec(c+d*x)^3*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(B*d^-1*tan(c+d*x)+1//3*B*d^-1*tan(c+d*x)^3+3//8*C*d^-1*arctanh(sin(c+d*x))+1//4*C*d^-1*sec(c+d*x)^3*tan(c+d*x)+3//8*C*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate(sec(c+d*x)^2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(C*d^-1*tan(c+d*x)+(1/2)*B*d^-1*arctanh(sin(c+d*x))+1//3*C*d^-1*tan(c+d*x)^3+(1/2)*B*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate((B*sec(c+d*x)+C*sec(c+d*x)^2)*sec(c+d*x), x) == :(B*d^-1*tan(c+d*x)+(1/2)*C*d^-1*arctanh(sin(c+d*x))+(1/2)*C*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate(B*sec(c+d*x)+C*sec(c+d*x)^2, x) == :(B*d^-1*arctanh(sin(c+d*x))+C*d^-1*tan(c+d*x))
@test integrate((B*sec(c+d*x)+C*sec(c+d*x)^2)*cos(c+d*x), x) == :(B*x+C*d^-1*arctanh(sin(c+d*x)))
@test integrate(cos(c+d*x)^2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(C*x+B*d^-1*sin(c+d*x))
@test integrate(cos(c+d*x)^3*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :((1/2)*B*x+C*d^-1*sin(c+d*x)+(1/2)*B*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^4*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :((1/2)*C*x+B*d^-1*sin(c+d*x)+-1//3*B*d^-1*sin(c+d*x)^3+(1/2)*C*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^5*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(3//8*B*x+C*d^-1*sin(c+d*x)+-1//3*C*d^-1*sin(c+d*x)^3+1//4*B*d^-1*cos(c+d*x)^3*sin(c+d*x)+3//8*B*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^6*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(3//8*C*x+B*d^-1*sin(c+d*x)+-2//3*B*d^-1*sin(c+d*x)^3+1//5*B*d^-1*sin(c+d*x)^5+1//4*C*d^-1*cos(c+d*x)^3*sin(c+d*x)+3//8*C*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//3*B*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//5*C*b^-1*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*C*b*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-6//5*C*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//3*B*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2*B*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+2//3*C*b^-1*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*C*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2)+-2*B*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(-2*C*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2*C*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+2*B*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2*B*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2*C*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//3*B*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2*C*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//3*B*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-7//2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//3*C*b^-3*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//5*B*b^-2*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+2//3*C*b^-4*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2)+6//5*B*b^-3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-9//2*(B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//5*C*b^-3*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+2//7*B*b^-2*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*B*b^-4*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+6//5*C*b^-4*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+10//21*B*b^-5*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate(sec(c+d*x)^4*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(1//5*d^-1*(4C+5A)*tan(c+d*x)+1//15*d^-1*tan(c+d*x)^3*(4C+5A)+3//8*B*d^-1*arctanh(sin(c+d*x))+1//4*B*d^-1*sec(c+d*x)^3*tan(c+d*x)+1//5*C*d^-1*sec(c+d*x)^4*tan(c+d*x)+3//8*B*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate(sec(c+d*x)^3*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(B*d^-1*tan(c+d*x)+1//3*B*d^-1*tan(c+d*x)^3+1//8*d^-1*(3C+4A)*arctanh(sin(c+d*x))+1//4*C*d^-1*sec(c+d*x)^3*tan(c+d*x)+1//8*d^-1*(3C+4A)*sec(c+d*x)*tan(c+d*x))
@test integrate(sec(c+d*x)^2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :((1/2)*B*d^-1*arctanh(sin(c+d*x))+1//3*d^-1*(2C+3A)*tan(c+d*x)+(1/2)*B*d^-1*sec(c+d*x)*tan(c+d*x)+1//3*C*d^-1*sec(c+d*x)^2*tan(c+d*x))
@test integrate((A+B*sec(c+d*x)+C*sec(c+d*x)^2)*sec(c+d*x), x) == :(B*d^-1*tan(c+d*x)+(1/2)*d^-1*(C+2A)*arctanh(sin(c+d*x))+(1/2)*C*d^-1*sec(c+d*x)*tan(c+d*x))
@test integrate(A+B*sec(c+d*x)+C*sec(c+d*x)^2, x) == :(A*x+B*d^-1*arctanh(sin(c+d*x))+C*d^-1*tan(c+d*x))
@test integrate((A+B*sec(c+d*x)+C*sec(c+d*x)^2)*cos(c+d*x), x) == :(B*x+A*d^-1*sin(c+d*x)+C*d^-1*arctanh(sin(c+d*x)))
@test integrate(cos(c+d*x)^2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(x*(C+(1/2)*A)+B*d^-1*sin(c+d*x)+(1/2)*A*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^3*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :((1/2)*B*x+d^-1*(A+C)*sin(c+d*x)+-1//3*A*d^-1*sin(c+d*x)^3+(1/2)*B*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^4*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(x*((1/2)*C+3//8*A)+B*d^-1*sin(c+d*x)+-1//3*B*d^-1*sin(c+d*x)^3+1//4*A*d^-1*cos(c+d*x)^3*sin(c+d*x)+1//8*d^-1*(3A+4C)*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^5*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(3//8*B*x+d^-1*(A+C)*sin(c+d*x)+-1//3*d^-1*sin(c+d*x)^3*(C+2A)+1//5*A*d^-1*sin(c+d*x)^5+1//4*B*d^-1*cos(c+d*x)^3*sin(c+d*x)+3//8*B*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(cos(c+d*x)^6*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(x*(3//8*C+5//16*A)+B*d^-1*sin(c+d*x)+-2//3*B*d^-1*sin(c+d*x)^3+1//5*B*d^-1*sin(c+d*x)^5+1//6*A*d^-1*cos(c+d*x)^5*sin(c+d*x)+1//16*d^-1*(5A+6C)*cos(c+d*x)*sin(c+d*x)+1//24*d^-1*cos(c+d*x)^3*(5A+6C)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//3*B*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//5*C*d^-1*(b*sec(c+d*x))^3//2*tan(c+d*x)+2//5*b*d^-1*(b*sec(c+d*x))^(1/2)*(3C+5A)*sin(c+d*x)+-2//5*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(3C+5A)*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//3*B*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2*B*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+2//3*C*d^-1*(b*sec(c+d*x))^(1/2)*tan(c+d*x)+1//3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(2C+6A)*Elliptic.F((1/2)*c+(1/2)*d*x,2)+-2*B*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2*C*d^-1*(b*sec(c+d*x))^-1//2*tan(c+d*x)+d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(-2C+2A)*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2*B*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//3*A*d^-1*(b*sec(c+d*x))^-3//2*tan(c+d*x)+2*B*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+1//3*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(2A+6C)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//5*A*d^-1*(b*sec(c+d*x))^-5//2*tan(c+d*x)+2//3*B*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+1//5*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*(6A+10C)*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//3*B*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-7//2*(A+B*sec(c+d*x)+C*sec(c+d*x)^2), x) == :(2//7*A*d^-1*(b*sec(c+d*x))^-7//2*tan(c+d*x)+1//21*b^-3*d^-1*(b*sec(c+d*x))^-1//2*(10A+14C)*sin(c+d*x)+2//5*B*b^-2*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+1//21*b^-4*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*(10A+14C)*Elliptic.F((1/2)*c+(1/2)*d*x,2)+6//5*B*b^-3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))

max_stars_repo_path: test/4 Trig functions/4.5 Secant/4.5.4.1 (a+b sec)^m (A+B sec+C sec^2).jl


text:
"""
ExtendedChebyshev{kind,T}()

is a quasi-matrix representing the Hilbert transform of ChebyshevU polynomials of the specified kind 2
on ℝ.

For an untransformed ExtendedChebyshevT, for x∈[-1,1]  ExtendedChebyshevT[x, n] = ChebyshevT()[x, n+1]
and outside the interval x∉[-1,1], there is an explicit formula for ExtendedChebyshevT[x, n].
"""

struct ExtendedChebyshev{kind,T} <: Basis{T} end
ExtendedChebyshev{kind}() where kind = ExtendedChebyshev{kind,Float64}()

axes(H::ExtendedChebyshev) = (Inclusion(ℝ), OneToInf())

==(a::ExtendedChebyshev{kind}, b::ExtendedChebyshev{kind}) where kind = true
==(a::ExtendedChebyshev, b::ExtendedChebyshev) = false

####
# Extended Chebyshev T
####

const ExtendedChebyshevT = ExtendedChebyshev{1}

function getindex(H::ExtendedChebyshevT{T}, x::Real, j::Int)::T where T
    x in ChebyshevInterval() && return ChebyshevT{T}()[x,j]

    j < 1 && throw_boundserror(H, [x,j])

    # ξ = inv(x + sqrtx2(x))
    # return transpose(ξ.^(j-1))
    if j == 1
        # return zero(T) # Experimental decaying support
        return one(T)
    else
        ξ = inv(x + sqrtx2(x))
        return transpose(ξ.^(j-1))
    end
end

"""
    extendedchebyshevt(n, z)

computes the `n`-th extended chebyshev T at `z`.
"""
extendedchebyshevt(n::Integer, z::Number) = Base.unsafe_getindex(ExtendedChebyshevT{typeof(z)}(), z, n+1)

# hilbertweightedchebyshevu() = HilbertWeightedChebyshevU()
# hilbertweightedchebyshevu(d::AbstractInterval{T}) where T = HilbertWeightedChebyshevU{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
# hilbertweightedchebyshevu(d::Inclusion) = hilbertweightedchebyshevu(d.domain)
# hilberweightedchebyshevu(S::AbstractQuasiMatrix) = hilberweightedchebyshevu(axes(S,1))

####
# Extended Chebyshev U
####
const ExtendedChebyshevU = ExtendedChebyshev{2}

function getindex(H::ExtendedChebyshevU{T}, x::Real, j::Int)::T where T

    j < 1 && throw_boundserror(H, [x,j])

    # eU_-1 = HilbertTransform[wT_1] - 1
    if j == 1
        x in ChebyshevInterval() && return zero(T)
        ξ = - sign(x) .* x ./ sqrt.(x.^2 .- one(T))
    # eU_0 = HilbertTransform[wT_0]
    elseif j == 2
        x in ChebyshevInterval() && return zero(T)
        ξ = - sign(x) ./ sqrt.(x.^2 .- one(T))
    elseif isodd(j)
        x in ChebyshevInterval() && return ChebyshevU{T}()[x,j-2]
        η = inv(x + sqrtx2(x))
        ξ = 2 .* sum(η .^(2:2:j-3)) + one(T) .- sign(x) .* x ./ sqrt.(x.^2 .- one(T))
        # ξ = 2 .* sum(η .^(2:2:j-3)) .- abs(x) ./ sqrt.(x.^2 .- one(T)) # decaying support
    else
        x in ChebyshevInterval() && return ChebyshevU{T}()[x,j-2]
        η = inv(x + sqrtx2(x))
        ξ = 2 .* sum(η .^(1:2:j-3)) .- sign(x) ./ sqrt.(x.^2 .- one(T))
    end

    return transpose(ξ)
end

struct ExtendedWeightedChebyshev{kind,T} <: Basis{T} end
ExtendedWeightedChebyshev{kind}() where kind = ExtendedWeightedChebyshev{kind,Float64}()

axes(H::ExtendedWeightedChebyshev) = (Inclusion(ℝ), OneToInf())

==(a::ExtendedWeightedChebyshev{kind}, b::ExtendedWeightedChebyshev{kind}) where kind = true
==(a::ExtendedWeightedChebyshev, b::ExtendedWeightedChebyshev) = false

###
# Extended Weighted Chebyshev T
###

const ExtendedWeightedChebyshevT = ExtendedWeightedChebyshev{1}

function getindex(H::ExtendedWeightedChebyshevT{T}, x::Real, j::Int)::T where T
    x in ChebyshevInterval() && return Weighted(ChebyshevT{T}())[x,j]
    return 0.
end

###
# Extended Weighted Chebyshev U
###

const ExtendedWeightedChebyshevU = ExtendedWeightedChebyshev{2}

function getindex(H::ExtendedWeightedChebyshevU{T}, x::Real, j::Int)::T where T
    x in ChebyshevInterval() && return Weighted(ChebyshevU{T}())[x,j]
    return 0.
end


text:
load("core.jl")

total_results = zeros(0, 6)

#env = TwoSeasonalArms(100)
env = FiveDriftingArms()

# for exp3_gamma in [0.1:0.1:0.5]
#   algo = Exp3(n_arms(env), exp3_gamma)
#   results = test_algorithm(algo, env, 5000, 400)
#   results = hcat(repmat([exp3_gamma], size(results, 1), 1), results)
#   total_results = vcat(total_results, results)
# end

for tau in [0.1:0.1:0.5]
  algo = Softmax(n_arms(env), tau)
  results = test_algorithm(algo, env, 5000, 400)
  results = hcat(repmat([tau], size(results, 1), 1), results)
  total_results = vcat(total_results, results)
end

csvwrite("env_demo.csv", total_results)

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
@doc doc"""
    SymmetricMatrices{n,F} <: Manifold

The [`Manifold`](@ref) $ \operatorname{Sym} (n)$ consisting of the real- or complex-valued
symmetric matrices of size $ n\times n$, i.e. the set

````math
\operatorname{Sym}(n) = \bigl\{A \in \mathbb F^{n\times n} \big| A^{\mathrm{H}} = A \bigr\},
````
where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed
and the field $\mathbb F \in \{ \mathbb R, \mathbb C\}$ is set by the
[`AbstractNumbers`](@ref) `F`.

Though it is slighty redundant, usually the matrices are safed as $n\times n$ arrays.

# Constructor

    SymmetricMatrices(n::Int, F::AbstractNumbers=ℝ)

Generate the manifold of $n\times n$ symmetric metrices.
"""
struct SymmetricMatrices{n,F} <: Manifold end

SymmetricMatrices(n::Int, F::AbstractNumbers = ℝ) = SymmetricMatrices{n,F}()

@doc doc"""
    check_manifold_point(M::SymmetricMatrices{n,F}, x; kwargs...)

Check whether `x` is a valid manifold point on the [`SymmetricMatrices`](@ref) `M`, i.e.
whether `x` is a symmetric matrix of size `(n,n)` with values from the corresponding
[`AbstractNumbers`](@ref) `F`.

The tolerance for the symmetry of `x` can be set using `kwargs...`.
"""
function check_manifold_point(M::SymmetricMatrices{n,F}, x; kwargs...) where {n,F}
    if (F === ℝ) && !(eltype(x) <: Real)
        return DomainError(
            eltype(x),
            "The matrix $(x) does not lie on $M, since its values are not real.",
        )
    end
    if (F === ℂ) && !(eltype(x) <: Real) && !(eltype(x) <: Complex)
        return DomainError(
            eltype(x),
            "The matrix $(x) does not lie on $M, since its values are not complex.",
        )
    end
    if size(x) != (n, n)
        return DomainError(
            size(x),
            "The point $(x) does not lie on $M since its size ($(size(x))) does not match the representation size ($(representation_size(M))).",
        )
    end
    if !isapprox(norm(x - transpose(x)), 0.0; kwargs...)
        return DomainError(
            norm(x - transpose(x)),
            "The point $(x) does not lie on $M, since it is not symmetric.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::SymmetricMatrices{n,F}, x, v; kwargs... )

Check whether `v` is a tangent vector to manifold point `x` on the
[`SymmetricMatrices`](@ref) `M`, i.e. `v` has to be a symmetric matrix of dimension `(n,n)`
and its values have to be from the correct [`AbstractNumbers`](@ref).

The tolerance for the symmetry of `x` and `v` can be set using `kwargs...`.
"""
function check_tangent_vector(M::SymmetricMatrices{n,F}, x, v; kwargs...) where {n,F}
    t = check_manifold_point(M, x; kwargs...)
    t === nothing || return t
    if (F === ℝ) && !(eltype(v) <: Real)
        return DomainError(
            eltype(v),
            "The matrix $(v) is not a tangent to a point on $M, since its values are not real.",
        )
    end
    if (F === ℂ) && !(eltype(v) <: Real) && !(eltype(v) <: Complex)
        return DomainError(
            eltype(v),
            "The matrix $(v) is not a tangent to a point on $M, since its values are not complex.",
        )
    end
    if size(v) != (n, n)
        return DomainError(
            size(v),
            "The vector $(v) is not a tangent to a point on $(M) since its size ($(size(v))) does not match the representation size ($(representation_size(M))).",
        )
    end
    if !isapprox(norm(v - transpose(v)), 0.0; kwargs...)
        return DomainError(
            norm(v - transpose(v)),
            "The vector $(v) is not a tangent vector to $(x) on $(M), since it is not symmetric.",
        )
    end
    return nothing
end

@doc doc"""
    distance(M::SymmetricMatrices, x, y)

Compute distance using the inherited metric, i.e. taking the Frobenius-norm of the
difference.
"""
distance(M::SymmetricMatrices, x, y) = norm(x - y)

@doc doc"""
    exp(M::SymmetricMatrices, x, v)

Compute the exponential map eminating from `x` in tangent direction `v` on the
[`SymmetricMatrices`](@ref) `M`, which reads

````math
\exp_xv = x + v.
````
"""
exp(::SymmetricMatrices, ::Any...)

exp!(M::SymmetricMatrices, y, x, v) = (y .= x .+ v)

@doc doc"""
    flat(M::SymmetricMatrices, x, w::FVector{TangentSpaceType})

Compute the [`flat`](@ref flat(M::Manifold, x, w::FVector)) isomorphism of the
[`SymmetricMatrices`](@ref) `M` on the manifold point `x` and tangent vector `w`.

Since `M` is already a vector space over $\mathbb R$, this returns just the vector `w`.
"""
flat(::SymmetricMatrices, ::Any...)

flat!(M::SymmetricMatrices, v::CoTFVector, x, w::TFVector) = copyto!(v, w)

function get_coordinates(
    M::SymmetricMatrices{N,ℝ},
    x,
    v,
    B::ArbitraryOrthonormalBasis{ℝ},
) where {N}
    dim = manifold_dimension(M)
    vout = similar(v, dim)
    @assert size(v) == (N, N)
    @assert dim == div(N * (N + 1), 2)
    k = 1
    for i = 1:N, j = i:N
        scale = ifelse(i == j, 1, sqrt(2))
        @inbounds vout[k] = v[i, j] * scale
        k += 1
    end
    return vout
end
function get_coordinates(
    M::SymmetricMatrices{N,ℂ},
    x,
    v,
    B::ArbitraryOrthonormalBasis{ℝ},
) where {N}
    dim = manifold_dimension(M)
    vout = similar(v, dim)
    @assert size(v) == (N, N)
    @assert dim == N * (N + 1)
    k = 1
    for i = 1:N, j = i:N
        scale = ifelse(i == j, 1, sqrt(2))
        @inbounds vout[k] = real(v[i, j]) * scale
        k += 1
        @inbounds vout[k] = imag(v[i, j]) * scale
        k += 1
    end
    return vout
end

function get_vector(
    M::SymmetricMatrices{N,ℝ},
    x,
    v,
    B::ArbitraryOrthonormalBasis{ℝ},
) where {N}
    dim = manifold_dimension(M)
    vout = allocate_result(M, get_vector, x)
    @assert size(v) == (div(N * (N + 1), 2),)
    @assert size(vout) == (N, N)
    k = 1
    for i = 1:N, j = i:N
        scale = ifelse(i == j, 1, 1 / sqrt(2))
        @inbounds vout[i, j] = v[k] * scale
        @inbounds vout[j, i] = v[k] * scale
        k += 1
    end
    return vout
end
function get_vector(
    M::SymmetricMatrices{N,ℂ},
    x,
    v,
    B::ArbitraryOrthonormalBasis{ℝ},
) where {N}
    dim = manifold_dimension(M)
    vout = allocate_result(M, get_vector, x, x .* 1im)
    @assert size(v) == (N * (N + 1),)
    @assert size(vout) == (N, N)
    k = 1
    for i = 1:N, j = i:N
        scale = ifelse(i == j, 1, 1 / sqrt(2))
        @inbounds vout[i, j] = Complex(v[k], v[k+1]) * scale
        @inbounds vout[j, i] = vout[i, j]
        k += 2
    end
    return vout
end

@doc doc"""
    inner(M::SymmetricMatrices, x, w, v)

Compute the inner product of the two tangent vectors `w`, `v` from the tangent
space at `x` on the [`SymmetricMatrices`](@ref) `M` using the restriction of the
metric from the embedding, i.e.
````math
(v,w)_x = \operatorname{tr}(v^{\mathrm{H}}w),
````
where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
@inline inner(M::SymmetricMatrices, x, w, v) = dot(w, v)

@doc doc"""
    log(M::SymmetricMatrices, x, y)
Compute the logarithmic map from `x` to `y` on the [`SymmetricMatrices`](@ref) `M`, which
reads

````math
\log_xy = y-x.
````
"""
log(::SymmetricMatrices, ::Any...)

log!(M::SymmetricMatrices, v, x, y) = (v .= y .- x)

@doc doc"""
manifold_dimension(M::SymmetricMatrices{n,𝔽})

Return the dimension of the [`SymmetricMatrices`](@ref) matrix `M` over the number system
`𝔽`, i.e.

````math
\dim \operatorname{Sym}(n,𝔽) = \frac{n(n+1)}{2} \dim_ℝ 𝔽,
````

where $\dim_ℝ 𝔽$ is the [`real_dimension`](@ref) of `𝔽`.
"""
function manifold_dimension(::SymmetricMatrices{N,𝔽}) where {N,𝔽}
    return div(N * (N + 1), 2) * real_dimension(𝔽)
end

@doc doc"""
    norm(M::SymmetricMatrices, x, v)

Compute the norm of the tangent vector `v` from the tangent space at `x` on the
[`SymmetricMatrices`](@ref) `M`, which is the norm from the embedding, i.e.

````math
\lVert v \rVert_x = \lVert v \rVert_2
````
"""
norm(M::SymmetricMatrices, x, v) = norm(v)

@doc doc"""
    project_point(M::SymmetricMatrices,x)

Projects `x` from the embedding onto the [`SymmetricMatrices`](@ref) `M`, i.e.

````math
\operatorname{proj}_{\operatorname{Sym}(n)}(x) = \frac{1}{2} \bigl( x + x^{\mathrm{H}} \bigr),
````

where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
project_point(::SymmetricMatrices, ::Any...)

project_point!(M::SymmetricMatrices, x) = (x .= (x + transpose(x)) ./ 2)

@doc doc"""
    project_tangent(M::SymmetricMatrices, x, v)

Project the matrix `v` onto the tangent space at `x` on the [`SymmetricMatrices`](@ref) `M`,

````math
\operatorname{proj}_x(v) = \frac{1}{2} \bigl( v + v^{\mathrm{H}} \bigr),
````

where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
project_tangent(::SymmetricMatrices, ::Any...)

project_tangent!(M::SymmetricMatrices, w, x, v) = (w .= (v .+ transpose(v)) ./ 2)

@doc doc"""
    representation_size(M::SymmetricMatrices)

Returns the size points on the [`SymmetricMatrices`](@ref) `M` are represented as, i.e.
for the $n\times n$ it's `(n,n)`.
"""
@generated representation_size(::SymmetricMatrices{N}) where {N} = (N, N)

@doc doc"""
    sharp(M::SymmetricMatrices, x, w::FVector{CotangentSpaceType})

Compute the [`sharp`](@ref sharp(M::Manifold, x, w::FVector)) isomorphism of the
[`SymmetricMatrices`](@ref) `M` on the manifold point `x` and cotangent vector `w`.

Since `M` is already a vector space over $\mathbb R$, this returns just the vector `w`.
"""
sharp(::SymmetricMatrices, ::Any...)

sharp!(M::SymmetricMatrices, v::TFVector, x, w::CoTFVector) = copyto!(v, w)

@doc doc"""
    vector_transport_to(M::SymmetricMatrices, x, v, y, ::ParallelTransport)

Compute the parallel
[`vector_transport_to`](@ref vector_transport_to(M::Manifold, x, v, y, ParallelTransport()))
of `v` from the tangent space at `x` on the [`SymmetricMatrices`](@ref) `M` to `y`.
Since the metric is inherited from the embedding space, this is just the identity, i.e.

````math
P_{y\gets x}(v) = v.
````
"""
vector_transport_to(::SymmetricMatrices, ::Any...)

function vector_transport_to!(M::SymmetricMatrices, vto, x, v, y, ::ParallelTransport)
    return copyto!(vto, v)
end

@doc doc"""
    zero_tangent_vector(M, x)

Return the zero tangent vector for the tangent space at `x` on the
[`SymmetricMatrices`](@ref) `M`, i.e. the zero matrix.
"""
zero_tangent_vector(::SymmetricMatrices, ::Any...)

zero_tangent_vector!(M::SymmetricMatrices, v, x) = fill!(v, 0)


text:
#! /usr/bin/julia
#
#	text_create.jl
#
#					Jul/04/2018
# --------------------------------------------------------------------
#
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
function data_prepare_proc()
	dict_aa = Dict()
#
	dict_aa = dict_append_proc(dict_aa,"t2381","名古屋",32761,"2013-5-30")
	dict_aa = dict_append_proc(dict_aa,"t2382","豊橋",87295,"2013-9-10")
	dict_aa = dict_append_proc(dict_aa,"t2383","岡崎",21674,"2013-2-14")
	dict_aa = dict_append_proc(dict_aa,"t2384","一宮",73612,"2013-3-9")
	dict_aa = dict_append_proc(dict_aa,"t2385","蒲郡",52391,"2013-8-4")
	dict_aa = dict_append_proc(dict_aa,"t2386","常滑",95187,"2013-1-21")
	dict_aa = dict_append_proc(dict_aa,"t2387","大府",31246,"2013-7-23")
	dict_aa = dict_append_proc(dict_aa,"t2388","瀬戸",85791,"2013-11-26")
	dict_aa = dict_append_proc(dict_aa,"t2389","犬山",24139,"2013-10-15")
#

	return	dict_aa
end
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_out = ARGS[1]
println(fname_out)
#
dict_aa = data_prepare_proc()
# dict_display_proc(dict_aa)
dict_to_text_proc(dict_aa,fname_out)
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
#####
##### Training hyperparameters
#####

Network = NetLib.ResNet

netparams = NetLib.ResNetHP(
  num_filters=128,
  num_blocks=5,
  conv_kernel_size=(3, 3),
  num_policy_head_filters=32,
  num_value_head_filters=32,
  batch_norm_momentum=0.1)

self_play = SelfPlayParams(
  sim=SimParams(
    num_games=1024,
    num_workers=32,
    batch_size=32,
    use_gpu=true,
    reset_every=2,
    flip_probability=0.,
    alternate_colors=false),
  mcts=MctsParams(
    num_iters_per_turn=200,
    cpuct=2.0,
    prior_temperature=1.0,
    temperature=PLSchedule([0, 20, 30], [1.0, 1.0, 0.3]),
    dirichlet_noise_ϵ=0.25,
    dirichlet_noise_α=0.3))

arena = ArenaParams(
  sim=SimParams(
    num_games=128,
    num_workers=32,
    batch_size=32,
    use_gpu=true,
    reset_every=2,
    flip_probability=0.0,
    alternate_colors=true),
  mcts=MctsParams(
    self_play.mcts,
    temperature=ConstSchedule(0.2),
    dirichlet_noise_ϵ=0.05),
  update_threshold=0.05)

learning = LearningParams(
  use_gpu=true,
  use_position_averaging=true,
  samples_weighing_policy=LOG_WEIGHT,
  batch_size=128,
  loss_computation_batch_size=128,
  optimiser=Adam(lr=2e-3),
  l2_regularization=1e-4,
  nonvalidity_penalty=1.,
  min_checkpoints_per_epoch=1,
  max_batches_per_checkpoint=2000,
  num_checkpoints=1)

params = Params(
  arena=arena,
  self_play=self_play,
  learning=learning,
  num_iters=45,
  ternary_rewards=true,
  use_symmetries=false,
  memory_analysis=nothing,
  mem_buffer_size=PLSchedule(
  [      0,        30],
  [50000, 180000]))

#####
##### Evaluation benchmark
#####

mcts_baseline =
  Benchmark.MctsRollouts(
    MctsParams(
      arena.mcts,
      num_iters_per_turn=1000,
      cpuct=1.))

minmax_baseline = Benchmark.MinMaxTS(
  depth=5,
  τ=0.2,
  amplify_rewards=true)

alphazero_player = Benchmark.Full(arena.mcts)

network_player = Benchmark.NetworkOnly(τ=0.5)

benchmark_sim = SimParams(
  arena.sim;
  num_games=128,
  num_workers=64,
  batch_size=64,
  alternate_colors=true)

benchmark = [
  Benchmark.Duel(alphazero_player, mcts_baseline,   benchmark_sim),
# Benchmark.Duel(alphazero_player, minmax_baseline, benchmark_sim),
  Benchmark.Duel(network_player,   mcts_baseline,   benchmark_sim),
# Benchmark.Duel(network_player,   minmax_baseline, benchmark_sim)
]

#####
##### Wrapping up in an experiment
#####

experiment = Experiment("amazon_full",
  GameSpec(), params, Network, netparams, benchmark)


text:
@doc raw"""
    GammaRationalKernel([α [,β [,γ]]])

The gamma-rational kernel is a generalization of the rational-quadratic kernel with an
additional shape parameter:

```math
\kappa(\mathbf{x},\mathbf{y})
= \left(1 +\alpha ||\mathbf{x},\mathbf{y}||^{\gamma}\right)^{-\beta}
\qquad \alpha > 0, \; \beta > 0, \; 0 < \gamma \leq 1
```

where ``\alpha`` is a scaling parameter and ``\beta`` and ``\gamma`` are shape parameters.
"""
struct GammaRationalKernel{T<:AbstractFloat} <: MercerKernel{T}
    alpha::HyperParameter{T}
    beta::HyperParameter{T}
    gamma::HyperParameter{T}
    GammaRationalKernel{T}(α::Real, β::Real, γ::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
        HyperParameter(convert(T,β), interval(OpenBound(zero(T)), nothing)),
        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
    )
end
function GammaRationalKernel(
        α::T1 = 1.0,
        β::T2 = one(T1),
        γ::T3 = one(floattype(T1,T2))
    ) where {T1<:Real,T2<:Real,T3<:Real}
    GammaRationalKernel{floattype(T1,T2,T3)}(α,β,γ)
end

@inline gammarationalkernel(z::T, α::T, β::T, γ::T) where {T<:AbstractFloat} = (1 + α*(z^γ))^(-β)

@inline basefunction(::GammaRationalKernel) = SquaredEuclidean()
@inline function kappa(κ::GammaRationalKernel{T}, z::T) where {T}
    gammarationalkernel(z, getvalue(κ.alpha), getvalue(κ.beta), getvalue(κ.gamma))
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7]);F = [484*x1^4+2420*x1^2*x2*x4+1496*x1^2*x2*x7-2288*x1^2*x3*x5-44*x1^2*x5*x7-3828*x1^2*x6*x7+3025*x2^2*x4^2+3740*x2^2*x4*x7+9025*x2^2*x6^2+12350*x2^2*x6*x7+5381*x2^2*x7^2-5720*x2*x3*x4*x5-3536*x2*x3*x5*x7-1900*x2*x3*x6^2-1300*x2*x3*x6*x7-110*x2*x4*x5*x7-7480*x2*x4*x6*x7+1430*x2*x4*x7^2-10640*x2*x5*x6*x7-7348*x2*x5*x7^2-5916*x2*x6*x7^2+2704*x3^2*x5^2+100*x3^2*x6^2-220*x3*x4*x6*x7+104*x3*x5^2*x7+10168*x3*x5*x6*x7+121*x4^2*x7^2-1232*x4*x5*x7^2+3137*x5^2*x7^2+174*x5*x6*x7^2+7569*x6^2*x7^2+14060*x2*x6+9620*x2*x7-1480*x3*x6+1628*x4*x7-8288*x5*x7+5476, 2420*x1^2*x4+1496*x1^2*x7+6050*x2*x4^2+7480*x2*x4*x7+18050*x2*x6^2+24700*x2*x6*x7+10762*x2*x7^2-5720*x3*x4*x5-3536*x3*x5*x7-1900*x3*x6^2-1300*x3*x6*x7-110*x4*x5*x7-7480*x4*x6*x7+1430*x4*x7^2-10640*x5*x6*x7-7348*x5*x7^2-5916*x6*x7^2+14060*x6+9620*x7, -2288*x1^2*x5-5720*x2*x4*x5-3536*x2*x5*x7-1900*x2*x6^2-1300*x2*x6*x7+5408*x3*x5^2+200*x3*x6^2-220*x4*x6*x7+104*x5^2*x7+10168*x5*x6*x7-1480*x6, 2420*x1^2*x2+6050*x2^2*x4+3740*x2^2*x7-5720*x2*x3*x5-110*x2*x5*x7-7480*x2*x6*x7+1430*x2*x7^2-220*x3*x6*x7+242*x4*x7^2-1232*x5*x7^2+1628*x7, -2288*x1^2*x3-44*x1^2*x7-5720*x2*x3*x4-3536*x2*x3*x7-110*x2*x4*x7-10640*x2*x6*x7-7348*x2*x7^2+5408*x3^2*x5+208*x3*x5*x7+10168*x3*x6*x7-1232*x4*x7^2+6274*x5*x7^2+174*x6*x7^2-8288*x7, -3828*x1^2*x7+18050*x2^2*x6+12350*x2^2*x7-3800*x2*x3*x6-1300*x2*x3*x7-7480*x2*x4*x7-10640*x2*x5*x7-5916*x2*x7^2+200*x3^2*x6-220*x3*x4*x7+10168*x3*x5*x7+174*x5*x7^2+15138*x6*x7^2+14060*x2-1480*x3, 1496*x1^2*x2-44*x1^2*x5-3828*x1^2*x6+3740*x2^2*x4+12350*x2^2*x6+10762*x2^2*x7-3536*x2*x3*x5-1300*x2*x3*x6-110*x2*x4*x5-7480*x2*x4*x6+2860*x2*x4*x7-10640*x2*x5*x6-14696*x2*x5*x7-11832*x2*x6*x7-220*x3*x4*x6+104*x3*x5^2+10168*x3*x5*x6+242*x4^2*x7-2464*x4*x5*x7+6274*x5^2*x7+348*x5*x6*x7+15138*x6^2*x7+9620*x2+1628*x4-8288*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _xyz{T<:AbstractFloat}(jd::T, equinox::T)
    t = (float(jd) - 15020.0)*inv(JULIANYEAR*100) # Reduced Julian century from 1900

    # NOTE: longitude arguments below are given in *equinox* of date.  Precess
    # these to equinox 1950 to give everything an even footing.  Compute
    # argument of precession from equinox of date back to 1950
    pp = (1.396041 + 0.000308*(t + 0.5))*(t-0.499998)
    # Compute mean solar longitude, precessed back to 1950
    el = 279.696678 + 36000.76892*t + 0.000303*t*t - pp

    # Compute Mean longitude of the Moon
    c = 270.434164 + 480960.0*t + 307.883142*t - 0.001133*t*t - pp

    # Compute longitude of Moon's ascending node
    n = 259.183275 - 1800.0*t - 134.142008*t + 0.002078*t*t - pp

    # Compute mean solar anomaly
    g = 358.475833 + 35999.04975*t - 0.00015*t*t

    # Compute the mean jupiter anomaly
    j = 225.444651 + 2880.0*t + 154.906654*t*t

    # Compute mean anomaly of Venus
    v = 212.603219 + 58320.0*t + 197.803875*t + 0.001286*t*t

    # Compute mean anomaly of Mars
    m = 319.529425 + 19080.0*t + 59.8585*t + 0.000181*t*t

    # Convert degrees to radians for trig functions
    el = deg2rad(el)
    g  = deg2rad(g)
    j  = deg2rad(j)
    c  = deg2rad(c)
    v  = deg2rad(v)
    n  = deg2rad(n)
    m  = deg2rad(m)

    # Calculate x, y, z using trigonometric series
    x = 0.999860*cos(el) -
        0.025127*cos(g - el) +
        0.008374*cos(g + el) +
        0.000105*cos(g + g + el) +
        0.000063*t*cos(g - el) +
        0.000035*cos(g + g - el) -
        0.000026*sin(g - el - j) -
        0.000021*t*cos(g + el) +
        0.000018*sin(2.0*g + el - 2.0*v) +
        0.000017*cos(c) -
        0.000014*cos(c - 2.0*el) +
        0.000012*cos(4.0*g + el - 8.0*m + 3.0*j) -
        0.000012*cos(4.0*g - el - 8.0*m + 3.0*j) -
        0.000012*cos(g + el - v) +
        0.000011*cos(2.0*g + el - 2.0*v) +
        0.000011*cos(2.0*g - el - 2.0*j)

    y = 0.917308*sin(el) +
        0.023053*sin(g - el) +
        0.007683*sin(g + el) +
        0.000097*sin(g + g + el) -
        0.000057*t*sin(g - el) -
        0.000032*sin(g + g - el) -
        0.000024*cos(g - el - j) -
        0.000019*t*sin(g + el) -
        0.000017*cos(2.0*g + el - 2.0*v) +
        0.000016*sin(c) +
        0.000013*sin(c - 2.0*el ) +
        0.000011*sin(4.0*g + el - 8.0*m + 3.0*j) +
        0.000011*sin(4.0*g - el - 8.0*m + 3.0*j) -
        0.000011*sin(g + el - v) +
        0.000010*sin(2.0*g + el - 2.0*v ) -
        0.000010*sin(2.0*g - el - 2.0*j )

    z = 0.397825*sin(el)     +
        0.009998*sin(g-el)   +
        0.003332*sin(g+el)   +
        0.000042*sin(g+g+el) -
        0.000025*t*sin(g-el) -
        0.000014*sin(g+g-el) -
        0.000010*cos(g-el-j)

    # Precess to new equator?  Avoid useless calculations.
    if isfinite(equinox) && equinox != 1950
        x, y, z = precess_xyz(x, y, z, 1950.0, equinox)
    end

    # Velocities
    xvel = -0.017200 * sin(el)   -
        0.000288*sin(g + el)     -
        0.000005*sin(2.0*g + el) -
        0.000004*sin(c)          +
        0.000003*sin(c - 2.0*el) +
        0.000001*t*sin(g + el)   -
        0.000001*sin(2.0*g - el)

    yvel = 0.015780*cos(el)       +
        0.000264*cos(g + el)      +
        0.000005*cos(2.0*g + el)  +
        0.000004*cos(c)           +
        0.000003*cos(c - 2.0*el)  -
        0.000001*t*cos(g + el)

    zvel = 0.006843*cos(el)      +
        0.000115*cos(g  + el)    +
        0.000002*cos(2.0*g + el) +
        0.000002*cos(c)          +
        0.000001*cos(c - 2.0*el)

    # Precess to new equator?  Avoid useless calculations.
    if isfinite(equinox) && equinox != 1950
        xvel, yvel, zvel = precess_xyz(xvel, yvel, zvel, 1950, equinox)
    end

    return x, y, z, xvel, yvel, zvel
end

"""
    xyz(jd[, equinox]) -> x, y, z, v_x, v_y, v_z

### Purpose ###

Calculate geocentric \$x\$, \$y\$, and \$z\$ and velocity coordinates of the
Sun.

### Explanation ###

Calculates geocentric \$x\$, \$y\$, and \$z\$ vectors and velocity coordinates
(\$dx\$, \$dy\$ and \$dz\$) of the Sun.  (The positive \$x\$ axis is directed
towards the equinox, the \$y\$-axis, towards the point on the equator at right
ascension 6h, and the \$z\$ axis toward the north pole of the equator).  Typical
position accuracy is \$<10^{-4}\$ AU (15000 km).

### Arguments ###

* `jd`: number of Reduced Julian Days for the wanted date.  It can be either a
  scalar or a vector.
* `equinox` (optional numeric argument): equinox of output. Default is 1950.

You can use `juldate` to get the number of Reduced Julian Days for the selected
dates.

### Output ###

The 6-tuple \$(x, y, z, v_x, v_y, v_z)\$, where

* \$x\$, \$y\$, \$z\$: scalars or vectors giving heliocentric rectangular
  coordinates (in AU) for each date supplied.  Note that \$\\sqrt{x^2 + y^2 + z^2}\$
  gives the Earth-Sun distance for the given date.
* \$v_x\$, \$v_y\$, \$v_z\$: velocity vectors corresponding to \$x\$, \$y\$, and
  \$z\$.

### Example ###

What were the rectangular coordinates and velocities of the Sun on
1999-01-22T00:00:00 (= JD 2451200.5) in J2000 coords?  Note: Astronomical
Almanac (AA) is in TDT, so add 64 seconds to UT to convert.

``` julia
jd = juldate(DateTime(1999, 1, 22))
# => 51200.5
xyz(jd + 64./86400., 2000)
# => (0.5145687092402946,-0.7696326261820777,-0.33376880143026394,0.014947267514081075,0.008314838205475709,0.003606857607574784)
```

Compare to Astronomical Almanac (1999 page C20)

                x  (AU)        y  (AU)     z (AU)
    xyz:      0.51456871   -0.76963263  -0.33376880
    AA:       0.51453130   -0.7697110   -0.3337152
    abs(err): 0.00003739    0.00007839   0.00005360
    abs(err)
        (km):   5609          11759         8040

NOTE: Velocities in AA are for Earth/Moon barycenter
      (a very minor offset) see AA 1999 page E3

               x vel (AU/day) y vel (AU/day)   z vel (AU/day)
    xyz:      -0.014947268   -0.0083148382    -0.0036068576
    AA:       -0.01494574    -0.00831185      -0.00360365
    abs(err):  0.000001583    0.0000029886     0.0000032076
    abs(err)
     (km/sec): 0.00265        0.00519          0.00557

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
xyz(jd::Real, equinox::Real=NaN) = _xyz(promote(float(jd), float(equinox))...)

# Can't use @vectorize_1arg because of the optional keyword.
function xyz{J<:Real}(jd::AbstractArray{J}, equinox::Real=NaN)
    typej = typeof(float(one(J)))
    x     = similar(jd, typej)
    y     = similar(jd, typej)
    z     = similar(jd, typej)
    xvel  = similar(jd, typej)
    yvel  = similar(jd, typej)
    zvel  = similar(jd, typej)
    for i in eachindex(jd)
        x[i], y[i], z[i], xvel[i], yvel[i], zvel[i] = xyz(jd[i], equinox)
    end
    return x, y, z, xvel, yvel, zvel
end


text:
module PotentialMap

using FindSteadyStates
using SteadyStateDiffEq
using LinearAlgebra
using ModelingToolkit
using ForwardDiff
using Optim
using Parameters
using Reexport

include("LeastAction/LeastAction.jl")

@reexport using .LeastAction

end