text:
"""
    TableDataset(table)
    TableDataset(path::AbstractString)

Wrap a Tables.jl-compatible `table` as a dataset container.
Alternatively, specify the `path` to a CSV file directly
to load it with CSV.jl + DataFrames.jl.
"""
struct TableDataset{T} <: AbstractDataContainer
    table::T

    # TableDatasets must implement the Tables.jl interface
    function TableDataset{T}(table::T) where {T}
        Tables.istable(table) ||
            throw(ArgumentError("TableDatasets must implement the Tabels.jl interface"))

        new{T}(table)
    end
end

TableDataset(table::T) where {T} = TableDataset{T}(table)
TableDataset(path::AbstractString) = TableDataset(DataFrame(CSV.File(path)))

# slow accesses based on Tables.jl
_getobs_row(x, i) = first(Iterators.peel(Iterators.drop(x, i - 1)))
function _getobs_column(x, i)
    colnames = Tuple(Tables.columnnames(x))
    rowvals = ntuple(j -> Tables.getcolumn(x, j)[i], length(colnames))

    return NamedTuple{colnames}(rowvals)
end
function Base.getindex(dataset::TableDataset, i)
    if Tables.rowaccess(dataset.table)
        return _getobs_row(Tables.rows(dataset.table), i)
    elseif Tables.columnaccess(dataset.table)
        return _getobs_column(dataset.table, i)
    else
        error("The Tables.jl implementation used should have either rowaccess or columnaccess.")
    end
end
function Base.length(dataset::TableDataset)
    if Tables.columnaccess(dataset.table)
        return length(Tables.getcolumn(dataset.table, 1))
    elseif Tables.rowaccess(dataset.table)
        # length might not be defined, but has to be for this to work.
        return length(Tables.rows(dataset.table))
    else
        error("The Tables.jl implementation used should have either rowaccess or columnaccess.")
    end
end

# fast access for DataFrame
Base.getindex(dataset::TableDataset{<:DataFrame}, i) = dataset.table[i, :]
Base.length(dataset::TableDataset{<:DataFrame}) = nrow(dataset.table)

# fast access for CSV.File
Base.getindex(dataset::TableDataset{<:CSV.File}, i) = dataset.table[i]
Base.length(dataset::TableDataset{<:CSV.File}) = length(dataset.table)

## Tables.jl interface

Tables.istable(::TableDataset) = true
for fn in (:rowaccess, :rows, :columnaccess, :columns, :schema, :materializer)
    @eval Tables.$fn(dataset::TableDataset) = Tables.$fn(dataset.table)
end


text:
function BinTreeCall(K, T, S, r, N, u, d)
        dt = T/N ## Time to Maturity/ Number of Steps
        p = (exp(r * dt) - d) / (u - d)
        disc = exp(-r * dt)

        St = zeros(N+1)

        St[1] = S * d^N

        for j in 2:(N+1)
                St[j] = St[j-1] * u/d
        end

        St = flipdim(St, 1)

        C = Array(Float64, (N+1, N+1))#linspace(1, N+1, N+1)

        for j in 1:(N+1)
                C[j, N+1] = max(0.0, St[j] - K)
        end

        for i = N:-1:1
                for j in 1:i
                      C[j, i] = disc * (p * C[j, i+1] + (1-p) * C[j+1, i+1])
                end
        end

C
end

BinTreeCall(100,1,100,0.06,3,1.1,0.9091)

ext: jl


text:
# Variable names have the following requirements:
# Starting with a letter
# Continuing with letters, digits, underscores, or exclamation marks

# In these examples, the Roman alphabet letters are used, but many other letters are accepted
# Best practice is to use lower case letters and digits with descriptive names

# Julia requires that the variable be on the left and the value be on the right

my_name = "My name is Viki"
println(my_name)

my_name1 = "My first name is Viktoria"
println(my_name1)

my_name2 = "My last name is Haghani"
println(my_name2)

vikis_pets = "I have four pets: Stella, Zoscia, Sunny, and Blueberry"
println(vikis_pets)

weird_l00king_var!abl3_name = "weird_l00king_var!abl3_name is a valid variable name in Julia"
println(weird_l00king_var!abl3_name)

# If a comma follows a valid variable name on the left of the assignment operator, the commma indicates that Julia should perform multiple assignment
# This assigns the variable to each character of the string in the order specified
stringexample = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
Hello, world, we, are, going, in, this, order = stringexample
println(Hello)
println(world)
println(we)
println(are)
println(going)
println(in)
println(this)
println(order)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
#= Interpolations.jl and ForwardDiff.jl are removed as dependencies in this version.
"""
    function_to_interpolated_function(f_xy; physwidth = 10.0m, physheight = 4.0m)

Input: f: (Q, Q)  → (Q, Q) Function taking a tuple quantities, outputs a tuple of quantities
    physwidth:  Q
    physheight: Q

Output: f: (Q, Q)  → (Q, Q) Function taking a tuple quantities, outputs a tuple of quantities

Convert the function f_xy into an interpolated function. The output function is linearly interpolated
between pixels, based on the currently defined sketch scale.

The output of the function can be otherwise; this is just the originally intended use.
"""
function function_to_interpolated_function(f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true)
    # Position iterators
    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
    tuplemat = [f_xy(x, y) for x in xs, y in reverse(ys)]
    fxy_inter = interpolate((xs, reverse(ys)),
        tuplemat,
        Gridded(Linear()));
    # Extend the domain, using the same values as on the domain border. Don't know how the closest point
    # on domain boundary is found.
    extrapolate(fxy_inter, Flat())
end

"""
    matrix_to_function(matrix::Array{Quantity{Complex{Float64}, D, U}, 2}) where {D, U}

Input: CQ²

Output: f: (Q, Q)  → (Q, Q) Function taking coordinates as a complex quantity, outputs a tuple of quantities

where
    Q = Quantity
    CQ² = Matrix of complex quantities

Convert the matrix with complex quantity elements into a vector-valued function.
The output function is linearly interpolated between pixels, based on the currently defined sketch scale.
"""
function matrix_to_function(matrix::Array{Quantity{Complex{Float64}, D, U}, 2}) where {D, U}
    # Image processing convention: a column has a horizontal line of pixels
    ny, nx = size(matrix)
    # We assume every element correspond to a position, unit of length
    physwidth = nx * scale_pt_to_unit(m)
    physheight = ny * scale_pt_to_unit(m)
    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
    # Adapt to Interpolations
    tuplemat = map( cmplx -> (real(cmplx), imag(cmplx)), transpose(matrix)[ : , end:-1:1])
    # Function that can interpolate between nodes.
    fxy_inter = interpolate((xs, reverse(ys)),
        tuplemat,
        Gridded(Linear()));
    # Extend the domain, using the same values as on the domain border. Don't know how the closest point
    # on domain boundary is found.
    extrapolate(fxy_inter, Flat())
end

"""
    matrix_to_function(matrix::Array{Float64,2})

Input: R²

Output: f: (R, R) → R function that interpolates between coordinates

where
    R = Real number
    R² = Matrix of real numbers

    (Q, Q) Tuple of quantities

Convert the matrix with real-valued elements into a real-valued function. The domain of the function
is quantities, i.e. two position coordinates. The size of the domain is taken from current sketch scale.
The output function is linearly interpolated between pixels, based on the currently defined sketch scale.
"""
function matrix_to_function(matrix)   # temp ::Array{Float64,2})
    # Image processing convention: a column has a horizontal line of pixels
    ny, nx = size(matrix)
    # We assume every element correspond to a position, unit of length
    physwidth = nx * scale_pt_to_unit(m)
    physheight = ny * scale_pt_to_unit(m)

    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
    # Adapt to Interpolations
    transposed_matrix = transpose(matrix)[ : , end:-1:1]
    # Function that can interpolate between nodes.
    fxy_inter = interpolate((xs, reverse(ys)),
        transposed_matrix,
        Gridded(Linear()));
    # Extend the domain, using the same values as on the domain border. Don't know how the closest point
    # on domain boundary is found.
    extrapolate(fxy_inter, Flat())
end

=#


text:
module HapiTests

using Hapi
using Hapi.FilePaths
using ReTest

@testset "Hapi.jl" begin
    rgx = build_regex("{thing1}_{thing2}.txt", (thing1=raw"\d+", thing2=raw"\w+"))
    @test rgx isa Regex
    @test occursin(rgx, "file1_a.txt")
    @test !occursin(rgx, "testa_1.txt")

    m1 = match(rgx, "file1_a.txt")
    @test m1[:thing1] == "1"
    @test m1[:thing2] == "a"

    glb = glob_pattern("projects/fileglob/", rgx)
    fd = first(glb)
    @test fd isa FileDependency
    @test path(fd) isa AbstractPath
    @test params(fd) isa NamedTuple
    @test collect(keys(params(fd))) == [:thing1, :thing2]
    @test fd[:thing1] == "1"
    @test fd[:thing2] == "1"

    @test length(glb) == 4
    @test all(f-> f isa FileDependency, glb)
end

end # module


text:
## Display string vector on fixed with screen
function show_string_vector(sV :: Vector{T1},  width :: T2 = 80;
    io :: IO = stdout)  where {T1 <: AbstractString,  T2 <: Integer}

    n = length(sV);
    if n < 1
        return nothing
    end

    iCol = 0;
    for s in sV
        len1 = length(s);
        if len1 > 0
            if iCol + len1 > width
                println(io, " ")
                iCol = 0;
            end
            print(io, s);
            iCol = iCol + len1;
            print(io, "    ");
            iCol = iCol + 4;
        end
    end
    println(io, " ")
end

# -------------


text:
module PlanarConvexHullsTest

using PlanarConvexHulls
using StaticArrays
using Test
using Random
using LinearAlgebra
using Statistics

using PlanarConvexHulls: vertex_order

const Point{T} = SVector{2, T}

@testset "is_ordered_and_strongly_convex" begin
    @test is_ordered_and_strongly_convex([Point(1, 2)], CCW)
    @test is_ordered_and_strongly_convex([Point(1, 2), Point(3, 4)], CCW)
    @test is_ordered_and_strongly_convex([Point(1, 2), Point(3, 4), Point(2, 4)], CCW)
    @test !is_ordered_and_strongly_convex([Point(1, 2), Point(3, 4), Point(2, 3)], CCW) # on a line

    v1 = Point(0, 0)
    v2 = Point(1, 0)
    v3 = Point(1, 1)
    v4 = Point(0, 1)
    vertices = [v1, v2, v3, v4]
    for i in 0 : 4
        shifted = circshift(vertices, i)
        @test is_ordered_and_strongly_convex(shifted, CCW)
        @test !is_ordered_and_strongly_convex(reverse(shifted), CCW)
        @test is_ordered_and_strongly_convex(reverse(shifted), CW)
    end

    for i in eachindex(vertices)
        for j in eachindex(vertices)
            if i != j
                vertices′ = copy(vertices)
                vertices′[i] = vertices[j]
                vertices′[j] = vertices[i]
                @test !is_ordered_and_strongly_convex(vertices′, CCW)
            end
        end
    end
end

@testset "area" begin
    @test area(ConvexHull{CCW}(SVector((SVector(1, 1),)))) == 0

    @test area(ConvexHull{CCW}(SVector(SVector(1, 1), SVector(2, 3)))) == 0

    triangle = ConvexHull{CCW}(SVector(Point(1, 1), Point(2, 1), Point(3, 3)))
    @test area(triangle) == 1.0

    square = ConvexHull{CCW}(SVector(Point(1, 1), Point(4, 1), Point(4, 3), Point(1, 3)))
    @test area(square) == 3 * 2
end

