text:
using GPUInspector
using Test
using CUDA
using InteractiveUtils: subtypes
using LinearAlgebra

if !GPUInspector.functional()
    error("Can't run testsuite since CUDA/GPU not present or not functional!")
end
if Threads.nthreads() == 1
    error("Can't run testsuite in single-threaded mode. Start Julia with multiple threads.")
end
if Threads.nthreads() < length(CUDA.devices()) + 1
    # we should have at least one thread per gpu + one monitoring thread
    @warn(
        "You should run the tests with at least $(length(CUDA.devices()) + 1) Julia threads.",
        Threads.nthreads(),
        length(CUDA.devices())
    )
end

@testset "GPUInspector.jl" begin
    @testset "UnitPrefixedBytes" begin
        # general stuff
        b = B(40_000_000)
        m = MiB(38.14697265625)
        m̃ = MB(40.0)
        @test typeof(b) == B
        @test typeof(m) == MiB
        @test typeof(m̃) == MB
        @test string(b) == "4.0e7 B"
        @test string(m) == "~38.15 MiB"
        @test string(m̃) == "40.0 MB"
        @test value(b) == 40_000_000
        @test value(m) ≈ 38.14697265625
        @test value(m̃) ≈ 40.0
        @test simplify(b) ≈ m
        @test simplify(m) ≈ m
        @test simplify(m̃) ≈ m
        @test simplify(b; base=10) ≈ m̃
        @test simplify(m; base=10) ≈ m̃
        @test simplify(m̃; base=10) ≈ m̃
        @test change_base(b) == b
        @test change_base(m) == m̃
        @test change_base(m̃) == m
        @test convert(MB, b) == m̃
        @test convert(MiB, b) == m
        @test convert(MB, m) == m̃
        @test convert(MiB, m̃) == m
        @test bytes(b) == 4.0e7
        @test bytes(m) == 4.0e7
        @test bytes(m̃) == 4.0e7
        @test bytes(4.0e7) ≈ m

        # GiB <-> GB conversion because it's particularly important
        @test convert(GiB, GB(300)) ≈ GiB(279.39677238464355)
        @test change_base(GB(300)) ≈ GiB(279.39677238464355)
        @test convert(GB, GiB(300)) ≈ GB(322.12254720000004)
        @test change_base(GiB(300)) ≈ GB(322.12254720000004)

        # ≈, ==, ===
        @test b == m
        @test b == m̃
        @test m == m̃
        @test b ≈ m
        @test b ≈ m̃
        @test m ≈ m̃
        @test b === b
        @test b !== m
        @test b !== m̃
        @test m !== m̃

        # basic arithmetics
        types = subtypes(UnitPrefixedBytes)
        for T in types
            @test T(1.2) + T(2.3) ≈ T(3.5)
            @test T(3.5) - T(2.3) ≈ T(1.2)
            @test 2 * T(1.234) ≈ T(2.468)
            @test T(10) / 2 == T(5)
        end
        for T in types, S in types
            @test bytes(T(1.2) + S(2.3)) ≈ bytes(T(1.2)) + bytes(S(2.3))
            @test abs(bytes(T(1.2) - S(2.3))) ≈ abs(bytes(T(1.2)) - bytes(S(2.3)))
        end
        @test B(40_000_000) + MB(3) - 2 * KiB(2) ≈ MB(42.995904)
    end

    @testset "toggle_tensorcoremath" begin
        @test isnothing(toggle_tensorcoremath(true; verbose=false))
        @test CUDA.math_mode() == CUDA.FAST_MATH
        @test isnothing(toggle_tensorcoremath(false; verbose=false))
        @test CUDA.math_mode() == CUDA.DEFAULT_MATH
        # test toggle
        @test isnothing(toggle_tensorcoremath(; verbose=false))
        @test CUDA.math_mode() == CUDA.FAST_MATH
        @test isnothing(toggle_tensorcoremath(; verbose=false))
        @test CUDA.math_mode() == CUDA.DEFAULT_MATH
    end

    @testset "gpu stresstest" begin
        @testset "different kinds" begin
            for dev in (device(), devices()) # single- and multi-gpu
                @test isnothing(stresstest(dev; duration=2, verbose=false))
                @test isnothing(stresstest(dev; enforced_duration=2, verbose=false))
                @test isnothing(stresstest(dev; approx_duration=2, verbose=false))
                @test isnothing(stresstest(dev; niter=100, verbose=false))
                @test isnothing(stresstest(dev; niter=100, verbose=false))
                @test isnothing(stresstest(dev; mem=0.2, verbose=false))
            end
        end

        @testset "keyword options" begin
            dev = device()
            @test isnothing(stresstest(dev; duration=2, verbose=false))
            @test isnothing(stresstest(dev; duration=2, size=3000, verbose=false))
            @test isnothing(stresstest(dev; duration=2, dtype=Float16, verbose=false))
            @test isnothing(stresstest(dev; duration=2, clearmem=true, verbose=false))
            @test isnothing(stresstest(dev; duration=2, batch_duration=1, verbose=false))
            # TODO: kwargs: threads, parallel
        end

        @testset "monitoring" begin
            @testset "automatic (monitoring=true)" begin
                for dev in (device(), devices()) # single- and multi-gpu
                    @test typeof(
                        stresstest(dev; duration=2, verbose=false, monitoring=true)
                    ) == MonitoringResults
                end
            end
            @testset "manual" begin
                for devs in ([device()], devices()) # single- and multi-gpu
                    @test isnothing(monitoring_start(; freq=1, devices=devs, verbose=false))
                    @test isnothing(
                        stresstest(devs; duration=2, verbose=false, monitoring=false)
                    )
                    @test typeof(monitoring_stop(; verbose=false)) == MonitoringResults
                end
                # thread kwarg
                # @test isnothing(monitoring_start(; freq=1, devices=[device()], verbose=false, thread=2))
                # sleep(0.5)
                # TODO: How to find the threadid() of a Task?
                # @test isnothing(monitoring_stop(; verbose=false))
            end
        end

        @testset "monitoring results" begin
            @testset "MonitoringResults" begin
                r = stresstest(device(); duration=2, verbose=false, monitoring=true)
                @test typeof(r) == MonitoringResults
                @test typeof(r.times) == Vector{Float64}
                @test typeof(r.devices) == Vector{Tuple{String,Base.UUID}}
                quant = first(keys(r.results))
                @test r.results[quant] == getproperty(r, quant)
            end
            @testset "save / load" begin
                d = Dict{Symbol,Vector{Vector{Float64}}}()
                ndevs = length(CUDA.devices())
                d[:asd] = [rand(ndevs) for _ in 1:5]
                d[:qwe] = [rand(ndevs) for _ in 1:5]
                d[:jkl] = [rand(ndevs) for _ in 1:5]
                devices = Tuple{String,Base.UUID}[
                    (GPUInspector._device2string(dev), uuid(dev)) for
                    dev in collect(CUDA.devices())
                ]
                r = MonitoringResults(rand(5), devices, d)
                cd(mktempdir()) do
                    save_monitoring_results("tmp.h5", r)
                    r2 = load_monitoring_results("tmp.h5")
                    @test r == r2
                end
            end
            @testset "plot" begin
                r = load_monitoring_results(joinpath(@__DIR__, "test.h5"))
                @test isnothing(plot_monitoring_results(r))
                @test isnothing(plot_monitoring_results(r, (:compute, :mem)))
            end

            @testset "savefig" begin
                using CairoMakie
                r = load_monitoring_results(joinpath(@__DIR__, "test.h5"))
                @test isnothing(savefig_monitoring_results(r))
                @test isnothing(
                    savefig_monitoring_results(r, (:compute, :mem))
                )
                @test isnothing(savefig_monitoring_results(r; ext=:png))
                @test isnothing(savefig_monitoring_results(r; ext=:pdf))
                rm.(filter(endswith(".pdf"), readdir())) # cleanup
            end
        end
    end

    @testset "gpuinfo / gpus" begin
        @test isnothing(gpus())
        @test isnothing(gpuinfo())
        @test isnothing(gpuinfo(0))
        @test isnothing(gpuinfo(device()))
    end

    @testset "p2p_bandwidth" begin
        @testset "unidirectional" begin
            # p2p_bandwidth
            @test typeof(p2p_bandwidth(; verbose=false)) == Float64
            @test 0 ≤ p2p_bandwidth(; verbose=false)
            # options
            @test typeof(p2p_bandwidth(MB(100); verbose=false)) == Float64
            @test typeof(
                p2p_bandwidth(; src=CuDevice(0), dst=CuDevice(1), verbose=false)
            ) == Float64
            @test typeof(p2p_bandwidth(; dtype=Float16, verbose=false)) == Float64
            @test typeof(p2p_bandwidth(; nbench=10, verbose=false)) == Float64
            @test typeof(p2p_bandwidth(; hist=true, verbose=true)) == Float64
            # p2p_bandwidth_all
            @test typeof(p2p_bandwidth_all(; verbose=false)) ==
                Matrix{Union{Nothing,Float64}}
            Mp2p = p2p_bandwidth_all(; verbose=false)
            @test all(isnothing, diag(Mp2p))
            @test all(
                !isnothing(Mp2p[i, j]) for i in axes(Mp2p, 1), j in axes(Mp2p, 2) if i != j
            )
        end
        @testset "bidirectional" begin
            # p2p_bandwidth_bidirectional
            @test typeof(p2p_bandwidth_bidirectional(; verbose=false)) == Float64
            @test 0 ≤ p2p_bandwidth_bidirectional(; verbose=false)
            # options
            @test typeof(p2p_bandwidth_bidirectional(MB(100); verbose=false)) == Float64
            @test typeof(p2p_bandwidth_bidirectional(; dtype=Float16, verbose=false)) ==
                Float64
            @test typeof(p2p_bandwidth_bidirectional(; nbench=10, verbose=false)) ==
                Float64
            @test typeof(p2p_bandwidth_bidirectional(; hist=true, verbose=true)) ==
                Float64
            # p2p_bandwidth_bidirectional_all
            @test typeof(p2p_bandwidth_bidirectional_all(; verbose=false)) ==
                Matrix{Union{Nothing,Float64}}
            Mp2p = p2p_bandwidth_bidirectional_all(; verbose=false)
            @test all(isnothing, diag(Mp2p))
            @test all(
                !isnothing(Mp2p[i, j]) for i in axes(Mp2p, 1), j in axes(Mp2p, 2) if i != j
            )
        end
    end

    @testset "host2device_bandwidth" begin
        @test isnothing(host2device_bandwidth())
        @test isnothing(host2device_bandwidth(MB(100)))
        @test isnothing(host2device_bandwidth(; dtype=Float16))
    end

    @testset "memory_bandwidth" begin
        @test typeof(memory_bandwidth()) == Float64
        @test typeof(memory_bandwidth(MiB(10))) == Float64
        @test typeof(memory_bandwidth(; dtype=Float32)) == Float64

        @test typeof(memory_bandwidth_saxpy()) == Float64
        @test typeof(memory_bandwidth_saxpy(; size=2^20*2)) == Float64
        @test typeof(memory_bandwidth_saxpy(; dtype=Float32)) == Float64
    end

    @testset "peakflops_gpu (CUDA cores)" begin
        @test typeof(peakflops_gpu(; verbose=false, tensorcores=false)) == Float64
        @test typeof(peakflops_gpu(; dtype=Float32, verbose=false, tensorcores=false)) == Float64
        @test typeof(peakflops_gpu(; dtype=Float64, verbose=false, tensorcores=false)) == Float64
    end

    @testset "peakflops_gpu (Tensor cores)" begin
        @test typeof(peakflops_gpu(; verbose=false, tensorcores=true)) == Float64
        @test typeof(peakflops_gpu(; dtype=Float16, verbose=false, tensorcores=true)) == Float64
    end

    @testset "peakflops_gpu_matmul / scaling" begin
        @test typeof(peakflops_gpu_matmul(; verbose=false)) == Float64
        @test typeof(peakflops_gpu_matmul(; size=1024, dtype=Float64, verbose=false)) == Float64
        @test typeof(peakflops_gpu_matmul(; nmatmuls=2, nbench=2, verbose=false)) == Float64
        @test typeof(peakflops_gpu_matmul_scaling(; verbose=false)) == Tuple{Vector{Int64}, Vector{Float64}}
    end
end


text:
using DataStructures

"""Node of the state tree to explore

It has to be mutable because nodes can be updated in the open set heap, and to do so we need to keep a handle of the node in the heap"""
mutable struct Node{TState, TCost <: Number}
  data::TState
  depth::Int32
  g::TCost
  h::TCost
  f::TCost
  parent::Union{Node{TState, TCost}, Nothing}
  heaphandle::Union{Int64,Nothing}
end

