text:
#######################################################################
#
# An example of creating Excel Line charts with XlsxWriter.
#
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_line.xlsx")
	ws = add_worksheet!(wb)
	bold = add_format!(wb, Dict("bold"=> 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Number", "Batch 1", "Batch 2"]
	data = [
		[2 3 4 5 6 7];
		[10 40 50 20 10 50];
		[30 60 70 50 40 30];
	]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", data[1,:])
	write_column!(ws, "B2", data[2,:])
	write_column!(ws, "C2", data[3,:])

	# Create a new chart object. In this case an embedded chart.
	chart1 = add_chart!(wb, Dict("type"=> "line"))

	# Configure the first series.
	add_series!(chart1, Dict(
		"name"=>       "=Sheet1!\$B\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second series. Note use of alternative syntax to define ranges.
	add_series!(chart1, Dict(
		"name"=>       ["Sheet1", 0, 2],
		"categories"=> ["Sheet1", 1, 0, 6, 0],
		"values"=>     ["Sheet1", 1, 2, 6, 2],
	))

	# Add a chart title and some axis labels.
	set_title!(chart1, Dict("name"=> "Results of sample analysis"))
	set_x_axis!(chart1, Dict("name"=> "Test number"))
	set_y_axis!(chart1, Dict("name"=> "Sample length (mm)"))

	# Set an Excel chart style. Colors with white outline and shadow.
	set_style!(chart1, 10)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D2", chart1, Dict("x_offset"=> 25, "y_offset"=> 10))

	close(wb)

	isfile("chart_line.xlsx")
end

test()


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"

set -ex
${JULIA} -e 'using Pkg; Pkg.instantiate()'

export JULIA_LOAD_PATH="@"
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

import BangBang
import JSON
import Literate
import LiterateTest
import LoadAllPackages
import OnlineStats
import Random
using Documenter
using Transducers

EXAMPLE_PAGES = [
    "Tutorial: Missing values" => "tutorials/tutorial_missings.md",
    "Tutorial: Parallelism" => "tutorials/tutorial_parallel.md",
    "Parallel word count" => "tutorials/words.md",
    "Upgrade to new `|>` of Transducers.jl 0.4.39" => "howto/upgrade-to-ixf.md",
    "Empty result handling" => "howto/empty_result_handling.md",
    "Writing transducers" => "howto/transducers.md",
    "Writing reducibles" => "howto/reducibles.md",
    "Useful patterns" => "howto/useful_patterns.md",
]

LoadAllPackages.loadall(joinpath((@__DIR__), "Project.toml"))

function transducers_literate(;
    inputbase = joinpath(@__DIR__, "..", "examples"),
    outputbase = joinpath(@__DIR__, "src"),
    examples = EXAMPLE_PAGES,
    kwargs...,
)
    for (_, outpath) in examples
        name, = splitext(basename(outpath))
        inputfile = joinpath(inputbase, "$name.jl")
        outputdir = joinpath(outputbase, dirname(outpath))
        if name == "words"
            config = Dict()
        else
            config = LiterateTest.config()
        end
        Literate.markdown(
            inputfile,
            outputdir;
            config = config,
            documenter = true,
            kwargs...,
        )
    end
end

function transducers_rm_duplicated_docs()
    shareddocs =
        Docs.Binding.(Ref(Transducers), [:whenstart, :whencomplete, :whencombine]) .=>
            Ref(Docs.Binding(Transducers, :wheninit))
    for (dup, canonical) in shareddocs
        @info "Simplifying docstring: $dup => $canonical"
        pop!(Docs.meta(dup.mod), dup, nothing)
        canstr = only(values(Docs.meta(canonical.mod)[canonical].docs))
        txt = "See [`$(canonical.var)`](@ref $canonical)"
        @eval dup.mod $Docs.@doc $txt $(dup.var)
    end
end

function transducers_redirection_mapping()
    mapping = [
        # old page => new page
        "manual" => "reference/manual",
        "interface" => "reference/interface",
    ]

    old_examples = [
        "tutorials/tutorial_missings",
        "tutorials/tutorial_parallel",
        "tutorials/words",
        "howto/empty_result_handling",
        "howto/transducers",
        "howto/reducibles",
    ]

    for e in old_examples
        push!(mapping, joinpath("examples", basename(e)) => e)
    end

    return mapping
end

function transducers_make_redirections(;
    mapping = transducers_redirection_mapping(),
    build = joinpath((@__DIR__), "build"),
)
    for (old, new) in mapping
        oldpath = joinpath(build, old)
        newpath = joinpath(build, new)
        relurl = relpath(newpath, oldpath)
        html = """<meta http-equiv="refresh" content="0; url=$relurl"/>"""
        mkpath(oldpath)
        write(joinpath(oldpath, "index.html"), html)
    end
end

function should_push_preview(event_path = get(ENV, "GITHUB_EVENT_PATH", nothing))
    event_path === nothing && return false
    event = JSON.parsefile(event_path)
    pull_request = get(event, "pull_request", nothing)
    pull_request === nothing && return false
    labels = [x["name"] for x in pull_request["labels"]]
    # https://developer.github.com/v3/activity/events/types/#pullrequestevent
    yes = "push_preview" in labels
    if yes
        @info "Trying to push preview as label `push_preview` is specified." labels
    else
        @info "Not pushing preview as label `push_preview` is not specified." labels
    end
    return yes
end

Random.seed!(1234)
transducers_rm_duplicated_docs()
transducers_literate()

examples = EXAMPLE_PAGES
strict = get(ENV, "CI", "false") == "true"
doctest = get(ENV, "CI", "false") == "true"

tutorials = filter(((_, path),) -> startswith(path, "tutorials/"), examples)
howto = filter(((_, path),) -> startswith(path, "howto/"), examples)
@assert issetequal(union(tutorials, howto), examples)
@info "`makedocs` with" strict doctest
makedocs(;
    modules = [Transducers],
    pages = [
        "Home" => "index.md",
        "Reference" =>
            ["Manual" => "reference/manual.md", "Interface" => "reference/interface.md"],
        "Tutorials" => tutorials,
        "How-to guides" => howto,
        "Explanation" => [
            "Parallelism" => "parallelism.md",  # TODO: merge this to index.md
            "Comparison to iterators" => "explanation/comparison_to_iterators.md",
            "Glossary" => "explanation/glossary.md",
            "State machines" => "explanation/state_machines.md",
            hide("Internals" => "explanation/internals.md"),
        ],
    ],
    repo = "https://github.com/JuliaFolds/Transducers.jl/blob/{commit}{path}#L{line}",
    sitename = "Transducers.jl",
    authors = "Takafumi Arakaki",
    strict = strict,
    doctest = doctest,
)

transducers_make_redirections()
deploydocs(;
    repo = "github.com/JuliaFolds/Transducers.jl",
    push_preview = should_push_preview(),
)

ext: jl


text:
@testset "brickgrids" begin
    TAs = ((Float64,  Array), (Float32,  Array), (BigFloat, Array))
    if CUDA.has_cuda_gpu()
        TAs = (TAs..., (Float32, CuArray))
    end

    for (T, A) in TAs
        cellsandcoordinates =
            (
             (LobattoCell{T, A}(5), (-one(T):1//2:one(T),)),
             (LobattoCell{T, A}(4), ((1, 2, 5, 8),)),
             (LobattoCell{T, A}(3, 5),
              (-one(T):1//2:one(T), zero(T):1//4:one(T))),
             (LobattoCell{T, A}(3, 5, 2),
              (-one(T):1//2:one(T), zero(T):1//4:one(T), 1:3))
            )
        for (cell, coordinates) in cellsandcoordinates
            grid = brickgrid(cell, coordinates)

            @test floattype(grid) == T
            @test arraytype(grid) <: A
            @test celltype(grid) == typeof(cell)
            @test ndims(grid) == length(coordinates)
            @test size(grid) == length.(coordinates) .- 1
            @test length(grid) == prod(length.(coordinates) .- 1)

            @test floattype(typeof(grid)) == T
            @test arraytype(typeof(grid)) <: A
            @test celltype(typeof(grid)) == typeof(cell)
            @test ndims(typeof(grid)) == length(coordinates)
            @test size(typeof(grid)) == length.(coordinates) .- 1
            @test length(typeof(grid)) == prod(length.(coordinates) .-1)

            cell = referencecell(grid)
            @test size(points(grid)) == (length(cell), length(grid))

            x = adapt(Array, points(grid))
            faceindices⁻, faceindices⁺ = adapt.(Array, faceindices(grid))
            @test isapprox(x[faceindices⁻], x[faceindices⁺], atol=100eps(T))
            matfaces = Bennu.materializefaces(cell, connectivity(grid))
            for (f, g) in zip(matfaces, faces(grid))
                @test size(f) == size(g)
            end
            for (b, d) in zip(adapt(Array,
                                    Bennu.materializeboundaryfaces(cell,
                                                                   matfaces)),
                              adapt(Array, boundaryfaces(grid)))
                @test (b ≠ 0 && d ≠ 0) || (b == d == 0)
            end

            Ds = derivatives(cell)

            for (i, D) in enumerate(Ds)
                otherdims = filter(!isequal(i), 1:length(Ds))
                D = adapt(Array, D)
                dx = reshape(D * x, (length(cell), size(grid)...))
                pdx = vec(minimum(dx, dims=(1, (otherdims .+ 1)...)))
                qdx = vec(maximum(dx, dims=(1, (otherdims .+ 1)...)))
                cdx = diff(collect(coordinates[i]))./2

                @test isapprox(pdx, qdx, atol=100eps(T))
                @test isapprox(getindex.(pdx, i), cdx, atol=100eps(T))
                for j in otherdims
                    @test isapprox(getindex.(pdx, j), zeros(T, length(pdx)),
                                   atol=100eps(T))
                end
            end

            if T != BigFloat
                mktempdir() do tmp
                    outfiles = vtk_grid(joinpath(tmp, "grid_$(T)_$(A)"), grid;
                                        append=false, ascii=true) do vtk
                        vtk["Time"] = 37.0
                    end
                    for file in outfiles
                        @test EzXML.readxml(file) isa EzXML.Document
                    end
                end
            end
        end
    end

    grid = brickgrid(LobattoCell(3), (-1:1,); periodic=(true,))
    @test faceindices(grid) == ([1 4; 3 6], [6 3; 4 1])
    @test boundaryfaces(grid) == [0  0;  0  0]

    grid = brickgrid(LobattoCell(3, 4), (-1:1, -1:1); periodic=(true, false))
    @test size(faces(grid)[1]) == (16, 10)
    @test size(faces(grid)[2]) == (16, 6)
    @test boundaryfaces(grid) == [0  0  0  0
                                  0  0  0  0
                                  3  3  0  0
                                  0  0  4  4]

    grid = brickgrid(LobattoCell(3, 4, 2), (-1:1, -1:1, -1:1);
                     periodic=(true, false, true))
    @test size(faces(grid)[1]) == (48, 28)
    @test size(faces(grid)[2]) == (96, 32)
    @test size(faces(grid)[3]) == (64, 12)
    @test boundaryfaces(grid) == [0  0  0  0  0  0  0  0
                                  0  0  0  0  0  0  0  0
                                  3  3  0  0  3  3  0  0
                                  0  0  4  4  0  0  4  4
                                  0  0  0  0  0  0  0  0
                                  0  0  0  0  0  0  0  0]
end

@testset "cubedshell" begin
    TAs = ((Float64,  Array), (Float32,  Array), (BigFloat, Array))
    if CUDA.has_cuda_gpu()
        TAs = (TAs..., (Float32, CuArray))
    end

    for (T, A) in TAs
        cell = LobattoCell{T, A}(3, 3)
        ncells_h_panel = 2

        (vert, conn) = Bennu.cubedshellconnectivity(cell, 1, ncells_h_panel);

        vert = adapt(Array, vert)
        conn = adapt(Array, conn)
        f1 = [[-1,  1, -1],
              [-1,  0, -1],
              [-1,  1,  0],
              [-1,  0,  0]]
        f2 = [[ 1, -1, -1],
              [ 1,  0, -1],
              [ 1, -1,  0],
              [ 1,  0,  0]]
        f3 = [[-1, -1, -1],
              [ 0, -1, -1],
              [-1, -1,  0],
              [ 0, -1,  0]]
        f4 = [[ 1,  1, -1],
              [ 0,  1, -1],
              [ 1,  1,  0],
              [ 0,  1,  0]]
        f5 = [[-1,  1, -1],
              [ 0,  1, -1],
              [-1,  0, -1],
              [ 0,  0, -1]]
        f6 = [[-1, -1,  1],
              [ 0, -1,  1],
              [-1,  0,  1],
              [ 0,  0,  1]]

        @test vert[[conn[1, 1, 1]...]] == f1
        @test vert[[conn[2, 1, 1]...]] == f1 .+ [[ 0, -1,  0]]
        @test vert[[conn[1, 2, 1]...]] == f1 .+ [[ 0,  0,  1]]
        @test vert[[conn[2, 2, 1]...]] == f1 .+ [[ 0, -1,  1]]

        @test vert[[conn[1, 1, 2]...]] == f2
        @test vert[[conn[2, 1, 2]...]] == f2 .+ [[ 0,  1,  0]]
        @test vert[[conn[1, 2, 2]...]] == f2 .+ [[ 0,  0,  1]]
        @test vert[[conn[2, 2, 2]...]] == f2 .+ [[ 0,  1,  1]]

        @test vert[[conn[1, 1, 3]...]] == f3
        @test vert[[conn[2, 1, 3]...]] == f3 .+ [[ 1,  0,  0]]
        @test vert[[conn[1, 2, 3]...]] == f3 .+ [[ 0,  0,  1]]
        @test vert[[conn[2, 2, 3]...]] == f3 .+ [[ 1,  0,  1]]

        @test vert[[conn[1, 1, 4]...]] == f4
        @test vert[[conn[2, 1, 4]...]] == f4 .+ [[-1,  0,  0]]
        @test vert[[conn[1, 2, 4]...]] == f4 .+ [[ 0,  0,  1]]
        @test vert[[conn[2, 2, 4]...]] == f4 .+ [[-1,  0,  1]]

        @test vert[[conn[1, 1, 5]...]] == f5
        @test vert[[conn[2, 1, 5]...]] == f5 .+ [[ 1,  0,  0]]
        @test vert[[conn[1, 2, 5]...]] == f5 .+ [[ 0, -1,  0]]
        @test vert[[conn[2, 2, 5]...]] == f5 .+ [[ 1, -1,  0]]

        @test vert[[conn[1, 1, 6]...]] == f6
        @test vert[[conn[2, 1, 6]...]] == f6 .+ [[ 1,  0,  0]]
        @test vert[[conn[1, 2, 6]...]] == f6 .+ [[ 0,  1,  0]]
        @test vert[[conn[2, 2, 6]...]] == f6 .+ [[ 1,  1,  0]]
    end
end

@testset "cubedsphere" begin
    TAs = ((Float64,  Array), (Float32,  Array), (BigFloat, Array))
    if CUDA.has_cuda_gpu()
        TAs = (TAs..., (Float32, CuArray))
    end

    for (T, A) in TAs
        cell = LobattoCell{T, A}(3, 3, 4)
        vert_coord = one(T):1//2:4one(T)
        ncells_h_panel = 2

        grid = cubedspheregrid(cell, vert_coord, ncells_h_panel)

        # Check the first level connectivity
        h_offset = 6 * (ncells_h_panel + 1)^2

        connectivity = adapt(Array, grid.connectivity)
        @test connectivity[1, 1, 1, 1][1:4] == connectivity[1, 1, 1, 1][5:8] .- h_offset == (1, 2, 4, 5)
        @test connectivity[1, 2, 1, 1][1:4] == connectivity[1, 2, 1, 1][5:8] .- h_offset == (2, 3, 5, 6)
        @test connectivity[1, 1, 2, 1][1:4] == connectivity[1, 1, 2, 1][5:8] .- h_offset == (4, 5, 7, 8)
        @test connectivity[1, 2, 2, 1][1:4] == connectivity[1, 2, 2, 1][5:8] .- h_offset == (5, 6, 8, 9)

        @test connectivity[1, 1, 1, 2][1:4] == connectivity[1, 1, 1, 2][5:8] .- h_offset == (10, 11, 13, 14)
        @test connectivity[1, 2, 1, 2][1:4] == connectivity[1, 2, 1, 2][5:8] .- h_offset == (11, 12, 14, 15)
        @test connectivity[1, 1, 2, 2][1:4] == connectivity[1, 1, 2, 2][5:8] .- h_offset == (13, 14, 16, 17)
        @test connectivity[1, 2, 2, 2][1:4] == connectivity[1, 2, 2, 2][5:8] .- h_offset == (14, 15, 17, 18)

        @test connectivity[1, 1, 1, 3][1:4] == connectivity[1, 1, 1, 3][5:8] .- h_offset == (3, 20, 6, 23)
        @test connectivity[1, 2, 1, 3][1:4] == connectivity[1, 2, 1, 3][5:8] .- h_offset == (20, 10, 23, 13)
        @test connectivity[1, 1, 2, 3][1:4] == connectivity[1, 1, 2, 3][5:8] .- h_offset == (6, 23, 9, 26)
        @test connectivity[1, 2, 2, 3][1:4] == connectivity[1, 2, 2, 3][5:8] .- h_offset == (23, 13, 26, 16)

        @test connectivity[1, 1, 1, 4][1:4] == connectivity[1, 1, 1, 4][5:8] .- h_offset == (12, 29, 15, 32)
        @test connectivity[1, 2, 1, 4][1:4] == connectivity[1, 2, 1, 4][5:8] .- h_offset == (29, 1, 32, 4)
        @test connectivity[1, 1, 2, 4][1:4] == connectivity[1, 1, 2, 4][5:8] .- h_offset == (15, 32, 18, 35)
        @test connectivity[1, 2, 2, 4][1:4] == connectivity[1, 2, 2, 4][5:8] .- h_offset == (32, 4, 35, 7)

        @test connectivity[1, 1, 1, 5][1:4] == connectivity[1, 1, 1, 5][5:8] .- h_offset == (1, 29, 2, 41)
        @test connectivity[1, 2, 1, 5][1:4] == connectivity[1, 2, 1, 5][5:8] .- h_offset == (29, 12, 41, 11)
        @test connectivity[1, 1, 2, 5][1:4] == connectivity[1, 1, 2, 5][5:8] .- h_offset == (2, 41, 3, 20)
        @test connectivity[1, 2, 2, 5][1:4] == connectivity[1, 2, 2, 5][5:8] .- h_offset == (41, 11, 20, 10)

        @test connectivity[1, 1, 1, 6][1:4] == connectivity[1, 1, 1, 6][5:8] .- h_offset == (9, 26, 8, 50)
        @test connectivity[1, 2, 1, 6][1:4] == connectivity[1, 2, 1, 6][5:8] .- h_offset == (26, 16, 50, 17)
        @test connectivity[1, 1, 2, 6][1:4] == connectivity[1, 1, 2, 6][5:8] .- h_offset == (8, 50, 7, 35)
        @test connectivity[1, 2, 2, 6][1:4] == connectivity[1, 2, 2, 6][5:8] .- h_offset == (50, 17, 35, 18)

        # Check the higher levels
        for z = 1:length(vert_coord) - 2
            for (a, b) in zip(connectivity[z, :, :, :],  connectivity[z + 1, :, :, :])
                @test a .+ h_offset == b
            end
        end

        # Check the vertices are at the right levels
        vertices = adapt(Array, grid.vertices)
        for z = 1:length(vert_coord)
            for p in vertices[:, :, :, z]
                @test norm(Bennu.cubespherewarp(p)) ≈ vert_coord[z]
            end
        end
    end
end


text:
using DotOverload
using Base.Test

type TT
  F::Int64
end

# basic assertions

@test macroexpand( :( @dotted a.f )) == :(getMember(a,:f))
@test macroexpand( :( @dotted a.b.c = b.g += e.f )) == :(setMember!(getMember(a,:b),:c,setMember!(b,:g,getMember(b,:g) + getMember(e,:f))))
@test macroexpand( :( @dotted (a).b )) == :(getMember(a,:b))

# there are some parsing changes between 0.4 - 0.5?
if (VERSION >= v"0.5.0-dev")
  @test macroexpand( :( @dotted a.(b) )) == :(getMember(a,(b,)))
  @test macroexpand( :( @dotted anything.(10) )) == :(getMember(anything,(10,)))
else
  @test macroexpand( :( @dotted a.(b) )) == :(getMember(a,b))
  @test macroexpand( :( @dotted anything.(10) )) == :(getMember(anything,10))
end

@test macroexpand( :( @dotted "weird".stuff )) == :(getMember("weird",:stuff))
@test macroexpand( :( @dotted "weird"."stuff" )) == :(getMember("weird","stuff"))

# basic usage

DotOverload.getMember(t, m::Symbol) = begin
  println("intercepted!")
  return getfield(t, m)
end

@dotted function dontDoThatAtHome()
  a = TT(42)
  return a.F
end

@test 42 == dontDoThatAtHome() #result: 42, prints: "intercepted"

# code equality

DotOverload.getMember(t, m::Symbol) = Base.getfield(t, m)
DotOverload.setMember!(t, m::Symbol, v) = Base.setfield!(t,m,v)

function test1()
  a = TT(3)
  return a.F
end

@dotted function test2()
  a = TT(3)
  return a.F
end

@test @code_lowered(test1()) != @code_lowered(test2())
@test @code_llvm(test1()) == @code_llvm(test2())

# lua table like Dict

@dotted function test3()
  getMember(t::Dict, k) = t[k]
  setMember!(t::Dict, k::Symbol, v) = t[k] = v

  d = Dict()
  d.key1 = 1
  d.key2 = 2

  return d.key1 + d.key2
end

@test test3() == 3

# more examples

@dotted function test4()
  function getMember(t, k)
    return getfield(t, k)
  end
  function setMember!(t, m, v)
    return Base.setfield!(t,m,v)
  end
  function getMember(t::Dict, k)
    v = nothing
    try
      v = t[k]
      try
        v = round(v*2)
      catch
      end
    catch
      t[k] = Dict()
      v = t[k]
    end
    return v
  end
  function setMember!(t::Dict, k::Symbol, v)
    t[k] = v*2
  end

  a = 10.5
  b = Dict()

  b.this.is.kinda.kool = a

  return b.this.is.kinda.kool
end

@test test4() == 42

# a module!

@dotted module m

using DotOverload
DotOverload.getMember(t::Dict, k::Symbol) = t[k]
DotOverload.setMember!(t::Dict, k::Symbol, v) = t[k] = v

a = Dict()
a[:stuff] = 25
a.stuff *= 2

foo() = a.stuff * 2

end

@test m.foo() == 100

# REPL example (0.4):
#
# julia> using DotOverload
#
# julia> v = @dotted function() a.b = 32; a.b end
# (anonymous function)
#
# julia> a = Dict()
# Dict{Any,Any} with 0 entries
#
# julia> v()
# ERROR: type Dict has no field b
#  in setMember! at C:\SN\prj\julia\DotOverload\src\DotOverload.jl:128
#  in anonymous at none:1
#
# julia> DotOverload.getMember(t::Dict, k) = t[k]
# getMember (generic function with 2 methods)
#
# julia> DotOverload.setMember!(t::Dict, k, v) = t[k] = v
# setMember! (generic function with 2 methods)
#
# julia> v()
# 32
#
# julia>


text:
using JGraphs
using Documenter

DocMeta.setdocmeta!(JGraphs, :DocTestSetup, :(using JGraphs); recursive = true)

makedocs(;
    modules = [JGraphs],
    authors = "Giovanni Puccetti",
    repo = "https://github.com/gpucce/JGraphs.jl/blob/{commit}{path}#{line}",
    sitename = "JGraphs.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://gpucce.github.io/JGraphs.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md",
        "References" => "references.md",
    ],
)

deploydocs(; repo = "github.com/gpucce/JGraphs.jl", devbranch = "main")


text:
using ProcessFlows
using Test

@testset "ProcessFlows.jl" begin
    # Write your tests here.
end

max_stars_repo_name: pbayer/pflow

max_issues_repo_name: pbayer/pflow

max_forks_repo_name: pbayer/pflow


text:
using CSV, Clustering, Plots

client_locations = CSV.read("clients.csv", header=true)
k_coordinates = (Matrix{Float64}(client_locations))'

init = scatter(client_locations.X, client_locations.Y, show=true)

result = kmeans(k_coordinates, 2, init=:kmpp, maxiter=1000, display=:iter)

res = scatter(client_locations.X, client_locations.Y, marker_z = result.assignments,
        color=:lightrainbow, legend=false)

a = assignments(result)
c = counts(result)
M = result.centers

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
using MultipleViewGeometry,Test, LinearAlgebra
using MultipleViewGeometry.ModuleTypes

𝐊 = Matrix{Float64}(I, 3, 3)
𝐑 = Matrix{Float64}(I, 3, 3)
𝐭 = [1.0, 1.0, 1.0]

@test construct(ProjectionMatrix(),𝐊,𝐑,𝐭) == [Matrix{Float64}(I, 3, 3) -ones(3)]

# 1. Construct a Fundamental matrix from Camera matrices.
# 2. Construct Projection matrices from the Fundamental matrix.
# 3. Construct a Fundamental matrix from the projection matrices.
# 4. The Fundamental matrices in step 2 and 3 should be equivalent up to sign
#    and scale.
𝐊₁ = Matrix{Float64}(I, 3, 3)
𝐑₁ = Matrix{Float64}(I, 3, 3)
𝐭₁ = [1.0, 1.0, 1.0]
𝐊₂ = Matrix{Float64}(I, 3, 3)
𝐑₂ = Matrix{Float64}(I, 3, 3)
𝐭₂ = [2.0, 2.0, 2.0]
𝐏₁ = construct(ProjectionMatrix(),𝐊₁,𝐑₁,𝐭₁)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂,𝐑₂,𝐭₂)
𝐅 = construct(FundamentalMatrix(),𝐊₁,𝐑₁,𝐭₁,𝐊₂,𝐑₂,𝐭₂)

𝐎₁, 𝐎₂ = construct(ProjectionMatrix(),𝐅)
𝐅₂ = construct(FundamentalMatrix(),𝐎₁, 𝐎₂)

# Ensure the matrices have the same scale and sign before comparing them.
𝐅 = 𝐅 / norm(𝐅)
𝐅 = 𝐅 / sign(𝐅[1,2])

𝐅₂ = 𝐅₂ / norm(𝐅₂)
𝐅₂ = 𝐅₂ / sign(𝐅₂[1,2])

@test 𝐅 ≈ 𝐅₂


text:
# canonical factor graph examples useful for development and learning.

export generateCanonicalFG_ZeroPose, generateCanonicalFG_Hexagonal, generateCanonicalFG_TwoPoseOdo, generateCanonicalFG_Circle
export buildGraphChain!
export warmUpSolverJIT

"""
    $SIGNATURES

Lower level utility function for simulated canonical graph generation, where a pose variable and factor is added using this function.

Notes
- Also adds the `refKey=:simulated` (default) PPE as in graph record of the simulated values.
- Use callback `postpose_cb(g::AbstractDFG,lastpose::Symbol)` to call user operations after each pose step.

Related

[`generateCanonicalFG_ZeroPose`](@ref), [`generateCanonicalFG_Helix2D!`](@ref)
"""
function _addPoseCanonical!(fg::AbstractDFG,
                            prevLabel::Symbol,
                            posecount::Int, # can be overriden with genLabel
                            factor::AbstractFactor;
                            poseRegex::Regex=r"x\d+",
                            genLabel = Symbol(match(r"[A-Za-z]+", poseRegex.pattern).match, posecount),
                            srcType::Type{<:InferenceVariable} = getVariableType(fg, prevLabel) |> typeof,
                            poseType::Type{<:InferenceVariable} = srcType, # control destination type TODO simplify
                            graphinit::Bool=false,
                            solvable::Integer=1,
                            inflation::Real=getSolverParams(fg).inflation,
                            variableTags::AbstractVector{Symbol}=Symbol[],
                            factorTags::AbstractVector{Symbol}=Symbol[],
                            doRef::Bool=true,
                            refKey::Symbol = :simulated,
                            overridePPE = doRef ? nothing : zeros(getDimension(poseType)),
                            postpose_cb::Function=(fg_,latestpose)->()  )
  #
  # calculate and add the reference value
  isAlready, simPPE, = IIF._checkVariableByReference(fg, prevLabel, poseRegex, poseType, factor, refKey=refKey, srcType=srcType, doRef=doRef, overridePPE=overridePPE)

  # dispatch on prior or binary factor
  _getlabels(fact::AbstractPrior) = [genLabel;]
  _getlabels(fact::AbstractRelative) = [prevLabel; genLabel]

  # add new pose variable
  v_n = addVariable!(fg, genLabel, poseType; solvable, tags=variableTags )
  addFactor!(fg, _getlabels(factor), factor; graphinit, solvable, tags=factorTags, inflation )

  # store simulated PPE for future use
  # ppe = DFG.MeanMaxPPE(refKey, simPPE, simPPE, simPPE)
  doRef ? setPPE!(v_n, refKey, typeof(simPPE), simPPE) : nothing

  # user callback in case something more needs to be passed down
  postpose_cb(fg, genLabel)

  # return the new variable
  return v_n
end

"""
    $SIGNATURES

Generate a canonical factor graph with a Pose2 `:x0` and MvNormal with covariance `P0`.

Notes
- Use e.g. `varType=Point2` to change from the default variable type `Pose2`.
- Use `priorArgs::Tuple` to override the default input arguments to `priorType`.
- Use callback `postpose_cb(g::AbstractDFG,lastpose::Symbol)` to call user operations after each pose step.
"""
function generateCanonicalFG_ZeroPose(; varType::Type{<:InferenceVariable}=Pose2,
                                        graphinit = nothing,
                                        solverParams::SolverParams=SolverParams(),
                                        dfg::AbstractDFG = LightDFG{SolverParams}(;solverParams),
                                        doRef::Bool=true,
                                        useMsgLikelihoods::Bool=getSolverParams(dfg).useMsgLikelihoods,
                                        label::Symbol=:x0,
                                        priorType::Type{<:AbstractPrior}=DFG._getPriorType(varType),
                                        μ0::AbstractVector{<:Real}= zeros(getDimension(varType)),
                                        Σ0::AbstractMatrix{<:Real}= diagm(0.01*ones(getDimension(varType))),
                                        priorArgs::Tuple = (MvNormal(μ0, Σ0),),
                                        solvable::Integer=1,
                                        variableTags::AbstractVector{Symbol}=Symbol[],
                                        factorTags::AbstractVector{Symbol}=Symbol[],
                                        postpose_cb::Function=(fg_,latestpose)->()  )
  #
  (graphinit isa Nothing) ? nothing : @error("generateCanonicalFG_ZeroPose keyword graphinit obsolete, use solverParams(graphinit=..) instead.")

  # only add the first variable if none others exist
  if !exists(dfg, label)
    # generate a default prior
    prpo = priorType(priorArgs...)
    # add the variable and prior with canonical helper function
    _addPoseCanonical!( dfg, label, 0, prpo; genLabel=label, solvable,
                        srcType=varType, variableTags, factorTags, graphinit = solverParams.graphinit,
                        doRef, postpose_cb )
    #
  else
    @warn "$label already exists in the factor graph, no new variables added."
  end

  return dfg
end

"""
    $SIGNATURES

Build a chain of variables with binary factors.

Notes:
- `fctData[1]` aligns with the first variable (likely `:x0` for default empty `dfg`), and
  - `:x0` has no previous variable/factor data.

DevNotes
- TODO Consolidate with [`IIF._buildGraphByFactorAndTypes!`](@ref) and wrap in RoME?
"""
function buildGraphChain!(  fctData::AbstractVector = [MvNormal([10;0;0.],diagm(0.1.*ones(3))) for _ in 1:3],
                            fctType::Type{<:AbstractRelative} = Pose2Pose2,
                            preFct_args_cb::Function = (fg_,data)->(data.currData,);
                            stopAfter::Integer=2^(Sys.WORD_SIZE-1)-1,
                            varType::IIF.InstanceType{InferenceVariable} = Pose2,
                            solverParams::SolverParams = SolverParams(),
                            solvable::Integer = 1,
                            fctKwargs::NamedTuple = (;),
                            varTags::AbstractVector{Symbol} = Symbol[],
                            fctTags::AbstractVector{Symbol} = Symbol[],
                            postpose_cb::Function = (fg_, lastpose) -> nothing,
                            doRef::Bool=true,
                            dfg::AbstractDFG = generateCanonicalFG_ZeroPose(; varType, solverParams, solvable, doRef, postpose_cb),
                            inflation_fct::Real = getSolverParams(dfg).inflation,
                            varRegex::Regex = r"x\d+",
                            varLast::Symbol = sortDFG(ls(dfg, varRegex))[end],
                            varCount::Integer = match(r"\d+", string(varLast)).match |> x->parse(Int,x),
                            varPrefix::Symbol = match(r"[a-zA-Z_]+", varRegex.pattern).match |> Symbol  )
  #

  for i in 1:length(fctData)
    i <= stopAfter ? nothing : break
    varCount += 1
    varCurr = Symbol(varPrefix, varCount)
    # prep small container for data
    prev_ = 1 <= i-1 ? fctData[i-1] : nothing
    next_ = i+1 <= length(fctData) ? fctData[i+1] : nothing
    dat_ = (;varCurr=varCurr, varLast=varLast, prevData=prev_, currData=fctData[i], nextData=next_)
    _addPoseCanonical!( dfg,
                        varLast,
                        varCount,
                        fctType(preFct_args_cb(dfg,dat_)...; fctKwargs...);
                        doRef=doRef,
                        poseRegex=varRegex,
                        genLabel=varCurr,
                        srcType=varType,
                        solvable=solvable,
                        inflation=inflation_fct,
                        variableTags=varTags,
                        factorTags=fctTags,
                        postpose_cb=postpose_cb )
    #
    varLast = varCurr
  end

  dfg
end

"""
    $SIGNATURES

Build a basic factor graph in Pose2 with two `Pose2` and one landmark `Point2` variables,
along with `PriorPose2` on `:x0` and `Pose2Pose2` to `:x1`.  Also a `Pose2Point2BearingRange`
to landmark `:l1`.

DevNotes
- TODO standardize to latest and greatest generator api pattern

See also: [`buildGraphChain!`](@ref)
"""
function generateCanonicalFG_TwoPoseOdo(;solverParams::SolverParams = SolverParams(),
                                        varType::Type{<:Pose2}=Pose2,
                                        solvable::Integer = 1,
                                        doRef::Bool=true,
                                        postpose_cb::Function = (fg_, lastpose) -> nothing,
                                        dfg::AbstractDFG=generateCanonicalFG_ZeroPose(; varType, solverParams, solvable, doRef, postpose_cb),
                                        addlandmark::Bool=true,
                                        fg=nothing,
                                        graphinit=nothing )
  #
  (graphinit isa Nothing) ? nothing : @error("generateCanonicalFG_TwoPoseOdo keyword graphinit obsolete, use solverParams=SolverParams(graphinit=..) instead.")
  (fg isa Nothing) ? nothing : @error("generateCanonicalFG_TwoPoseOdo keyword fg obsolete, use dfg= instead.")

  buildGraphChain!( [MvNormal([10.0;0;0.0],diagm([1.0;1.0;0.01]));];
                    dfg,
                    varType )

  if addlandmark
    addVariable!(dfg, :l1, Point2)
    addFactor!( dfg, [:x1;:l1], Pose2Point2BearingRange(Normal(0.0,0.01), Normal(20.0, 1.0)), graphinit=solverParams.graphinit )
  end

  # return the built graph
  return dfg
end

#


text:
using Surrogates

#1D

obj_1D = x-> 2*x+1
a = 0.0
b = 10.0
x = sample(5,a,b,SobolSample())
y = obj_1D.(x)
my_svm_1D = SVMSurrogate(x,y,a,b)
val = my_svm_1D(5.0)
add_point!(my_svm_1D,3.1,7.2)
add_point!(my_svm_1D,[3.2,3.5],[7.4,8.0])

#ND
obj_N = x -> x[1]^2*x[2]
lb = [0.0,0.0]
ub = [10.0,10.0]
x = sample(100,lb,ub,UniformSample())
y = obj_N.(x)
my_svm_ND = SVMSurrogate(x,y,lb,ub)
val = my_svm_ND((5.0,1.2))
add_point!(my_svm_ND,(1.0,1.0),1.0)
add_point!(my_svm_ND,[(1.2,1.2),(1.5, 1.5)],[1.728,3.375])


text:
"""
Logs metrics to the console output. By default it will only log at the end of an epoch
and log the epoch, step, loss and validation loss. This is also the default configuration
when you run fit! without specifying your own meters.

Example:

    meter = ConsoleMeter([:loss, :val_accuracy, :accuracy]; epochOnly=false)
    fit!(workout, data, val_data; epochs=5, meters=[meter])

"""
mutable struct ConsoleMeter <: Meter
    throttle::Float64
    next::Float64
    metricnames::Vector{Symbol}
    epochOnly

    function ConsoleMeter(metrics=[:loss, :val_loss]; throttle=1.0, epochOnly=true)
        new(throttle, 0.0, metrics, epochOnly)
    end
end

function (meter::ConsoleMeter)(workout::Workout, phase::Symbol)
    meter.epochOnly && phase == :train && return
    now = time()
    if now > meter.next || phase == :valid
        result = ""
        for metricname in meter.metricnames
            getmetricvalue(workout, metricname) do value
                    result *= @sprintf  " - %s: %1.4f" metricname value
            end
        end
        if result != ""
            s = @sprintf "[%4d:%7d]" workout.epochs workout.steps
            print("\r", s, result)
            meter.next = now + meter.throttle
        end
        phase == :valid && println()
    end
end

"""
Logs metrics to a TensorBoard file so it can be viewed with TensorBoard. By default
it will log the metrics loss and val_loss at the end of each training step
and end of the validation phase.

This meter depends on the TensorBoardLogger to be installed. So if you didn't do so
already, please run:

    import Pkg; Pkg.add("TensorBoardLogger")
"""
mutable struct TensorBoardMeter <: Meter
    logger
    path::String
    metrics
    last_processed::IdDict{Symbol,Int}

    function TensorBoardMeter(path="./tensorboard_logs/runs", metrics=[:loss, :val_loss])
        try
            @eval import TensorBoardLogger
        catch
            @warn "Package TensorBoardLogger not installed"
        end
        new(nothing, path, metrics, IdDict())
    end
end

function (meter::TensorBoardMeter)(workout::Workout, phase::Symbol)
    if meter.logger === nothing
        meter.logger = TensorBoardLogger.TBLogger(meter.path)
    end

    for metric in meter.metrics
        getmetricvalue(workout, metric) do value
            last = get(meter.last_processed, metric, -1)
            workout.steps > last && TensorBoardLogger.log_value(meter.logger,
                string(metric), value, step=workout.steps)
            meter.last_processed[metric] = workout.steps
        end
    end
end

"""
Logs metrics to a delimeted text file. By default it will log the metrics loss and
val_loss at the end of each training step and end of the validation phase.

"""
mutable struct FileMeter <: Meter
    fileio
    filename::String
    metrics
    last_processed::IdDict{Symbol,Int}
    epochOnly::Bool

    function FileMeter(filename="traininglog.txt", metrics=[:loss, :val_loss], epochOnly=true)
        try
            @eval import DelimitedFiles
        catch
            @warn "Package DelimitedFiles not installed"
        end
        io = open(filename, "w")
        new(io, filename, metrics, IdDict(), epochOnly)
    end
end

function (meter::FileMeter)(workout::Workout, phase::Symbol)
    meter.epochOnly && phase == :train && return
    for metric in meter.metrics
        getmetricvalue(workout, metric) do value
            last = get(meter.last_processed, metric, -1)
            workout.steps > last && DelimitedFiles.writedlm(meter.fileio,
                [(metric, value, workout.steps)])
            meter.last_processed[metric] = workout.steps
        end
    end
    flush(meter.fileio)
end

"""
Plot metrics using the Plots module. At the end of each epoch the plot will
be updated with the values of the metrics. This works especially nice if you
are prototyping some model in an IDE like Juno.
"""
struct PlotMeter <: Meter
    plt
    metrics::Vector{Symbol}
    Plots::Module

    function PlotMeter(Plots::Module, metrics=[:loss, :val_loss])
        plt = Plots.plot(length(metrics),xlabel = "steps", ylabel="values", label=metrics)
        Plots.display(plt)
        new(plt, metrics, Plots)
    end
end

function (m::PlotMeter)(workout::Workout, phase::Symbol)
    phase == :train && return
    for (i, metric) in enumerate(m.metrics)
        if hasmetric(workout, metric)
            m.plt[i] = history(workout, metric)
        end
    end
    m.Plots.display(m.plt)
end


text:
filename() = "day07.txt"
read_data() = split(String(read(filename())), "\n")

# part 1

struct BagRule
    color::String
    children::Vector{Pair{String,Int}}
end

function parse_rule(s)
    m = match(r"^(.*) bags contain (.*)\.", s)
    bag = m.captures[1]
    return if startswith(m.captures[2], "no other")
        [(bag = bag, qty = 0, inner = "")]
    else
        inner = split(m.captures[2], r", *")
        [let m = match(r"^(\d+) (.*)$", replace(bag_with_number, r" bags*" => ""))
            (bag = bag, qty = parse(Int, m.captures[1]), inner = m.captures[2])
        end for bag_with_number in inner]
    end
end

function build_bags_dict(lines)
    dct = Dict()
    for line in lines
        rules = parse_rule(line)
        color = rules[1].bag
        rule = BagRule(color, Pair{String,Int}[])
        for r in rules
            r.qty > 0 && push!(rule.children, r.inner => r.qty)
        end
        dct[color] = rule
    end
    return dct
end

function part1(dct, target)
    function dfs(color)
        !haskey(dct, color) || isempty(dct[color].children) && return 0
        color == target && return 1
        return sum(dfs(c[1]) for c in dct[color].children)
    end
    return count(dfs(k) > 0 for k in keys(dct) if k != target)
end

function part1()
    dct = build_bags_dict(read_data())
    part1(dct, "shiny gold")
end

# part 2

function weight(dct, color)
    helper(color) = 1 + mapreduce(c -> c[2] * helper(c[1]), +, dct[color].children; init = 0)
    return helper(color) - 1  # exclude myself
end

function part2()
    dct = build_bags_dict(read_data())
    return weight(dct, "shiny gold")
end

#------------------------------------------
# LightGraphs

using LightGraphs, SimpleWeightedGraphs

bag_index(rules, color::AbstractString) = findfirst(r -> r[1].bag == color, rules)
bag_color(rules, idx::Int) = rules[idx][1].bag

function build_digraph(lines)
    graph = SimpleWeightedDiGraph(length(lines))
    rules = parse_rule.(lines)  # array of array of named tuples
    for r in rules
        outer_bag_index = bag_index(rules, r[1].bag)
        for directed in r
            if directed.qty > 0
                inner_bag_index = bag_index(rules, directed.inner)
                add_edge!(graph, outer_bag_index, inner_bag_index, directed.qty)
            end
        end
    end
    return (; graph, rules)
end

# e.g. part1_graph("shiny gold")
function part1_graph(target)
    graph, rules = build_digraph(readlines(filename()))
    target_index = bag_index(rules, target)
    return count([has_path(graph, i, target_index) for i in 1:nv(graph) if i !== target_index])
end

# --------------------------------------------------
# Visualization

using GraphPlot
using Cairo, Compose, Colors

function plot_pdf()
    lines = readlines(filename())
    graph, rules = build_digraph(lines)

    nodelabel = [bag_color(rules, i) for i in 1:length(rules)]
    nodecolor = [bag_color(rules, i) == "shiny gold" ? colorant"firebrick1" : colorant"honeydew"
        for i in 1:length(rules)]

    draw(PDF("day07.pdf", 100cm, 100cm),
        gplot(graph, nodelabel=nodelabel, nodefillc=nodecolor, arrowlengthfrac=0.01))

    return gplot(graph, nodelabel=nodelabel)
end


text:
immutable Dirichlet <: ContinuousMultivariateDistribution
    alpha::Vector{Float64}
    alpha0::Float64

    function Dirichlet{T <: Real}(alpha::Vector{T})
        alpha0::Float64 = 0.0
        for i in 1:length(alpha)
            if alpha[i] < 0.0
                error("Dirichlet: elements of alpha must be non-negative")
            else
                alpha0 += alpha[i]
            end
        end
        new(float64(alpha), alpha0)
    end
end

Dirichlet(d::Integer, alpha::Real) = Dirichlet(fill(alpha, d))

Dirichlet(dim::Integer) = Dirichlet(ones(dim))

function entropy(d::Dirichlet)
    k = length(d.alpha)
    en = lmnB(d)
    en += (d.alpha0 - k) * digamma(d.alpha0)
    for j in 1:k
        en -= (d.alpha[j] - 1.0) * digamma(d.alpha[j])
    end
    return en
end

function insupport{T <: Real}(d::Dirichlet, x::Vector{T})
    n = length(x)
    if length(d.alpha) != n
        return false
    end
    s = 0.0
    for i in 1:n
        if x[i] < 0.0
            return false
        end
        s += x[i]
    end
    if abs(s - 1.0) > 1e-8
        return false
    end
    return true
end

mean(d::Dirichlet) = d.alpha ./ d.alpha0

function modes(d::Dirichlet)
    k = length(d.alpha)
    x = Array(Float64, k)
    s = d.alpha0 - k
    for i in 1:k
        if d.alpha[i] <= 1.0
            error("modes only defined when alpha[i] > 1 for all i")
        end
        x[i] = (d.alpha[i] - 1.0) / s
    end
    return [x]
end

pdf{T <: Real}(d::Dirichlet, x::Vector{T}) = exp(logpdf(d, x))

function logpdf{T <: Real}(d::Dirichlet, x::Vector{T})
    return dot((d.alpha - 1.0), log(x)) - lmnB(d)
end

function logpdf!{T <: Real}(r::AbstractArray, d::Dirichlet, x::Matrix{T})
    if size(x, 1) != length(d.alpha)
        throw(ArgumentError("Inconsistent argument dimensions."))
    end

    n = size(x, 2)
    if length(r) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end
    b::Float64 = lmnB(d)
    At_mul_B(r, log(x), d.alpha - 1.0)
    for i in 1:n
        r[i] -= b
    end
    return
end

function rand!(d::Dirichlet, x::Vector)
    s = 0.0
    n = length(x)
    for i in 1:n
        tmp = rand(Gamma(d.alpha[i]))
        x[i] = tmp
        s += tmp
    end
    for i in 1:n
        x[i] /= s
    end
    return x
end

function rand(d::Dirichlet)
    x = Array(Float64, length(d.alpha))
    return rand!(d, x)
end

function rand!(d::Dirichlet, X::Matrix)
    m, n = size(X)
    for j in 1:n
        for i in 1:m
            X[i, j] = rand(Gamma(d.alpha[i]))
        end
    end
    for j in 1:n
        isum = 0.0
        for i in 1:m
            isum += X[i, j]
        end
        for i in 1:m
            X[i, j] /= isum
        end
    end
    return X
end

function var(d::Dirichlet)
    n = length(d.alpha)
    tmp = d.alpha0^2 * (d.alpha0 + 1.0)
    S = Array(Float64, n, n)
    for j in 1:n
        for i in 1:n
            if i == j
                S[i, j] = d.alpha[i] * (d.alpha0 - d.alpha[i]) / tmp
            else
                S[i, j] = -d.alpha[i] * d.alpha[j] / tmp
            end
        end
    end
    return S
end

# Log multinomial beta
function lmnB(d::Dirichlet)
    s = 0.0
    for i in 1:length(d.alpha)
        s += lgamma(d.alpha[i])
    end
    return s - lgamma(d.alpha0)
end

mnB(d::Dirichlet) = exp(lmnB(d))


text:
#Demo 4: Multi-dimensional datasets (advanced usage)
#-------------------------------------------------------------------------------

using CMDimData
using CMDimData.MDDatasets
using CMDimData.EasyPlot
import Printf: @sprintf

#==Attributes
===============================================================================#
LBL_AXIS_TIME = "Time (s)"
LBL_AXIS_AMPLITUDE = "Amplitude"

#==Input data
===============================================================================#
tfund = 1e-9 #Fundamental
osr = 20 #(fundamental) oversampling ratio
nfund = 20 #cycles of the fundamental

#==Computations
===============================================================================#
t = DataF1(0:(tfund/osr):(nfund*tfund))
tmax = maximum(t)

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("period", [1,2,3]*tfund)
	PSweep("slope", [-1,0,1]*(1/tmax))
	PSweep("offset", [-.4, 0, .9])
]

#Generate data:
lines = DataHR{DataF1}(sweeplist) #Create empty dataset
tones = DataHR{DataF1}(sweeplist) #Create empty dataset
for inds in subscripts(lines)
	(T, m, b) = coordinates(lines, inds)
	lines.elem[inds...] = DataF1(t.x, (x)->m*x+b)
	tones.elem[inds...] = sin(t*(2pi/T))
end

tones += lines #Create shifted dataset
#Filter on 2T harmonic:
tones_2T = getsubarray(tones, period=2tfund)
#Filter on 3T harmonic:
tones_3T = getsubarray(tones, period=3tfund)
#Filter on increasing slope (3rd index):
tones_incm = getsubarray(tones, :,3,:)

#==Generate plot
===============================================================================#
strns(T) = @sprintf("%.1f ns", T/1e-9) #Generate string with # of ns from supplied value

plot = cons(:plot, nstrips=4, legend=false, #Legend too busy with GracePlot
	ystrip1 = set(axislabel=LBL_AXIS_AMPLITUDE, striplabel="Tones"),
	ystrip2 = set(axislabel=LBL_AXIS_AMPLITUDE, striplabel="Tones ($(strns(2tfund)))"),
	ystrip3 = set(axislabel=LBL_AXIS_AMPLITUDE, striplabel="Tones ($(strns(3tfund)))"),
	ystrip4 = set(axislabel=LBL_AXIS_AMPLITUDE, striplabel="Tones (increasing slope)"),
	xaxis = set(label=LBL_AXIS_TIME),
)
push!(plot,
	cons(:wfrm, tones, strip=1),
	cons(:wfrm, tones_2T, strip=2),
	cons(:wfrm, tones_3T, strip=3),
	cons(:wfrm, tones_incm, strip=4),
)

pcoll = cons(:plot_collection, title="Mulit-Dataset Tests: Subarrays", ncolumns=1)
	push!(pcoll, plot)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
# ---
# title: 710. Random Pick with Blacklist
# id: problem710
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Binary Search, Sort, Random
# link: <https://leetcode.com/problems/random-pick-with-blacklist/description/>
# hidden: true
# ---
#
# Given a blacklist `B` containing unique integers from `[0, N)`, write a
# function to return a uniform random integer from `[0, N)` which is **NOT**  in
# `B`.
#
# Optimize it such that it minimizes the call to system's `Math.random()`.
#
# **Note:**
#
#   1. `1 <= N <= 1000000000`
#   2. `0 <= B.length < min(100000, N)`
#   3. `[0, N)` does NOT include N. See [interval notation](https://en.wikipedia.org/wiki/Interval_\(mathematics\)).
#
# **Example 1:**
#
#
#
#     Input: ["Solution","pick","pick","pick"]
#     [[1,[]],[],[],[]]
#     Output: [null,0,0,0]
#
#
# **Example 2:**
#
#
#
#     Input: ["Solution","pick","pick","pick"]
#     [[2,[]],[],[],[]]
#     Output: [null,1,1,1]
#
#
# **Example 3:**
#
#
#
#     Input: ["Solution","pick","pick","pick"]
#     [[3,[1]],[],[],[]]
#     Output: [null,0,0,2]
#
#
# **Example 4:**
#
#
#
#     Input: ["Solution","pick","pick","pick"]
#     [[4,[2]],[],[],[]]
#     Output: [null,1,3,1]
#
#
# **Explanation of Input Syntax:**
#
# The input is two lists: the subroutines called and their arguments.
# `Solution`'s constructor has two arguments, `N` and the blacklist `B`. `pick`
# has no arguments. Arguments are always wrapped with a list, even if there
# aren't any.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Pythia

export sglasso, glasso, lasso

include("sglasso.jl")

end


text:
"""
3. Veloctiy and stress only (ρ_0 = 1.0)
f₀(x) = f_eq(ρ₀(x), u₀(x), p₀(x))
ρ₀ = 1.0
u₀ = u₀(x)
p₀ = unkown
"""
struct ConstantDensity <: InitializationStrategy end

function initial_condition(
    ::ConstantDensity,
    q::Quadrature,
    problem::FluidFlowProblem,
    x::T,
    y::T,
) where {T <: Real}
    u = lattice_velocity(q, problem, x, y)

    return equilibrium(q, one(T), u, one(T))
end


text:
"""
Package: DataWrangler

outlie\\[!\\]([x], y, σ = 2; q)

Replace or detect outlier values in vector `x` describing a time series

# Parameters
- `x`: Optional vector containing the support for `y` (no missing value allowed)
- `y`: Vector with outliers to be dealt with
- `σ`: Number of sigmas used to identify outliers, it defaults to four.
- `q`: Number of closest vector values to the imputation to be considered, it defaults to 3*length(y)÷4).

# Description

When using `outlie!` the function replaces outliers with a `missing` value and return the index of the values replaced, if `outlie` is used isntead only the outliers index poistion are returned.

Outliers are indentified by measureing its distance in number of sigmas estimated from a loess fitting with window `q`. When datasets with noise close to zero all outliers will still be detected but non-outliers might be indentified as such, a way to prevent non-outliers to be indentified in this situtaion is by introducing a small amount of noise ϵ in `y`, e.g. `outlie!(x,y+ϵ)`.

# Examples
```julia
n = 1000
x = sort(rand(n))*2*pi;
y = Array{Union{Missing,Float64}}(undef,n);
y[:] = sin.(x).+randn(n)/10
mid = vcat(100,300,600,950);
y[mid] .= x[mid] .+ 2*(randn(length(mid)).+1)
y[500] = x[500] - 2*(randn(1)[1]+1)

outlie(x,y)

julia> outlie(x,y)
5-element Vector{Int64}:
 100
 300
 500
 600
 950
```
"""
function outlie(y::AbstractVector{<:Union{Missing,T}}, nσ::Real=4.0;
                q = 3*length(y)÷4) where  {R<:Real,T<:Real}
    cy = copy(y)
    outlie!(1:length(y),cy,nσ;q)
end

function outlie(x::AbstractVector{R}, y::AbstractVector{<:Union{Missing,T}}, nσ::Real=4.0;
                 q = 3*length(y)÷4) where  {R<:Real,T<:Real}
    cy = copy(y)
    outlie!(x,cy,nσ;q)
end

function outlie!(y::AbstractVector{<:Union{Missing,T}}, nσ::Real=4.0;
                q = 3*length(y)÷4) where  {R<:Real,T<:Real}
    outlie!(1:length(y),y,nσ;q)
end

function outlie!(x::AbstractVector{R}, y::AbstractVector{<:Union{Missing,T}}, nσ::Real=4.0;
                 q = 3*length(y)÷4) where  {R<:Real,T<:Real}
    allout = Int64[]
    eo = Int64[0]
    while !isempty(eo)
        eo = outlie1!(x,y,nσ;q)
        push!(allout,eo...)
    end
    sort(allout)
end

function outlie1!(x::AbstractVector{R}, y::AbstractVector{<:Union{Missing,T}}, nσ::Real=4.0;
                  q = 3*length(y)÷4) where  {R<:Real,T<:Real}

    @assert nσ > 1 "`nσ` should be larger than one"
    @assert !any(ismissing, x) "No missing values allowed in `x`"
    @assert issorted(x) "`x` values must be sorted"
    @assert length(x) == length(y) "`x` and `y` must have the same length"

    loq = loess(x,y; q)
    loqy = loq(x)

    dy = y-loqy
    dy = dy .- median(skipmissing(dy))
    sdy = collect(skipmissing(sort(dy)))
    lsdy = length(sdy)
    sdy = sdy[lsdy÷10:end-lsdy÷10]  # robust sd
    dσ = std(skipmissing(sdy))*T(sqrt(2))

    outlist = findall(v-> abs(v) > nσ*dσ, skipmissing(dy))

    y[outlist] .= missing
    outlist

end


text:
struct ODEInterpolatingAdjointSensitivityFunction{C<:AdjointDiffCache,Alg<:InterpolatingAdjoint,
                                                  uType,SType,CPS,fType<:DiffEqBase.AbstractDiffEqFunction,CV} <: SensitivityFunction
  diffcache::C
  sensealg::Alg
  discrete::Bool
  y::uType
  sol::SType
  checkpoint_sol::CPS
  f::fType
  colorvec::CV
end

mutable struct CheckpointSolution{S,I,T}
  cpsol::S # solution in a checkpoint interval
  intervals::I # checkpoint intervals
  cursor::Int # sol.prob.tspan = intervals[cursor]
  tols::T
end

function ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,checkpoints,colorvec,tols)
  tspan = reverse(sol.prob.tspan)
  checkpointing = ischeckpointing(sensealg, sol)
  (checkpointing && checkpoints === nothing) && error("checkpoints must be passed when checkpointing is enabled.")

  checkpoint_sol = if checkpointing
    intervals = map(tuple, @view(checkpoints[1:end-1]), @view(checkpoints[2:end]))
    interval_end = intervals[end][end]
    tspan[1] > interval_end && push!(intervals, (interval_end, tspan[1]))
    cursor = lastindex(intervals)
    interval = intervals[cursor]
    cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1])), sol.alg; tols...)
    CheckpointSolution(cpsol, intervals, cursor, tols)
  else
    nothing
  end

  diffcache, y = adjointdiffcache(g,sensealg,discrete,sol,dg,sol.prob.f)

  return ODEInterpolatingAdjointSensitivityFunction(diffcache,sensealg,
                                                    discrete,y,sol,
                                                    checkpoint_sol,sol.prob.f,colorvec)
end

function findcursor(intervals, t)
  # equivalent with `findfirst(x->x[1] <= t <= x[2], intervals)`
  lt(x, t) = <(x[2], t)
  return searchsortedfirst(intervals, t, lt=lt)
end

# u = λ'
# add tstop on all the checkpoints
function (S::ODEInterpolatingAdjointSensitivityFunction)(du,u,p,t)
  @unpack sol, y, checkpoint_sol, discrete = S
  idx = length(y)
  f = sol.prob.f

  if checkpoint_sol === nothing
    if typeof(t) <: ForwardDiff.Dual && eltype(S.y) <: AbstractFloat
      y = sol(t)
    else
      sol(y,t)
    end
  else
    intervals = checkpoint_sol.intervals
    interval = intervals[checkpoint_sol.cursor]
    if !(interval[1] <= t <= interval[2])
      cursor′ = findcursor(intervals, t)
      interval = intervals[cursor′]
      cpsol_t = checkpoint_sol.cpsol.t
      if typeof(t) <: ForwardDiff.Dual && eltype(S.y) <: AbstractFloat
        y = sol(interval[1])
      else
        sol(y, interval[1])
      end
      prob′ = remake(sol.prob, tspan=intervals[cursor′], u0=y)
      cpsol′ = solve(prob′, sol.alg; dt=abs(cpsol_t[end] - cpsol_t[end-1]), checkpoint_sol.tols...)
      checkpoint_sol.cpsol = cpsol′
      checkpoint_sol.cursor = cursor′
    end
    checkpoint_sol.cpsol(y, t)
  end

  λ     = @view u[1:idx]
  grad  = @view u[idx+1:end]
  dλ    = @view du[1:idx]
  dgrad = @view du[idx+1:end]

  vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad)

  dλ .*= -one(eltype(λ))

  discrete || accumulate_dgdu!(dλ, y, p, t, S)
  return nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,sensealg::InterpolatingAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     reltol=nothing, abstol=nothing,
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  tspan = reverse(tspan)
  discrete = t != nothing

  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")
  numstates = length(u0)
  numparams = length(p)

  len = numstates+numparams

  λ = similar(p, len)
  λ .= false
  sense = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,
                                                     checkpoints,f.colorvec,
                                                     (reltol=reltol,abstol=abstol))

  init_cb = t !== nothing && tspan[1] == t[end]
  cb = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  original_mm = sol.prob.f.mass_matrix
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    mm = zeros(len, len)
    copyto!(@view(mm[1:numstates, 1:numstates]), sol.prob.f.mass_matrix')
    copyto!(@view(mm[numstates+1:end, numstates+1:end]), I)
  end
  odefun = ODEFunction(sense, mass_matrix=mm)
  return ODEProblem(odefun,z0,tspan,p,callback=cb)
end


text:
"""
    AsymmetricLaplace(μ,λ,κ)

The *Asymmetric Laplace distribution* with location `μ` scale `λ` and asymmetry `κ` has probability density function

```math
f(x; \\mu, \\lambda, \\kappa) = \\left( \\frac{\\lambda}{\\kappa + 1/\\kappa} \\right) \\exp - (x - \\mu) \\lambda \\sgn(x - \\mu) \\kappa^\\sgn(x - \\mu)
```

```julia
AsymmetricLaplace()          # Asymmetric Laplace distribution with zero location, unit scale and asymmetry 1, i.e. AsymmetricLaplace(0, 1, 1)
AsymmetricLaplace(μ)         # Asymmetric Laplace distribution with location μ, unit scale and asymmetry 1, i.e. AsymmetricLaplace(μ, 1, 1)
AsymmetricLaplace(μ, λ)      # Asymmetric Laplace distribution with location μ, scale λ and asymmetry 1
AsymmetricLaplace(μ, λ, κ)   # Asymmetric Laplace distribution with location μ, scale λ and asymmetry κ

params(d)                    # Get the parameters, i.e., (μ, λ, κ)
location(d)                  # Get the location parameter, i.e. μ
scale(d)                     # Get the scale parameter, i.e. λ
asymmetry(d)                 # Get the asymmetry parameter, i.e. κ
```

External links

* [Asymmetric Laplace distribution on Wikipedia](http://en.wikipedia.org/wiki/Asymmetric_Laplace_distribution)

"""
struct AsymmetricLaplace{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    λ::T
    κ::T
    AsymmetricLaplace{T}(µ::T, λ::T, κ::T) where {T} = new{T}(µ, λ, κ)
end

function AsymmetricLaplace(μ::T, λ::T, κ::T; check_args=true) where {T <: Real}
    check_args && @check_args(AsymmetricLaplace, λ > zero(λ) && κ > zero(κ))
    return AsymmetricLaplace{T}(μ, λ, κ)
end

AsymmetricLaplace(μ::Real, λ::Real, κ::Real) = AsymmetricLaplace(promote(μ, λ, κ)...)
AsymmetricLaplace(μ::Integer, λ::Integer, κ::Integer) = AsymmetricLaplace(float(μ), float(λ), float(κ))
AsymmetricLaplace(μ::T) where {T <: Real} = AsymmetricLaplace(μ, one(T), one(T))
AsymmetricLaplace(μ::T, λ::T) where {T <: Real} = AsymmetricLaplace(μ, λ, one(T))
AsymmetricLaplace() = AsymmetricLaplace(0.0, 1.0, 1.0, check_args=false)

@distr_support AsymmetricLaplace -Inf Inf

#### Conversions
function convert(::Type{AsymmetricLaplace{T}}, μ::S, λ::S, κ::S) where {T <: Real, S <: Real}
    AsymmetricLaplace(T(μ), T(λ), T(κ))
end
function convert(::Type{AsymmetricLaplace{T}}, d::AsymmetricLaplace{S}) where {T <: Real, S <: Real}
    AsymmetricLaplace(T(d.μ), T(d.λ), T(d.κ), check_args=false)
end

#### Parameters

location(d::AsymmetricLaplace) = d.μ
scale(d::AsymmetricLaplace) = d.λ
asymmetry(d::AsymmetricLaplace) = d.κ
params(d::AsymmetricLaplace) = (d.μ, d.λ, d.κ)
@inline partype(d::AsymmetricLaplace{T}) where {T<:Real} = T

#### Statistics

mean(d::AsymmetricLaplace) = d.μ + (1 - d.κ^2) / (d.λ * d.κ)
function median(d::AsymmetricLaplace)
    μ, λ, κ = params(d)
    if κ > 1
        μ + κ/λ * log((1 + κ^2)/(2κ^2))
    else
        μ - 1/(κ * λ) * log((1 + κ^2)/2)
    end
end
mode(d::AsymmetricLaplace) = d.μ

var(d::AsymmetricLaplace) = (1 + d.κ^4) / (d.λ^2 * d.κ^2)
std(d::AsymmetricLaplace) = sqrt(1 + d.κ^4) / (d.λ * d.κ)
skewness(d::AsymmetricLaplace) = 2(1 - d.κ^6) / (d.κ^4 + 1)^(3/2)
kurtosis(d::AsymmetricLaplace) = 6(1 + d.κ^8) / (1 + d.κ^4)^2

entropy(d::AsymmetricLaplace) = log(ℯ * (1 + d.κ) / (d.λ * d.κ))

#### Evaluations

function pdf(d::AsymmetricLaplace, x::Real)
    μ, λ, κ = params(d)
    if x ≤ μ
        exp((x - μ) / (λ * κ)) / λ / (κ + 1 / κ)
    else
        exp(- κ / λ * (x - μ)) / λ / (κ + 1 / κ)
    end
end
function logpdf(d::AsymmetricLaplace, x::Real)
    μ, λ, κ = params(d)
    if x ≤ μ
        (x - μ) / (λ * κ) - log(λ) - log(κ + 1 / κ)
    else
        - κ * (x - μ) / λ - log(λ) - log(κ + 1 / κ)
    end
end

function cdf(d::AsymmetricLaplace, x::Real)
    μ, λ, κ = params(d)
    if x ≤ μ
        κ^2 / (1 + κ^2) * exp((x - μ) / (λ * κ))
    else
        1 - 1 / (1 + κ^2) * exp(- κ / λ * (x - μ))
    end
end
ccdf(d::AsymmetricLaplace, x::Real) = 1 - cdf(d, x)
function logcdf(d::AsymmetricLaplace, x::Real)
    μ, λ, κ = params(d)
    if x ≤ μ
        2 * log(κ) - log(1 + κ^2) + (x - μ) / (λ * κ)
    else
        log(1 - 1 / (1 + κ^2) * exp(- κ / λ * (x - μ)))
    end
end
function logccdf(d::AsymmetricLaplace, x::Real)
    μ, λ, κ = params(d)
    if x ≤ μ
        log(1 - κ^2 / (1 + κ^2) * exp((x - μ) / (λ * κ)))
    else
        -(κ * (x - μ) / λ  + log(1 + κ^2))
    end
end

function quantile(d::AsymmetricLaplace, p::Real)
    μ, λ, κ = params(d)
    if p ≤ cdf(d, d.μ)
        μ + log((1 + κ^2) / κ^2 * p) * (λ * κ)
    else
        μ - λ / κ * log((1 - p) * (1 + κ^2))
    end
end
cquantile(d::AsymmetricLaplace, p::Real) = 1 - quantile(d, p)
function invlogcdf(d::AsymmetricLaplace, lp::Real)
    μ, λ, κ = params(d)
    if lp ≤ 0
        println("branch less than")
        μ - (λ / κ) * (log(1 - exp(lp)) + log(1 + κ^2))
    else
        println("branch larger than")
        μ + λ * κ * (log(lp) + log(1 + κ^2) - 2 * log(κ))
    end
end
function invlogccdf(d::AsymmetricLaplace, lp::Real)
    μ, λ, κ = params(d)
    if lp ≤ log(cdf(d, μ))
        μ - (λ / κ) * (lp + log(1 + κ^2))
    else
        μ + λ * κ * log((1 + κ^2) * (1 - exp(lp)) / κ^2)
    end
end

function gradlogpdf(d::AsymmetricLaplace, x::Real)
    x == d.μ && error("Gradient is undefined at the location point")
    x > d.μ ? - d.κ / d.λ : 1 / (d.λ * d.κ)
end

function mgf(d::AsymmetricLaplace, t::Real)
    exp(t * d.μ) / ((1 + t * d.κ / d.λ) * (1 - t / (d.κ * d.λ)))
end
function cf(d::AsymmetricLaplace, t::Real)
    cis(t * d.μ) / ((1 + t * d.κ * im / d.λ) * (1 - t * im / (d.κ * d.λ)))
end

#### Sampling

rand(rng::AbstractRNG, d::AsymmetricLaplace) = quantile(d, rand(rng))


text:
using Distributed

addprocs(
    readlines(ENV["PBS_NODEFILE"]);
    enable_threaded_blas = true,
    topology = :master_worker,
    exeflags = [
        "--project=$(Base.active_project())",
        "--threads=$(Threads.nthreads())",
        "--optimize",
    ],
)

# Load packages
@everywhere using TestProj
@everywhere MASTER_ARGS = ARGS

pmap(1:10) do id
    println(hello(join(MASTER_ARGS, ", ", ", and ")))
    println("PID:          $(getpid())")
    println("Host name:    $(gethostname())")
    println("Num. threads: $(Threads.nthreads())")
    sleep(0.5 + rand())
end

nothing


text:
export HyperPlanesIntersection, LinesHull, affinehull, linespace

# It is easy to go from H-rep of affine space to V-rep of affine space by computing the kernel of a matrix using RowEchelon
# However, it is really worth it since Base.in for an HRepElement in HyperPlanesIntersection and Base.in for an VRepElement in LinesHull are false already.
function remproj(x::RepElement{<:Integer}, l::RepElement{<:Integer})
    # generates large numbers but keeps the type integer
    x * dot(coord(l), coord(l)) - l * dot(coord(x), coord(l))
end
function remproj(x, l)
    ll = dot(coord(l), coord(l))
    if isapproxzero(ll)
        x
    else
        simplify(x - l * ((dot(coord(x), coord(l)) / ll)))
    end
end

# H-representation

"""
    hrep(hyperplanes::HyperPlaneIt; d::FullDim)

Creates an affine space of full dimension `d` from the list of hyperplanes
`hyperplanes`.

### Examples
```julia
hrep([HyperPlane([0, 1, 0], 1), HyperPlane([0, 0, 1], 0)])
```
creates the 1-dimensional affine subspace containing all the points ``(x_1, 0, 0)``, i.e. the ``x_1``-axis.

```julia
hrep([HyperPlane([1, 1], 1), HyperPlane([1, 0], 0)])
```
creates the 0-dimensional affine subspace only containing the point ``(0, 1)``.
"""
function hrep(hyperplanes::HyperPlaneIt; d::FullDim = FullDim_rec(hyperplanes))
    return HyperPlanesIntersection(d, hyperplanes)
end

# Representation of an affine space as the intersection of hyperplanes.
# Also called affine set, affine manifold, affine variety, linear variety or flat.
# An affine space L satisfies:
# λx + (1-λ)y ∈ L, ∀x, y ∈ L, ∀ λ ∈ R
# Note that λ is not rhyperplaneuired to be between 0 and 1 as in convex sets.
mutable struct HyperPlanesIntersection{T, AT, D<:FullDim} <: HAffineSpace{T}
    d::D
    # HyperPlanes whose intersection is the affine space
    hyperplanes::Vector{HyperPlane{T, AT}}
    orthogonalized::Bool
    function HyperPlanesIntersection{T, AT, D}(d::FullDim,
                                               hps::HyperPlaneIt{T}) where {T, AT, D}
        new{T, AT, D}(FullDim_convert(D, d), lazy_collect(hps), false)
    end
end
function HyperPlanesIntersection(d::FullDim,
                                 it::ElemIt{HyperPlane{T, AT}}) where {T, AT}
    return HyperPlanesIntersection{T, AT, typeof(d)}(d, it)
end
HyperPlanesIntersection{T, AT}(d::FullDim) where {T, AT} = HyperPlanesIntersection(d, HyperPlane{T, AT}[])

FullDim(h::HyperPlanesIntersection) = h.d
hvectortype(L::Type{<:HyperPlanesIntersection{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:HyperPlanesIntersection}, d::FullDim, ::Type{T}) where T = HyperPlanesIntersection{T, similar_type(hvectortype(PT), d, T), typeof(d)}

function Base.empty!(L::HyperPlanesIntersection{T, AT}) where {T, AT}
    # We want the dimension to be `-1` and the hyperplanes to make
    # an infeasible system so we build the system `x = 0 & 0 = 1`.
    empty!(L.hyperplanes)
    for i in 1:fulldim(L)
        push!(L.hyperplanes, HyperPlane(basis(AT, L.d, i), zero(T)))
    end
    push!(L.hyperplanes, HyperPlane(origin(AT, fulldim(L)), one(T)))
    return L
end
function Base.intersect!(L::HyperPlanesIntersection, h::HyperPlane)
    if L.orthogonalized
        h_p = remproj(h, L)
        if isapproxzero(h_p.a)
            if !isapproxzero(h_p.β)
                # `L` is empty
                empty!(L)
            end
        else
            push!(L.hyperplanes, h_p)
        end
    else
        push!(L.hyperplanes, h)
    end
    return L
end

@vecrepelem HyperPlanesIntersection HyperPlane hyperplanes

hreptype(::Type{<:HyperPlanesIntersection{T, AT}}) where {T, AT} = Intersection{T, AT}

# Returns an HyperPlanesIntersection representing the affine hull of p.
# The affine hull is defined as
# {λx + (1-λ)y | x, y ∈ p, λ ∈ R}
function affinehull(h::HRep, current=false)
    if !current
        detecthlinearity!(h)
    end
    HyperPlanesIntersection(FullDim(h), hyperplanes(h))
end

function detecthlinearity!(L::HyperPlanesIntersection, solver = nothing)
    if !L.orthogonalized
        V = detecthlinearity(L)
        L.hyperplanes = V.hyperplanes
        L.orthogonalized = true
    end
    return L
end
function detecthlinearity(L::HyperPlanesIntersection{T, AT}, solver=nothing) where {T, AT}
    H = HyperPlanesIntersection{T, AT}(FullDim(L))
    for h in hyperplanes(L)
        hp = remproj(h, H)
        if !(h in H)
            intersect!(H, h)
        end
    end
    H
end

# V-representation
struct VEmptySpace{T, AT, D <: FullDim} <: VLinearSpace{T}
    d::D
end
FullDim(v::VEmptySpace) = v.d
function emptyspace(v::VRep{T}) where {T}
    d = FullDim(v)
    return VEmptySpace{T, vvectortype(typeof(v)), typeof(d)}(d)
end

Base.in(v::VRepElement, L::VEmptySpace) = isapproxzero(v)

"""
    vrep(lines::LineIt; d::FullDim)

Creates an affine space of full dimension `d` from the list of lines `lines`.

### Examples
```julia
vrep([Line([1, 0, 0]), Line([0, 1, 0])])
```
creates the 2-dimensional affine subspace containing all the points ``(x_1, x_2, 0)``, i.e. the ``x_1````x_2``-plane.
"""
vrep(lines::LineIt; d::FullDim = FullDim_rec(lines)) = LinesHull(d, lines)

# Representation of an affine space containing the origin by the minkowsky sum of lines
mutable struct LinesHull{T, AT, D<:FullDim} <: VLinearSpace{T}
    d::D
    lines::Vector{Line{T, AT}}
    orthogonalized::Bool
    function LinesHull{T, AT, D}(d::FullDim, lines::LineIt{T}) where {T, AT, D}
        new{T, AT, D}(FullDim_convert(D, d), lazy_collect(lines), false)
    end
end
function LinesHull(d::FullDim, it::ElemIt{Line{T, AT}}) where {T, AT}
    return LinesHull{T, AT, typeof(d)}(d, it)
end
LinesHull{T, AT}(d::FullDim) where {T, AT} = LinesHull(d, Line{T, AT}[])

FullDim(v::LinesHull) = v.d
vvectortype(::Type{<:LinesHull{T, AT}}) where {T, AT} = AT
similar_type(PT::Type{<:LinesHull}, d::FullDim, ::Type{T}) where T = LinesHull{T, similar_type(vvectortype(PT), d, T), typeof(d)}

function convexhull!(L::LinesHull, l::Line)
    if L.orthogonalized
        l_p = remproj(l, L)
        if !isapproxzero(l_p)
            push!(L.lines, l_p)
        end
    else
        push!(L.lines, l)
    end
    return L
end

@vecrepelem LinesHull Line lines

conetype(::Type{LinesHull{T, AT, D}}) where {T, AT, D} = RaysHull{T, AT, D}
vreptype(::Type{LinesHull{T, AT, D}}) where {T, AT, D} = Hull{T, AT, D}

# Returns a LinesHull representing the following set (TODO does it have a name?, does someone has a reference talking about it ?)
# {x | y in p => x + y in p}
function linespace(v::VRep, current=false)
    if !current
        detectvlinearity!(v)
    end
    LinesHull(FullDim(v), lines(v))
end

function detectvlinearity!(L::LinesHull, solver = nothing)
    if !L.orthogonalized
        V = detectvlinearity(L)
        L.lines = V.lines
        L.orthogonalized = true
    end
    return L
end
function detectvlinearity(L::LinesHull{T, AT}, solver = nothing) where {T, AT}
    V = LinesHull{T, AT}(FullDim(L))
    V.orthogonalized = true
    for l in lines(L)
        convexhull!(V, l)
    end
    return V
end

detectlinearity!(h::HRepresentation) = detecthlinearity!(h)
detectlinearity!(v::VRepresentation) = detectvlinearity!(v)
_linearity(h::HRepresentation) = hyperplanes(h)
_linearity(v::VRepresentation) = lines(v)
function remproj(el::RepElement, L::Union{HyperPlanesIntersection, LinesHull})
    detectlinearity!(L)
    # Modified Gram-Schmidt process, it is more stable than
    # the classical Gram-Schmidt process.
    # see https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process#Numerical_stability
    for lin in _linearity(L)
        el = remproj(el, lin)
    end
    return el
end
function Base.in(el::HRepElement, L::HyperPlanesIntersection)
    return isapproxzero(remproj(el, L))
end
function Base.in(el::VRepElement, L::LinesHull)
    return isapproxzero(remproj(el, L))
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
using quantumoptics.sparsematrix
using Base.LinAlg.BLAS

N = 600

x1 = complex(sprand(N, N, 0.1))
x2 = complex(sprand(N, N, 0.1))
y1 = quantumoptics.sparsematrix.SparseMatrix(full(x1))
y2 = quantumoptics.sparsematrix.SparseMatrix(full(x2))

x1*x2
v1 = ones(Complex128, N)
v2 = zeros(Complex128, N)

function f_mv(alpha::Complex128, M::SparseMatrixCSC{Complex128}, v::Vector{Complex128}, beta::Complex128, result::Vector{Complex128})
    result[:] = beta*result + alpha*(M*v)
end

function f_mv(alpha::Complex128, v::Vector{Complex128}, M::SparseMatrixCSC{Complex128}, beta::Complex128, result::Vector{Complex128})
    result[:] = beta*result + alpha*transpose(transpose(v)*M)
end

function f_mm(alpha::Complex128, M::SparseMatrixCSC{Complex128}, B::Matrix{Complex128}, beta::Complex128, result::Matrix{Complex128})
    result[:] = beta*result + (alpha*M)*B
end

function f_mm(alpha::Complex128, B::Matrix{Complex128}, M::SparseMatrixCSC{Complex128}, beta::Complex128, result::Matrix{Complex128})
    result[:] = beta*result + B*(alpha*M)
end

function mv(alpha::Complex128, M::SparseMatrixCSC{Complex128}, v::Vector{Complex128}, beta::Complex128, result::Vector{Complex128})
    @inbounds for i=1:length(result)
        result[i] *= beta
    end
    @inbounds for colindex = 1:M.n
        vj = alpha*v[colindex]
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            result[M.rowval[i]] += M.nzval[i]*vj
        end
    end
    nothing
end

function mv(alpha::Complex128, v::Vector{Complex128}, M::SparseMatrixCSC{Complex128}, beta::Complex128, result::Vector{Complex128})
    @inbounds for i=1:length(result)
        result[i] *= beta
    end
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            result[colindex] += M.nzval[i]*alpha*v[M.rowval[i]]
        end
    end
    nothing
end

function mm(alpha::Complex128, M::SparseMatrixCSC{Complex128}, M2::Matrix{Complex128}, beta::Complex128, result::Matrix{Complex128})
    @inbounds for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for j=1:M.m
        for colindex = 1:M.n
            m2 = alpha*M2[colindex, j]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i], j] += M.nzval[i]*m2
            end
        end
    end
    nothing
end

function mm(alpha::Complex128, v::Matrix{Complex128}, M::SparseMatrixCSC{Complex128}, beta::Complex128, result::Matrix{Complex128})
    @inbounds for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            mi = M.nzval[i]*alpha
            mrowvali = M.rowval[i]
            for j=1:size(result,1)
                result[j,colindex] += mi*v[j,mrowvali]
            end
        end
    end
    nothing
end

M = full(x2)
result = full(x2)

println("spMatrix-Vector")
@time f_mv(Complex(1.), x1, v1, Complex(0.), v2)
@time f_mv(Complex(1.), x1, v1, Complex(0.), v2)

@time mv(Complex(1.), x1, v1, Complex(0.), v2)
@time mv(Complex(1.), x1, v1, Complex(0.), v2)

@time quantumoptics.sparsematrix.gemv!(Complex(1.), y1, v1, Complex(0.), v2)
@time quantumoptics.sparsematrix.gemv!(Complex(1.), y1, v1, Complex(0.), v2)

println("Vector-spMatrix")
@time f_mv(Complex(1.), v1, x1, Complex(0.), v2)
@time f_mv(Complex(1.), v1, x1, Complex(0.), v2)

@time mv(Complex(1.), v1, x1, Complex(0.), v2)
@time mv(Complex(1.), v1, x1, Complex(0.), v2)

# @time quantumoptics.sparsematrix.gemv!(Complex(1.), v1, y1, Complex(0.), v2)
# @time quantumoptics.sparsematrix.gemv!(Complex(1.), v1, y1, Complex(0.), v2)

println("spMatrix-Matrix")
@time f_mm(Complex(1.), x1, M, Complex(0.), result)
@time f_mm(Complex(1.), x1, M, Complex(0.), result)

@time mm(Complex(1.), x1, M, Complex(0.), result)
@time mm(Complex(1.), x1, M, Complex(0.), result)

@time quantumoptics.sparsematrix.gemm!(Complex(1.), y1, M, Complex(0.), result)
@time quantumoptics.sparsematrix.gemm!(Complex(1.), y1, M, Complex(0.), result)

println("Matrix-spMatrix")
@time f_mm(Complex(1.), M, x1, Complex(0.), result)
@time f_mm(Complex(1.), M, x1, Complex(0.), result)

@time mm(Complex(1.), M, x1, Complex(0.), result)
@time mm(Complex(1.), M, x1, Complex(0.), result)

@time quantumoptics.sparsematrix.gemm!(Complex(1.), M, y1, Complex(0.), result)
@time quantumoptics.sparsematrix.gemm!(Complex(1.), M, y1, Complex(0.), result)


text:
module RayTraceEllipsoid

using CoordinateTransformations, StaticArrays, UnitfulAngles
import LinearAlgebra: normalize, ⋅

export Vec, Ray, Ellipsoid, Interface, OpticUnit, raytrace!, advance!, bend!

"""
    Vec = SVector{3,Float64}

Points and directions are just a point in 3D space best described as a static vector, `SVector`. `Vec` is an alias for that.
"""
const Vec = SVector{3,Float64}

"""
    Ray(orig::Vec, dir::Vec)

The main Ray type with a ray origin, `orig`, and direction, `dir`. The ray's direction gets normalized.
"""
mutable struct Ray
    orig::Vec
    dir::Vec
    Ray(o::Vec, d::Vec) = new(o, normalize(d))
end

Ray() = Ray(Vec(0,0,0), Vec(1,0,0))

"""
    Ellipsoid(c::Vec, r::Vec, dir::Vec, open::Float64)

    An ellipsoid with a center, `c`, and radii, `r`, as well as a direction (gets automatically normalized), `dir`, and an opening angle, `open`, creating a dome (or window), that the ellipsoid is defined in. Note that `open` is the angle between the dome's edge and the direction of the dome (so actually half the opening angle) and is defined in **some angular units** (using UnitfulAngles, for example: u"°").

`Ellipsoid` has 6 additional fields all relating to various spatial transformations that convert the ellipsoid to a unit-sphere and back again. These are all `CoordinateTransformations`.

# Examples
For an ellipsoid upwards-pointing hemisphere with a center at (1,2,3), and radii (4,5,6):
```jldoctest
julia> Ellipsoid(Vec(1,2,3), Vec(4,5,6), Vec(0,0,1), 0.0)
Rayden.Ellipsoid([1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [0.0, 0.0, 1.0], 0.0, Translation(-1.0, -2.0, -3.0), LinearMap([0.25 0.0 0.0; 0.0 0.2 0.0; 0.0 0.0 0.166667]), AffineMap([0.25 0.0 0.0; 0.0 0.2 0.0; 0.0 0.0 0.166667], [-0.25, -0.4, -0.5]), Translation(1.0, 2.0, 3.0), LinearMap([4.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 6.0]), AffineMap([4.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 6.0], [1.0, 2.0, 3.0]))
```
"""
struct Ellipsoid
    c::Vec
    r::Vec
    dir::Vec
    open::Float64 # cos(α), where α is half the opening angle of the dome
    # all of the following are transformations
    center::Translation{Vec} # translate the ellipsoid to zero
    scale::LinearMap{SDiagonal{3,Float64}}
    center_scale::AffineMap{SDiagonal{3,Float64},Vec} # translate and scale to a unit-sphere
    uncenter::Translation{Vec}
    unscale::LinearMap{SDiagonal{3,Float64}}
    uncenter_unscale::AffineMap{SDiagonal{3,Float64},Vec}
    function Ellipsoid(c::Vec, r::Vec, dir::Vec, α)
        uncenter = Translation(c)
        unscale = LinearMap(SDiagonal(r))
        center = inv(uncenter)
        scale = inv(unscale)
        center_scale = scale∘center
        uncenter_unscale = inv(center_scale)
        new(c, r, normalize(dir), cos(α), center, scale, center_scale, uncenter, unscale, uncenter_unscale)
    end
end

"""
    distance(orig::Vec, dir::Vec)

Return the two distances between a point with origin, `orig`, and direction, `dir`, and the two (potentially identical) intersection points with a unit-sphere.
"""
function distance(orig::Vec, dir::Vec)
    b = -orig⋅dir
    disc = b^2 - orig⋅orig + 1
    if disc ≥ 0
        d = sqrt(disc)
        t2 = b + d
        if t2 ≥ 0
            t1 = b - d
            return t1 > 0 ? (t1, t2) : (Inf, t2)
        end
    end
    return (Inf, Inf)
end

"""
    advance!(r::Ray, s::Ellipsoid)

Find the shortest point of intersection that is within the ellipsoid's dome and reassign the origin of the ray. Returns failure of the intersection.
"""
function advance!(r::Ray, s::Ellipsoid)
    # move the ray's origin to the intersection point that is within the ellipsoid's dome, and return failure
    orig = s.center_scale(r.orig) # transform the ray's origin according to the ellipsoid
    dir = normalize(s.scale(r.dir)) # scale the ray's direction as well
    ls = distance(orig, dir)
    for l in ls
        if !isinf(l)
            o = orig + l*dir
            p = s.unscale(o)
            cosα = s.dir⋅normalize(p)
            if cosα > s.open
                r.orig = s.uncenter(p)
                return false
            end
        end
    end
    return true
end

"""
    Interface(normal::AffineMap, n::Float64)

Build an optical interface from a AffineMap that transforms a point on the ellipsoid of the interface to the normal at that point, `normal`, and the refractive index ratio between the inside and the outside of the ellipsoid, `n`.
"""
struct Interface
    normal::AffineMap{SDiagonal{3, Float64}, SVector{3,Float64}}
    n::Float64
    n2::Float64
    Interface(normal::AffineMap{SDiagonal{3,Float64},SVector{3,Float64}}, n::Float64) = new(normal, n, n^2)
end

"""
    OpticUnit(body::Ellipsoid, interface::Interface, register::Bool, name::String)

Build an optical unit. `register` indicates whether this unit should register intersection points (thus functions as a retina).
"""
struct OpticUnit
    body::Ellipsoid
    interface::Interface
    register::Bool
    name::String
end

"""
    OpticUnit(body::Ellipsoid, pointin::Bool, n::Float64, register::Bool, name::String)

Convenience function to build optical units. Build an optical unit depending on if the normal should be pointing in or out, `pointin`. The correct AffineMap transformation will be automatically calculated.
"""
function OpticUnit(body::Ellipsoid, pointin::Bool, n::Float64, register::Bool, name::String)
    i = ifelse(pointin, -1.0, 1.0)
    dir = LinearMap(SDiagonal{3, Float64}(i, i, i))
    normal = dir∘body.scale∘body.scale∘body.center
    interface = Interface(normal, n)
    return OpticUnit(body, interface, register, name)
end

"""
    bend!(r::Ray, i::Interface)

Refract or reflect a ray with an interface. Update the direction of the ray. Returns if event failed, which is always `false` (see `raytrace!` for details why that is so).
"""
function bend!(r::Ray, i::Interface)
    i.n == 1 && return false # we're done cause the refractive indices are the same
    N = normalize(i.normal(r.orig))
    a = -r.dir⋅N
    b = i.n2*(1 - a^2)
    dir = if b ≤ 1 # refract or reflect? refract!
        i.n*r.dir + (i.n*a - sqrt(1 - b))*N
    else # reflect!
        r.dir + 2a*N
    end
    r.dir = normalize(dir)
    return false
end

"""
    raytrace!(r::Ray, c::OpticUnit)

Advance a ray to the intersection point with an ellipsoid. If the intersection was successful, bend the ray. Updates the ray accordingly. Returns intersection failure.
"""
raytrace!(r::Ray, c::OpticUnit) = advance!(r, c.body) || bend!(r, c.interface)

end #module


text:
export getReadyForAuto!

auto_default_specs = Dict(
# Dimension of the system of equations as specified in the user-supplied routine FUNC
"NDIM" => 0,
# The number of boundary conditions as specified in the user-supplied routine BCND.
"NBC" => 0,
# The number of integral conditions as specified in the user-supplied routine ICND.
"NINT" => 0,
# Maximum parameter number that can be used in all user-supplied routines.
"NPAR" => 0,
# Used to indicate whether derivatives are supplied by the user or to be obtained by differencing
"JAC" => 0,
# The number of mesh intervals used for discretization.
"NTST" => 50,
# The number of Gauss collocation points per mesh interval, (2 ≤ NCOL ≤ 7).
"NCOL" => 4,
# This constant controls the mesh adaption
"IAD" => 3,
# Relative convergence criterion for equation parameters in the Newton/Chord method.
"EPSL" => 10e-7,
# Relative convergence criterion for solution components in the Newton/Chord method
"EPSU" =>10e-7,
# Relative arclength convergence criterion for the detection of special solutions.
"EPSS" =>10e-5,
# The maximum number of iterations allowed in the accurate location of special solutions,
# such as bifurcations, folds, and user output points, by Müller’s method with bracketing.
"ITMX" => 8,
# After NWTN Newton iterations the Jacobian is frozen, i.e., AUTO uses full Newton for the first
# NWTN iterations and the Chord method for iterations NWTN+1 to ITNW.
"NWTN" => 3,
# The maximum number of combined Newton-Chord iterations
"ITNW" =>5,
# AUTO uses pseudo-arclength continuation for following solution families. The pseudo-arclength
# stepsize is the distance between the current solution and the next solution on a family. By
# default, this distance includes all state variables (or state functions) and all free parameters.
# The constant DS defines the pseudo-arclength stepsize to be used for the first attempted step
# along any family. (Note that if IADS>0 then DS will automatically be adapted for subsequent
# steps and for failed steps.) DS may be chosen positive or negative; changing its sign reverses the
# direction of computation. The relation DSMIN ≤ | DS | ≤ DSMAX must be satisfied. The precise
# choice of DS is problem-dependent; the demos use a value that was found appropriate after some
# experimentation.
"DS" => 0.01,
# This is minimum allowable absolute value of the pseudo-arclength stepsize
"DSMIN" => 0.001,
# The maximum allowable absolute value of the pseudo-arclength stepsize
"DSMAX" =>0.1,
# This constant controls the frequency of adaption of the pseudo-arclength stepsize
"IADS" => 1,
# By default, the pseudo-arclength stepsize includes all state variables (or state functions) and all
# free parameters. Under certain circumstances one may want to modify the weight accorded to
# individual parameters in the definition of stepsize. For this purpose, THL defines the parameters
# whose weight is to be modified. If THL={} then all weights will have default value 1.0, else one
# must enter pairs, {Parameter Index : Weight, ...} .
"THL" => "",
# Under certain circumstances one may want to modify the weight accorded to individual state
# variables (or state functions) in the definition of stepsize. For this purpose, THU defines the
# number of states whose weight is to be modified. If THU={} then all weights will have default
# value 1.0, else one must enter pairs, {State Index : Weight, ...} . At present none of the demos
# use this option.
"THU" => "",
# This constant adds stopping conditions. It is specified as a list of bifurcation type strings fol-
# lowed by a number n greater than zero. These strings specify that the contination should stop
# as soon as the nth bifurcation of the associated type has been reached. Example:
# STOP=[’HB3’,’UZ3’] Stop at the third Hopf bifurcation or third user defined point (see Sec-
# tion 10.9.11), whichever comes first.
"STOP" => "[]",
#The maximum number of steps to be taken along any family.
"NMX" =>200,
# The lower bound on the principal continuation parameter. (This is the parameter which appears
# first in the ICP list; see Section 10.7.1.)
"RLO" =>0,
# The upper bound on the principal continuation parameter
"RL1" =>0,
# The lower bound on the principal solution measure.
"A0" =>0,
# The upper bound on the principal solution measure.
"A1" =>0,
)

function floatFromTerminal()
    return(float(split(readline(STDIN), "\n")[1]))
end

function stringFromTerminal()
    return(split(readline(STDIN), "\n")[1])
end

function getReadyForAuto!(M::Model; specs = auto_default_specs, interactive = true, sol = None)
    M.auto_specs = auto_default_specs
    M.auto_specs["NDIM"] = length(M.odes)
    M.auto_specs["NPAR"] = minimum([length(M.pars), 36])
    if interactive == true
        println("Which parameters do you want to vary?")
        par = stringFromTerminal()
        println("Minimum for this parameter?")
        M.auto_specs["RL0"] = floatFromTerminal()
        println("Maximum for this parameter?")
        M.auto_specs["RL1"] = floatFromTerminal()
        println("What's your prinicpal solution measure?")
        var = stringFromTerminal()
        println("Minimum for this variable?")
        M.auto_specs["A0"] = floatFromTerminal()
        println("Maximum for this variable?")
        M.auto_specs["A1"] = floatFromTerminal()
        println("stepsize (DS)?")
        M.auto_specs["DS"] = floatFromTerminal()
        println("Minimum stepsize (DSMin, absolute value)?")
        M.auto_specs["DSMIN"] = floatFromTerminal()
        println("Maximum stepsize (DSMAX, absolute value)?")
        M.auto_specs["DSMAX"] = floatFromTerminal()
    end
    if sol != None
        parseToAUTO(M, sol)
    end
    return(M)
end

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
using Documenter, Temporal

makedocs(modules=[Temporal],
         #doctest=false,
         #clean=true,
         #format=:html,
         authors="Jacob Amos",
         pages=Any["Home" => "index.md",
                   "TS Object" => "ts.md",
                   "Operations" => "operations.md",
                   "Data Access" => "io.md",
                   "Data Manipulation" => Any["indexing.md",
                                              "combining.md",
                                              "aggregation.md"]])

deploydocs(deps=Deps.pip("mkdocs", "python-markdown-math"),
           repo="github.com/dysonance/Temporal.jl.git",
           julia="0.6")


text:
import Bio.Align.BAM

#count = 0
#reader = open(BAM.Reader, "E:\\bam_test\\PD001\\1870157\\g_1870157_chr8_48711871.bam")

#for record in reader
#    if BAM.ismapped(record) && BAM.refname(record)=="chr8"
#        if BAM.position(record) > 48711900 && BAM.position(record) <=48711950
#            global count += 1
#            println(string(BAM.tempname(record), "\t", string(BAM.alignlength(record))))
#        end
#    end
#end

#println(count)

#close(reader)

function get_AlignLen(file_bam)
    count = 0
    reader = open(BAM.Reader, file_bam)

    for record in reader
        if BAM.ismapped(record) && BAM.refname(record)=="chr8"
            if BAM.position(record) > 48711900 && BAM.position(record) <=48711950
                count += 1
                println(string(BAM.tempname(record), "\t", string(BAM.alignlength(record))))
            end
        end
    end
    close(reader)
    return count
end

AlignNum =  get_AlignLen("E:\\bam_test\\X001\\1870157\\g_1870157_chr8_48711871.bam")
println(AlignNum)


text:
# to_drop = ["Population",
#         "General Expenditure, by Function:",
#         "Personal income",
#         "Total Expenditure - General Expenditure - Intergovernmental General Expenditure",
#         "Insurance Trust Expenditure - Unemployment Compensation Systems",
#         "Insurance Trust Expenditure - Workers' Compensation Systems",
#         "Insurance Trust Expenditure - State-Administered Pension Systems",
#         "Insurance Trust Expenditure - Other Insurance Trust Systems"]

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
import Dagger: @par, @spawn, spawn

@everywhere begin
    checkwid() = myid()==1
    function dynamic_fib(n)
        n <= 1 && return n
        t = Dagger.spawn(dynamic_fib, n-1)
        y = dynamic_fib(n-2)
        return (fetch(t)::Int) + y
    end
end

@testset "@par" begin
    @testset "per-call" begin
        x = 2
        a = @par x + x
        @test a isa Dagger.Thunk
        b = @par sum([x,1,2])
        c = @par a * b
        @test collect(c) == 20
    end
    @testset "block" begin
        c = @par begin
            x = 2
            a = x + x
            b = sum([x,1,2])
            c = a * b
        end
        @test x isa Int
        @test a isa Dagger.Thunk
        @test c isa Dagger.Thunk
        @test collect(c) == 20
    end
end

@testset "@spawn" begin
    @test Dagger.Sch.eager_context() === nothing
    @testset "per-call" begin
        x = 2
        a = @spawn x + x
        @test a isa Dagger.EagerThunk
        b = @spawn sum([x,1,2])
        c = spawn(*, a, b)
        @test c isa Dagger.EagerThunk
        @test fetch(a) == 4
        @test fetch(b) == 5
        @test fetch(c) == 20
    end
    @test Dagger.Sch.eager_context() isa Context
    @testset "waiting" begin
        a = @spawn sleep(1)
        @test !isready(a)
        wait(a)
        @test isready(a)
    end
    @testset "options" begin
        s = 1
        m = true
        a = @spawn single=s checkwid()
        b = @spawn meta=m ((_a)->_a isa Dagger.Chunk)(a)
        @test fetch(a)
        @test fetch(b)
    end
    @testset "errors" begin
        @testset "independent" begin
            a = @spawn error("Test")
            wait(a)
            @test isready(a)
            @test_throws_unwrap Dagger.ThunkFailedException fetch(a)
            b = @spawn 1+2
            @test fetch(b) == 3
        end
        @testset "direct vs indirect" begin
            a = @spawn error("Test")
            b = @spawn a+1

            ex = try
                fetch(a)
            catch err
                err
            end
            ex = Dagger.Sch.unwrap_nested_exception(ex)
            ex_str = sprint(io->Base.showerror(io,ex))
            @test occursin(r"^Thunk.*failure:", ex_str)
            @test occursin("Test", ex_str)
            @test !occursin("due to a failure in", ex_str)

            ex = try
                fetch(b)
            catch err
                err
            end
            ex = Dagger.Sch.unwrap_nested_exception(ex)
            ex_str = sprint(io->Base.showerror(io,ex))
            @test occursin(r"^Thunk.*failure due to a failure in", ex_str)
            @test occursin("Test", ex_str)
        end
        @testset "single dependent" begin
            a = @spawn error("Test")
            b = @spawn a+2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(a)
        end
        @testset "multi dependent" begin
            a = @spawn error("Test")
            b = @spawn a+2
            c = @spawn a*2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(b)
            @test_throws_unwrap Dagger.ThunkFailedException fetch(c)
        end
        @testset "dependent chain" begin
            a = @spawn error("Test")
            @test_throws_unwrap Dagger.ThunkFailedException fetch(a)
            b = @spawn a+1
            @test_throws_unwrap Dagger.ThunkFailedException fetch(b)
            c = @spawn b+2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(c)
        end
        @testset "single input" begin
            a = @spawn 1+1
            b = @spawn (a->error("Test"))(a)
            @test fetch(a) == 2
            @test_throws_unwrap Dagger.ThunkFailedException fetch(b)
        end
        @testset "multi input" begin
            a = @spawn 1+1
            b = @spawn 2*2
            c = @spawn ((a,b)->error("Test"))(a,b)
            @test fetch(a) == 2
            @test fetch(b) == 4
            @test_throws_unwrap Dagger.ThunkFailedException fetch(c)
        end
        @testset "diamond" begin
            a = @spawn 1+1
            b = @spawn a+1
            c = @spawn a*2
            d = @spawn ((b,c)->error("Test"))(b,c)
            @test fetch(a) == 2
            @test fetch(b) == 3
            @test fetch(c) == 4
            @test_throws_unwrap Dagger.ThunkFailedException fetch(d)
        end
    end
    @testset "remote spawn" begin
        a = fetch(Distributed.@spawnat 2 Dagger.spawn(+, 1, 2))
        @test Dagger.Sch.EAGER_INIT[]
        @test fetch(Distributed.@spawnat 2 !(Dagger.Sch.EAGER_INIT[]))
        @test a isa Dagger.EagerThunk
        @test fetch(a) == 3

        # Mild stress-test
        @test dynamic_fib(10) == 55
    end
end


text:
using Coexist, PyCall
Coexist.pyCoexist_setup()

pydir = joinpath(@__DIR__,"..", "coexist_python", "model_COVID_testing-ODEoutput.py")
py"""
import os
exec(open(os.path.join($pydir)).read())
"""

state = 50*ones(9*8*4*4)
dstate = similar(state)

function solveSystem(
	state,
	timeSpan = (0.0, 80.0),
	p=nothing
	)
	prob = ODEProblem(Coexist.dydt,state,(0.0,80.0),p=nothing)
	sol = solve(prob,Tsit5(),reltol=1e-3,abstol=1e-3)
	sol = convert(Array, sol)
	return reshape(sol, (4,4,8,9, Int64(prod(size(sol))/(4*4*8*9))))
end

dydt_benchmark = @benchmark Coexist.dydt(dstate,state,0,nothing)
solveSystem_benchmark = @benchmark solveSystem(state)
py_dydt_benchmark = @benchmark py"solveSystem(state, 80, **paramDict_current)"
py_solveSystem_benchmark = @benchmark py"solveSystem(state, 80, **paramDict_current)"

println("BENCHMARK OF dydt_Complete")
display(dydt_benchmark)
println()
println("BENCHMARK OF solveSystem")
display(solveSystem_benchmark)
println()
println("BENCHMARK OF py_dydt_Complete")
display(py_dydt_benchmark)
println()
println("BENCHMARK OF py_solveSystem")
display(py_solveSystem_benchmark)
println()

max_stars_repo_name: vollmersj/coexist-julia

max_forks_repo_name: vollmersj/coexist-julia


text:
# ---
# title: 1004. Max Consecutive Ones III
# id: problem1004
# author: Indigo
# date: 2020-01-30
# difficulty: Medium
# categories: Two Pointers, Sliding Window
# link: <https://leetcode.com/problems/max-consecutive-ones-iii/description/>
# hidden: true
# ---
#
# Given an array `A` of 0s and 1s, we may change up to `K` values from 0 to 1.
#
# Return the length of the longest (contiguous) subarray that contains only 1s.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1,1,1,0,0,0,1,1,1,1,0], K = 2
#     Output: 6
#     Explanation:
#     [1,1,1,0,0, _ **1** ,1,1,1,1, **1**_ ]
#     Bolded numbers were flipped from 0 to 1.  The longest subarray is underlined.
#
# **Example 2:**
#
#
#
#     Input: A = [0,0,1,1,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1], K = 3
#     Output: 10
#     Explanation:
#     [0,0, _1,1, **1** , **1** ,1,1,1, **1** ,1,1_,0,0,0,1,1,1,1]
#     Bolded numbers were flipped from 0 to 1.  The longest subarray is underlined.
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 20000`
#   2. `0 <= K <= A.length`
#   3. `A[i]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

function longest_ones(A::Vector{Int}, K::Int)
    l = r = 1
    len = length(A)
    cnt = res = 0
    while r <= len
        (A[r] == 0) && (cnt += 1)
        r += 1
        while cnt > K
            (A[l] == 0) && (cnt -= 1)
            l += 1
        end
        res = max(res, r - l)
    end
    return res
end
## @lc code=end


text:
function babylonian_sqrt(S)
    # Your code goes here!
    return 0
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
cd(@__DIR__); include("setups/rgc100.jl")
gr(dpi = 200)

## construct full HGLET and Lapped-HGLET dictionaries
@time HGLET_dic = HGLET_dictionary(GP, G_Sig; gltype = :L)
@time LPHGLET_dic = LPHGLET_dictionary(GP, G_Sig; gltype = :L, ϵ = 0.3)

## generate figures
# pre-selected (j, k, l)s
T = [(2, 1, 1), (2, 1, 12), (3, 1, 3), (3, 1, 5),
    (4, 5, 6), (6, 23, 3), (7, 5, 2), (8, 15, 1)]

XLs = [
    [-50, 50], [-50, 50], [-110, 40], [-110, 40],
    [-20, 50], [50, 120], [-80, -20], [-40, -10]
]
YLs = [
    [-50, 50], [-50, 50], [-50, 100], [-50, 100],
    [-50, 20], [-70, 0], [-90, -30], [-160, -130]
]

for i in 1:length(T)
    (j, k, l) = T[i]
    𝚽H = HGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
    𝚽lH = LPHGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
    p1 = gplot(W, X; width = 1)
        scatter_gplot!(X; marker = 𝚽H[:, l], plotOrder = :s2l, ms = 5)
        plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = false)
        plot!(left_margin = 3mm, bottom_margin = 3mm, clims = (-0.05, 0.05))
    p2 = gplot(W, X; width = 1)
        scatter_gplot!(X; marker = 𝚽lH[:, l], plotOrder = :s2l, ms = 5)
        plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = false)
        plot!(left_margin = 3mm, bottom_margin = 3mm, clims = (-0.05, 0.05))
    plt = plot(p1, p2, size = (640, 300), layout = Plots.grid(1, 2),
        xlims = XLs[i], ylims = YLs[i])
        # xlims = [-120, 120], ylims = [-150, 150])
    # gplot(W, X; width = 1)
    # ord = :s2l
    # if i == 8
    #     ord = :l2s
    # end
    # scatter_gplot!(X; marker = 𝚽lH[:, l] - 𝚽H[:, l], plotOrder = ord, ms = 4)
    # plot!(xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = true)
    # plot!(left_margin = 3mm, bottom_margin = 3mm, right_margin = 3mm)
    # plt = plot!(size = (500, 500))
    savefig(plt, "../figs/RGC100_HGLET_LPHGLET_j$(j-1)k$(k-1)l$(l-1).png")
end