@testset "in" begin
    @testset "point" begin
        rng = MersenneTwister(1)
        p = SVector(1, 1)
        C = ConvexHull{CCW}(SVector((p,)))
        @test p ∈ C
        @test Float64.(p) ∈ C
        for i in 1 : 10
            @test p + SVector(randn(rng), randn(rng)) ∉ C
        end
    end

    @testset "line segment" begin
        p1 = SVector(1, 1)
        p2 = SVector(3, 5)
        linesegment = ConvexHull{CCW}([p1, p2])
        @test p1 ∈ linesegment
        @test p2 ∈ linesegment
        @test div.(p1 + p2, 2) ∈ linesegment
        @test p1 + 2 * (p2 - p1) ∉ linesegment
    end

    @testset "triangle" begin
        rng = MersenneTwister(1)
        triangle = ConvexHull{CCW}(SVector(Point(1, 1), Point(2, 1), Point(3, 3)))
        for p in vertices(triangle)
            @test p ∈ triangle
        end
        for i = 1 : 100_000
            weights = normalize(rand(rng, SVector{3}), 1)
            p = reduce(+, vertices(triangle) .* weights)
            @test p ∈ triangle
        end
    end

    @testset "rectangle" begin
        rng = MersenneTwister(1)
        width = 4
        height = 3
        origin = Point(2, 4)
        rectangle = ConvexHull{CCW}(map(x -> x + origin, SVector(Point(0, 0), Point(width, 0), SVector(width, height), SVector(0, height))))
        for p in vertices(rectangle)
            @test p ∈ rectangle
        end
        for i = 1 : 100_000
            p = origin + SVector(width * rand(rng), height * rand(rng))
            @test p ∈ rectangle
        end
        for i = 1 : 10
            p = origin + SVector(width * rand(rng), height * rand(rng))
            @test setindex(p, origin[1] + width + rand(rng), 1) ∉ rectangle
            @test setindex(p, origin[1] - rand(rng), 1) ∉ rectangle
            @test setindex(p, origin[2] + height + rand(rng), 2) ∉ rectangle
            @test setindex(p, origin[2] - rand(rng), 2) ∉ rectangle
        end
    end
end

@testset "centroid" begin
    @testset "point" begin
        p = Point(1, 2)
        @test centroid(ConvexHull{CCW}([p])) === Float64.(p)
    end

    @testset "line segment" begin
        p1 = SVector(1, 1)
        p2 = SVector(3, 5)
        linesegment = ConvexHull{CCW}([p1, p2])
        @test centroid(linesegment) == Point(2.0, 3.0)
    end

    @testset "triangle" begin
        triangle = ConvexHull{CCW}(SVector(Point(1, 1), Point(2, 1), Point(3, 3)))
        @test centroid(triangle) ≈ mean(vertices(triangle)) atol=1e-15
    end
end

function convex_hull_alg_test(hull_alg!)
    @testset "random $order" for order in [CCW, CW]
        hull = ConvexHull{order, Float64}()
        rng = MersenneTwister(2)
        for n = 1 : 10
            sizehint!(hull, n) # just to get code coverage
            for _ = 1 : 10_000
                points = [rand(rng, Point{Float64}) for i = 1 : n]
                hull_alg!(hull, points)
                @test is_ordered_and_strongly_convex(vertices(hull), order)
            end
        end
    end

    @testset "collinear input $order" for order in [CCW, CW]
        hull = ConvexHull{order, Float64}()
        points = [Point(0., 0.), Point(0., 1.), Point(0., 2.), Point(1., 0.), Point(1., 1.), Point(1., 2.)]
        for i = 1 : 10
            shuffle!(points)
            hull_alg!(hull, points)
            @test is_ordered_and_strongly_convex(vertices(hull), order)
            @test isempty(symdiff(vertices(hull), [Point(0., 0.), Point(1., 0.), Point(1., 2.), Point(0., 2.)]))
        end
    end

    @testset "numerical problem $order" for order in [CCW, CW]
        hull = ConvexHull{order, Float64}()
        points = Point{Float64}[
            [-0.17559141363793285, 0.1777893970092983],
            [-0.17574071124767324, 0.04921003844438379],
            [0.08480367424885946, 0.17745906024313485],
            [0.08465437663911904, 0.04887970167822034],
            [-0.13059226250542266, 0.17773231116261107],
            [-0.13074156011516305, 0.04915295259769656],
            [-0.1753504285942606, -0.1772960226863406],
            [-0.1756647405301049, -0.04875022362480264],
            [0.08503864201675085, -0.17672638211747324],
            [0.08472433008090652, -0.04818058305593528],
            [-0.13035231731355584, -0.17719758249586814],
            [-0.13066662924940015, -0.04865178343433017]]
        hull_alg!(hull, points)
        @test is_ordered_and_strongly_convex(vertices(hull), order)
    end
end

@testset "jarvis_march!" begin
    convex_hull_alg_test(jarvis_march!)
end

@testset "closest_point $order" for order in [CCW, CW]
    hull = ConvexHull{order, Float64}()
    rng = MersenneTwister(3)
    for n = 1 : 10
        for _ = 1 : 100
            points = [rand(rng, Point{Float64}) for i = 1 : n]
            jarvis_march!(hull, points)

            for point in vertices(hull)
                closest = closest_point(point, hull)
                @test closest == point
            end

            for _ = 1 : 100
                point = rand(rng, Point{Float64})
                closest = closest_point(point, hull)
                if point ∈ hull
                    @test closest == point
                else
                    closest_to_closest = closest_point(closest, hull)
                    @test closest_to_closest ≈ closest atol=1e-14
                end
            end
        end
    end

    # visualize = true
    # if visualize
    #     flush(stdout)
    #     hull = ConvexHull{CCW, Float64}()
    #     rng = MersenneTwister(3)
    #     n = 5
    #     points = [rand(rng, Point{Float64}) for i = 1 : n]
    #     jarvis_march!(hull, points)
    #     projections = map(1 : 1_000_000) do _
    #         closest_point(rand(rng, Point{Float64}), hull)
    #     end
    #     plt = scatterplot(getindex.(projections, 1), getindex.(projections, 2))
    #     scatterplot!(plt, getindex.(vertices(hull), 1), getindex.(vertices(hull), 2), color = :red)
    #     display(plt)
    # end
end

function make_oversized(A::AbstractMatrix, b::AbstractVector)
    n = length(b)
    n_oversized = n + 1
    similar(A, (n_oversized, 2)), similar(b, n_oversized)
end

function make_oversized(A::SMatrix, b::SVector{n}) where n
    n_oversized = n + 1
    similar(A, Size(n_oversized, 2)), similar(b, Size(n_oversized))
end

function hreptest(hull::H, rng) where {H<:ConvexHull}
    A, b = hrep(hull)
    A_oversized, b_oversized = make_oversized(A, b)
    hrep!(A_oversized, b_oversized, hull)
    for _ = 1 : 100
        testpoint = rand(rng, Point{Float64})
        @test (testpoint ∈ hull) == all(A * testpoint .<= b) == all(A_oversized * testpoint .<= b_oversized)
    end
    if Length(vertices(hull)) !== Length{StaticArrays.Dynamic()}()
        @test(@allocated(hrep(hull)) == 0)
    end
end

@testset "hrep $order" for order in [CCW, CW]
    dynamichull = ConvexHull{order, Float64}()
    rng = MersenneTwister(4)
    for n = 2 : 10
        for _ = 1 : 100
            points = [rand(rng, Point{Float64}) for i = 1 : n]
            jarvis_march!(dynamichull, points)
            hreptest(dynamichull, rng)

            h = num_vertices(dynamichull)
            statichull = ConvexHull{order}(SVector{h}(vertices(dynamichull)))
            hreptest(statichull, rng)
        end
    end
end

struct CustomPoint{T} <: FieldVector{2, T}
    x::T
    y::T
end

@testset "Constructors" begin
    width = 4
    height = 3
    origin = Point(2, 4)
    verts = map(x -> x + origin, [Point(0, 0), Point(width, 0), SVector(width, height), SVector(0, height)])
    rectangle = ConvexHull{CW}(reverse(verts))

    srectangle = SConvexHull{4, Float64}(rectangle)
    drectangle = DConvexHull{Float64}(rectangle)
    @test vertices(srectangle) isa SVector{4}
    @test vertices(drectangle) isa Vector
    for hull in [srectangle, drectangle]
        @test vertex_order(hull) == CCW
        @test eltype(hull) == Float64
        @test is_ordered_and_strongly_convex(vertices(hull), vertex_order(hull))
    end

    let hull = ConvexHull{CCW, Float64}(verts)
        @test vertices(hull) == verts
        @test eltype(hull) == Float64
    end
    let hull = ConvexHull{CCW, Float64, CustomPoint{Float64}}(verts)
        @test vertices(hull) == verts
        @test eltype(hull) == Float64
        @test typeof(vertices(hull)) == Vector{CustomPoint{Float64}}
    end
    let hull = ConvexHull{CCW, Float64, Point{Float64}, SVector{4, Point{Float64}}}(verts)
        @test vertices(hull) == verts
        @test eltype(hull) == Float64
        @test typeof(vertices(hull)) == SVector{4, SVector{2, Float64}}
    end

    @test isempty(vertices(DConvexHull{Float64}()))
end

@testset "benchmarks" begin
    include(joinpath("..", "perf", "runbenchmarks.jl"))
end

end # module


text:
using SparseArrays

# I am using a separate type here because I want to preserve my ability to load
# Histogram types using JLD2. Oh boy...

export SparseHistogram, to_sparse, set_threshold!, drop_below_threshold!,
  sparsity

# Construct a sparse vector based on the input data
function to_sparse(x; threshold)
  ret = spzeros(length(x))
  for (ii, val) in enumerate(x)
    if val > threshold
      push!(ret.nzind, ii)
      push!(ret.nzval, val)
    end
  end

  ret
end

mutable struct SparseHistogram{RangeType <: Real}
  range::Tuple{Vector{RangeType},Vector{RangeType}}
  data::SparseVector{Float64, Int64}

  buffer::SparseVector{Float64, Int64}

  # Threshold defines the minimum value to keep in the filter
  threshold::Float64

  # The default constructor specializes for when we can infer the type of the
  # range from the signature. (Alternatively, see the outer constructor).
  function SparseHistogram(range::Tuple{Vector{D}, Vector{D}},
                           data::SparseVector;
                           threshold) where D <: Real
    new{D}(range, data, spzeros(length(data)), threshold)
  end
end

AnyHistogram = Union{SparseHistogram, Histogram}

# Constructor for when the input data is not sparse
function SparseHistogram(range, data; threshold)
  SparseHistogram(range, to_sparse(data, threshold=threshold);
                  threshold=threshold)
end

# Outer constructor that defers determination of the histogram type until after
# pulling the ranges
function SparseHistogram(range, data::SparseVector; threshold)
  SparseHistogram(map(collect, range), data; threshold=threshold)
end

# Copy constructor. Note that this only duplicates the data.
# We assume that nobody is crazy enough to modify the range.
duplicate(x::SparseHistogram) = SparseHistogram(x)
# Note, this will fail if no threshold is given when converting a Histogram
function SparseHistogram(x::AnyHistogram; threshold = x.threshold)
  SparseHistogram(get_range(x), to_sparse(get_values(x), threshold=threshold),
                  threshold=threshold)
end

set_threshold!(x::SparseHistogram; threshold) = (x.threshold = threshold)

# Returns the data matrix
size(x::SparseHistogram) = map(length, x.range)
get_data(x::SparseHistogram) = reshape(x.data, size(x))

#
# Sparse helper methods
#

function clear!(x::SparseVector)
  resize!(x.nzind, 0)
  resize!(x.nzval, 0)
end

function copy_sparse!(out, in, threshold)
  nzind = in.nzind
  nzval = in.nzval
  n = length(nzind)

  out_ind = out.nzind
  out_val = out.nzval

  # Presume that we will copy most of the input
  resize!(out_ind, n)
  resize!(out_val, n)

  # Copy values above the threshold
  end_ind = 1
  @inbounds for ii in 1:n
    if nzval[ii] > threshold
      out_ind[end_ind] = nzind[ii]
      out_val[end_ind] = nzval[ii]

      end_ind += 1
    end
  end

  resize!(out_ind, end_ind-1)
  resize!(out_val, end_ind-1)

  nothing