"Results structure"
struct AStarResult{TState, TCost <: Number}
  status::Symbol
  path::Array{TState,1}
  cost::TCost
  closedsetsize::Int64
  opensetsize::Int64
end

"order by f = g + h"
nodeorderingkey(n::Node) = n.f

"By default every transition from a state to each neighbour costs 1"
defaultcost(s1, s2) = one(Int64)

"By default, the herustic returns 0 (Breadth First Search)"
defaultheuristic(state, goal) = zero(Int64)

defaultisgoal(state, goal) = state == goal

"reconstruct the path of states up to the found final node"
function reconstructpath(n::Node)
  res = [n.data]
  while !isnothing(n.parent)
    n = n.parent
    push!(res, n.data)
  end
  return reverse(res)
end

"""
  astar(neighbours, start, goal;
        heuristic=defaultheuristic, cost=defaultcost, isgoal=defaultisgoal, hashfn=hash, timeout=Inf, maxcost=Inf, maxdepth=Inf)

Execute the A* algorithm to get the best path from the start state to reach a goal condition.
Only the first 3 arguments are mandatory, all the others are optional.

It returns a structure in which the `status` field is a Symbol that can be either:
- `:success`: the algorithm found a path from start to goal
- `:timeout`: the algorithm timed out, a partial path to the best state is returned in the `path` field
- `:nopath`: the algorithm didn't find any path to a goal, the path to the best state is still returned

The other fields are:
- `path`: an array of states from the start state to the goal or the best found state
- `cost`: the cost of the returned path
- `closedsetsize`: how many states the algorithm tested if they were a goal (size of the closed set)
- `opensetsize`: how many states were still in the open set when the algorithm ended

# Arguments
- `neighbours`: a function that takes a state and returns the neighbour states as an array (or iterable)
- `start`: the starting state, the type of the state is completely unrestricted
- `goal`: the goal state, the type is unrestricted, usually it's the same as the start
- `heuristic`: a function that given a state and the goal returns an estimate of the cost to reach goal. This estimate should be optimistic if you want to be sure to get the best path. Notice that the best path could be very expensive to find, so if you want a good but not guaranteed optimal path, you could multiply your heuristic by a constant, the algorithm will usually be much faster
- `cost`: a function that takes the current state and a neighbour and returns the cost to do that state transition. By default all transitions cost 1
- `isgoal`: a function that takes a state and the goal and evaluates if the goal is reached (by default ==)
- `hashfn`: a function that takes a state and returns a compact representation to use as dictionary key (usually one of UInt, Int, String), by default it is the `hash` function. This is a very important field for composite states in order to avoid duplications.
- `timeout`: timeout in number of seconds after which the algorithm stops returning the best partial path to the state with the lowest heuristic, by default it is unrestricted. Please notice that the algorithm wil run _AT LEAST_ the specified time.
- `maxcost`: a maximum bound of the accumulated cost of the path, this can result in a :nopath result even if a path to the goal (with a greater cost) exists. By default it is Inf
- `maxdepth`: the maximum depth the algorithm is allowed to go down while expanding the search state, the same considerations as the `maxcost` parameter apply. By default it is Inf
"""
function astar(neighbours, start, goal;
               heuristic=defaultheuristic, cost=defaultcost, isgoal=defaultisgoal, hashfn=hash, timeout=Inf, maxcost=Inf, maxdepth=Inf)
  starttime = time()
  startheuristic = heuristic(start, goal)
  startcost = zero(cost(start, start))
  nodetype = typeof(start)
  costtype = typeof(startcost)
  startnode = Node{nodetype, costtype}(start, zero(Int32), startcost, startheuristic, startheuristic, nothing, nothing)
  bestnode = startnode
  starthash = hashfn(start)

  closedset = Set{typeof(starthash)}()
  openheap = MutableBinaryHeap(Base.By(nodeorderingkey), Node{nodetype, costtype}[])
  starthandle = push!(openheap, startnode)
  startnode.heaphandle = starthandle
  opennodedict = Dict(starthash=>startnode)

  while !isempty(openheap)
    if time() - starttime > timeout
      return AStarResult{nodetype, costtype}(:timeout, reconstructpath(bestnode), bestnode.g, length(closedset), length(openheap))
    end

    node = pop!(openheap)
    nodehash = hashfn(node.data)
    delete!(opennodedict, nodehash)

    if isgoal(node.data, goal)
      return AStarResult{nodetype, costtype}(:success, reconstructpath(node), node.g, length(closedset), length(openheap))
    end

    push!(closedset, nodehash)

    if node.h < bestnode.h
      bestnode = node
    end

    neighbour_states = neighbours(node.data)

    for neighbour in neighbour_states
      neighbourhash = hashfn(neighbour)
      if neighbourhash in closedset
        continue
      end

      gfromthisnode = node.g + cost(node.data, neighbour)

      if gfromthisnode > maxcost
        continue
      end

      if node.depth > maxdepth
        continue
      end

      if neighbourhash in keys(opennodedict)
        neighbournode = opennodedict[neighbourhash]
        if gfromthisnode < neighbournode.g
          neighbournode.depth = node.depth + one(Int32)
          neighbournode.g = gfromthisnode
          neighbournode.parent = node
          neighbournode.f = gfromthisnode + neighbournode.h
          update!(openheap, neighbournode.heaphandle, neighbournode)
        end
      else
        neighbourheuristic = heuristic(neighbour, goal)
        neighbournode = Node{nodetype, costtype}(neighbour, node.depth + one(Int32), gfromthisnode, neighbourheuristic, gfromthisnode + neighbourheuristic, node, nothing)
        neighbourhandle = push!(openheap, neighbournode)
        neighbournode.heaphandle = neighbourhandle
        push!(opennodedict, neighbourhash => neighbournode)
      end
    end
  end

  return AStarResult{nodetype, costtype}(:nopath, reconstructpath(bestnode), bestnode.g, length(closedset), length(openheap))
end

ext: jl


text:
"""
Internal implementation detail.
"""
module Common

include("bytes.jl")
include("errors.jl")
include("Utils.jl")
include("BufferParser.jl")
include("CDObject.jl")

end


text:
using DelayDiffEq
using Test

h(p, t) = 1.0
const tspan = (2.0, 0.0)

f(u, h, p,t) = h(p, t + 1)

# solution on [0,2]
f_analytic(u₀, ::typeof(h), p, t) = t < 1 ? (t^2 - 1)/2 : t - 1

const dde_f = DDEFunction(f, analytic = f_analytic)

@testset "Without lags" begin
    sol = solve(DDEProblem(dde_f, h, tspan), MethodOfSteps(RK4()))
    @test sol.errors[:l∞] < 1.2e-12 # 1.2e-16
end

@testset "Constant lags" begin
    @testset "incorrect" begin
        prob = DDEProblem(dde_f, h, tspan; constant_lags = [1.0])
        @test_throws ErrorException solve(prob, MethodOfSteps(Tsit5()))
    end

    @testset "correct" begin
        prob = DDEProblem(dde_f, h, tspan; constant_lags = [-1.0])

        dde_int = init(prob, MethodOfSteps(Tsit5()))
        @test dde_int.opts.d_discontinuities.valtree == [Discontinuity(-1.0, 2)]

        sol = solve!(dde_int)
        @test sol.errors[:l∞] < 3.9e-12 # 3.9e-15
        @test dde_int.tracked_discontinuities == [Discontinuity(-2.0, 1)]
    end
end

@testset "Dependent lags" begin
    prob = DDEProblem(dde_f, h, tspan; dependent_lags = ((u, p, t) -> -1.0,))

    dde_int = init(prob, MethodOfSteps(Tsit5()))
    @test isempty(dde_int.opts.d_discontinuities)

    sol = solve!(dde_int)
    @test sol.errors[:l∞] < 3.9e-12 # 3.9e-15
    @test dde_int.tracked_discontinuities == [Discontinuity(-2.0, 1)]
end

@testset "dt and dtmax" begin
    prob = DDEProblem(dde_f, h, tspan)

    dde_int = init(prob, MethodOfSteps(RK4()); dt = 0.1, dtmax = 0.5)
    @test dde_int.dt == -0.1
    @test dde_int.opts.dtmax == -0.5

    sol1 = solve!(dde_int)
    sol2 = solve(prob, MethodOfSteps(RK4()); dt = -0.1, dtmax = -0.5)
    @test sol1.t == sol2.t
    @test sol1.u == sol2.u
end


text:
using Submodular, Convex
using FactCheck

Tol = 1e-3

facts("Prox") do

  context("prox of associated polyhedra") do
    S = SetVariable(4)
    p(z) = -0.5*z^2 + 4*z + 0.5 * z
    perm_func = compose(p, card(S))
    P = SubmodPoly(perm_func)

    x = Variable(4)
    w₁ = [1, 2, 5, 3]
    w₂ = [-2, -2, -5, -3]
    @fact prox(x in P, w₁) --> roughly([1, 2, 4, 3], Tol)
    @fact prox(x in P, w₂) --> roughly([-2, -2, -5, -3], Tol)
  end

  context("prox of convex functions") do
    w = rand(4)
    x = Variable(4)
    f = 0.5 * norm(x)^2
    prox(f, w)
    @fact evaluate(0.5 * w - x) --> roughly(zeros(4, 1), Tol)
  end

end


text:
@testset "FiniteDiff comparison" begin
    q = [1,2,3,4.0]
    q = Quaternion(q ./ norm(q)...,true)
    @test norm(Dojo.dmrpdq(Dojo.vector(q)) -
        ForwardDiff.jacobian(Dojo.mrp, Dojo.vector(q)), Inf) < 1.0e-8
    @test norm(Dojo.daxisdq(Dojo.vector(q)) -
        ForwardDiff.jacobian(Dojo.axis, Dojo.vector(q))) < 1.0e-8
    @test norm(Dojo.drotation_vectordq(Dojo.vector(q)) -
        ForwardDiff.jacobian(Dojo.rotation_vector, Dojo.vector(q))) < 1.0e-8
    q = one(Quaternion{Float64})
    @test norm(Dojo.dmrpdq(Dojo.vector(q)) -
        ForwardDiff.jacobian(Dojo.mrp, Dojo.vector(q)), Inf) < 1.0e-8
    # TODO For zero rotation, ForwardDiff runs into a numeric singularity
    @test norm(Dojo.drotation_vectordq(Dojo.vector(q)) -
        FiniteDiff.finite_difference_jacobian(Dojo.rotation_vector, Dojo.vector(q))) < 1.0e-5
end

@testset "axes pair to quaternion" begin
    n1 = [0,0,0.0]
    n2 = [0,0,0.0]
    q = Dojo.axes_pair_to_quaternion(n1, n2)
    @test norm(n2 - Dojo.vector_rotate(n1, q), Inf) < 1e-8

    n1 = rand(3)
    n1 /= norm(n1)
    n2 = n1
    q = Dojo.axes_pair_to_quaternion(n1, n2)
    @test norm(n2 - Dojo.vector_rotate(n1, q), Inf) < 1e-8

    n1 = rand(3)
    n1 /= norm(n1)
    n2 = -n1
    q = Dojo.axes_pair_to_quaternion(n1, n2)
    @test norm(n2 - Dojo.vector_rotate(n1, q), Inf) < 1e-5

    n1 = rand(3)
    n1 /= norm(n1)
    n2 = rand(3)
    n2 /= norm(n2)
    q = Dojo.axes_pair_to_quaternion(n1, n2)
    @test norm(n2 - Dojo.vector_rotate(n1, q), Inf) < 1e-8
end


text:
function main()
    println("hello from Julia")
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function test_vKITTI_Conv(model,predict,w; batchsize=32,saveImages=true,
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testMode="even",
    experiment="First_Second"
    )
    @time begin
        printlnFlush("Loading weights from $(w)")
        w = loadNetwork(w)
        totalloss = totalcorrect = totalcount = 0
        for world in worlds
            for variation in testVariations
                x,y,colors = loadData_vKITTI([world],[variation],model,testMode)
                batchCount = Int64(ceil(size(y,4)/batchsize))
                varloss = imcounter = varcorrect = varcount = 0
                for batch=1:batchCount
                    printUpdate(@sprintf("Batch: %5d/%-5d",batch,batchCount))
                    #Calculate current indices
                    startIndex = 1+(batch-1)*batchsize
                    endIndex = min(batch*batchsize,size(y,4))
                    #Minibatch
                    batchx = KnetArray(x[:,:,:,startIndex:endIndex])
                    batchy = y[:,:,:,startIndex:endIndex]
                    #Accuracy
                    ypred = Array(predict(w,batchx))
                    varloss += pixelwiseSoftloss(ypred,batchy)
                    testAcc = pixelwiseAccuracy(ypred,batchy)
                    varcorrect += testAcc[1]
                    varcount += testAcc[2]
                    #Image conversion
                    if saveImages
                        for im=1:size(ypred,4)
                            savePrediction(ypred[:,:,:,im],colors,outputPath(@sprintf("vKITTI/%s/Conv/%s/%s/%s/%05d.png",experiment,model,world,variation,imcounter)))
                            imcounter += 1
                        end
                    end
                end
                printlnFlush(@sprintf("World: %s\tVariation: %s\tLoss: %0.10f\tAccuracy: %0.4f%%",world,variation,varloss/batchCount,varcorrect/varcount*100))
                totalloss += varloss/batchCount
                totalcorrect += varcorrect
                totalcount += varcount
            end
        end
        printlnFlush(@sprintf("Total Loss: %0.10f\tTotal Accuracy: %0.4f%%",totalloss,totalcorrect/totalcount*100))
    end