end

# Copies filter data, keeping nonzeros
function copy_filter!(x::SparseHistogram; out::SparseHistogram)
  copy_sparse!(out.data, x.data, x.threshold)
  clear!(out.buffer)

  out.threshold = x.threshold
end

# Remove values from the histogram below a given threshold
#
# Trim resizes the matrix. We will generally continue to operate in place so we
# will end up using the empty space
function drop_below_threshold!(x::SparseHistogram;
                               threshold = x.threshold,
                               trim = true
                              )
  droptol!(x.data, threshold, trim = trim)
  droptol!(x.buffer, threshold, trim = trim)
end
drop_below_threshold!(x::Histogram; kwargs...) = nothing

sparsity(x::SparseHistogram) =  1.0 - nnz(get_values(x)) / length(get_values(x))

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
## Solve the PBCC objective for a range of parameter settings on a suite of small graphs
include("../src/PBCC.jl")
using MAT
using SparseArrays

graphs = ["Cities_H", "Newsgroups100", "Zoo_H", "Amazon_Gift_Cards", "Amazon_Appliances", "Amazon_Amazon_Fashion"]

## Beta experiments
for graph = graphs
    M = matread("data/"*graph*".mat")
    H = M["H"]
    # For each dataset, solve for a range of beta values, and store X
    betas = 0.05:0.05:0.95
    mu = 0.0
    for beta = betas
        outputmat = "Output/"*graph*"/beta_$beta"*"_mu_$mu.mat"
        println(outputmat)
        start = time()
        X, LPbound = PBCC_LP(H,mu,beta,false,false)
        runtime = time()-start
        delta = 2/5
        C, Labels, obj = ManyRoundX(H,X,delta,mu,beta,5)
        x,y = Biclustering(H,Labels)
        ratio = obj/LPbound
        matwrite(outputmat,Dict("X"=>X,"LPbound"=>LPbound,"C"=>C,"Labels"=>Labels,
        "row_labels"=>x, "col_labels"=>y, "obj"=>obj,"ratio"=>ratio,"runtime"=>runtime))

    end

end

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using Random, OrdinaryDiffEq
import GLMakie

function FitzHugh_Nagumo_ODE(du,u,p,t)
    a, b, d, epsilon, hsq6 = p
    uu = @view u[1,:,:]
    vv = @view u[2,:,:]

    # This code can be optimized a lot if time allows
    # Create padded u matrix to incorporate Newman boundary conditions - 9-point stencil
    uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ]
    diff_term = d .* (
        4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+
            uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+
            uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+
            uuu[1:end-2,1:end-2] .- 20 .*uu
    ) ./ hsq6

    du[1,:,:] = a .* uu .*(1 .-uu).*(uu.-b) .- vv .+ diff_term
    du[2,:,:] = epsilon .* (uu .- vv )
end

# --------------------------  main program  -----------------------------

@time begin

# parameters
L = 300   # size of domain
N = 150   # no. of grid points
h = L/N # spatial steps size
hsq6 = 6*h*h

# FHN parameters
a = 3.
b = 0.2
d = 1.   #  diffusion
epsilon = 0.01

# open figure
kplt = 0  # no. of subplot

allsols = []
allts = []

for ic = 1:2    # initial conditions ( = rows of the plot )

println("ic=",ic)

uu = zeros(N,N)
vv = zeros(N,N)
Random.seed!(13371)
uu = uu .+ 0.01*randn(N,N)    # random perturbation

# local stimulius = concentric waves
if ic == 1
   uu[40:41,40:41] .= 0.999
   uu[75:76,75:76] .= 0.999
end

# local stimulius = spiral waves
if ic == 2
   uu[1:end,10:11] .= 0.999
end

# initial values
u0 = zeros(2,N,N)
u0[1,:,:] = uu
u0[2,:,:] = vv

p = a, b, d, epsilon, hsq6

if ic == 1
    saveat = 0.0:2:600
    tspan = (saveat[1], saveat[end])

    prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
    sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat)
    tvec = sol.t
    uout = [u[1, :, :] for u in sol.u]
    push!(allsols, uout)
    push!(allts, tvec)

elseif ic == 2
    saveat = 0.0:2.0:300.0
    tspan = (saveat[1], saveat[end])
    prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
    sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat)

    uout = [u[1, :, :] for u in sol.u]
    # second perturbation line
    u = sol.u[end]
    u[1, 55:73,1:70] .= 0.999
    saveat = 300.0:2:900
    tspan = (saveat[1], saveat[end])
    prob = ODEProblem(FitzHugh_Nagumo_ODE, u, tspan, p)
    sol = solve(prob, Tsit5(); reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat)

    tvec = sol.t
    uout2 = [u[1, :, :] for u in sol.u]
    append!(uout, uout2)
    push!(allsols, uout)
    push!(allts, tvec)

end
end # ic loop
end # cputime

# %% Animate
names = ("concentric", "plane_to_spiral")

for i in (1,2)
    name = names[i]
    uout = allsols[i]
    tvec = allts[i]

    heatobs = GLMakie.Observable(uout[1])
    tobs = GLMakie.Observable(0.0)
    titobs = GLMakie.lift(t -> "t = $(t)", tobs)

    fig = GLMakie.Figure(resolution = (600, 550))
    ax = fig[1,1] = GLMakie.Axis(fig; title = titobs)
    hmap = GLMakie.heatmap!(ax, heatobs; colormap = :tokyo, colorrange = (-0.2, 1))
    cb = GLMakie.Colorbar(fig[1, 2], hmap; width = 20)
    display(fig)

    GLMakie.record(
        fig, projectdir("animations", "11", "fitzhugh_$(name).mp4"),
        1:length(tvec); framerate = 15
    ) do i

        tobs[] = tvec[i]
        heatobs[] = uout[i]
    end

end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
facts("fast exponentiation methods") do

  context("odd") do

    @fact odd(0) => false
    @fact odd(1) => true
    @fact odd(2) => false
    @fact odd(3) => true

    @fact odd(2^100) => false

    @fact odd(BigInt(2)^2011) => false

  end

  @fact fast_exponentiation(10, 0, 2) => 1

  @fact fast_exponentiation(1, 2, 3) => 1

  @fact fast_exponentiation(2, 2011, 15) => int(mod(BigFloat(2)^2011, 15))

end


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

function _sphdist{T<:AbstractFloat}(long1::T, lat1::T, long2::T, lat2::T,
                                    degrees::Bool)
    # Convert both points to rectangular coordinates.
    rxy, z1 = polrec(1.0, lat1,  degrees=degrees)
    x1, y1  = polrec(rxy, long1, degrees=degrees)
    rxy, z2 = polrec(1.0, lat2,  degrees=degrees)
    x2, y2  = polrec(rxy, long2, degrees=degrees)
    # Compute vector dot product for both points.
    cs = x1*x2 + y1*y2 + z1*z2
    # Compute the vector cross product for both points.
    xc = y1*z2 - z1*y2
    yc = z1*x2 - x1*z2
    zc = x1*y2 - y1*x2
    sn = vecnorm((xc, yc, zc))
    # Convert to polar coordinates.
    radius, angle = recpol(cs, sn)
    return degrees ? rad2deg(angle) : angle
end

"""
    sphdist(long1, lat1, long2, lat2[, degrees=true]) -> angular_distance

### Purpose ###

Angular distance between points on a sphere.

### Arguments ###

* `long1`:  longitude of point 1
* `lat1`: latitude of point 1
* `long2`: longitude of point 2
* `lat2`: latitude of point 2
* `degrees` (optional boolean keyword): if `true`, all angles, including the
  output distance, are assumed to be in degrees, otherwise they are all in
  radians.  It defaults to `false`.

### Output ###

Angular distance on a sphere between points 1 and 2, as an `AbstractFloat`.  It
is expressed in radians unless `degrees` keyword is set to `true`.

### Example ###

``` julia
sphdist(120, -43, 175, +22)
# => 1.5904422616007134
```

### Notes ###

* `gcirc` function is similar to `sphdist`, but may be more suitable for
  astronomical applications.
* If `long1`, `lat1` are scalars, and `long2`, `lat2` are vectors, then the
 output is a vector giving the distance of each element of `long2`, `lat2` to
 `long1`, `lat1`.  Similarly, if `long1`,`de1` are vectors, and `long2`,` lat2`
 are scalars, then the output is a vector giving the distance of each element of
 `long1`, `lat1` to `long2`, `lat2`.  If both `long1`, `lat1` and `long2`,
 `lat2` are vectors then the output is a vector giving the distance of each
 element of `long1`, `lat1` to the corresponding element of `long2`, `lat2`.

Code of this function is based on IDL Astronomy User's Library.
"""
sphdist(long1::Real, lat1::Real,
        long2::Real, lat2::Real; degrees::Bool=false) =
            _sphdist(promote(float(long1), float(lat1),
                             float(long2), float(lat2))... , degrees)

function sphdist{LO1<:Real, LA1<:Real}(long1::AbstractArray{LO1},
                                       lat1::AbstractArray{LA1},
                                       long2::Real,
                                       lat2::Real;
                                       degrees::Bool=false)
    @assert length(long1) == length(lat1)
    dist = similar(long1, typeof(float(one(LO1))))
    for i in eachindex(long1)
        dist[i] = sphdist(long1[i], lat1[i], long2, lat2, degrees=degrees)
    end
    return dist
end

function sphdist{LO2<:Real, LA2<:Real}(long1::Real,
                                       lat1::Real,
                                       long2::AbstractArray{LO2},
                                       lat2::AbstractArray{LA2};
                                       degrees::Bool=false)
    @assert length(long2) == length(lat2)
    dist = similar(long2, typeof(float(one(LO2))))
    for i in eachindex(long2)
        dist[i] = sphdist(long1, lat1, long2[i], lat2[i], degrees=degrees)
    end
    return dist
end

function sphdist{LO1<:Real, LA1<:Real, LO2<:Real, LA2<:Real}(long1::AbstractArray{LO1},
                                                             lat1::AbstractArray{LA1},
                                                             long2::AbstractArray{LO2},
                                                             lat2::AbstractArray{LA2};
                                                             degrees::Bool=false)
    @assert length(long1) == length(lat1) == length(long2) == length(lat2)
    dist = similar(long1, typeof(float(one(LO1))))
    for i in eachindex(long1)
        dist[i] = sphdist(long1[i], lat1[i], long2[i], lat2[i], degrees=degrees)
    end
    return dist
end


text:
module CUDAatomics

using CUDAnative, LLVM
using LLVM.Interop

let type_map = Dict( Int32 => ("u32", "r", "b32", "s32"), UInt32 => ("u32", "r", "b32", "u32"), Int64 => ("u64", "l", "b64", "s64"), UInt64 => ("u64", "l", "b64", "s64"), Float32 => ("f32", "f", "b32", "f32"), Float64 => ("f64", "d", "b64", "f64") ), atomiclist = [("add",2,1), ("exch",2,3), ("min",2,4), ("max",2,4), ("inc",2,1), ("dec",2,1), ("cas",3,3), ("and",2,3), ("or",2,3), ("xor",2,3)]

global @generated function cvt( a::Type{T}, b::U ) where {T,U}
    round_str = ""
    if U==Float32
        if !(T in [Float32, Float64])
            round_str = "rzi."
        end
    elseif U==Float64
        if T==Float32
            round_str = "rn."
        elseif T!=Float64
            round_str = "rzi."
        end
    else
        if T in [Float32,Float64]
            round_str = "rn."
        end
    end

    call_string = string("cvt.",round_str,type_map[T][1],".",type_map[U][1], " \$0, \$1;")
    ex = :(@asmcall)
    append!(ex.args, [call_string, string("=",type_map[T][2],",",type_map[U][2]), false, Symbol(T), :(Tuple{$U}), :b])
    return :(Base.@_inline_meta; $ex)
end

function atomicexpression(instr::String, nargs, typeindex)
    fex = :(@generated function $(Symbol("atomic"*instr))(a) end)
    for i=1:(nargs-1)
        push!(fex.args[3].args[1].args, Symbol("a"*"$i"))
    end
    push!(fex.args[3].args[1].args, Expr(:kw, :index, 1))
    push!(fex.args[3].args[1].args, Expr(:kw, :field, Val(nothing)))
    fargs = fex.args[3].args[2].args
    append!(fargs, (quote
        type_map = $type_map
        fieldsym = field.parameters[1]
        if fieldsym == nothing
            base_type = a.parameters[1]
            offset = 0
        else
            field_index = findfirst(fieldnames(a.parameters[1]) .== fieldsym)
            base_type = a.parameters[1].types[field_index]
            offset = (cumsum(sizeof.(a.parameters[1].types)) .- sizeof.(a.parameters[1].types))[field_index]
        end
        ASstr = a.parameters[3] == CUDAnative.AS.Shared ? "shared" : "global"
    end).args)

    append!(fargs, (quote
        call_string = string( "cvta.to.",ASstr,".u64 %rd1, \$1;\natom.",ASstr,".",$instr,".", type_map[base_type][$typeindex], " \$0, [%rd1],")
        call_string = string(call_string, string(string(string.(" \$", collect(2:$nargs), ",")...)[1:end-1],";"))
    end).args)

    for i=1:(nargs-1)
        varsym = Symbol("a"*"$i")
        subex = :($(Symbol("a"*"$i"*"val")) = $varsym==base_type ? $(QuoteNode(:($varsym))) : :(cvt($base_type, $varsym)))
        subex.args[end].args[end].args[end].args[end] = varsym
        push!(fargs, subex)
    end

    append!(fargs, (quote
        ex = :(@asmcall)
        constraint_list = Array{String}(["=", type_map[base_type][2], ",l"])
        argtype_expr = :(Tuple{UInt64})
        for i=2:$nargs
            push!(constraint_list,string(","*type_map[base_type][2]))
            push!(argtype_expr.args,Symbol(base_type))
        end
        append!(ex.args, [call_string, string(constraint_list...), true, Symbol(base_type), argtype_expr, :( UInt64(pointer(a)) + $offset + $(sizeof(base_type))*(index-1)), a1val])
        return :(Base.@_inline_meta; $ex)
    end).args)

    for i=2:(nargs-1)
        push!(fargs[end-2].args[end].args, Symbol("a"*"$i"*"val"))
    end

    return fex
end

for atomicfunc in atomiclist
    eval(atomicexpression(atomicfunc[1], atomicfunc[2], atomicfunc[3]))
    eval(:(export $(Symbol("atomic"*atomicfunc[1]))))
end

end

function atomicsub(a::CuDeviceArray{Int32, N, A}, b, index=1, field=Val(nothing)) where {N,A}
    atomicadd(a, -b, index, field)
end
export atomicsub

export cvt

end


text:
### State type

# just an ADM scene
const OCState = Records.Frame{Records.Entity{AutomotiveDrivingModels.VehicleState,AutomotiveDrivingModels.VehicleDef,Int64}}

#### Observaton type

# # for LidarFeatures
# const OCObs = LidarSensor

# for vector features
const OCObs = Array{Float64, 1}
#### Action type

struct OCAction
    acc::Float64
end

function Base.copyto!(a::OCAction, b::OCAction)
    a.acc = b.acc
end

function Base.hash(a::OCAction, h::UInt64 = zero(UInt64))
    return hash(a.acc, h)
end

function Base.:(==)(a::OCAction, b::OCAction)
    return a.acc == b.acc
end

#### POMDP type

@with_kw mutable struct OCPOMDP <: POMDP{OCState, OCAction, OCObs}
    env::CrosswalkEnv = CrosswalkEnv(CrosswalkParams(obstacles_visible=true))
    sensor::AbstractSensor = GaussianSensor()
    models::Dict{Int64, DriverModel} = Dict{Int64, DriverModel}(1=>EgoDriver(OCAction(0.)))
    ego_type::VehicleDef =  VehicleDef()
    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)
    max_acc::Float64 = 2.0
    ego_start::Float64 = 5.0
    ped_start::Float64 = -env.params.crosswalk_length/4
    ego_goal::Float64 = env.params.roadway_length/2 + 2*env.params.crosswalk_width
    ped_goal::Float64 = env.params.crosswalk_length/4
    ΔT::Float64 = 0.5 # decision frequency
    p_birth::Float64 = 0.3 # probability of appearance
    max_peds::Int64 = 10
    no_ped_prob::Float64 = 0.3 # probability that there is no pedestrian (for training)
    no_ped::Bool = false
    v_noise::Float64 = 1.0
    pos_obs_noise::Float64 = 0.5
    vel_obs_noise::Float64 = 0.5
    collision_cost::Float64 = -1.0
    action_cost::Float64 = 0.0
    goal_reward::Float64 = 1.0
    γ::Float64 = 0.95 # discount factor
end

# function OCPOMDP(; env::CrosswalkEnv = CrosswalkEnv(),
#                    sensor::LidarSensor = LidarSensor(50, max_range=30., angle_spread=float(pi)),
#                    models::Dict{Int64, DriverModel} = Dict{Int64, DriverModel}(1=>EgoDriver(OCAction(0.))),
#                    ego_type::VehicleDef = VehicleDef(),
#                    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0),#,env.params.crosswalk_width),
#                    max_acc::Float64 = 2.0,
#                    ego_start::Float64 = 5.,
#                    ped_start::Float64 = -env.params.crosswalk_length/4,
#                    ego_goal::Float64 = env.params.roadway_length/2 + 2*env.params.crosswalk_width,
#                    ped_goal::Float64 = env.params.crosswalk_length/4,
#                    ΔT::Float64  = 0.5,
#                    p_birth::Float64 = 0.3,
#                    max_peds::Int64=10,
#                    no_ped_prob::Float64 = 0.3,
#                    no_ped::Bool = false,
#                    v_noise::Float64 = 1.,
#                    pos_obs_noise::Float64 = 0.5,
#                    vel_obs_noise::Float64 = 0.5,
#                    collision_cost::Float64 = -1.,
#                    action_cost::Float64 = 0.0,
#                    goal_reward::Float64 = 1.,
#                    γ::Float64  = 0.95)
#     return OCPOMDP(env,
#                    sensor,
#                    models,
#                    ego_type,
#                    ped_type,
#                    max_acc,
#                    ego_start,
#                    ped_start,
#                    ego_goal,
#                    ped_goal,
#                    ΔT,
#                    p_birth,
#                    max_peds,
#                    no_ped_prob,
#                    no_ped,
#                    v_noise,
#                    pos_obs_noise,
#                    vel_obs_noise,
#                    collision_cost,
#                    action_cost,
#                    goal_reward,
#                    γ)
# end

### HELPERS

function POMDPs.discount(pomdp::OCPOMDP)
    return pomdp.γ
end

POMDPs.actions(pomdp::OCPOMDP) = [OCAction(-4.0), OCAction(-2.0), OCAction(0.0), OCAction(2.0)]

function POMDPs.actionindex(pomdp::OCPOMDP, action::OCAction)
    if action.acc == -4.0
        return 1
    elseif action.acc == -2.0
        return 2
    elseif action.acc == 0.
        return 3
    else
        return 4
    end
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
@testset "Tide files (vav & eterna)" begin
	file_vav = joinpath(dirname(@__DIR__),"test","input","vav_analysis.dat");
	file_tsoft = joinpath(dirname(@__DIR__),"test","output","vav2tsoft.txt");
	isfile(file_tsoft) ? rm(file_tsoft) : nothing
	vav2tsoft(file_vav,file_tsoft,site="Cantlay",name="test");
	@test isfile(file_tsoft)
	file_eterna = joinpath(dirname(@__DIR__),"test","input","eterna_analysis.prn");
	file_tsoft = joinpath(dirname(@__DIR__),"test","output","eterna2tsoft.txt")
	isfile(file_tsoft) ? rm(file_tsoft) : nothing
	eterna2tsoft(file_eterna,file_tsoft,site="Cantlay",name="test");
	@test isfile(file_tsoft)
end


text:
using PyCall, PyPlot

@testset "plots" begin
    data = load_arviz_data("centered_eight")
    data2 = load_arviz_data("non_centered_eight")

    rng = MersenneTwister(42)
    arr1 = randn(rng, 4, 100)
    arr2 = randn(rng, 4, 100)
    arr3 = randn(rng, 100)

    @testset "$(f)" for f in (plot_trace, plot_pair, plot_joint)
        f(data; var_names = ["tau", "mu"])
        close(gcf())
        f((x = arr1, y = arr2); var_names = ["x", "y"])
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test f((x = arr1, y = arr2); var_names = ["x", "y"], backend = :bokeh) isa
                  ArviZ.BokehPlot
        end
    end

    @testset "$(f)" for f in
                        (plot_autocorr, plot_ess, plot_mcse, plot_posterior, plot_violin)
        f(data; var_names = ["tau", "mu"])
        close(gcf())
        f(arr1)
        close(gcf())
        f((x = arr1, y = arr2); var_names = ["x", "y"])
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test f((x = arr1, y = arr2); var_names = ["x", "y"], backend = :bokeh) isa
                  ArviZ.BokehPlot
        end
    end

    @testset "$(f)" for f in (plot_energy, plot_parallel)
        f(data)
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test f(data; backend = :bokeh) isa ArviZ.BokehPlot
        end
    end

    @testset "$(f)" for f in (plot_density, plot_forest)
        f(data; var_names = ["tau", "mu"])
        close(gcf())
        f([(x = arr1,), (x = arr2,)], var_names = ["x"])
        close(gcf())
        f(arr3)
        close(gcf())
        f((x = arr1, y = arr2); var_names = ["x", "y"])
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test f([(x = arr1,), (x = arr2,)]; var_names = ["x"], backend = :bokeh) isa
                  ArviZ.BokehPlot
        end
    end

    @testset "plot_rank" begin
        plot_rank(data; var_names = ["tau", "mu"])
        close(gcf())
        plot_rank(arr1)
        close(gcf())
        plot_rank((x = arr1, y = arr2); var_names = ["x", "y"])
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test isa(
                plot_rank((x = arr1, y = arr2); var_names = ["x", "y"], backend = :bokeh),
                ArviZ.BokehPlot,
            )
        end
    end

    @testset "plot_compare" begin
        df = compare(Dict("a" => data, "b" => data2))
        plot_compare(df)
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test plot_compare(df; backend = :bokeh) isa ArviZ.BokehPlot
        end
    end

    @testset "$(f)" for f in (plot_dist, plot_kde)
        f(arr1)
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test f(arr1; backend = :bokeh) isa ArviZ.BokehPlot
        end
    end

    @testset "$f" for f in (plot_hpd, plot_hdi)
        x_data = randn(rng, 100)
        y_data = 2 .+ x_data .* 0.5
        y_data_rep = 0.5 .* randn(rng, 200, 100) .+ transpose(y_data)
        f(x_data, y_data_rep)
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test f(x_data, y_data_rep; backend = :bokeh) isa ArviZ.BokehPlot
        end
    end

    @testset "plot_elpd" begin
        plot_elpd(Dict("a" => data, "b" => data2))
        close(gcf())
        plot_elpd(Dict(
            "a" => loo(data; pointwise = true),
            "b" => loo(data2; pointwise = true),
        ))
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test plot_elpd(Dict("a" => data, "b" => data2); backend = :bokeh) isa
                  ArviZ.BokehPlot
        end
    end

    @testset "plot_khat" begin
        l = loo(data; pointwise = true)
        plot_khat(l)
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test plot_khat(l; backend = :bokeh) isa ArviZ.BokehPlot
        end
    end

    @testset "plot_loo_pit" begin
        plot_loo_pit(data; y = "obs")
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test plot_loo_pit(data; y = "obs", backend = :bokeh) isa ArviZ.BokehPlot
        end
    end

    @testset "plot_loo_pit" begin
        plot_ppc(data)
        close(gcf())
        ispynull(ArviZ.bokeh) || @testset "bokeh" begin
            @test plot_ppc(data; backend = :bokeh) isa ArviZ.BokehPlot
        end
    end