end

function test_vKITTI_RNN(model,predict,w; saveImages=true,
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testMode="even",
    experiment="First_Second"
    )
    @time begin
        printlnFlush("Loading weights from $(w)")
        @time w = loadNetwork(w)
        totalloss = totalcorrect = totalcount = 0
        for world in worlds
            for variation in testVariations
                varloss = imcounter = varcorrect = varcount = 0
                x,y,colors = loadData_vKITTI([world],[variation],model,testMode)
                batchCount = size(y,4)
                h = KnetArray(zeros(Float32,28,28,512,1))
                for t=1:batchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,batchCount))
                    #Accuracy
                    ypred = Array(predict(w,h,KnetArray(x[:,:,:,t:t])))
                    varloss += pixelwiseSoftloss(ypred,y[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,y[:,:,:,t:t])
                    varcorrect += testAcc[1]
                    varcount += testAcc[2]
                    #Image conversion
                    if saveImages
                        savePrediction(ypred[:,:,:,1],colors,outputPath(@sprintf("vKITTI/%s/RNN/%s/%s/%s/%05d.png",experiment,model,world,variation,imcounter)))
                        imcounter += 1
                    end
                end
                printlnFlush(@sprintf("World: %s\tVariation: %s\tLoss: %0.10f\tAccuracy: %0.4f%%",world,variation,varloss/batchCount,varcorrect/varcount*100))
                totalloss += varloss/batchCount
                totalcorrect += varcorrect
                totalcount += varcount
            end
        end
        printlnFlush(@sprintf("Total Loss: %0.10f\tTotal Accuracy: %0.4f%%",totalloss,totalcorrect/totalcount*100))
    end
end

function test_vKITTI_LSTM(model,predict,w; saveImages=true,
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testMode="even",
    experiment="First_Second"
    )
    @time begin
        printlnFlush("Loading weights from $(w)")
        @time w = loadNetwork(w)
        totalloss = totalcorrect = totalcount = 0
        for world in worlds
            for variation in testVariations
                varloss = imcounter = varcorrect = varcount = 0
                x,y,colors = loadData_vKITTI([world],[variation],model,testMode)
                batchCount = size(y,4)
                h = KnetArray(zeros(Float32,28,28,512,1))
                c = KnetArray(zeros(Float32,28,28,512,1))
                for t=1:batchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,batchCount))
                    #Accuracy
                    ypred = Array(predict(w,h,c,KnetArray(x[:,:,:,t:t])))
                    varloss += pixelwiseSoftloss(ypred,y[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,y[:,:,:,t:t])
                    varcorrect += testAcc[1]
                    varcount += testAcc[2]
                    #Image conversion
                    if saveImages
                        savePrediction(ypred[:,:,:,1],colors,outputPath(@sprintf("vKITTI/%s/LSTM/%s/%s/%s/%05d.png",experiment,model,world,variation,imcounter)))
                        imcounter += 1
                    end
                end
                printlnFlush(@sprintf("World: %s\tVariation: %s\tLoss: %0.10f\tAccuracy: %0.4f%%",world,variation,varloss/batchCount,varcorrect/varcount*100))
                totalloss += varloss/batchCount
                totalcorrect += varcorrect
                totalcount += varcount
            end
        end
        printlnFlush(@sprintf("Total Loss: %0.10f\tTotal Accuracy: %0.4f%%",totalloss,totalcorrect/totalcount*100))
    end
end

function test_vKITTI_GRU(model,predict,w; saveImages=true,
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testMode="even",
    experiment="First_Second"
    )
    @time begin
        printlnFlush("Loading weights from $(w)")
        @time w = loadNetwork(w)
        totalloss = totalcorrect = totalcount = 0
        for world in worlds
            for variation in testVariations
                varloss = imcounter = varcorrect = varcount = 0
                x,y,colors = loadData_vKITTI([world],[variation],model,testMode)
                batchCount = size(y,4)
                h = KnetArray(zeros(Float32,28,28,512,1))
                for t=1:batchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,batchCount))
                    #Accuracy
                    ypred = Array(predict(w,h,KnetArray(x[:,:,:,t:t])))
                    varloss += pixelwiseSoftloss(ypred,y[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,y[:,:,:,t:t])
                    varcorrect += testAcc[1]
                    varcount += testAcc[2]
                    #Image conversion
                    if saveImages
                        savePrediction(ypred[:,:,:,1],colors,outputPath(@sprintf("vKITTI/%s/GRU/%s/%s/%s/%05d.png",experiment,model,world,variation,imcounter)))
                        imcounter += 1
                    end
                end
                printlnFlush(@sprintf("World: %s\tVariation: %s\tLoss: %0.10f\tAccuracy: %0.4f%%",world,variation,varloss/batchCount,varcorrect/varcount*100))
                totalloss += varloss/batchCount
                totalcorrect += varcorrect
                totalcount += varcount
            end
        end
        printlnFlush(@sprintf("Total Loss: %0.10f\tTotal Accuracy: %0.4f%%",totalloss,totalcorrect/totalcount*100))
    end
end

#E1 : Degrees
#E2 : Half Train/Test

#   Conv
function run_vKITTI_Conv_RGBD_E2_test()
    test_vKITTI_Conv("RGBD_VGG",vKITTI_Conv_RGBD,"vKITTI/First_Second/Conv_RGBD.jld"; batchsize=32,saveImages=true,
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testMode="secondHalf",
        experiment="First_Second")
end

#   LSTM
function run_vKITTI_LSTM_RGBD_E2_test()
    test_vKITTI_LSTM("RGBD_VGG",vKITTI_LSTM_RGBD,"vKITTI/First_Second/LSTM_RGBD.jld"; saveImages=true,
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testMode="secondHalf",
        experiment="First_Second")
end

export run_vKITTI_Conv_RGBD_E2_test
export run_vKITTI_LSTM_RGBD_E2_test

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
module Uploaders

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export uploader

register_normal_element("q__uploader", context = @__MODULE__)

function uploader(args...;
              wrap::Function = StippleUI.DEFAULT_WRAPPER,
              kwargs...)
  wrap() do
    q__uploader(args...; kwargs...)
  end
end

end


text:
# port form musl library, commit e216951f509b71da193da2fc63e25b998740d58b
# see http://git.musl-libc.org/cgit/musl/tree/src/math/erf.c for implementation details
# origin: FreeBSD /usr/src/lib/msun/src/s_erf.c Copyright (c) 1993: Sun Microsystems, Inc.

let
global erf
global erfc

erx  = 8.45062911510467529297e-01 # 0x3FEB0AC1, 0x60000000

# Coefficients for approximation to  erf on [0,0.84375]
efx8 =  1.02703333676410069053e+00 # 0x3FF06EBA, 0x8214DB69
pp0  =  1.28379167095512558561e-01 # 0x3FC06EBA, 0x8214DB68
pp1  = -3.25042107247001499370e-01 # 0xBFD4CD7D, 0x691CB913
pp2  = -2.84817495755985104766e-02 # 0xBF9D2A51, 0xDBD7194F
pp3  = -5.77027029648944159157e-03 # 0xBF77A291, 0x236668E4
pp4  = -2.37630166566501626084e-05 # 0xBEF8EAD6, 0x120016AC
qq1  =  3.97917223959155352819e-01 # 0x3FD97779, 0xCDDADC09
qq2  =  6.50222499887672944485e-02 # 0x3FB0A54C, 0x5536CEBA
qq3  =  5.08130628187576562776e-03 # 0x3F74D022, 0xC4D36B0F
qq4  =  1.32494738004321644526e-04 # 0x3F215DC9, 0x221C1A10
qq5  = -3.96022827877536812320e-06 # 0xBED09C43, 0x42A26120

# Coefficients for approximation to  erf  in [0.84375,1.25]
pa0  = -2.36211856075265944077e-03 # 0xBF6359B8, 0xBEF77538
pa1  =  4.14856118683748331666e-01 # 0x3FDA8D00, 0xAD92B34D
pa2  = -3.72207876035701323847e-01 # 0xBFD7D240, 0xFBB8C3F1
pa3  =  3.18346619901161753674e-01 # 0x3FD45FCA, 0x805120E4
pa4  = -1.10894694282396677476e-01 # 0xBFBC6398, 0x3D3E28EC
pa5  =  3.54783043256182359371e-02 # 0x3FA22A36, 0x599795EB
pa6  = -2.16637559486879084300e-03 # 0xBF61BF38, 0x0A96073F
qa1  =  1.06420880400844228286e-01 # 0x3FBB3E66, 0x18EEE323
qa2  =  5.40397917702171048937e-01 # 0x3FE14AF0, 0x92EB6F33
qa3  =  7.18286544141962662868e-02 # 0x3FB2635C, 0xD99FE9A7
qa4  =  1.26171219808761642112e-01 # 0x3FC02660, 0xE763351F
qa5  =  1.36370839120290507362e-02 # 0x3F8BEDC2, 0x6B51DD1C
qa6  =  1.19844998467991074170e-02 # 0x3F888B54, 0x5735151D

# Coefficients for approximation to  erfc in [1.25,1/0.35]
ra0  = -9.86494403484714822705e-03 # 0xBF843412, 0x600D6435
ra1  = -6.93858572707181764372e-01 # 0xBFE63416, 0xE4BA7360
ra2  = -1.05586262253232909814e+01 # 0xC0251E04, 0x41B0E726
ra3  = -6.23753324503260060396e+01 # 0xC04F300A, 0xE4CBA38D
ra4  = -1.62396669462573470355e+02 # 0xC0644CB1, 0x84282266
ra5  = -1.84605092906711035994e+02 # 0xC067135C, 0xEBCCABB2
ra6  = -8.12874355063065934246e+01 # 0xC0545265, 0x57E4D2F2
ra7  = -9.81432934416914548592e+00 # 0xC023A0EF, 0xC69AC25C
sa1  =  1.96512716674392571292e+01 # 0x4033A6B9, 0xBD707687
sa2  =  1.37657754143519042600e+02 # 0x4061350C, 0x526AE721
sa3  =  4.34565877475229228821e+02 # 0x407B290D, 0xD58A1A71
sa4  =  6.45387271733267880336e+02 # 0x40842B19, 0x21EC2868
sa5  =  4.29008140027567833386e+02 # 0x407AD021, 0x57700314
sa6  =  1.08635005541779435134e+02 # 0x405B28A3, 0xEE48AE2C
sa7  =  6.57024977031928170135e+00 # 0x401A47EF, 0x8E484A93
sa8  = -6.04244152148580987438e-02 # 0xBFAEEFF2, 0xEE749A62

# Coefficients for approximation to  erfc in [1/.35,28]
rb0  = -9.86494292470009928597e-03 # 0xBF843412, 0x39E86F4A
rb1  = -7.99283237680523006574e-01 # 0xBFE993BA, 0x70C285DE
rb2  = -1.77579549177547519889e+01 # 0xC031C209, 0x555F995A
rb3  = -1.60636384855821916062e+02 # 0xC064145D, 0x43C5ED98
rb4  = -6.37566443368389627722e+02 # 0xC083EC88, 0x1375F228
rb5  = -1.02509513161107724954e+03 # 0xC0900461, 0x6A2E5992
rb6  = -4.83519191608651397019e+02 # 0xC07E384E, 0x9BDC383F
sb1  =  3.03380607434824582924e+01 # 0x403E568B, 0x261D5190
sb2  =  3.25792512996573918826e+02 # 0x40745CAE, 0x221B9F0A
sb3  =  1.53672958608443695994e+03 # 0x409802EB, 0x189D5118
sb4  =  3.19985821950859553908e+03 # 0x40A8FFB7, 0x688C246A
sb5  =  2.55305040643316442583e+03 # 0x40A3F219, 0xCEDF3BE6
sb6  =  4.74528541206955367215e+02 # 0x407DA874, 0xE79FE763
sb7  = -2.24409524465858183362e+01 # 0xC03670E2, 0x42712D62

function erfc1(x::T) where T<:FloatTypes
    s = abs(x) - 1
    P = @horner_oftype s pa0 pa1 pa2 pa3 pa4 pa5 pa6
    Q = @horner_oftype s 1.0 qa1 qa2 qa3 qa4 qa5 qa6
    return 1 - T(erx) - P/Q
end

function erfc2(ix::UInt32, x::T) where T<:FloatTypes
    if ix < highword(T(1.25))
        # 0.84375 <= |x| < 1.25
        return erfc1(x)
    end
    # 1.25 <= |x| < 28
    x = abs(x)
    s = 1/(x*x)
    if ix < highword(T(1/0.35000001))
        # 1.25 <= |x| < 1/.35 ~ 2.85714
        R = @horner_oftype s ra0 ra1 ra2 ra3 ra4 ra5 ra6 ra7
        S = @horner_oftype s 1.0 sa1 sa2 sa3 sa4 sa5 sa6 sa7 sa8
    else
        # 1/.35 <= |x| < 28
        R = @horner_oftype s rb0 rb1 rb2 rb3 rb4 rb5 rb6
        S = @horner_oftype s 1.0 sb1 sb2 sb3 sb4 sb5 sb6 sb7
    end
    z = trunclo(x)
    return exp(-z*z-T(0.5625))*exp((z-x)*(z+x)+R/S)/x
end

function erf(x::T) where T<:FloatTypes
    ix = highword(x)
    sign = (ix>>31) % Bool
    ix &= 0x7fffffff
    if ix >= highword(T(Inf)) # erf(nan)=nan, erf(+-inf)=+-1
        return 1-2*sign + 1/x
    end
    if ix < highword(T(0.84375)) # |x| < 0.84375
        if ix < highword(T(2)^-28) #|x| < 2**-28  avoid underflow
            return (8*x +T(efx8)*x)/8
        end
        z = x*x
        r = @horner_oftype z pp0 pp1 pp2 pp3 pp4
        s = @horner_oftype z 1 qq1 qq2 qq3 qq4 qq5
        y = r/s
        return x + x*y
    end
    if ix < highword(T(6)) # 0.84375 <= |x| < 6
        y = 1 - erfc2(ix,x)
    else
        y = 1 - realmin(T)
    end
    return sign ? -y : y
end

function erfc(x::T) where T<:FloatTypes
    ix = highword(x)
    sign = (ix>>31) % Bool
    ix &= 0x7fffffff
    if ix >= highword(T(Inf)) # erfc(nan)=nan, erfc(+-inf)=0,2
        return 2*sign + 1/x
    end
    if ix < highword(T(0.84375)) # |x| < 0.84375
        if ix < highword(T(2)^-56)  # |x| < 2**-56
            return 1 - x
        end
        z = x*x
        r = @horner_oftype z pp0 pp1 pp2 pp3 pp4
        s = @horner_oftype z 1 qq1 qq2 qq3 qq4 qq5
        y = r/s
        if sign || ix < highword(T(0.25)) # x < 1/4
            return 1 - (x+x*y)
        end
        return T(0.5) - (x - T(0.5) + x*y)
    end
    if ix < highword(T(28)) # 0.84375 <= |x| < 28
        return sign ? 2 - erfc2(ix,x) : erfc2(ix,x)
    end
    return sign ? 2 - realmin(T) : realmin(T)*realmin(T)
end

end


text:
function RandomGraph_ER_noise_experiments(testing=true)

    if testing

        default_n = [20]
        default_p = [0.05]

        alphas = [.5]
        betas = [1.0]

        trial_count = 1
        parameter_sets = [
            ([50,60],default_p),
            (default_n,[0.0,0.05]),
        ]
        result_path = all_results_path*"/RG_ERNoise/test/"

    else
        default_n = [250]
        default_p = [0.05]

        alphas = [.5,1.0]
        betas = [0.0,1.0,1e1,1e2]

        trial_count = 20
        parameter_sets = [
            (default_n,[0.01,0.05,.1,0.15,.2,0.25,0.3,0.4]),
            ([100,250,500,1000,1250,1500],default_p),
        ]
        result_path = all_results_path*"/RG_ERNoise/"
    end
    seeded=true

    if !isdir(result_path)
        mkpath(result_path)
    end

    for (n_sizes, p_remove) in parameter_sets
        #
        #    Lambda TAME
        #
        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=KlauAlgo(),method=ΛTAME_M(),alphas,betas,
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),n_sizes, p_remove,matchingMethod=ΛTAME_GramMatching(),profile=true)

        filename=filter(x-> !isspace(x),"LambdaTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:KlauAlgo_trialcount:$(trial_count).json") #Bug: trialcount shortened because filename is too long
        save_results(result_path*filename,X)

        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=LocalSearch(),method=ΛTAME_M(),alphas,betas,
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),n_sizes, p_remove,matchingMethod=ΛTAME_GramMatching(),profile=true)
        filename=filter(x-> !isspace(x),"LambdaTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:LocalSearch_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)

        #
        #    LowRankTAME
        #
        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=LocalSearch(),method=LowRankTAME_M(),alphas,betas,
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),n_sizes, p_remove,profile=true)
        filename=filter(x-> !isspace(x),"LowRankTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:LocalSearch_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)

        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=KlauAlgo(),method=LowRankTAME_M(),alphas,betas,
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),n_sizes, p_remove,profile=true)
        filename=filter(x-> !isspace(x),"LowRankTAME_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:KlauAlgo_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)
        #= low rank matching routines
        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=KlauAlgo(),method=LowRankTAME_M(),alphas,betas,
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1),n_sizes, p_remove,profile=true,low_rank_matching=true)
        filename=filter(x-> !isspace(x),"LowRankTAME-lrm_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:KlauAlgo_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)
        =#

        #
        #    LowRankEigenAlign
        #
        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=LocalSearch(),method=LowRankEigenAlign_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, p_remove,profile=true)
        filename=filter(x-> !isspace(x),"LowRankEigenAlign_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:LocalSearch_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)

        X = distributed_random_trials(trial_count,ErdosRenyiNoise(),seeded;postProcessing=KlauAlgo(),method=LowRankEigenAlign_M(),
                                    graph=RandomGeometric(),degreedist=LogNormal(log(5),1), n_sizes, p_remove,profile=true)
        filename=filter(x-> !isspace(x),"LowRankEigenAlign_graphType:RG_alphas:$(alphas)_betas:$(betas)_degdist:LogNormal-log5_n:$(n_sizes)_noiseModel:ErdosReyni_p:$(p_remove)_postProcess:KlauAlgo_trialcount:$(trial_count).json")
        save_results(result_path*filename,X)

    end

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
push!(LOAD_PATH,"../src/")

using Documenter, ClosedGroupFunctions

makedocs(
    sitename="ClosedGroupFunctions.jl",
    authors = "Felix Höddinghaus",
    pages = [
        "index.md",
        "Generating and Labelling" => "gens.md",
        "Conjugacy Classes" => "conj.md",
        "storage.md",
        "other.md",
    ],
)


text:
using LabelledArrays, OrdinaryDiffEq, Test

LorenzVector = @SLArray (3,) (:x,:y,:z)
LorenzParameterVector = @SLArray (3,) (:σ,:ρ,:β)

function f(u,p,t)
  x = p.σ*(u.y-u.x)
  y = u.x*(p.ρ-u.z) - u.y
  z = u.x*u.y - p.β*u.z
  LorenzVector(x,y,z)
end

u0 = LorenzVector(1.0,0.0,0.0)
p = LorenzParameterVector(10.0,28.0,8/3)
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Rosenbrock23())
@test sol.retcode === :Success
sol = solve(prob,Tsit5())
@test prob.u0 === sol.u[1] === u0
@test prob.p === p
@test sol[10].x > 0