end


text:
# Abstract Types
abstract type AbstractInstrument end
abstract type AbstractEchelleObservationOrSet end
abstract type AbstractEchelleObservation <: AbstractEchelleObservationOrSet end
abstract type AbstractEchelleObservationSet <: AbstractEchelleObservationOrSet end
#=
abstract type AbstractBlazeType end
abstract type AbstractBlaze1DType <: AbstractBlazeType end
abstract type AbstractBlaze2DType <: AbstractBlazeType end
=#
abstract type AbstractPSF end
abstract type AbstractSpectralOrderModel end
abstract type AbstractSpectralModel end
abstract type AbstractObservationModel end
abstract type AbstractObsSeriesModel end
abstract type AbstractManifestFormatEntry end
export  AbstractEchelleObservationOrSet, AbstractEchelleObservation, AbstractEchelleObservationSet
export AbstractBlazeType, AbstractBlaze1DType, AbstractBlaze2DType
export AbstractPSF
export AbstractSpectralOrderModel, AbstractSpectralModel
export AbstractObservationModel, AbstractObsSeriesModel
export AbstractManifestFormatEntry

# Type Aliases
#const AbstractEchelleObservationOrSet = Union{AbstractEchelleObservation, AbstractEchelleObservationSet}
const OrdersType = Union{Int64,UnitRange{Int64},AbstractVector{Int64}}
const Wavelength1DType{T1} = Union{AbstractArray{T1,1}, AbstractRange{T1} } where T1<:Real
const Wavelength2DType{T1} = Union{AbstractArray{T1,2}, AbstractArray{AbstractRange{T1},1} } where T1<:Real
export OrdersType, Wavelength1DType, Wavelength2DType


text:
module CategoricalConversions

export Continuous, HasRefPool, Categorical,
    catetorical_labels, categorical_range, categorical_trait, categorical_positions

using DataAPI: DataAPI
using AbstractPlotting: Automatic

# Trait for categorical values
struct Categorical end
struct Continuous end
struct HasRefPool end # the better Categorical

categorical_trait(x::AbstractVector) = !isnothing(DataAPI.refpool(x)) ? HasRefPool() : Categorical()
categorical_trait(::AbstractVector{<: Number}) = Continuous()

categorical_labels(xs) = categorical_labels(categorical_trait(xs), xs)
categorical_labels(::Categorical, xs) = unique(xs)
categorical_labels(::Continuous,  _)  = Automatic() # we let them be automatic
categorical_labels(::HasRefPool,  xs) = DataAPI.levels(xs) # could also use values(DataAPI.refpool(xs))

categorical_range(xs) = categorical_range(categorical_trait(xs), xs)
categorical_range(::Categorical, xs) = 1:length(categorical_labels(xs))
categorical_range(::Continuous,  _)  = Automatic() # we let them be automatic
categorical_range(::HasRefPool,  xs) = keys(DataAPI.refpool(xs))

categorical_position(x, xs) = categorical_position(categorical_trait(xs), x, xs)
categorical_position(::Categorical, x, xs, labels = categorical_labels(xs)) = findfirst(l -> l == x, labels)
categorical_position(::Continuous,  x, _)  = x
categorical_position(::HasRefPool,  x, xs) = DataAPI.invrefpool(xs)[x]

categorical_positions(xs) = categorical_positions(categorical_trait(xs), xs)
categorical_positions(::Continuous, xs) = xs
function categorical_positions(t::Categorical, xs)
    labels = categorical_labels(xs)
    categorical_position.(Ref(t), xs, Ref(xs), Ref(labels))
end
categorical_positions(::HasRefPool, xs) = DataAPI.refarray(xs)

end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
@testset "Bitstream" begin
    x = SBitstream(0.5)
    xbit = pop!(x)

    @test (push!(x, xbit); true)
    @test observe(x) == xbit
    @test length(x) == 1
    @test pop!(x) == xbit
end


text:
export load_model, get_device, seisdata2torch, slide, detect, trigger_onset
export get_picks, load_test_data

"""
  load_model(phase;device=nothing)

Load generalized-phase-detection model onto device.

# Arguments
- `phase::String`: Seismic phase to load - "P" or "S".

# Keywords
- `device::Union{PyObject,Nothing}`: Pytorch device - specify "gpu" or "cpu".
"""
function load_model(phase::String;device::Union{PyObject,Nothing}=nothing)
    # check phase input
    if uppercase(phase) ∉ ["P","S"]
        error("Only available models are 'P' and 'S'.")
    end

    # get correct model
    model_dir = joinpath(dirname(@__FILE__),"..","deps","downloads")
    modelP = joinpath(model_dir,"model_P.pt")
    modelS = joinpath(model_dir,"model_S.pt")
    model_file = uppercase(phase) == "P" ? modelP : modelS

    # load model
    if isa(device,Nothing)
        device = get_device()
    end
    model = UNet(3,1).to(device)
    checkpoint = torch.load(model_file, map_location=device)

    # load state dictionary
    state_dict = Dict()
    for key in keys(checkpoint["model_state_dict"])
        if occursin("tracked",key)
            continue
        end

        state_dict[key] = checkpoint["model_state_dict"][key]
    end
    model.load_state_dict(state_dict)
    model.eval()
    return model
end

"""
  get_device()

Returns current PyTorch device.
"""
get_device() = torch.device(torch.cuda.is_available() ? "cuda" : "cpu")

"""
  seisdata2torch(S)

Convert 3 component `SeisData` struct `S` into Torch tensor.

# Arguments
- `S::SeisData`: 3-component `SeisData` for detection.
- `window_samples::Int`: Number of samples per detection window.
"""
function seisdata2torch(S::SeisData,window_samples::Int)
    # check only 3 channels in SeisChannel
    @assert S.n == 3
    # create sliding windows
    X = cat([slide(S[ii].x,window_samples) for ii = 1:S.n]...,dims=3)
    X ./= maximum(abs.(X),dims=(1,3))
    return torch.from_numpy(X).float().permute(1,2,0)
end

"""
  slide(A,window_samples)

Split vector `A` into non-overlapping windows of length `window_samples`.

# Arguments
- `S::SeisData`: 3-component `SeisData` for detection.
- `window_samples::Int`: Number of samples per detection window.
"""
function slide(A::AbstractVector, window_samples::Int)
    N = size(A,1)
    if N % window_samples == 0
        return Array(reshape(A,window_samples,N ÷ window_samples))
    else
        return Array(reshape(A[1 : N - N % window_samples], window_samples, N ÷ window_samples))
    end
end

"""
  detect(X,model,batch_size;device=nothing)

Run generalized-phase-detection on tensor `X` in batches of size `batch_size`.

# Arguments
- `X::PyObject`: Torch tensor.
- `model::PyObject`: Torch `UNet` model.
- `batch_size::Int`: Number of windows per detection. Specifying too large a
    batch size will cause the GPU to run out of VRAM.

# Keywords
- `device::Union{PyObject,Nothing}`: Pytorch device - specify "gpu" or "cpu".
"""
function detect(X,model,batch_size::Int;device::Union{PyObject,Nothing}=nothing)
    # use default device
    if isa(device,Nothing)
        device = get_device()
    end

    # run predictions in batches
    Y_pred = zeros(Float32, X.size(2),X.size(0))
    for ii in 1:batch_size:size(Y_pred,2)
        i_start = ii
        i_stop = min(size(Y_pred,2),ii+batch_size-1)
        X_test = get(X,i_start-1:i_stop-1)
        X_test = X_test.to(device)
        @pywith torch.no_grad() begin
            out = model(X_test)
            Y_pred[:,i_start:i_stop] .= out.cpu().permute(2,0,1).numpy()[:,:,1]
        end
        GC.gc(false)
    end
    return Y_pred
end

"""
  trigger_onset(charfct, thresh1, thresh2; max_len=1e99, max_len_delete=false)

Calculate trigger on and off times.

Given `thresh1` and `thresh2` calculate trigger on and off times from
characteristic function.

# Arguments
- `charfct::AbstractArray`: Characteristic function of e.g. STA/LTA trigger.
- `thresh1::Real`: Value above which trigger (of characteristic function)
    is activated (higher threshold).
- `thresh2::Real`: Value below which trigger (of characteristic function)
    is deactivated (lower threshold).

# Keywords
- `max_len::Int`: Maximum length of triggered event in samples. A new
    event will be triggered as soon as the signal reaches again above thresh1.
- `max_len_delete::Bool`: Do not write events longer than max_len into report file.
"""
function trigger_onset(charfct::AbstractArray, thresh1::Real, thresh2::Real;
    max_len::Int=10^10, max_len_delete::Bool=false,
)
    ind1 = findall(charfct .> thresh1)
    if length(ind1) == 0
        return []
    end
    ind2 = findall(charfct .> thresh2)

    on = [ind1[1]]
    off = [-1]

    # determine the indices where charfct falls below off-threshold
    ind2_ = Array{Bool}(undef,length(ind2))
    ind2_[1:end-1] .= diff(ind2) .> 1
    # last occurence is missed by diff, add it manually
    ind2_[end] = true
    append!(off,ind2[ind2_])
    append!(on,ind1[findall(diff(ind1) .> 1) .+ 1])
    # include last pick if trigger is on or drop it
    if max_len_delete
        # drop it
        append!(off,max_len)
        append!(on,on[end])
    else
        # include it
        append!(off,ind2[end])
    end

    pick = []
    while on[end] > off[1]
        while on[1] <= off[1]
            deleteat!(on,1)
        end
        while off[1] < on[1]
            deleteat!(off,1)
        end
        if off[1] - on[1] > max_len
            if max_len_delete
                deleteat!(on,1)
                continue
            end
            prepend!(off,on[1] + max_len)
        end
        push!(pick,[on[1],off[1]])
    end
    return permutedims(hcat(pick...))
end

"""
  get_picks(S,probs,thresh1,thresh2;min_trig_dur=0.)

Get picks from detection probabilities `probs` from SeisData `S`.

# Arguments
- `S::SeisData`: 3-component `SeisData` for detection.
- `probs::AbstractArray`: Array of probabilities of phase detection per window.
- `thresh1::Real`: Probability (0-1) above which trigger is activated (higher threshold).
- `thresh2::Real`: Probability (0-1) below which trigger is deactivated (lower threshold).

# Keywords
`min_trig_dur::Real`: Minimum duration (in seconds) of trigger.

"""
function get_picks(
    S::SeisData,
    probs::AbstractArray,
    thresh1::Real,
    thresh2::Real;
    min_trig_dur::Real=0.,
)
    @assert 0. < thresh2 < thresh1 <= 1.
    tt = collect(0:length(probs)-1) ./ S.fs[1]
    picks = Array{DateTime}(undef,0)
    trigs = trigger_onset(vcat(probs...),thresh1,thresh2)
    for ii in 1:size(trigs,1)
        if trigs[ii,1] == trigs[ii,2]
            continue
        end

        # check trigger duration
        trig_dur = tt[trigs[ii,2]] - tt[trigs[ii,1]]
        if trig_dur < min_trig_dur
            continue
        end

        pick = argmax(probs[trigs[ii,1]:trigs[ii,2]]) + trigs[ii,1]

        tstamp = u2d(S.t[1][1,2] * 1e-6+ tt[pick])
        push!(picks,tstamp)
    end
    return picks, tt
end

"""
  load_test_data()

Load 3-component seismic data from the 2016 Mw 5.2 Anza, California sequence.
"""
function load_test_data()
    sacfiles = joinpath(dirname(@__FILE__),"..","deps","downloads","*.SAC")
    return read_data("sac",sacfiles)
end


text:
module CollabUtils2JournalTrigger

using ..Ahorn, Maple

@mapdef Trigger "CollabUtils2/JournalTrigger" JournalTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    levelset::String="Celeste")

const placements = Ahorn.PlacementDict(
    "Journal Trigger (Collab Utils 2 / READ DOCS)" => Ahorn.EntityPlacement(
        JournalTrigger,
        "rectangle",
    ),
)

function Ahorn.editingOptions(trigger::JournalTrigger)
    return Dict{String, Any}(
        "positionMode" => Maple.trigger_position_modes
    )
end

function Ahorn.nodeLimits(trigger::JournalTrigger)
    return 0, 1
end

end


text:
filename() = "day18.txt"

using MacroTools

function read_data()
    return readlines(filename())
end

#=
Replace + operators with / to equalize precedence with *.
Then, leverage the Julia parser. Finally, restore the +
operators using MacroTools.postwallk.
=#
function part1()
    function evaluate(line)
        revised = replace(line, "+" => "/")
        expr = Meta.parse(revised)
        expr = MacroTools.postwalk(x -> x == :(/) ? :(+) : x, expr)
        return eval(expr)
    end
    return sum(evaluate.(read_data()))
end

#=
Same trick as above to implement the reversed precedence rules.
=#
function part2()
    function evaluate(line)
        revised = replace(replace(line, "*" => "-"), "+" => "/")
        expr = Meta.parse(revised)
        expr = MacroTools.postwalk(x -> ifelse(x == :(/), :(+), ifelse(x == :(-), :(*), x)), expr)
        return eval(expr)
    end
    return sum(evaluate.(read_data()))
end


text:
using Pkg; Pkg.activate("..")
using Turing
using Turkie
using GLMakie # You could also use CairoMakie or another backend
using CairoMakie
GLMakie.activate!()
CairoMakie.activate!()
Turing.@model function demo(x) # Some random Turing model
    m0 ~ Normal(0, 2)
    s ~ InverseGamma(2, 3)
    m ~ Normal(m0, √s)
    for i in eachindex(x)
        x[i] ~ Normal(m, √s)
    end
end

xs = randn(100) .+ 1;
m = demo(xs);
cb = TurkieCallback(m) # Create a callback function to be given to the sample function
chain = sample(m, NUTS(0.65), 30; callback = cb)

record(cb.figure, joinpath(@__DIR__, "video.gif")) do io
    addIO!(cb, io)
    sample(m,  NUTS(0.65), 50; callback = cb)
end

## Let's test if Soss work as well!
using Soss
using Random

sossdemo = Soss.@model x begin
    m0 ~ Normal(0, 2)
    s ~ InverseGamma(2, 3)
    m ~ Normal(m0, √s)
    x ~ For(eachindex(x)) do i
        Normal(m, √s)
    end
end

cb = TurkieCallback(Dict(:m0 => [:trace, :mean],
                        :s => [:autocov, :var]))

advancedHMC(sossdemo(), (x=xs,), 100; callback = cb)

## Test for array of parameters
using LinearAlgebra
D = 1
N = 20
Turing.@model function vectordemo(x, y, σ)
    m ~ Normal(0, 10)
    β ~ MvNormal(m * ones(D + 1), ones(D + 1))
    for i in eachindex(y)
        y[i] ~ Normal(dot(β, vcat(1, x[i])), σ)
    end
end

x = [rand(D) for _ in 1:N]
β = randn(D + 1) * 2
σ = 0.1
y = dot.(Ref(β), vcat.(1, x)) .+ σ * randn(N)

m = vectordemo(x, y, σ)

cb = TurkieCallback(m) # Create a callback function to be given to the sample function

chain = sample(m, NUTS(0.65), 200; callback = cb)


text:
using Test, YaoBlocks, LinearAlgebra

@testset "block map" begin
    @test BlockMap(ComplexF64, X) isa BlockMap{ComplexF64, typeof(X)}

    @test_throws ErrorException BlockMap(X) * rand(2)
    @test content(BlockMap(X)) === X

    st = rand(ComplexF64, 2)
    @test BlockMap(X) * st ≈ mat(X) * st
    @test issymmetric(BlockMap(X))
    @test size(BlockMap(X)) == (2, 2)
end


text:
using Test
using ReadVTK

# Commit in the example file repository for which the test files will be downloaded
# Note: The purpose of using a specific commit hash (instead of `main`) is to be able to tie a given
#       version of ReadVTK to a specific version of the test file repository. This way, also tests
#       for older ReadVTK releases should continue to work.
TEST_EXAMPLES_COMMIT = "92b44ef4666cae5aa5ffe1c5c35b5823c0073c31"

# Local folder to store downloaded example files. If you change this, also adapt `../.gitignore`!
TEST_EXAMPLES_DIR = "examples"

get_test_example_file(filename) = get_example_file(filename, head=TEST_EXAMPLES_COMMIT,
                                                   output_directory=TEST_EXAMPLES_DIR)

# Start with a clean environment: remove example file directory if it exists
isdir(TEST_EXAMPLES_DIR) && rm(TEST_EXAMPLES_DIR, recursive=true)
mkpath(TEST_EXAMPLES_DIR)

@time @testset "ReadVTK" begin
  @testset "basic tests" begin
    @testset "VTKFile" begin
      @test VTKFile(get_test_example_file("celldata_inline_binary_uncompressed.vtu")) isa VTKFile

      mktemp() do path, io
        write(io, "# vtk DataFile Version")
        flush(io)

        @test_throws ErrorException VTKFile(path)
      end
    end

    vtk_file = VTKFile(get_test_example_file("celldata_inline_binary_uncompressed.vtu"))

    @testset "get_cell_data" begin
      @test get_cell_data(vtk_file) isa ReadVTK.VTKData
    end

    cell_data = get_cell_data(vtk_file)

    @testset "VTKData auxiliary functions" begin
      @test firstindex(cell_data) == "cell_ids"
      @test lastindex(cell_data) == "indicator_shock_capturing"
      @test length(cell_data) == 5
      @test size(cell_data) == (5,)
      @test keys(cell_data) == ("cell_ids", "element_ids", "levels",
                                "indicator_amr", "indicator_shock_capturing")
      @test iterate(cell_data) == ("cell_ids" => cell_data["cell_ids"], 2)
      @test_throws KeyError cell_data["does_not_exist"]
      @test eltype(cell_data) == Pair{String, ReadVTK.VTKDataArray}
    end

    @testset "extract VTKDataArray from VTKData" begin
      @test cell_data["cell_ids"] isa ReadVTK.VTKDataArray
    end

    cell_ids = cell_data["cell_ids"]

    @testset "get_data" begin
      @test get_data(cell_ids) isa Base.ReinterpretArray{Int, 1}
    end

    data = get_data(cell_ids)

    @testset "validate data" begin
      @test length(data) == vtk_file.n_cells
      @test first(data) == 4
      @test last(data) == 4113
      @test sum(data) == 6357314
    end

    @testset "get_points" begin
      @test size(get_points(vtk_file)) == (3, 4434)
    end

    points = get_points(vtk_file)

    @testset "validate points" begin
      @test points[:, 1] ≈ [-64.0, -64.0, 0.0]
      @test sum(points) ≈ -6442.5
    end

    @testset "get_cells" begin
      @test get_cells(vtk_file) isa ReadVTK.VTKCells
    end

    cells = get_cells(vtk_file)

    @testset "validate cells" begin
      @test size(cells) == (3085,)
      @test div(Int(sum(cells.types)), 8) == 3085
      @test sum(cells.offsets) == 19040620
      @test cells.connectivity[1000] == 421
    end

    @testset "show" begin
      @test isnothing(show(devnull, vtk_file))
      @test isnothing(show(devnull, cell_data))
      @test isnothing(show(devnull, cell_ids))
    end
  end

  @testset "binary compressed file" begin
    vtk_file = VTKFile(get_test_example_file("celldata_inline_binary_compressed.vtu"))
    cell_data = get_cell_data(vtk_file)
    cell_ids = cell_data["cell_ids"]
    data = get_data(cell_ids)

    @testset "validate data" begin
      @test length(data) == vtk_file.n_cells
      @test first(data) == 4
      @test last(data) == 4113
      @test sum(data) == 6357314
    end
  end

  @testset "appended uncompressed file" begin
    vtk_file = VTKFile(get_test_example_file("celldata_appended_binary_uncompressed.vtu"))
    cell_data = get_cell_data(vtk_file)
    cell_ids = cell_data["cell_ids"]
    data = get_data(cell_ids)

    @testset "validate data" begin
      @test length(data) == vtk_file.n_cells
      @test first(data) == 4
      @test last(data) == 4113
      @test sum(data) == 6357314
    end
  end

  @testset "appended compressed file" begin
    vtk_file = VTKFile(get_test_example_file("celldata_appended_binary_compressed.vtu"))
    cell_data = get_cell_data(vtk_file)
    cell_ids = cell_data["cell_ids"]
    data = get_data(cell_ids)

    @testset "validate data" begin
      @test length(data) == vtk_file.n_cells
      @test first(data) == 4
      @test last(data) == 4113
      @test sum(data) == 6357314
    end
  end

  @testset "point data" begin
    vtk_file = VTKFile(get_test_example_file("pointdata_appended_binary_compressed.vtu"))

    @testset "get_point_data" begin
      @test get_point_data(vtk_file) isa ReadVTK.VTKData
    end

    point_data = get_point_data(vtk_file)
    pressure = point_data["p"]
    data = get_data(pressure)

    @testset "validate data" begin
      @test length(data) == vtk_file.n_points
      @test first(data) ≈ 0.7999332810225936
      @test last(data) ≈ 0.8004962389182811
      @test sum(data) ≈ 192.1204941112099
    end
  end
end

# Clean up afterwards: delete example file directory
@test_nowarn rm(TEST_EXAMPLES_DIR, recursive=true)


text:
import Base: eltype

# Gregorian calendar-ic globals
const second = 1.0
const minute = 60second
const hour = 60minute
const day = 24hour
const year = 365day
const stellaryear = 23hour + 56minute + 4.098903691
const Ω = 2π/stellaryear

function pressenter()
    println("\nPress enter to continue.")
    chomp(readline())
end

macro zeros(T, dims, vars...)
    expr = Expr(:block)
    append!(expr.args, [:($(esc(var)) = zeros($(esc(T)), $(esc(dims))); ) for var in vars])
    expr
end

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

function build_homogeneous_solution(names, typeprefix=:CellField; prefix=Symbol(""))
    nfields = length(names)
    fieldtype = Expr(:curly, typeprefix, :A, :G, :T)
    exprs = [ :( $(names[i]) :: $fieldtype ) for i = 1:nfields ]

    fullname = Symbol(prefix, :Solution)
    signature = Expr(:curly, fullname, :A, :G, :T)

    return quote
        #import StaticArrays: FieldVector

        struct $signature <: AbstractSolution{$(nfields), $fieldtype}
            $(exprs...)
        end
    end
end

function build_accessories(names; prefix=Symbol(""))
    nfields = length(names)
    bcsfields = [ :( $(names[i]) :: FieldBoundaryConditions ) for i = 1:nfields ]
    accfields = [ :( $(names[i]) :: T ) for i = 1:nfields ]
    lhsfields = [ :( $(names[i]) :: Tridiagonal{T, A} ) for i = 1:nfields ]

    bcs_signature = Symbol(prefix, :BoundaryConditions)
    acc_signature = Expr(:curly, Symbol(prefix, :Accessory), :T)
    lhs_signature = Expr(:curly, Symbol(prefix, :LeftHandSide), :T, :A)

    return quote
        import StaticArrays: FieldVector
        import LinearAlgebra: Tridiagonal
        import OceanTurb: build_lhs

        struct $bcs_signature <: FieldVector{$(nfields), FieldBoundaryConditions}
            $(bcsfields...)
        end

        struct $acc_signature <: FieldVector{$(nfields), T}
            $(accfields...)
        end

        struct $lhs_signature <: FieldVector{$(nfields), Tridiagonal{T, A}}
            $(lhsfields...)
            function LeftHandSide(solution::AbstractSolution{1, AbstractField})
                lhs = build_lhs(solution)[1]
                new{eltype(solution[1]), arraytype(solution[1])}(lhs)
            end
        end

        function LeftHandSide(solution::AbstractSolution{N, AbstractField}) where N
            lhs = build_lhs(solution)
            LeftHandSide{eltype(solution[1]), arraytype(solution[1])}(lhs...)
        end
    end