function iip_f(du,u,p,t)
  du.x = p.σ*(u.y-u.x)
  du.y = u.x*(p.ρ-u.z) - u.y
  du.z = u.x*u.y - p.β*u.z
end

u0 = @LArray [1.0,0.0,0.0] (:x,:y,:z)
p = LorenzParameterVector(10.0,28.0,8/3)
tspan = (0.0,10.0)
prob = ODEProblem(iip_f,u0,tspan,p)

@test similar(u0) isa LArray
@test zero(u0) isa LArray

sol = solve(prob,Rosenbrock23())
sol = solve(prob,Tsit5())
@test typeof(prob.u0) == eltype(sol.u) == typeof(u0)
@test prob.p === p
@test sol[10].x > 0


text:
using Test, LinearAlgebra
using BellScenario: stirling2, stirling2_partitions, stirling2_matrices, permutation_matrices, base_n_val, n_choose_k_matrices

@testset "./src/combinatorics.jl" begin

@testset "stirling2()" begin
    @testset "domain errors" begin
        @test_throws DomainError stirling2(1,0)
        @test_throws DomainError stirling2(1,2)
    end

    @test stirling2(1,1) == 1

    @testset "n = 2" begin
        @test stirling2(2,1) == 1
        @test stirling2(2,2) == 1
    end

    @testset "n = 3" begin
        @test stirling2(3,1) == 1
        @test stirling2(3,2) == 3
        @test stirling2(3,3) == 1
    end

    @testset "n = 4" begin
        @test stirling2(4,1) == 1
        @test stirling2(4,2) == 7
        @test stirling2(4,3) == 6
        @test stirling2(4,4) == 1
    end

    @testset "n = 5" begin
        @test stirling2(5,1) == 1
        @test stirling2(5,2) == 15
        @test stirling2(5,3) == 25
        @test stirling2(5,4) == 10
        @test stirling2(5,5) == 1
    end

    @testset "n = 6" begin
        @test stirling2(6,1) == 1
        @test stirling2(6,2) == 31
        @test stirling2(6,3) == 90
        @test stirling2(6,4) == 65
        @test stirling2(6,5) == 15
        @test stirling2(6,6) == 1
    end

    @testset "n = 10" begin
        @test stirling2(10,1) == 1
        @test stirling2(10,2) == 511
        @test stirling2(10,3) == 9330
        @test stirling2(10,4) == 34105
        @test stirling2(10,5) == 42525
        @test stirling2(10,6) == 22827
        @test stirling2(10,7) == 5880
        @test stirling2(10,8) == 750
        @test stirling2(10,9) == 45
        @test stirling2(10,10) == 1
    end
end

@testset "stirling2_partitions()" begin
    @testset "base cases" begin
        @test stirling2_partitions(0,0) == []
        @test stirling2_partitions(1,0) == []
        @test stirling2_partitions(1,2) == []
        @test stirling2_partitions(5,1) == [[[1,2,3,4,5]]]
    end

    @testset "simple cases" begin
        @test stirling2_partitions(3,2) == [[[1,2],[3]],[[2],[1,3]],[[1],[2,3]]]

        @test stirling2_partitions(4,2) == [
            [[1, 2, 3], [4]],
            [[3], [1, 2, 4]],
            [[1, 2], [3, 4]],
            [[1, 3], [2, 4]],
            [[2], [1, 3, 4]],
            [[2, 3], [1, 4]],
            [[1], [2, 3, 4]],
        ]

        @test stirling2_partitions(4,4) == [[[1],[2],[3],[4]]]
    end

    @testset "general partitions" begin
        for n in 2:7
            for k in 1:n
                partitions = stirling2_partitions(n,k)
                @test length(partitions) == stirling2(n,k)
                @test length(unique(partitions)) == stirling2(n,k)
                @test all(s -> length(s) == k, partitions)
                @test all(s -> sort(cat(s...,dims=1)) == [1:n...], partitions)
            end
        end
    end
end

@testset "stirling2_matrices()" begin
    @testset "domain errors" begin
        @test_throws DomainError n_choose_k_matrices(2,3)
        @test_throws DomainError n_choose_k_matrices(2,0)
    end

    @testset "simple cases" begin
        @test stirling2_matrices(4,1) == [
            [1 1 1 1]
        ]
        @test stirling2_matrices(3,2) == [
            [1 1 0;0 0 1],[0 1 0;1 0 1],[1 0 0;0 1 1]
        ]
    end

    @testset "general cases" begin
        for n in 2:7
            for k in 1:n
                matrices = stirling2_matrices(n,k)

                @test length(matrices) == stirling2(n,k)
                @test length(unique(matrices)) == stirling2(n,k)
            end
        end
    end
end

@testset "permutation_matrices()" begin
    @testset "dim = 2" begin
        @test permutation_matrices(2) == [[1 0; 0 1], [0 1; 1 0]]
    end

    @testset "dim = 3" begin
        @test issetequal(permutation_matrices(3), [
            [1 0 0; 0 1 0; 0 0 1],
            [1 0 0; 0 0 1; 0 1 0],
            [0 1 0; 1 0 0; 0 0 1],
            [0 1 0; 0 0 1; 1 0 0],
            [0 0 1; 1 0 0; 0 1 0],
            [0 0 1; 0 1 0; 1 0 0],
        ])
    end

    @testset "permtutation map validity" for N in 4:5
        perms = permutation_matrices(N)

        @testset "maps are unique" begin
            @test size(perms) == (factorial(N),)
            @test size(perms) == size(unique(perms))
        end

        @testset "matrix dimension, determinant, and invertibility" for perm in perms
            # dimension
            @test size(perm) == (N,N)

            # parity
            parity = det(perm)
            @test (parity == -1) || (parity == 1)

            # invertibility
            @test perm*perm' == Matrix{Int64}(I,N,N)
        end
    end