end

macro solution(names...)
    esc(quote
        $(build_homogeneous_solution(names))
        $(build_accessories(names))
    end
    )
end

macro prefixed_solution(prefix, names...)
    esc(quote
        $(build_homogeneous_solution(names; prefix=prefix))
        $(build_accessories(names; prefix=prefix))
    end
    )
end

@def add_standard_model_fields begin
    clock       :: Clock{T}
    grid        :: G
    timestepper :: TS
    solution    :: Solution
    bcs         :: BoundaryConditions
end

@def add_clock_grid_timestepper begin
    clock       :: Clock{T}
    grid        :: G
    timestepper :: TS
end

function run_until!(model, dt, tfinal)
    nt = floor(Int, (tfinal - time(model))/dt)
    iterate!(model, dt, nt)

    last_dt = tfinal - time(model)
    iterate!(model, last_dt)

    return nothing
end

"""
    diffusive_flux!(flux, fieldname, model)

Calculates the diffusive flux associated with `fieldname` in `model.solution.fieldname`i
and stores the result in the `FaceField` `flux`.
"""
function diffusive_flux!(flux, fieldname, model)
    field = getproperty(model.solution, fieldname)
    K = getproperty(model.timestepper.eqn.K, fieldname)

    for i in interiorindices(flux)
        @inbounds flux[i] = - K(model, i) * ∂z(field, i)
    end

    return nothing
end

"""
    diffusive_flux(fieldname, model)

Returns `flux::FaceField` with the diffusive flux associated with `fieldname`
in `model.solution.fieldname`.
"""
function diffusive_flux(fieldname, model)
    flux = FaceField(model.grid)
    diffusive_flux!(flux, fieldname, model)
    return flux
end


text:
"""
    prep_data_psi(data, seglen, verbose)

Checks and prepares data shape.

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels.
  - `seglen::Integer`: segment length (determines the frequency resolution).
  - `verbose::Bool`: if `true`, warnings and info logs would be echoed.

### Returns

  - `data::AbstractArray`: NxM array for N data points in M channels.
"""
function prep_data_psi(data::AbstractArray, seglen::Integer, verbose::Bool)

    if ndims(data) != 2  # data dimension
        data = squeeze(data)
        ndims(data) != 2 && throw(DimensionMismatch("Data must be a 2D-array!"))
        verbose && @info "Data is squeezed to a 2D-array)"
    end
    if size(data, 1) < size(data, 2)  # should be NxM array for N data points in M channels
        verbose && @info "Data is transposed to (#samples, #channels)"
        data = reshape(data, size(data, 2), size(data, 1))
    end
    if size(data, 1) < seglen  # seglen must be smaller than number of samples
        throw(DimensionMismatch("seglen must be smaller than number of samples!"))
    end
    return data
end

"""
    prep_freq(freqlist, seglen, verbose)

Checks and prep frequency list and bands

### Arguments

  - `freqlist::AbstractArray`: a UnitRange or 2D-Array where each column is a frequency band
  - `seglen::Integer`: segment length (determines the frequency resolution).
  - `verbose::Bool`: if `true`, warnings and info logs would be echoed.

"""
function prep_freq(freqlist::AbstractArray{Int}, seglen::Integer, verbose::Bool)
    # size(freqlist) = (freqs, nfbands)
    if length(freqlist) == 0
        freqlist = reshape(1:((int(seglen / 2))-1), :, 1)
    elseif ndims(freqlist) == 1
        freqlist = reshape(freqlist, :, 1)
    elseif ndims(freqlist) > 2
        throw("freqlist must be a UnitRange or a 2D-Array!")
    elseif size(freqlist, 1) < size(freqlist, 2)
        freqlist = freqlist'
        verbose && @info "freqlist is transposed to (#freq, #nfbands)"
    end
    if maximum(freqlist) >= int(seglen / 2)
        throw("Maximum frequency for freqlist is larger than the Nyquist frequency!")
    elseif minimum(freqlist) < 1
        throw("Minimum frequency for freqlist is 0 or negative!")
    end
    maxfreq = maximum(freqlist)  # max frequency of all frequency bands
    nfbands = size(freqlist, 2)  # number of frequency bands
    return freqlist, maxfreq, nfbands
end

"""
    make_eposeg(data, seglen, nep, nseg, nchan, segshift)

Partitioning data into epochs and segments

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels.
  - `seglen::Integer`: segment length
  - `nep::Integer`: number of epochs
  - `nseg::Integer`: number of segments per epoch
  - `nchan::Integer`: number of channels
  - `segshift::Integer`: number of bins by which neighboring segments are shifted.

### Returns

  - `epseg::AbstractArray`: partitioned data into shape `(seglen, nep, nseg, nchan)`

**Note**: returned Array may have more data entries than input data.
"""
function make_eposeg(
    data::AbstractArray,
    seglen::Integer,
    eplen::Integer,
    nep::Integer,
    nseg::Integer,
    nchan::Integer,
    segshift::Integer,
)::AbstractArray

    # preallocation
    epseg = Array{Float64}(undef, seglen, nep, nseg, nchan)

    for (i, e) in zip(1:nep, 1:eplen:(nep*eplen))
        for (j, s) in zip(1:nseg, 1:segshift:(nseg*seglen))
            @views epseg[:, i, j, :] = data[e:(e+eplen-1), :][s:(s+seglen-1), :]
        end
    end

    return epseg
end

"""
    cs2ps(cs)

Cross Spectra to Phase Slope as defined in Eq. 3 [Nolte et al. 2008]:

```math
\\tilde\\Psi_{ij}=\\mathfrak{I}\\left(\\sum_{f \\in F} C_{ij}^*(f)~C_{ij}(f + \\delta f)\\right)
```

where:

```math
C_{ij}(f) = \\frac{S_{ij}(f)}{\\sqrt{S_{ii}(f)~S_{jj}(f)}}
```

is the complex coherency, and \$S\$ is the cross spectral matrix
(returned by `data2cs` function), \$\\delta f\$ is the frequency resolution,
and \$\\mathfrak{I}\$ denotes taking the imaginary part.

### Arguments

  - `cs::AbstractArray`: Cross Spectral array with shape `(seglen, nchan, nchan)`

### Returns

  - phase slope index (::AbstractArray) with shape `(nchan, nchan)`

**Note**: frequency resolution is assumed to be the resolution of freq. band!
"""
function cs2ps(cs::AbstractArray)

    # size(cs) = (seglen, nchan, nchan)
    # we don't use df but we slice the fband

    # complex coherency
    @einsum coh[f, i, j] := cs[f, i, j] / sqrt(cs[f, i, i] * cs[f, j, j])

    # phase slope (Eq. 3)
    @views imag.(sum(conj(coh[1:(end-1), :, :]) .* coh[2:end, :, :]; dims = 1))
end

"""
    data2cs(data)

Segmented data to Cross Spectra based on Eq. 2 [Nolte et al. 2008]:

```math
S_{ij}(f) = \\langle \\hat{y}_i(f) \\hat{y}_i^*(f)\\rangle .
```

### Arguments

  - `data::AbstractArray`: Segmented data of shape `(maxfreq, nep, nseg, nchan)`

### Return

  - `cs::AbstractArray{Complex}`: Cross Spectral of shape `(maxfreq, nep, nseg, nchan, nchan)`
"""
function data2cs(data::AbstractArray)
    # cs: cross-spectral matrix

    # Eq. 2 size(csepseg) = (maxfreq, nep, nseg, nchan, nchan)
    @einsum cs[f, e, s, i, j] := data[f, e, s, i] * conj(data[f, e, s, j])
end

"""
    _cs_ = cs2cs_(data, cs, fband, nep, segave, subave, method)

preparing Cross Spectra for Phase Slope by segment averaging and subtraction

### Arguments

  - `data::AbstractArray`: Fourier-transformed detrended epoched segmented data.
  - `cs::AbstractArray`: Cross Spectra of data
  - `fband::AbstractArray`: 1D-Array of frequency range for PSI calculation
  - `nep::Integer`: number of epochs
  - `segave::Bool`: if `true`, averages across CS segments
  - `subave::Bool`: if `true`, subtract average across CS segments
  - `method::String`: standard deviation estimation method

### Returns

  - `_cs_::AbstractArray`: segment averaged and subtracted Cross Spectra
"""
function cs2cs_(
    data::AbstractArray,
    cs::AbstractArray,
    fband::AbstractArray,
    nep::Integer,
    segave::Bool,
    subave::Bool,
    method::String,
)
    if segave
        if method == "bootstrap"
            randboot = rand(1:nep, nep)
            cs_ = dropmean(view(cs, :, randboot, :, :, :), (2, 3))
            av_ = dropmean(view(data, fband, randboot, :, :), (2, 3))
        elseif method == "psi"
            cs_ = dropmean(cs, (2, 3))
            av_ = dropmean(view(data, fband, :, :, :), (2, 3))
        elseif method == "jackknife"
            cs_ = dropmean(cs, 3)
            av_ = dropmean(view(data, fband, :, :, :), 3)
        end
    else
        if method == "bootstrap"
            randboot = rand(1:nep, nep)
            cs_ = dropmean(view(cs, :, randboot, 1, :, :), 2)
            av_ = dropmean(view(data, fband, randboot, 1, :), 2)
        elseif method == "psi"
            cs_ = dropmean(view(cs, :, :, 1, :, :), 2)
            av_ = dropmean(view(data, fband, :, 1, :), 2)
        elseif method == "jackknife"
            cs_ = view(cs, :, :, 1, :, :)
            av_ = view(data, fband, :, 1, :)
        end
    end

    if subave
        if method == "bootstrap" || method == "psi"
            @einsum av_ij[f, i, j] := av_[f, i] * conj(av_[f, j])
        elseif method == "jackknife"
            @einsum av_ij[f, e, i, j] := av_[f, e, i] * conj(av_[f, e, j])
        end
        _cs_ = squeeze(cs_ - av_ij)
    else
        if method == "bootstrap" || method == "psi" || method == "jackknife"
            _cs_ = squeeze(cs_)
        end
    end

    return _cs_
end

"""
    psi_est(data, seglen ; segshift, eplen, freqlist, method,
            nboot, segave, subave, detrend)

calculates phase slope index (PSI)

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels
  - `seglen::Integer`: segment length (determines the frequency resolution).
    If defining `freqlist`, `seglen` must be the same as sampling frequency.

*optional arguments*

  - `segshift::Integer`: number of bins by which neighboring segments are shifted
    (default is `seglen / 2`)
  - `eplen::Integer`: length of epochs (default is number of samples)
  - `freqlist::AbstractArray`: a UnitRange or 2D-Array where each column is a frequency band
    (default is full range). Note that the DC component (0th frequency of FFT)
    is discarded, and `freqlist` shall only include integer values.
  - `method::String`: standard deviation estimation method (default is `"jackknife"`)
  - `subave::Bool`: if `true`, subtract average across Cross Spectra segments (default is `false`)
  - `segave::Bool`: if `true`, average across Cross Spectra segments (default is `true`)
  - `nboot::Integer`: number of bootstrap resampling iterations (default is `100`)
  - `detrend::Bool`: if `true`, performs a 0th-order detrend across raw segments (default is `false`)
  - `window::Function`: window function with interval length as sole necessary argument (default is `Hanning`)
  - `verbose::Bool`: if `true`, warnings and info logs would be echoed. (default is `false`)

### Returns

  - `psi::AbstractArray`: Phase Slope Index with shape `(channel, channel, frequency bands)`
  - `psi_std::AbstractArray`: PSI estimated standard deviation with shape `(channel, channel, frequency bands)`
"""
function psi_est(
    data::AbstractArray,
    seglen::Integer;
    segshift::Integer = 0,
    eplen::Integer = 0,
    freqlist::AbstractArray{Int} = Int[],
    method::String = "jackknife",
    subave::Bool = false,
    segave::Bool = true,
    nboot::Integer = 100,
    detrend::Bool = false,
    window::Function = hanning_fun,
    verbose::Bool = false,
)

    # check and reshape data if necessary and possible
    data = prep_data_psi(data, seglen, verbose)
    nsamples, nchan = size(data)  # number of samples per channel and number of channels

    # method shall always be lowercase
    method = lowercase(method)
    if eplen == 0
        verbose && @warn "Epoch length = 0 => No estimation of standard deviation."
        method = "none"
    end

    # if eplen = nsamples: continuous recording
    eplen == 0 && (eplen = nsamples)
    nep = int(nsamples / eplen)  # number of epochs

    if nep == 1 && subave == true
        subave = false
        verbose && @warn "subave is set to false for continuous data (nep = 1)"
    end

    # segment shift
    segshift == 0 && (segshift = int(seglen / 2))
    nseg = int((eplen - seglen) / segshift) + 1

    # frequency list and bands
    freqlist, maxfreq, nfbands = prep_freq(freqlist, seglen, verbose)

    # rolling window segment data
    eposeg = make_eposeg(data, seglen, eplen, nep, nseg, nchan, segshift)

    # detrending data if required
    if detrend
        detrend!(eposeg, 0)
    end

    # applying window function (Hanning window is default)
    eposeg .*= window(seglen)

    # Fast Fourier Transforming the data
    eposeg = view(fft(eposeg, 1), 2:(maxfreq+1), :, :, :)

    # preallocation of Phase-Slope-index arrays
    psi = Array{Float64}(undef, nchan, nchan, nfbands)
    if method == "jackknife"
        psi_est = Array{Float64}(undef, nchan, nchan, nfbands, nep)
    elseif method == "bootstrap"
        psi_est = Array{Float64}(undef, nchan, nchan, nfbands, nboot)
    end

    # calculating Cross-Spectra and PSI
    for (f, fband) in enumerate(eachrow(freqlist'))
        cs_full = data2cs(view(eposeg, fband, :, :, :))

        cs_psi = cs2cs_(eposeg, cs_full, fband, nep, segave, subave, "psi")
        psi[:, :, f] = cs2ps(cs_psi)

        if method == "jackknife"
            cs_jack = cs2cs_(eposeg, cs_full, fband, nep, segave, subave, "jackknife")
            for e = 1:nep
                cs_jack_se = (nep * cs_psi - view(cs_jack, :, e, :, :)) / (nep + 1)
                psi_est[:, :, f, e] = cs2ps(cs_jack_se)
            end
        elseif method == "bootstrap"
            for e = 1:nboot
                cs_boot = cs2cs_(eposeg, cs_full, fband, nep, segave, subave, "bootstrap")
                psi_est[:, :, f, e] = cs2ps(cs_boot)
            end
        end
    end

    # calculating the standard error
    if method == "jackknife"
        psi_std = sqrt(nep) * squeeze(std(psi_est; corrected = true, dims = 4))
    elseif method == "bootstrap"
        psi_std = squeeze(std(psi_est; corrected = true, dims = 4))
    else
        psi_std = fill(NaN, (nchan, nchan, nfbands))
    end

    psi = squeeze(psi)
    psi_std = squeeze(psi_std)

    return psi, psi_std
end


text:
function Base.convert(::Type{SisoZpk{T,TR}}, f::SisoRational{T2}) where {T<:Number, TR<:Number, T2<:Number}
    if length(f.num) == 0
        return SisoZpk(T[],TR[],0)
    elseif all(f.den == zero(f.den))
        error("Zero denominator, this should not be possible") # NOTE: Should we keep this?
    end

    return SisoZpk{T,TR}(roots(f.num), roots(f.den), f.num[end]/f.den[end]) # NOTE: polynomials are stored with increasing coefficient orders
end

function Base.convert(::Type{<:SisoZpk}, f::SisoRational{T}) where {T<:Number}
    TR = complex(Base.promote_op(/, T, T)) # Type of roots(f.z)
    Base.convert(SisoZpk{T,TR}, f)
end

function Base.convert(::Type{SisoZpk{T,TR}}, f::SisoZpk) where {T<:Number, TR<:Number}
    return SisoZpk{T,TR}(f.z, f.p, f.k)
end

function Base.convert(::Type{SisoRational{T}}, f::SisoRational) where {T<:Number}
    return SisoRational{T}(convert(Poly{T}, f.num), convert(Poly{T}, f.den))
end

function Base.convert(::Type{SisoRational{T1}}, f::SisoZpk{T2,TR}) where {T1<:Number, T2<:Number,TR<:Number}
    if T2 <: Real
        num = prod(roots2real_poly_factors(f.z))*f.k
        den = prod(roots2real_poly_factors(f.p))
    else
        num = prod(roots2poly_factors(f.z))*f.k
        den = prod(roots2poly_factors(f.p))
    end
    return SisoRational{T1}(num, den)
end

function Base.convert(::Type{SisoRational}, f::SisoZpk{T1,TR}) where {T1<:Number, TR<:Number}
    if T1 <: Real
        T = promote_type(T1,typeof(real(one(TR))))
        # Alternative version if T1 rich enough
        # T = T1
        convert(SisoRational{T}, f)
    else
        convert(SisoRational{promote_type(T1,TR)}, f)
    end
end

Base.convert(::Type{S}, b::T2) where {T1, S <: SisoTf{T1}, T2<:Number} = convert(T1, b)*one(S)


text:
# Sum of the largest k components

# export SumLargest

# TODO: SumLargest(r) is (the postcomposition of) the conjugate of
#  (1) ind{0 <= x <= 1 : sum(x) = r},
# where instead IndSimplex(r) corresponds to
#  (2) ind{0 <= x : sum(x) = r}.
# Therefore SumLargest(r) now is only correct when r = 1.
# To make SumLargest correct we should extend IndSimplex to allow for that
# additional bound in its definition, by adding a second argument to the
# constructor. Then in the following line we should replace IndSimplex(k)
# with IndSimplex(k, 1.0). Note that (1) is proper only if x ∈ Rⁿ for n ⩾ r.

"""
    SumLargest(k::Integer=1, λ::Real=1.0)

Returns the function `g(x) = λ⋅sum(x_[1], ..., x_[k])`, for an integer k ⩾ 1 and `λ ⩾ 0`.
"""
SumLargest(k::I=1, lambda::R=1.0) where {I <: Integer, R <: Real} = Postcompose(Conjugate(IndSimplex(k)), lambda)

function (f::Conjugate{IndSimplex{I}})(x::AbstractArray{S}) where {I <: Integer, S <: Real}
    if f.f.a == 1
        return maximum(x)
    end
    v = zero(S)
    if ndims(x) == 1
        p = partialsortperm(x, 1:f.f.a, rev=true)
        v = sum(x[p])
    else
        p = partialsortperm(x[:], 1:f.f.a, rev=true)
        v = sum(x[p])
    end
    return v
end

fun_name(f::Postcompose{Conjugate{IndSimplex{I}}, R}) where {I <: Integer, R <: Real} = "sum of k largest components"
fun_expr(f::Postcompose{Conjugate{IndSimplex{I}}, R}) where {I <: Integer, R <: Real} = "x ↦ λ⋅sum(x_[1], ..., x_[k])"
fun_params(f::Postcompose{Conjugate{IndSimplex{I}}, R}) where {I <: Integer, R <: Real} = "k = $(f.f.f.a), λ = $(f.a)"


text:
function activepowervariables(m::JuMP.AbstractModel, devices_netinjection::A, devices::Array{T,1}, time_periods::Int64) where {A <: JumpExpressionMatrix, T <: PowerSystems.HydroGen}

    on_set = [d.name for d in devices if d.available == true]

    t = 1:time_periods

    phy = @variable(m, phy[on_set,t]) # Power output of generators

    devices_netinjection = varnetinjectiterate!(devices_netinjection,  phy, t, devices)

    return phy, devices_netinjection

end


text:
using HierarchialPerformanceTest
using Test

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


text:
"""
    visualize_tree(q, n, t, root_state)

"""
function visualize_tree(q::Dict, n::Dict, t::Dict, root_state; title="MCTS Tree", init_expand=10, init_duration=0, kwargs...)
    S = typeof(root_state)
    total_n = Dict{S, Int}()
    for ((s, a), nn) in n
        total_n[s] = get(total_n, s, 0) + nn
    end

    maxq = maximum(values(q))
    minq = minimum(values(q))-0.00001

    children = [Int[]]
    text = ["$root_state\nN: $(total_n[root_state])"]
    style = [""]
    link_style = [""]
    tooltip = [""]

    # create sa nodes
    state_children = Dict{S, Vector{Int}}()
    sa_nodeindices = Dict{keytype(q), Int}()
    for (s, a) in keys(q)
        txt = "$a\nN: $(n[(s, a)])\nQ: $(@sprintf("%8.4g", q[(s,a)]))"
        push!(children, Int[])
        push!(text, txt)
        rel_q = (q[(s,a)]-minq)/(maxq-minq)
        color = weighted_color_mean(rel_q, colorant"green", colorant"red")
        push!(style, "stroke:#$(hex(color))")
        w = max(20.0*sqrt(n[(s, a)]/total_n[s]), 1)
        push!(link_style, "stroke-width:$(w)px")
        push!(tooltip, txt)
        sa_nodeindices[(s, a)] = length(children)
        if !haskey(state_children, s)
            state_children[s] = Int[]
        end
        push!(state_children[s], length(children))
    end

    children[1] = state_children[root_state]

    # create sasp nodes and push children to sa
    for (s, a, sp) in keys(t)
        txt = "$sp\nN: $(get(total_n, sp, 0))"
        push!(children, get(state_children, sp, Int[]))
        push!(text, txt)
        push!(style, "")
        w = 20.0*sqrt(t[(s,a,sp)]/n[(s, a)])
        push!(link_style, "stroke-width:$(w)px")
        push!(tooltip, txt)
        push!(children[sa_nodeindices[(s, a)]], length(children))
    end

    return D3Tree(children,
                  text=text,
                  style=style,
                  tooltip=tooltip,
                  link_style=link_style,
                  init_expand=init_expand,
                  init_duration=init_duration,
                  kwargs...
                 )
end


text:
import Pkg; Pkg.activate(".")
import CSV
using DataFrames
using StatsBase
using Statistics

##
# posterior samples for beta binomial model
betab_posterior = CSV.read(joinpath("data", "posterior_distributions", "beta_binomial_posterior.csv"))
# posterior samples from previous models
lssl_posterior = CSV.read(joinpath("data", "posterior_distributions", "lssl.csv"))
const_posterior = CSV.read(joinpath("data", "posterior_distributions", "const_posterior.csv"))
powerlaw_posterior = CSV.read(joinpath("data", "posterior_distributions", "powerlaw_posterior.csv"))

## BetaBinomial
betab_posterior[:mu] |> mean
betab_posterior[:mu] |> std

betab_posterior[:phi] .|> exp |> mean
betab_posterior[:phi] .|> exp |> std

## Power Law
powerlaw_posterior[:a] .|> exp |> mean
powerlaw_posterior[:a] .|> exp |> std

powerlaw_posterior[:b] |> mean
powerlaw_posterior[:b] |> std

powerlaw_posterior[:phi] .|> exp |> mean
powerlaw_posterior[:phi] .|> exp |> std

## Constant Connectance
const_posterior[:a] |> mean
const_posterior[:a] |> std

const_posterior[:phi] .|> exp |> mean
const_posterior[:phi] .|> exp |> std

## LSSL
lssl_posterior[:a] .|> exp |> mean
lssl_posterior[:a] .|> exp |> std

lssl_posterior[:phi] .|> exp |> mean
lssl_posterior[:phi] .|> exp |> std