end

@testset "n_choose_k_matrices()" begin
    @testset "Domain Errors" begin
        @test_throws DomainError n_choose_k_matrices(2,3)
        @test_throws DomainError n_choose_k_matrices(2,0)
    end

    @testset "simple testcases" begin
        @test n_choose_k_matrices(3,3) == [[1 0 0;0 1 0;0 0 1]]

        @test n_choose_k_matrices(3,2) == [
            [1 0;0 1;0 0],[1 0;0 0;0 1],[0 0;1 0;0 1]
        ]

        @test n_choose_k_matrices(4,2) == [
            [1 0;0 1;0 0;0 0],[1 0;0 0;0 1;0 0],[1 0;0 0;0 0;0 1],
            [0 0;1 0;0 1;0 0],[0 0;1 0;0 0;0 1],[0 0;0 0;1 0;0 1]
        ]
    end

    @testset "general tests" begin
        for n in 2:5
            for k in 1:n
                matrices = n_choose_k_matrices(n,k)

                @test length(matrices) == binomial(n,k)
                @test length(unique(matrices)) == binomial(n,k)
            end
        end
    end
end

@testset "base_n_val()" begin
    @testset "little endian" begin
        @test 1 == base_n_val(digits(1,base=2,pad=2), 2, big_endian = false)
        @test 57 == base_n_val(digits(57,base=3,pad=5), 3, big_endian = false)
    end

    @testset "big endian" begin
        @test 12 == base_n_val([1,1,0,0],2)
        @test 974 == base_n_val([1,2,3,4,4],5)
    end

    @testset "argument errors" begin
        @test_throws ArgumentError base_n_val([1,2,0,5],5)
        @test_throws ArgumentError base_n_val([1,2,0,5],4)
        @test_throws ArgumentError base_n_val([-1,0,1],2)
    end
end

end


text:
using Plugins
import Plugins.symbol, Plugins.setup!, Plugins.setup!
using Test
using Random

# To run only selected tests, use e.g.:
#
#   Pkg.test("Plugins", test_args=["deps"])
#
enabled_tests = lowercase.(ARGS)
function addtests(fname)
    key = lowercase(splitext(fname)[1])
    if isempty(enabled_tests) || key in enabled_tests
        Random.seed!(42)
        include(fname)
    end
end

addtests("basics.jl")
addtests("deps.jl")
addtests("assembled_types.jl")


text:
struct ∇SemWLS{Vt <: Union{AbstractArray, UniformScaling{Bool}},
        St <: AbstractArray, Lt <: AbstractArray} <: DiffFunction
    V::Vt
    s::St
    L::Lt
end

function ∇SemWLS(semwls::A, nobs) where {A <: SemWLS}
    L = elimination_matrix(nobs)
    return ∇SemWLS(semwls.V, semwls.s, L)
end

function (diff::∇SemWLS)(par, grad::AbstractVector, model::Sem{O, I, L, D}) where
        {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    model.imply.imp_fun(model.imply.imp_cov, par)
    model.imply.gradient_fun(model.imply.∇Σ, par)
    grad .= (-2*(diff.s-model.imply.imp_cov)'*diff.V*model.imply.∇Σ)'
end


text:
function score(x, y)
    distance = hypot(x, y)
    distance > 10 && return 0
    distance > 5 && distance <= 10 && return 1
    distance > 1 && distance <= 5  && return 5
    distance <= 1 && return 10
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
ProjDir = dirname(@__FILE__)
cd(ProjDir) do

A = Float64[
  1/2 1/3 1/4 1/5 1/6;
  1/3 1/4 1/5 1/6 1/7;
  1/4 1/5 1/6 1/7 1/8;
  1/5 1/6 1/7 1/8 1/9;
  1/6 1/7 1/8 1/9 1/10;
]
b = ones(5)

(Al, Au, Ap) = lu(A)

@show y = Al \ b[Ap]
println()

@show x = Au \ y
println()

@show As = sparse(A)
println()

@show As\b
println()

@show F = lu(As)
println()

@show F.L * F.U == F.Rs .* A[F.p, :]
println()

@show F.L*F.U
println()

@show sparse(F.Rs.* A[F.p, F.q])
println()

@show lu(A)
println()

end


text:
"""Test whether the interface for your `geom` has been implemented correctly."""
function test_interface_for_geom(geom)
    try
        type = geomtype(geom)
        ncoord(geom)
        ngeom(geom)
        getgeom(geom, 1)
        if type == Point()
            getcoord(geom, 1)
        end
    catch e
        println("You're missing an implementation: $e")
        return false
    end
    return true
end


text:
#=
This file contains all function used to parse the data and return the data to an
appropriate format to be compute.

Link of the data : http://www.roadef.org/challenge/2005/fr/sujet.php
The data should be in ../data/Instances_Type/name_instance/
There is 3 types of instances labeled A, B and X, so folders
should be named Instances_A/, Instances_B/, and Instances_X/.
=#

# Constant of the problem
const OPTIMISATION_FILE_NAME = "optimization_objectives.txt"
const PAINT_FILE_NAME = "paint_batch_limit.txt"
const RATIO_FILE_NAME = "ratios.txt"
const VEHICLES_FILE_NAME = "vehicles.txt"

# An Instance structure that is used to format data as we want.
struct Instance
    # objectif function
    HPRC_rank::Int
    LPRC_rank::Int # If there is no LPRC, LPRC_rank = -1
    PCB_rank::Int
    # paint limitation
    nb_paint_limitation::Int
    # ratio constraint
    RC_p::Array{Int, 1}
    RC_q::Array{Int, 1}
    nb_HPRC::Int
    nb_LPRC::Int
    # sequence vehicle data
    RC_flag::Array{Bool, 2}
    RC_keys::Array{Int, 1}
    HPRC_keys::Array{Int, 1}
    LPRC_keys::Array{Int, 1}
    same_RC::Dict{Int, Array{Int, 1}}
    same_HPRC::Dict{Int, Array{Int, 1}}
    same_LPRC::Dict{Int, Array{Int, 1}}
    color_code::Array{Int, 1}
    same_color::Dict{Int, Array{Int, 1}}
    # Number of vehicles that were not build the precedent day.
    nb_late_prec_day::Int # Usage 1:nb_late_prec_day give the list of index of those vehicles.
    nb_cars::Int
end

"""
    parser(instance_name::String, instance_type::String, path_folder::String=string(@__DIR__)*"/../data/Instances_")

Returns return the instance of type instance_type parsed form the files in path_folder.
"""
function parser(instance_name::String, instance_type::String, path_folder::String=string(@__DIR__)*"/../data/Instances_")
    path = path_folder * instance_type * "/" * instance_name * "/"

    #Lecture  OPTIMISATION_FILE_NAME
    HPRC_rank = 0
    LPRC_rank = -1
    PCB_rank = 0
    open(path * OPTIMISATION_FILE_NAME) do f
        lines = readlines(f)
        for i in 2:length(lines)
            if lines[i] != ""
                if lines[i][3] == 'h'
                    HPRC_rank = i-1
                elseif lines[i][3] == 'l'
                    LPRC_rank = i-1
                elseif lines[i][3] == 'p'
                    PCB_rank = i-1
                end
            end
        end
    end
    @assert HPRC_rank != 0
    @assert PCB_rank != 0

    #Lecture PAINT_FILE_NAME
    nb_paint_limitation = 0
    open(path * PAINT_FILE_NAME) do f
        lines = readlines(f)
        nb_paint_limitation = parse(Int, lines[2][1:(end-1)])
    end

    #Lecture RATIO_FILE_NAME
    nb_HPRC = 0
    nb_LPRC = 0
    RC_p = Array{Int,1}()
    RC_q = Array{Int,1}()
    open(path * RATIO_FILE_NAME) do f
        lines = readlines(f)
        for line in lines[2:end]
            if line != ""
                values = split(line, ";", limit=2)
                ratio = split(values[1], "/")
                push!(RC_p, parse(Int, ratio[1]))
                push!(RC_q, parse(Int, ratio[2]))
                if values[2][1] == '1'
                    nb_HPRC += 1
                elseif values[2][1] == '0'
                    nb_LPRC += 1
                end
            end
        end
    end

    # TODO:
    # string(Int(RC_flag[car, option])) * RC_value
    # become RC_value *string(Int(RC_flag[car, option]))
    # then:
    # two cars -> difference in RC value is less than 2^nb_LPRC
    # one array not needed but more time to compute difference later
    #Lecture VEHICLES_FILE_NAME
    day = 0
    indice = 0
    nb_RC = nb_HPRC + nb_LPRC
    color_code = Array{Int, 1}()
    RC_keys = Array{Int, 1}()
    HPRC_keys = Array{Int, 1}()
    LPRC_keys = Array{Int, 1}()
    RC_flag = Array{Bool, 2}(undef, 0, nb_RC)
    same_RC = Dict{Int, Array{Int, 1}}()
    same_HPRC = Dict{Int, Array{Int, 1}}()
    same_LPRC = Dict{Int, Array{Int, 1}}()
    same_color = Dict{Int, Array{Int, 1}}()
    nb_cars = 0
    nb_cars_total = 0
    nb_late_prec_day = 0
    open(path * VEHICLES_FILE_NAME) do f
        lines = readlines(f)
        if lines[end] == ""
            nb_cars_total = length(lines) - 2
        else
            nb_cars_total = length(lines) - 1
        end
        RC_flag = Array{Bool, 2}(undef, nb_cars_total, nb_RC)
        color_code = Array{Int, 1}(undef, nb_cars_total)
        RC_keys = Array{Int, 1}(undef, nb_cars_total)
        HPRC_keys = Array{Int, 1}(undef, nb_cars_total)
        LPRC_keys = Array{Int, 1}(undef, nb_cars_total)
        for line in lines[2:end]
            if line != ""
                values = split(line, ";", keepempty = false)
                nb_cars += 1
                car_RC_value = "0"
                for rc in 5:length(values)
                    RC_flag[nb_cars, rc - 4] = (values[rc] == "1")
                    car_RC_value = values[rc] * car_RC_value
                end
                RC_keys[nb_cars] = parse(Int, car_RC_value, base = 2)
                HPRC_keys[nb_cars] = parse(Int, SubString(car_RC_value, nb_LPRC+1), base = 2)
                LPRC_keys[nb_cars] = parse(Int, "0"*SubString(car_RC_value, 1:nb_LPRC), base = 2)
                color_code[nb_cars] = parse(Int, values[4])
                if nb_late_prec_day == 0
                    date = split(values[1], " ")
                    if day == 0
                        day = date[3]
                    elseif date[3] != day
                        day = date[3]
                        nb_late_prec_day = nb_cars-1
                        same_RC[RC_keys[nb_cars]] = Array{Int, 1}([nb_cars])
                        same_HPRC[HPRC_keys[nb_cars]] = Array{Int, 1}([nb_cars])
                        same_LPRC[LPRC_keys[nb_cars]] = Array{Int, 1}([nb_cars])
                        same_color[color_code[nb_cars]] = Array{Int, 1}([nb_cars])
                    end
                else
                    if !haskey(same_RC, RC_keys[nb_cars])
                        same_RC[RC_keys[nb_cars]] = Array{Int, 1}([nb_cars])
                        if !haskey(same_HPRC, HPRC_keys[nb_cars])
                            same_HPRC[HPRC_keys[nb_cars]] = Array{Int, 1}([nb_cars])
                        else
                            push!(same_HPRC[HPRC_keys[nb_cars]], nb_cars)
                        end
                        if !haskey(same_LPRC, LPRC_keys[nb_cars])
                            same_LPRC[LPRC_keys[nb_cars]] = Array{Int, 1}([nb_cars])
                        else
                            push!(same_LPRC[LPRC_keys[nb_cars]], nb_cars)
                        end
                    else
                        push!(same_RC[RC_keys[nb_cars]], nb_cars)
                        push!(same_HPRC[HPRC_keys[nb_cars]], nb_cars)
                        push!(same_LPRC[LPRC_keys[nb_cars]], nb_cars)
                    end
                    if !(color_code[nb_cars] in keys(same_color))
                        same_color[color_code[nb_cars]] = Array{Int, 1}([nb_cars])
                    else
                        push!(same_color[color_code[nb_cars]], nb_cars)
                    end
                end
            end
        end
    end
    @assert nb_cars_total == nb_cars

    # for i in (nb_late_prec_day+1):nb_cars
    #     @assert i in same_RC[RC_keys[i]]
    #     @assert i in same_LPRC[LPRC_keys[i]]
    #     @assert i in same_HPRC[HPRC_keys[i]]
    # end

    return Instance(
            HPRC_rank, LPRC_rank, PCB_rank,                            # objectives file
            nb_paint_limitation,                                       # paint file
            RC_p, RC_q, nb_HPRC, nb_LPRC,                              # ratio file
            RC_flag, RC_keys, HPRC_keys, LPRC_keys,
            same_RC, same_HPRC, same_LPRC,
            color_code, same_color, nb_late_prec_day, nb_cars
        )
end


text:
using JSON
using Test
using ComparativeAutograder

const ARITHMETIC_SOLUTION = """
function h(x, y, z)
    return x*y, y*z, x*z, x+y+z
end
"""

@testset "Arithmetic" begin

    cases = [TestCase((randn(), randn(), randn())) for _ in 1:100]
    suite = TestSuite("h", cases)

    soln_path, soln_io = mktemp()
    sub_path, sub_io = mktemp()
    sub_stdout, sub_stderr = nothing, nothing

    try
        write(soln_io, ARITHMETIC_SOLUTION)
        close(soln_io)
        write(sub_io, """h(x, y, z) = x*y, y*z, x*z, x+y+z""")
        close(sub_io)
        sub_stdout, sub_stderr = @capture_stdstreams runtestsuite(
            suite, soln_path, sub_path
        )
    finally
        rm.([soln_path, sub_path])
    end

    @test sub_stdout != nothing
    sub_data = JSON.parse(sub_stdout)
    @test sub_data["passed"] == true
end


text:
# Generate Task and TaskGDS types with names e.g. Attractor and AttractorGDS
# and no other fields
macro standard_task_type(name)
    esc(quote
        struct $name{F} <: Task{F}
            task_map::F
        end

        struct $(Symbol("$name",:GDS)){F} <: TaskGDS{F}
            task_map::F
        end
    end)
end

@standard_task_type BlankTask
@standard_task_type Damping
@standard_task_type Attractor
@standard_task_type CollisionAvoidance
@standard_task_type VelocityLimit

macro ext_standard_task_type(name)
    esc(quote
        struct $name{F} <: PBDS.Task{F}
            task_map::F
        end

        struct $(Symbol("$name",:GDS)){F} <: PBDS.TaskGDS{F}
            task_map::F
        end
    end)
end


text:
export SvmRanker

@with_kw mutable struct SvmRanker <: BaseEstimator
    svm::Vector{UInt8} = UInt8[]
    c::Float32 = 0.01
    p::Int = 1
    o::Int = 2
end

is_classifier(::SvmRanker) = false
is_ranker(::SvmRanker) = true

function gridparams(m::SvmRanker)
    grid = ["c" => [0.01, 0.1, 1, 10, 100]]
    params = gridparams(grid)
end

function fit!(m::SvmRanker, x, y, w = nothing; group = nothing, columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    if isnothing(group) && ndims(x) == 3
        group = repeat([size(x, 2)], size(x, 3))
    end
    @unpack c, p, o = m
    dst = to_svm(x, y, w, group = group)
    run(`$SVM_RANK_LEARN -c $c -p $p -o $o $dst svm`)
    m.svm = read("svm")
    foreach(rm, ["svm", dst])
    BSON.bson("model.bson", model = m)
    return m
end

function predict(m::SvmRanker, x)
    @unpack svm = m
    dst = to_svm(x)
    write("svm", svm)
    run(`$SVM_RANK_CLASSIFY $dst svm preds`)
    ŷ = readdlm("preds", Float32)
    foreach(rm, ["svm", "preds", dst])
    return ŷ
end


text:
function build_indvar_dict(ex,depvar)
  indvar_dict = OrderedDict{Symbol,Int}()
  cur_sym = 0
  for i in 2:2:length(ex.args) #Every odd line is line number
    arg = ex.args[i].args[1] #Get the first thing, should be dsomething
    firstarg = Symbol(first(string(arg))) # Check for d
    if firstarg == :d
      nodarg = Symbol(join(Base.Iterators.drop(string(arg), 1)))
      if nodarg == depvar
        warn("$depvar is fixed as the independent variable but is also used as a dependent variable. Results my be incorrect.")
      end
      if !haskey(indvar_dict,nodarg)
        cur_sym += 1
        indvar_dict[nodarg] = cur_sym
      else
        error("The derivative term for $nodarg is repeated. This is not allowed.")
      end
    end
  end
  syms = indvar_dict.keys
  indvar_dict,syms
end

function build_param_list(params)
  param_list = Vector{Symbol}();
  for i in 1:length(params)
    if typeof(params[i]) <: Symbol
      push!(param_list,params[i])
    elseif params[i].head == :call || params[i].head == :(=)
      warn("p=>val and p=val are deprecated. Simply list the parameters. See the DifferentialEquations.jl documentation for more information on the syntax change.")
    end
  end
  param_dict
end


text:
### AbstractSciMLOperator Interface

#=
1. Function call and multiplication: L(du, u, p, t) for inplace and du = L(u, p, t) for
   out-of-place, meaning L*u and mul!(du, L, u).
2. If the operator is not a constant, update it with (u,p,t). A mutating form, i.e.
   update_coefficients!(A,u,p,t) that changes the internal coefficients, and a
   out-of-place form B = update_coefficients(A,u,p,t).
3. isconstant(A) trait for whether the operator is constant or not.
4. islinear(A) trait for whether the operator is linear or not.
=#

Base.eltype(::Type{AbstractSciMLOperator{T}}) where T = T
Base.eltype(::AbstractSciMLOperator{T}) where T = T
update_coefficients!(L,u,p,t) = nothing
update_coefficients(L,u,p,t) = L

# Traits
isconstant(::AbstractSciMLOperator) = false
islinear(::AbstractSciMLOperator) = false
has_expmv!(L::AbstractSciMLOperator) = false # expmv!(v, L, t, u)
has_expmv(L::AbstractSciMLOperator) = false # v = exp(L, t, u)
has_exp(L::AbstractSciMLOperator) = false # v = exp(L, t)*u
has_mul(L::AbstractSciMLOperator) = true # du = L*u
has_mul!(L::AbstractSciMLOperator) = false # mul!(du, L, u)
has_ldiv(L::AbstractSciMLOperator) = false # du = L\u
has_ldiv!(L::AbstractSciMLOperator) = false # ldiv!(du, L, u)

### AbstractDiffEqLinearOperator Interface

#=
1. AbstractDiffEqLinearOperator <: AbstractSciMLOperator
2. Can absorb under multiplication by a scalar. In all algorithms things like
   dt*L show up all the time, so the linear operator must be able to absorb
   such constants.
4. isconstant(A) trait for whether the operator is constant or not.
5. Optional: diagonal, symmetric, etc traits from LinearMaps.jl.
6. Optional: exp(A). Required for simple exponential integration.
7. Optional: expmv(A,u,p,t) = exp(t*A)*u and expmv!(v,A::DiffEqOperator,u,p,t)
   Required for sparse-saving exponential integration.
8. Optional: factorizations. A_ldiv_B, factorize et. al. This is only required
   for algorithms which use the factorization of the operator (Crank-Nicholson),
   and only for when the default linear solve is used.
=#

# Extra standard assumptions
isconstant(::AbstractDiffEqLinearOperator) = true
islinear(o::AbstractDiffEqLinearOperator) = isconstant(o)

# Other ones from LinearMaps.jl
# Generic fallbacks
LinearAlgebra.exp(L::AbstractDiffEqLinearOperator,t) = exp(t*L)
has_exp(L::AbstractDiffEqLinearOperator) = true
expmv(L::AbstractDiffEqLinearOperator,u,p,t) = exp(L,t)*u
expmv!(v,L::AbstractDiffEqLinearOperator,u,p,t) = mul!(v,exp(L,t),u)
# Factorizations have no fallback and just error


text:
using Test, Astro

# references are to a book by Meeus

print("Degrees to Degrees Minutes Seconds")
@test(d_to_dms(45.5) == (45.0,30.0,0.0))
println(" passed")

print("7.a Convert Gregorian date to Julian day number")
jd = cal_to_jd(1957, 10, 4.81)
@test isapprox(jd, 2436116.31, atol=0.01)
println(" passed")

print("7.b Convert Julian date to Julian day number")
jd = cal_to_jd(333, 1, 27.5, false)
@test isapprox(jd, 1842713.0, atol=0.01)
println(" passed")

print("7.c Convert Julian day number to Gregorian date")
yr, mo, day = jd_to_cal(2436116.31)
@test(yr == 1957)
@test(mo == 10)
@test isapprox(day, 4.81, atol=0.01)
println(" passed")

print("7.c(1) Convert Julian day number to Julian date")
yr, mo, day = jd_to_cal(1842713.0, false)
@test(yr == 333)
@test(mo == 1)
@test isapprox(day, 27.5, atol=0.01)
println(" passed")

print("7.c(2) Convert Julian day number to Julian date")
yr, mo, day = jd_to_cal(1507900.13, false)
@test(yr == -584)
@test(mo == 5)
@test isapprox(day, 28.63, atol=0.01)
println(" passed")

print("7.d Time interval in days")
jd0 = cal_to_jd(1910, 4, 20.0)
jd1 = cal_to_jd(1986, 2, 9.0)
@test(jd1 - jd0 == 27689)
println(" passed")

print("7.d(1) Time interval in days")
jd = cal_to_jd(1991, 7, 11)
jd = jd + 10000
yr, mo, day = jd_to_cal(jd)
@test(yr == 2018)
@test(mo == 11)
@test(day == 26)
println(" passed")

print("7.e Day of the week")
jd = cal_to_jd(1954, 6, 30.0)
@test(jd == 2434923.5)
dow = jd_to_day_of_week(jd)
@test(dow == 3)
println(" passed")

print("7.f Day of the year")
N = cal_to_day_of_year(1978, 11, 14)
@test(N == 318)
println(" passed")

print("7.g Day of the year")
N = cal_to_day_of_year(1988, 4, 22)
@test(N == 113)
println(" passed")

print("7(pg 66-1) Day of the year to calendar")
mo, day = day_of_year_to_cal(1978, 318)
@test(mo == 11)
@test(day == 14)
println(" passed")

print("7(pg 66-2) Day of the year to calendar")
mo, day = day_of_year_to_cal(1988, 113)
@test(mo == 4)
@test(day == 22)
println(" passed")

print("difference between angles")
diffr = diff_angle(pi/3, pi/5)
@test(diffr == 5.8643062867009474)
println(" passed")

print("difference between angles 2")
# TODO utils-test.jl check these results more
@test isapprox(diff_angle(deg2rad(10), deg2rad(0))   |> rad2deg, 350, atol=.1)
@test isapprox(diff_angle(deg2rad(0), deg2rad(10))   |> rad2deg, 10,  atol=.1)
@test isapprox(diff_angle(deg2rad(-10), deg2rad(10)) |> rad2deg, 20,  atol=.1)
println(" passed")

print("degree/minute/second to decimal degrees")
@test dms_to_d(10,45,0) == 10.75
@test dms_to_d(-10,45,0) == -10.75
println(" passed")

print("hangle to d")
@test hangle_to_dec_deg(5,0,0) == 75
@test hangle_to_dec_deg(12,0,0) == 180.0
@test isapprox(hangle_to_dec_deg(23,55,55), 358.9791666666667, atol=0.0001)
println(" passed")

print("fractional day and hms day conversions")
@test fday_to_hms(0) == (0.0,0.0,0.0)
@test fday_to_hms(0.1) == (2.0,24.0,0.0)
@test fday_to_hms(0.2) == (4.0,48.0,0.0)
@test fday_to_hms(0.25) == (6.0,0.0,0.0)
@test fday_to_hms(0.5) == (12.0,0.0,0.0)
@test fday_to_hms(0.75) == (18.0,0.0,0.0)
@test fday_to_hms(0.99) == (23.0,45.0,36.0)
@test hms_to_fday(23.0, 45, 36) == 0.99
@test hms_to_fday(18.0, 0, 0) == 0.75
println(" passed")

print("3.1 Interpolation")
y = interpolate3(0.18125, [0.884226, 0.877366, 0.870531])
@test isapprox(y, 0.876125, atol=1e-6)
@test isapprox(interpolate_angle3(1, [0, 30,60]), 34.86725, atol=0.01)
@test isapprox(interpolate_angle3(0, [0, 30,60]), 30.0,     atol=0.01)
@test isapprox(interpolate_angle3(-1, [0, 30,60]), 25.13,   atol=0.01)
println(" passed")

print("polynomials")
@test @polynomial_horner(18, 7, 4, 1) == 403
@test polynomial([7, 4, 1], 18) == 403
@test isapprox(@polynomial_horner(18, 9, 13, 7.4, 4, 1), 130944.59, atol=.1)
@test isapprox(polynomial((9, 13, 7.4, 4, 1), 18), 130944.6, atol=.1)
println(" passed")


text:
module TestDay03

import AOC2021.Day03.part1, AOC2021.Day03.part2
using Test

test_input_raw = """
00100
11110
10110
10111
10101
01111
00111
11100
10000
11001
00010
01010"""

test_input = [string(x) for x in split(test_input_raw, "\n")]

@testset "Day 03" begin
    @testset "part 1" begin
        @test part1(test_input) == 198
    end
    @testset "part 2" begin
        @test part2(test_input) == 230
    end
end

end # module


text:
# Copied and slightly modified from
# https://github.com/JuliaWeb/HTTP.jl/blob/master/src/Handlers.jl

using HTTP: Stream
# using HTTP.Handlers: Node, Leaf, find, segment, insert!, match
using URIs
using  AbstractTrees

export register!, match_middleware

# tree-based router handler
mutable struct Variable
    # would it be better to use symbol
    # so captured variable can easily be converted to namedtuples
    name::String
    pattern::Union{Nothing,Regex}
end

Base.show(io::IO, x::Variable) = print(io, "Variable(name=:$(x.name), pattern=$(x.pattern)")

const VARREGEX = r"^{([[:alnum:]]+):?(.*)}$"

function Variable(pattern)
    re = match(VARREGEX, pattern)
    if re === nothing
        error("problem parsing path variable for route: `$pattern`")
    end
    pat = re.captures[2]
    return Variable(re.captures[1], pat == "" ? nothing : Regex(pat))
end

struct Leaf
    method::String # should be a lowercase symbol
    variables::Vector{Tuple{Int,String}} # [(path_segement_idx, value)]
    path::String
    handler::Any
end

Base.show(io::IO, x::Leaf) = print(io, "Leaf($(x.method))")

mutable struct Node
    segment::Union{String,Variable}
    exact::Vector{Node} # sorted alphabetically, all x.segment are String
    conditional::Vector{Node} # unsorted; will be applied in source-order; all x.segment are Regex
    wildcard::Union{Node,Nothing} # unconditional variable or wildcard
    doublestar::Union{Node,Nothing} # /** to match any length of path; must be final segment

    # should this be a dict instead
    methods::Vector{Leaf}
end

function Base.show(io::IO, x::Node)
    print(io, "Node($(x.segment))")
end

isvariable(x) = startswith(x, "{") && endswith(x, "}")
segment(x) = segment == "*" ? String(segment) : isvariable(x) ? Variable(x) : String(x)

Node(x) = Node(x, Node[], Node[], nothing, nothing, Leaf[])
Node() = Node("*")

function AbstractTrees.children(node::Node)
    l = []
    # the order they are return in the list will
    # determin which order they are iterated in
    if !isempty(node.exact)
        push!(l, node.exact...)
    end
    if !isempty(node.conditional)
        push!(l, node.conditional...)
    end
    if !isnothing(node.wildcard)
        push!(l, node.wildcard)
    end
    if !isnothing(node.doublestar)
        push!(l, node.wildcard)
    end
    filter!(x -> !isnothing(x), l)
    return l
end

AbstractTrees.nodetype(::Type{Node}) = Node
# AbstractTrees.nodevalue(n::Node) = n.methods

#   nodevalue(tree, node_index)

function AbstractTrees.printnode(io, x)
     print(io, "Node($(x.segment))")
end

function find(y, itr; by=identity, eq=(==))
    for (i, x) in enumerate(itr)
        eq(by(x), y) && return i
    end
    return nothing
end

function Base.insert!(node::Node, leaf, segments, i)
    if i > length(segments)
        # time to insert leaf method match node
        j = find(leaf.method, node.methods; by=x -> x.method, eq=(x, y) -> x == "*" || x == y)
        if j === nothing
            push!(node.methods, leaf)
        else
            # hmmm, we've seen this route before, warn that we're replacing
            @warn "replacing existing registered route; $(node.methods[j].method) => \"$(node.methods[j].path)\" route with new path = \"$(leaf.path)\""
            node.methods[j] = leaf
        end
        return
    end
    segment = segments[i]
    # @show segment, segment isa Variable
    if segment isa Variable
        # if we're inserting a variable segment, add variable name to leaf vars array
        push!(leaf.variables, (i, segment.name))
    end
    # figure out which kind of node this segment is
    if segment == "*" || (segment isa Variable && segment.pattern === nothing)
        # wildcard node
        if node.wildcard === nothing
            node.wildcard = Node(segment)
        end
        return Base.insert!(node.wildcard, leaf, segments, i + 1)
    elseif segment == "**"
        # double-star node
        if node.doublestar === nothing
            node.doublestar = Node(segment)
        end
        if i < length(segments)
            error("/** double wildcard must be last segment in path")
        end
        return Base.insert!(node.doublestar, leaf, segments, i + 1)
    elseif segment isa Variable
        # conditional node
        # check if we've seen this exact conditional segment before
        j = find(segment.pattern, node.conditional; by=x -> x.segment.pattern)
        if j === nothing
            # new pattern
            n = Node(segment)
            push!(node.conditional, n)
        else
            n = node.conditional[j]
        end
        return Base.insert!(n, leaf, segments, i + 1)
    else
        # exact node
        @assert segment isa String
        j = find(segment, node.exact; by=x -> x.segment)
        if j === nothing
            # new exact match segment
            n = Node(segment)
            push!(node.exact, n)
            sort!(node.exact; by=x -> x.segment)
            return Base.insert!(n, leaf, segments, i + 1)
        else
            # existing exact match segment
            return Base.insert!(node.exact[j], leaf, segments, i + 1)
        end
    end
end

function Base.match(node::Node, method::String, s::AbstractString)
    params = Dict()
    m = Base.match(node, params,  method, split_route(s), 1)
    return m, Dict(Symbol(k) => v for (k,v) in params)
end

function Base.match(node::Node, params, method, segments, i)
    # @show node.segment, i, segments
    if i > length(segments)
        if isempty(node.methods)
            return nothing
        end
        j = find(method, node.methods; by=x -> x.method, eq=(x, y) -> x == "*" || x == y)
        if j === nothing
            # we return missing here so we can return a 405 instead of 404
            # i.e. we matched the route, but there wasn't a matching method
            return missing
        else
            leaf = node.methods[j]
            # @show leaf.variables, segments
            if !isempty(leaf.variables)
                # we have variables to fill in
                for (i, v) in leaf.variables
                    params[v] = segments[i]
                end
            end
            return leaf.handler
        end
    end
    segment = segments[i]
    anymissing = false
    # first check for exact matches
    j = find(segment, node.exact; by=x -> x.segment)
    if j !== nothing
        # found an exact match, recurse
        m = Base.match(node.exact[j], params, method, segments, i + 1)
        anymissing = m === missing
        m = coalesce(m, nothing)
        # @show :exact, m
        if m !== nothing
            return m
        end
    end
    # check for conditional matches
    for node in node.conditional
        # @show node.segment.pattern, segment
        if match(node.segment.pattern, segment) !== nothing
            # matched a conditional node, recurse
            m = Base.match(node, params, method, segments, i + 1)
            anymissing = m === missing
            m = coalesce(m, nothing)
            if m !== nothing
                return m
            end
        end
    end
    if node.wildcard !== nothing
        m = Base.match(node.wildcard, params, method, segments, i + 1)
        anymissing = m === missing
        m = coalesce(m, nothing)
        if m !== nothing
            return m
        end
    end
    if node.doublestar !== nothing
        m = Base.match(node.doublestar, params, method, segments, length(segments) + 1)
        anymissing = m === missing
        m = coalesce(m, nothing)
        if m !== nothing
            return m
        end
    end
    return anymissing ? missing : nothing
end

function match_middleware(node::Node, params, method, segments, i)
    matches = []
    # @info "Node"
    # @show node.segment, i, segments
    if i > length(segments)
        if isempty(node.methods)
            return nothing
        end
        j = find(method, node.methods; by=x -> x.method, eq=(x, y) -> x == "*" || x == y)
        if j === nothing
            # we return missing here so we can return a 405 instead of 404
            # i.e. we matched the route, but there wasn't a matching method
            return missing
        else
            leaf = node.methods[j]
            # @show leaf.variables, segments
            if !isempty(leaf.variables)
                # we have variables to fill in
                for (i, v) in leaf.variables
                    params[v] = segments[i]
                end
            end
            push!(matches, leaf.handler)
            return matches
        end
    end
    segment = segments[i]
    anymissing = false
    # first check for exact matches
    j = find(segment, node.exact; by=x -> x.segment)
    # @info "Exact"
    if j !== nothing
        # found an exact match, recurse
        m = match_middleware(node.exact[j], params, method, segments, i + 1)
        anymissing = m === missing
        m = coalesce(m, nothing)
        # @show :exact, m
        if m !== nothing
            push!(matches, m...)
            # return m
        end
    end
    # @info "Conditional" node=node.conditional
    # check for conditional matches
    for node in node.conditional
        # @show node.segment.pattern, segment
        if match(node.segment.pattern, segment) !== nothing
            # matched a conditional node, recurse
            m = match_middleware(node, params, method, segments, i + 1)
            anymissing = m === missing
            m = coalesce(m, nothing)
            if m !== nothing
                push!(matches, m...)
                # return m
            end
        end
    end
    # @info "Wildcard" wildcard=node.wildcard
    if node.wildcard !== nothing
        m = match_middleware(node.wildcard, params, method, segments, i + 1)
        anymissing = m === missing
        m = coalesce(m, nothing)
        # @show :wildcard, m
        if m !== nothing
            push!(matches, m...)
            # return m
        end
    end
    # @info "Double Star"
    if node.doublestar !== nothing
        m = match_middleware(node.doublestar, params, method, segments, length(segments) + 1)
        anymissing = m === missing
        m = coalesce(m, nothing)
        if m !== nothing
            push!(matches, m...)
            # return m
        end
    end
    return matches
end

"""
    HTTP.register!(r::Node, [method,] path, handler)

Register a handler function that should be called when an incoming request matches `path`
and the optionally provided `method` (if not provided, any method is allowed). Can be used
to dynamically register routes.
The following path types are allowed for matching:
  * `/api/widgets`: exact match of static strings
  * `/api/*/owner`: single `*` to wildcard match any string for a single segment
  * `/api/widget/{id}`: Define a path variable `id` that matches any valued provided for this segment; path variables are available in the request context like `req.context[:params]["id"]`
  * `/api/widget/{id:[0-9]+}`: Define a path variable `id` that only matches integers for this segment
  * `/api/**`: double wildcard matches any number of trailing segments in the request path; must be the last segment in the path
"""
function register! end

function register!(n::Node, method, path, handler)
    method = uppercase(String(method))
    segments = if path == "/"
        ["/"]
    else
        map(segment, split(path, '/'; keepempty=false))
    end
    Base.insert!(n, Leaf(method, Tuple{Int,String}[], path, handler), segments, 1)
    return
end

function safe_which(fn, args)
    try
        which(fn, args)
    catch
        nothing
    end
end

function wrap_handler(handler)
    if !isnothing(safe_which(handler, Tuple{Any}))
        return HttpHandler(handler)
    elseif !isnothing(safe_which(handler, Tuple{Any, Any}))
        return Middleware(handler)
    else
        error("Unable to infer correct handler type")
    end
end

function split_route(s)
    if s == "/"
        ["/"]
    else
        split(s, '/'; keepempty=false)
    end
end

function Base.match(app, req::Request)
    url = URI(req.target)

    handler, params = match(app.paths, req.method, url.path)
    segments = split_route(url.path)
    # todo clean up this matching logic so it uses the 3 args version
    middelware = match_middleware(app.middleware, params, req.method, segments, 1)

    # needed for HTTP 0.9 compat
    if hasfield(Request, :context)
        req.context[:params] = params
    end
    # handler and be nothing or missing
    # nothing - didn't match a registered route
    # missing - matched the path, but method not supported

    if ismissing(middelware) || isnothing(middelware)
        middelware = Middleware[]
    end

    return handler, middelware
end

Base.match(app, stream::Stream) = Base.match(app, stream.message)


text:
function smolyak_grid_full(node_type::Function,d::S,mu::S) where {S<:Integer}

  T = typeof(1.0)

  multi_index        = generate_multi_index(d,mu)
  unique_multi_index = sort(unique(multi_index))
  unique_node_number = m_i(unique_multi_index)

  # Create base nodes to be used in the sparse grid

  base_nodes   = Array{Array{T,1},1}(undef,length(unique_node_number))
  base_weights = Array{Array{T,1},1}(undef,length(unique_node_number))
  @inbounds for i = 1:length(unique_node_number)
    base_nodes[i], base_weights[i] = node_type(unique_node_number[i])
  end

  # Select the relevant polynomials from the multi index

  ii = (sum(multi_index,dims=2) .>= max(d,mu+1)).*(sum(multi_index,dims=2) .<= d+mu)
  multi_index_full = zeros(S,sum(ii),size(multi_index,2))
  j = 1
  @inbounds for i = 1:size(multi_index,1)
    if ii[i] == true
      multi_index_full[j,:] = multi_index[i,:]
      j += 1
    end
  end

  # Construct the sparse grid from the nodes

  nodes = Array{T,2}(undef,determine_grid_size_full(multi_index_full))
  l = 1
  for j = 1:size(multi_index_full,1)
    new_nodes = base_nodes[multi_index_full[j,1]]  # Here new_nodes is a 1d array
    for i = 2:d
      new_nodes = combine_nodes(new_nodes,base_nodes[multi_index_full[j,i]])  # Here new_nodes becomes a 2d array
    end
    m = size(new_nodes,1)
    nodes[l:l+m-1,:] = new_nodes
    l += m
  end

  if d == 1
    nodes = nodes[:]
  end

  return nodes, multi_index_full

end

function smolyak_grid_full(node_type::Function,d::S,mu::S,domain::Union{Array{T,1},Array{T,2}}) where {S<:Integer, T<:AbstractFloat}

  if size(domain,2) != d
    error("domain is inconsistent with the number of dimensions")
  end

  nodes, multi_index_full = smolyak_grid_full(node_type,d,mu)

  nodes = scale_nodes(nodes,domain)

  return nodes, multi_index_full

end

function determine_grid_size_full(mi)

  temp = similar(mi)

  @inbounds for i = 1:size(mi,1)
    @inbounds for j = 1:size(mi,2)
      if mi[i,j] == 1
        temp[i,j] = 1
      else
        temp[i,j] = 2^(mi[i,j]-1)+1
      end
    end
  end

  s = 0
  @inbounds for i = 1:size(mi,1)
    t = 1
    @inbounds for j = 1:size(mi,2)
      t *= temp[i,j]
    end
    s += t
  end

  return (s, size(mi,2))

end

function master_index(multi_index::Array{S,2}) where {S<:Integer}

  temp_ind   = similar(multi_index)
  master_ind = zeros(S,size(multi_index,1),2)

  @inbounds for i in eachindex(multi_index)
    if multi_index[i] == 1
      temp_ind[i] = 1
    else
      temp_ind[i] = m_i(multi_index[i])
    end
  end

  master_ind[1,1] = 1
  master_ind[1,2] = prod(temp_ind[1,:])

  @inbounds for i = 2:size(master_ind,1)
    master_ind[i,1] = master_ind[i-1,1] + master_ind[i-1,2]
    master_ind[i,2] = prod(temp_ind[i,:])
  end

  return master_ind

end

function cheb_poly(order::S,x::T) where {S<:Integer, T<:AbstractFloat}

  p  = 1.0
  p1 = 0.0
  p2 = 0.0

  for i = 2:order+1
    if i == 2
	  p1, p = p, x
    else
      p2, p1 = p1, p
  	  p  = 2*x*p1-p2
    end
  end

  return p

end

function prod_cjs(max_grid::Union{Array{T,1},Array{T,2}},min_grid::Union{Array{T,1},Array{T,2}},poly_grid::Array{T,2}) where {T<:AbstractFloat}

  cjs = ones(size(poly_grid))

  @inbounds for i = 1:size(poly_grid,1)
    @inbounds for j = 1:size(poly_grid,2)
      if poly_grid[i,j] == max_grid[j] || poly_grid[i,j] == min_grid[j]
        cjs[i,j] *= 2
      end
    end
  end

  return prod(cjs,dims=2)

end

function compute_scale_factor(multi_index::Array{S,1}) where {S<:Integer}

  scale_factor = 1.0

  @inbounds for i = 1:length(multi_index)
    if multi_index[i] > 1
      scale_factor *= 2.0/(m_i(multi_index[i])-1)
    end
  end

  return scale_factor

end

function smolyak_weights_full(y_f::Array{T,1},grid::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {S<:Integer, T<:AbstractFloat}

  mi = sum(multi_index,dims=2)
  d  = size(multi_index,2)
  mu = maximum(mi)-d

  max_grid = maximum(grid,dims=1)
  min_grid = minimum(grid,dims=1)

  weights = Array{Array{T,1},1}(undef,size(multi_index,1))
  g_ind = master_index(multi_index)

  @inbounds for i = 1:size(g_ind,1) # This loops over the number of polynomials
    ws = zeros(g_ind[i,2])
    poly_grid = grid[g_ind[i,1]:g_ind[i,1]+g_ind[i,2]-1,:]
    poly_y    = y_f[g_ind[i,1]:g_ind[i,1]+g_ind[i,2]-1]
	  if size(grid,1) == 1 # This is to accommodate the mu = 0 case
	    cjs_prod = ones(size(poly_grid,1))
	  else
      cjs_prod = prod_cjs(max_grid,min_grid,poly_grid)
    end
    @inbounds for l = 1:g_ind[i,2] # This loops over the weights
      ll = CartesianIndices(Tuple(m_i(multi_index[i,:])))[l]
      @inbounds for j = 1:g_ind[i,2] # This loops over the nodes
        rhs_term = cheb_poly(ll[1]-1,poly_grid[j,1])*poly_y[j]
        @inbounds for k = 2:size(poly_grid,2) # This loops over the polynomial terms in the product
          rhs_term *= cheb_poly(ll[k]-1,poly_grid[j,k])
        end
        ws[l] += rhs_term/cjs_prod[j]
      end
      scale_factor = compute_scale_factor(multi_index[i,:])
      ws[l] = scale_factor*(1/cjs_prod[l])*ws[l]
    end
    weights[i] = (-1)^(d+mu-mi[i])*factorial(d-1)/(factorial(d+mu-mi[i])*factorial(-1-mu+mi[i]))*ws
  end

  return weights

end

function smolyak_weights_full(y_f::Array{T,1},grid::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Array{T,2}) where {S<:Integer, T<:AbstractFloat}

  d = size(multi_index,2)
  grid = copy(grid)
  for i = 1:d
    grid[:,i] = normalize_node(grid[:,i],domain[:,i])
  end

  weights = smolyak_weights_full(y_f,grid,multi_index)

  return weights

end

function smolyak_evaluate_full(weights::Array{Array{T,1},1},point::Array{T,1},multi_index::Array{S,2}) where {S<:Integer, T<:AbstractFloat}

  d = size(multi_index,2)

  evaluated_polynomials = zeros(size(multi_index,1))

  @inbounds for i = 1:size(multi_index,1) # This loops over the number of polynomials
  	@inbounds for l = 1:length(weights[i])
      ll = CartesianIndices(Tuple(m_i(multi_index[i:i,:])))[l]
      temp = weights[i][l]*cheb_poly(ll[1]-1,point[1])
      @inbounds for k = 2:d
  	    temp *= cheb_poly(ll[k]-1,point[k])
      end
      evaluated_polynomials[i] += temp
    end
    #evaluated_polynomials[i] *= (-1)^(d+mu-mi[i])*factorial(d-1)/(factorial(d+mu-mi[i])*factorial(-1-mu+mi[i]))
  end

  return sum(evaluated_polynomials)

end

function smolyak_evaluate_full(weights::Array{Array{T,1},1},point::Array{T,1},multi_index::Array{S,2},domain::Array{T,2}) where {S<:Integer, T<:AbstractFloat}

  d = size(multi_index,2)
  point = copy(point)
  for i = 1:d
    point[i] = normalize_node(point[i],domain[:,i])
  end

  estimate = smolyak_evaluate_full(weights,point,multi_index)

  return estimate

end

function deriv_cheb_poly(order::S,x::T) where {S<:Integer, T<:AbstractFloat}

  p0 = one(T)
  p1 = zero(T)
  p2 = zero(T)
  pd0 = zero(T)
  pd1 = zero(T)
  pd2 = zero(T)

  for i = 2:order+1
    if i == 2
	  p1, p0 = p0, x
	  pd1, pd0 = pd0, one(T)
    else
      p2, p1 = p1, p0
      pd2, pd1 = pd1, pd0
  	  p0  = 2*x*p1-p2
      pd0 = 2*p1+2*x*pd1-pd2
    end
  end

  return pd0

end

function smolyak_derivative_full(weights::Array{Array{T,1},1},point::Array{T,1},multi_index::Array{S,2},pos::S) where {S<:Integer, T<:AbstractFloat}

  mi = sum(multi_index,dims=2)
  d  = size(multi_index,2)
  mu = maximum(mi)-d

  evaluated_polynomials = zeros(size(multi_index,1))

  for i = 1:size(multi_index,1) # This loops over the number of polynomials
    for l = 1:length(weights[i])
      ll = CartesianIndices(Tuple(m_i(multi_index[i:i,:])))[l]
      if pos == 1
        temp = weights[i][l]*deriv_cheb_poly(ll[1]-1,point[1])
      else
        temp = weights[i][l]*cheb_poly(ll[1]-1,point[1])
      end
      for k = 2:d
        if k == pos
          temp *= deriv_cheb_poly(ll[k]-1,point[k])
        else
          temp *= cheb_poly(ll[k]-1,point[k])
        end
	    end
	    evaluated_polynomials[i] += temp
    end
	  #evaluated_polynomials[i] *= (-1)^(d+mu-mi[i])*factorial(d-1)/(factorial(d+mu-mi[i])*factorial(-1-mu+mi[i]))
  end

  evaluated_derivative = sum(evaluated_polynomials)

  return evaluated_derivative

end

function smolyak_derivative_full(weights::Array{Array{T,1},1},point::Array{T,1},multi_index::Array{S,2},domain::Array{T,2},pos::S) where {S<:Integer, T<:AbstractFloat}

  d = size(multi_index,2)
  point = copy(point)
  for i = 1:d
    point[i] = normalize_node(point[i],domain[:,i])
  end

  evaluated_derivative = smolyak_derivative_full(weights,point,multi_index,pos)

  return evaluated_derivatives

end

function smolyak_gradient_full(weights::Array{Array{T,1},1},point::Array{T,1},multi_index::Array{S,2}) where {S<:Integer, T<:AbstractFloat}

  d = size(multi_index,2)

  gradient = Array{T,2}(undef,1,d)
  for i = 1:d
    gradient[i] = smolyak_derivative_full(weights,point,multi_index,i)
  end

  return gradient

end

function smolyak_gradient_full(weights::Array{Array{T,1},1},point::Array{T,1},multi_index::Array{S,2},domain::Array{T,2}) where {S<:Integer, T<:AbstractFloat}

  d = size(multi_index,2)

  gradient = Array{T,2}(undef,1,d)
  for i = 1:d
    gradient[i] = smolyak_derivative_full(weights,point,multi_index,domain,i)
  end

  return gradient

end

###########################################################################

function smolyak_evaluate_full(weights::Array{T,1},multi_index_full::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_evaluate_full(weights,x,multi_index_full)

  end

  return goo

end

function smolyak_evaluate_full(weights::Array{T,1},multi_index_full::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_evaluate_full(weights,x,multi_index_full,domain)

  end

  return goo

end

function smolyak_derivative_full(weights::Array{T,1},multi_index_full::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_derivative_full(weights,x,multi_index_full)

  end

  return goo

end

function smolyak_derivative_full(weights::Array{T,1},multi_index_full::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_derivative_full(weights,x,multi_index_full,domain)

  end

  return goo

end

function smolyak_derivative_full(weights::Array{T,1},multi_index_full::Array{S,2},pos::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_derivative_full(weights,x,multi_index_full,pos)

  end

  return goo

end

function smolyak_derivative_full(weights::Array{T,1},multi_index_full::Array{S,2},domain::Union{Array{T,1},Array{T,2}},pos::Array{S,1}) where {T<:AbstractFloat,S<:Integer}

  function goo(x::Array{T,1}) where {T <: AbstractFloat}

    return smolyak_derivative_full(weights,x,multi_index_full,domain,pos)

  end

  return goo

end


text:
"""
    convert_WOAunits_to_unitful(s)

Converts the string from WOA's `units` attribute to a Unitful.jl unit.
"""
convert_to_Unitful(s) = UNITFUL_WOA[s]

#=====================================
Dictionary to translate
units from WOA
to units of Unitful
=====================================#

const UNITFUL_WOA = Dict(
    "micromoles_per_liter"    => u"μmol/l",
    "micromoles_per_kilogram" => u"1.025μmol/l", # 1 l ≈ 1.025 kg from ocean.ices.dk/Tools/UnitConversion.aspx
    "percent" => u"percent"
)