text:
"""
Allows finding and enumerating all exact covers for a given collection of sets.

Given some set X and a collection S of subsets of X, an exact cover of X is
a selection S* of subsets from S such that each element of X is contained
in S* exactly once. [See this.](https://en.wikipedia.org/wiki/Exact_cover)

This module defines a type `CoverMatrix` that represents such a collection of
subsets along with methods to query it for exact covers. Instances can be
created from a `Dict{K, Set{V}}`, where the values `Set{V}` represent the
subsets and the universe X is implicitly defined by their union. The keys `K`
are used to identify the subsets in the exact covers.

For example, to find a simple exact cover for the sets {{1, 2}, {1, 2, 3}, {3}}:

    subsets = CoverMatrix(Dict(
      :a => Set([1, 2]),
      :b => Set([1, 2, 3]),
      :c => Set([3])
    ))

    find_exact_cover(subsets) # gives either `Set([:b])` or `Set([:a, :c])`

Of course there might be no exact cover for the given subsets, in which case
`find_exact_cover` will return `nothing`. If there are multiple exact covers,
it may return any of them. If you need all of them, you can call
`exact_cover_producer`, which returns a `Task` that produces each of them
individually. For instance:

    for cover in exact_cover_producer(subsets)
      println(cover)
    end

But note that the `CoverMatrix` is modified when enumerating covers; it cannot
be used otherwise while the enumeration is in progress and will only be
restored to its original state once the consumer finishes.

If you are just interested in the number of covers, you can use
`exact_cover_exists` or `exact_cover_count`.

The implemented algorithm for finding the exact covers is Knuth's
["Algorithm X" with Dancing Links (DLX)]
(http://www-cs-faculty.stanford.edu/~uno/papers/dancing-color.ps.gz).

"""
module ExactCovers

export
  CoverMatrix,
  exact_cover_exists,
  exact_cover_count,
  exact_cover_producer,
  find_exact_cover,
  cover_subset!,
  uncover_subset!

"""
The node type for the (sparse) subset matrix.

It has three specializations: the root (`CoverMatrix`), the column/element
headers (`ElementNode`) and the actual cells representing the entries in the
incidence matrix (`Cell`), where each has a different "payload".

All nodes are members of two cyclic lists, one horizontally and one vertically.
The `ElementNode`s and `CoverMatrix` act as list sentinels and are only linked
to themselves on creation.

There are some general functions on `Node`s to simplify code further down. These
are:
  * accessors for their neighbors,
  * functions to temporarily detach and reattach them from their containing
    lists, and
  * functions to initially insert them into the lists at specified positions.
"""
type Node{T}
  payload::T
  left::Node
  right::Node
  up::Node
  down::Node

  function Node(data::T)
    self = new(data)
    self.left = self.right = self.up = self.down = self
    self
  end
end

left(node::Node) = node.left
right(node::Node) = node.right
up(node::Node) = node.up
down(node::Node) = node.down

function detach_horizontally!(node::Node)
  node.left.right = node.right
  node.right.left = node.left
end

function reattach_horizontally!(node::Node)
  node.left.right = node
  node.right.left = node
end

function detach_vertically!(node::Node)
  node.up.down = node.down
  node.down.up = node.up
end

function reattach_vertically!(node::Node)
  node.up.down = node
  node.down.up = node
end

function append_down!(position::Node, node::Node)
  node.up = position
  node.down = position.down
  reattach_vertically!(node)

  node
end

function append_right!(position::Node, node::Node)
  node.left = position
  node.right = position.right
  reattach_horizontally!(node)

  node
end

type ElementHeader
  element
  count::Int
  isprimary::Bool
end

typealias ElementNode Node{ElementHeader}

type CellData
  element_node::ElementNode
  subset_key
end

typealias Cell Node{CellData}

type CoverMatrixHeader
  elements::Dict{Any, ElementNode}
  subsets::Dict{Any, Cell}
  visible_primary_element_count::Int

  CoverMatrixHeader() = new(Dict{Any, ElementNode}(), Dict{Any, Cell}(), 0)
end

typealias CoverMatrix Node{CoverMatrixHeader}

"Constructs a `CoverMatrix` from the given `Dict`. (See the module description)"
function CoverMatrix{K,V}(data::Dict{K, Set{V}}; optionals::Set{V} = Set{V}())
  result = CoverMatrix(CoverMatrixHeader())

  # Constructs the `CoverMatrix` row-by-row.
  for (key, subset) in data
    rest = copy(subset)
    # For each row (subset), the `ElementNode`s (colums) that are already
    # present in the matrix are iterated.
    for element_node in element_nodes(result)
      # If the corresponding element is contained in the current subset, a new
      # cell is inserted at the bottom of that column. If there are already
      # cells present in that row, the new cell must also be linked to them.
      # `append_cell!` takes care of this; it will always put the new cell at
      # the right end of the (cyclically linked) row list.
      element = element_node.payload.element
      if in(element, rest)
        delete!(rest, element)
        append_cell!(result, element_node, key)
      end
    end

    # For all subset elements that do not have columns in the matrix yet, we
    # need to create one (with `append_element!`). After that, it's the same as
    # above.
    for element in rest
      isprimary = !in(element, optionals)
      element_node = append_element!(result, element, isprimary)
      append_cell!(result, element_node, key)
    end
  end

  result
end

"""
Creates a new `ElementNode`, attaches it at the end of the element list, and
then returns it.

This function is a helper for the construction of the matrix. It must be called
at most once for each element.
"""
function append_element!(matrix::CoverMatrix, element, isprimary::Bool)
  matrix.payload.visible_primary_element_count += isprimary
  # Note: `matrix` is the sentinel for the (linked cyclic) element list, hence
  #   `matrix.left` is the last `ElementNode`.
  matrix.payload.elements[element] = append_right!(
    matrix.left,
    ElementNode(ElementHeader(element, 0, isprimary))
  )
end

"""
Inserts a new cell into the incidence matrix.

The cell is linked so that it is in `element_node`'s column and in the row
corresponding to the subset identified by `key`. This must be called in
left-to-right, top-down order.
"""
function append_cell!(matrix::CoverMatrix, element_node::ElementNode, key)
  cell = append_down!(element_node.up, Cell(CellData(element_node, key)))
  element_node.payload.count += 1
  if haskey(matrix.payload.subsets, key)
    first_cell = matrix.payload.subsets[key]
    append_right!(first_cell.left, cell)
  else
    matrix.payload.subsets[key] = cell
  end
end

Base.show(io::IO, matrix::CoverMatrix) =
  print(io, "CoverMatrix#$(object_id(matrix))")

Base.size(matrix::CoverMatrix) =
  (length(matrix.payload.subsets), length(matrix.payload.elements))

"""
Creates a `Task` that produces `Node`s from the (cyclic) linked lists.

It starts at some node `root` of the subset matrix and moves along neighboring
nodes until it arrives back at `root`, producing all nodes along the way. Which
direction is taken at each step is determined by `next`. `skip_first` specified
whether `root` should be produced as well.

This function is used as a primitive to simplify iteration.
"""
node_producer(root::Node, next::Function; skip_first = false) = @task let
  if !skip_first
    produce(root)
  end

  current = root |> next
  while current != root
    produce(current)
    current = current |> next
  end
end

element_nodes(node::CoverMatrix) = node_producer(node, right, skip_first = true)
cells(node::ElementNode) = node_producer(node, down, skip_first = true)
cells_reverse(node::ElementNode) = node_producer(node, up, skip_first = true)
cells_in_subset(node::Cell) = node_producer(node, right)
cells_in_subset_reverse(node::Cell) = node_producer(node.left, left)
other_cells_in_subset(node::Cell) =
  node_producer(node, right, skip_first = true)
other_cells_in_subset_reverse(node::Cell) =
  node_producer(node, left, skip_first = true)

"""
Enumerates exact covers and returns their count.

This function implements the DLX algorithm from the referenced paper, but is
slightly extended to allow aborting the search after the first hit as well as
counting the exact covers found.

If the actual covers are relevant (as opposed to simply the count),
`produce_traces` can be enabled. This will record the covered nodes in `stack`,
and cause the stack to be `produce`d when an exact cover is found. At this
point, `stack` contains all choices that led to the exact cover; it contains
the keys of all subsets covered at this point.
"""
function enumerate_exact_covers!(
  matrix::CoverMatrix;
  stack = Cell[],
  only_first = false,
  produce_traces = false
)
  if is_completely_covered(matrix)
    produce_traces && produce(stack)
    return 1
  end

  found = 0
  element_node = find_minimal_element_node(matrix)

  cover_element_node!(matrix, element_node)
  for element_cell in cells(element_node)
    produce_traces && push!(stack, element_cell)
    for cell in other_cells_in_subset(element_cell)
      cover_element_node!(matrix, cell.payload.element_node)
    end

    found += enumerate_exact_covers!(
      matrix,
      stack = stack,
      only_first = only_first,
      produce_traces = produce_traces
    )

    for cell in other_cells_in_subset_reverse(element_cell)
      uncover_element_node!(matrix, cell.payload.element_node)
    end
    produce_traces && pop!(stack)

    if only_first && found > 0
      break
    end
  end
  uncover_element_node!(matrix, element_node)

  found
end

"Determines whether at least one uncovered primary element exists."
is_completely_covered(matrix::CoverMatrix) =
  matrix.payload.visible_primary_element_count == 0

"""
Finds the uncovered primary element node that is contained in the fewest uncovered
subsets.

Precondition: At least one uncovered primary element must exist.
"""
function find_minimal_element_node(matrix::CoverMatrix)
  nodes = element_nodes(matrix)

  minimal_node = consume(nodes)
  for node in nodes
    if !minimal_node.payload.isprimary
      minimal_node = node
    elseif node.payload.isprimary && node.payload.count < minimal_node.payload.count
      minimal_node = node
    end
  end

  minimal_node
end

"Implements the cover operation from the referenced paper."
function cover_element_node!(matrix::CoverMatrix, element_node::ElementNode)
  detach_horizontally!(element_node)
  matrix.payload.visible_primary_element_count -= element_node.payload.isprimary
  for element_cell in cells(element_node)
    for cell in other_cells_in_subset(element_cell)
      detach_vertically!(cell)
      cell.payload.element_node.payload.count -= 1
    end
  end
end

"""
Looks up the `ElementNode` for `element` and then executes an cover operation
on it.
"""
function cover_element!(matrix::CoverMatrix, element)
  cover_element_node!(matrix, matrix.payload.elements[element])
end

"Implements the uncover operation from the referenced paper."
function uncover_element_node!(matrix::CoverMatrix, element_node::ElementNode)
  for element_cell in cells_reverse(element_node)
    for cell in other_cells_in_subset_reverse(element_cell)
      reattach_vertically!(cell)
      cell.payload.element_node.payload.count += 1
    end
  end
  reattach_horizontally!(element_node)
  matrix.payload.visible_primary_element_count += element_node.payload.isprimary
end

"""
Looks up the `ElementNode` for `element` and then executes an uncover operation
on it.
"""
function uncover_element!(matrix::CoverMatrix, element)
  uncover_element_node!(matrix, matrix.payload.elements[element])
end

"""
Similar to `cover_element!`, but covers all elements from the subset identified
by `key`.

This was not in the paper, and is also not used in the DLX algorithm. Rather,
it is public API that can be used to constrain the exact cover search so that
only covers which contain the given subset are considered.

This operation can be undone by `uncover_subset!`, which means that a single
`CoverMatrix` can be reused. TODO sudoku reference
"""
function cover_subset!(matrix::CoverMatrix, key)
  for cell = cells_in_subset(matrix.payload.subsets[key])
    detach_vertically!(cell)
    cover_element_node!(matrix, cell.payload.element_node)
  end
end

"""
Reverses `cover_subset!` for the subset identified by `key`.

This *must not* be called without a preceding call to `cover_subset!` on the
same `key`. `cover_subset!` calls *must* be undone *in reverse order*. Failure
to heed this advice will result in chaos and madness.
"""
function uncover_subset!(matrix::CoverMatrix, key)
  for cell = cells_in_subset_reverse(matrix.payload.subsets[key])
    uncover_element_node!(matrix, cell.payload.element_node)
    reattach_vertically!(cell)
  end
end

"Counts the possible exact covers of `matrix`."
exact_cover_count(matrix::CoverMatrix) = enumerate_exact_covers!(matrix)

"Checks whether an exact cover for `matrix` exists."
exact_cover_exists(matrix::CoverMatrix) =
  0 != enumerate_exact_covers!(matrix, only_first = true)

"""
Returns a `Task` that produces the possible exact covers of `matrix`.

The produced covers are `Set`s of keys that identify the actual subsets,
as defined when `matrix` was constructed. If `only_first` is set, this function
aborts after the first found exact cover.

During enumeration (i.e. when a value is consumed) the `matrix` is in a "dirty"
state; it should then not be used to execute another search or be modified
by `cover_subset!` or `uncover_subset!`.
"""
exact_cover_producer(matrix::CoverMatrix; only_first = false) = @task let
  covers = @task enumerate_exact_covers!(
    matrix,
    only_first = only_first,
    produce_traces = true
  )

  for cover in covers
    keys = map(cell -> cell.payload.subset_key, cover)
    set = Set(keys)
    produce(set)
  end
end

"""
Searches for the first exact cover of `matrix` and returns it.

If none is found, `nothing` is returned.
"""
function find_exact_cover(matrix::CoverMatrix)
  producer = exact_cover_producer(matrix, only_first = true)
  cover = consume(producer)
  consume(producer) # <- We need to do this because the CoverMatrix would not be
                    #    restored to its original state otherwise, as the
                    #    producing task would simply 'hang' after the cover was
                    #    found. This forces the producer to finish.
  cover
end

end


text:
using Documenter, MDToolbox

makedocs(sitename="MDToolbox.jl",
        pages = [
        "MDToolbox.jl" => "index.md",
        "Installation" => "installation.md",
        "Getting started" => [
            "getting_started01.md",
            "getting_started02.md",
            "getting_started03.md"
        ],
        "Examples" => [
            "superimpose_rmsd.md",
            "free_energy_surface.md",
            "wham_example.md"
        ],
        "References" => [
            "fileio.md",
            "structure.md",
            "reduction.md",
            "wham.md",
            "mbar.md",
            "msm.md"],
        "Workflow for developers" => "workflow_for_developers.md"
            ]
)

deploydocs(
    repo = "github.com/matsunagalab/MDToolbox.jl.git",
)

ext: jl


text:
# This file was generated, do not modify it. # hide
pipe.transformed_target_model.model.K = 2
pipe.one_hot_encoder.drop_last = true;

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function error_all_different(x::V; param = nothing, dom_size = 0
) where {T <: Number, V <: AbstractVector{T}}
    acc = Dictionary{T, Int}()
    foreach(y -> insert_or_inc(acc, y), x)
    return Float64(sum(acc .- 1))
end

concept_all_different(x) = allunique(x)

"""
    _all_different
Global constraint ensuring that all the values of a given configuration are unique.
"""
const all_different = Constraint(
    concept = concept_all_different,
    error = make_error(:all_different),
    syms = Set([:permutable]),
)


text:
include("attributes.jl")
include("hyperplane.jl")
include("channels.jl")
include("AbstractLShaped.jl")
include("feasibility.jl")
include("subproblem.jl")


text:
using CovarianceEstimation
using Statistics
using LinearAlgebra
using Test
using Random
using DelimitedFiles
using StatsBase

include("ref_lw_lshrink.jl")
include("ref_lw_nlshrink.jl")
include("legacy.jl")
Random.seed!(1234)

const CE = CovarianceEstimation

const X   = randn(3, 8)
const Z   = [2 -1 2; -1 2 -1; -1 -1 -1]
const X2s = [[1 0; 0 1; -1 0; 0 -1], [1 0; 0 1; 0 -1; -1 0]]
const Xa  = randn(5, 5)
const Xb  = randn(8, 3)
const Xc  = randn(15, 20)
const Xd  = randn(30, 50)
const Xe  = randn(50, 30)

const test_matrices = [X, Xa, Xb, Xc, Xd, Xe, Z]

function testTransposition(ce::CovarianceEstimator, X)
    @test cov(ce, X; dims=1) ≈ cov(ce, transpose(X); dims=2)
    @test cov(ce, X; dims=2) ≈ cov(ce, transpose(X); dims=1)
end

function testDims(ce::CovarianceEstimator, X)
    @test_throws ArgumentError cov(ce, X, dims=0)
    @test_throws ArgumentError cov(ce, X, dims=3)
end

function testUncorrelated(ce::CovarianceEstimator)
    for X2 ∈ X2s
        @test isdiag(cov(ce, X2))
    end

end

function testTranslation(ce::CovarianceEstimator, X)
    C1 = cov(ce, X)
    C2 = cov(ce, X .+ randn(1, size(X, 2)))
    @test C1 ≈ C2 atol = 1e-12 rtol = 1e-16
    C1t = cov(ce, X')
    C2t = cov(ce, X' .+ randn(1, size(X, 1)))
    @test C1t ≈ C2t atol = 1e-12 rtol = 1e-16
end

include("test_linearshrinkage.jl")
include("test_nonlinearshrinkage.jl")


text:
# SPDX-License-Identifier: X11
# 2020-12-02
# Peaks (300pt)

function main()
  N,M = parse.(Int, (split ∘ readline)())
  H = parse.(Int, (split ∘ readline)())
  MA = zeros(N)

  for i ∈ 1:M
    A,B = parse.(Int, (split ∘ readline)())
    MA[A] = max(MA[A], H[B])
    MA[B] = max(MA[B], H[A])
  end

  ans = 0
  for i ∈ 1:N
    H[i] > MA[i] && (ans += 1)
  end

  println(ans)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# This code is derived from the Julia SubArray implementation (MIT license)

# L is true if the view itself supports fast linear indexing
"""
    SSubArray

Statically sized SubArray.
"""
struct SSubArray{S,T,N,P,I,L} <: StaticArray{S, T, N}
    parent::P
    indices::I
    offset1::Int       # for linear indexing and pointer, only valid when L==true
    stride1::Int       # used only for linear indexing
    function SSubArray{S,T,N,P,I,L}(parent, indices, offset1, stride1) where {S,T,N,P,I,L}
        Base.@_inline_meta
        Base.check_parent_index_match(parent, indices)
        new(parent, indices, offset1, stride1)
    end
end
# Compute the linear indexability of the indices, and combine it with the linear indexing of the parent
function SSubArray(parent::Union{StaticArray{S}, HybridArray{S}}, indices::Tuple) where S
    Base.@_inline_meta
    SSubArray(IndexStyle(Base.viewindexing(indices), IndexStyle(parent)), all_dynamic_fixed_val(S, indices...), parent, Base.ensure_indexable(indices), Base.index_dimsum(indices...))
end
function SSubArray(::IndexCartesian, ::Val{:dynamic_fixed_true}, parent::Union{StaticArray{S}, HybridArray{S}}, indices::I, ::NTuple{N,Any}) where {S,I,N}
    Base.@_inline_meta
    newsize = new_out_size(S, indices...)
    SSubArray{newsize, eltype(typeof(parent)), N, typeof(parent), I, false}(parent, indices, 0, 0)
end
function SSubArray(::IndexLinear, ::Val{:dynamic_fixed_true}, parent::Union{StaticArray{S}, HybridArray{S}}, indices::I, ::NTuple{N,Any}) where {S,I,N}
    Base.@_inline_meta
    # Compute the stride and offset
    stride1 = Base.compute_stride1(parent, indices)
    newsize = new_out_size(S, indices...)
    SSubArray{newsize, eltype(typeof(parent)), N, typeof(parent), I, true}(parent, indices, Base.compute_offset1(parent, stride1, indices), stride1)
end
function SSubArray(indexing::Any, ::Val{:dynamic_fixed_false}, parent::Union{StaticArray{S}, HybridArray{S}}, indices::I, idsum::NTuple{N,Any}) where {S,I,N}
    return SubArray(indexing, parent, indices, idsum)
end

if !isdefined(Base, :__subarray_throw_boundserror)
    @noinline __subarray_throw_boundserror(::Type{T}, parent, indices, offset1, stride1, I) where {T} =
    throw(BoundsError(T(parent, indices, offset1, stride1), I))
else
    using Base: __subarray_throw_boundserror
end

# This makes it possible to elide view allocation in cases where the
# view is indexed with a boundscheck but otherwise all its uses
# are inlined
@inline Base.throw_boundserror(A::SSubArray, I) =
    __subarray_throw_boundserror(typeof(A), A.parent, A.indices, A.offset1, A.stride1, I)

# Simple utilities
size(V::SSubArray) = (Base.@_inline_meta; map(n->Int(Base.unsafe_length(n)), axes(V)))

similar(V::SSubArray, T::Type, dims::Dims) = similar(V.parent, T, dims)

Base.sizeof(V::SSubArray) = length(V) * sizeof(eltype(V))

copy(V::SSubArray) = V.parent[V.indices...]

Base.parent(V::SSubArray) = V.parent
Base.parentindices(V::SSubArray) = V.indices

## Aliasing detection
dataids(A::SSubArray) = (dataids(A.parent)..., Base._splatmap(dataids, A.indices)...)
Base.unaliascopy(A::SSubArray) = typeof(A)(Base.unaliascopy(A.parent), map(Base.unaliascopy, A.indices), A.offset1, A.stride1)

# When the parent is an Array we can trim the size down a bit. In the future this
# could possibly be extended to any mutable array.
function Base.unaliascopy(V::SSubArray{S,T,N,A,I,LD}) where {S,T,N,A<:Array,I<:Tuple{Vararg{Union{Real,AbstractRange,Array}}},LD}
    dest = Array{T}(undef, index_lengths(V.indices...))
    copyto!(dest, V)
    SSubArray{S,T,N,A,I,LD}(dest, map(_trimmedindex, V.indices), 0, Int(LD))
end

function Base.unsafe_view(A::Union{StaticArray, HybridArray}, I::Vararg{Base.ViewIndex,<:Any})
    Base.@_inline_meta
    SSubArray(A, I)
end

Base.unsafe_view(V::SSubArray, I::Vararg{Base.ViewIndex,N}) where {N} =
    (Base.@_inline_meta; Base._maybe_reindex(V, I))
Base._maybe_reindex(V::SSubArray, I) = (Base.@_inline_meta; Base._maybe_reindex(V, I, I))
Base._maybe_reindex(V::SSubArray, I, ::Tuple{AbstractArray{<:Base.AbstractCartesianIndex}, Vararg{Any}}) =
    (Base.@_inline_meta; SSubArray(V, I))
# But allow arrays of CartesianIndex{1}; they behave just like arrays of Ints
Base._maybe_reindex(V::SSubArray, I, A::Tuple{AbstractArray{<:Base.AbstractCartesianIndex{1}}, Vararg{Any}}) =
    (Base.@_inline_meta; Base._maybe_reindex(V, I, Base.tail(A)))
Base._maybe_reindex(V::SSubArray, I, A::Tuple{Any, Vararg{Any}}) = (Base.@_inline_meta; Base._maybe_reindex(V, I, Base.tail(A)))
function Base._maybe_reindex(V::SSubArray, I, ::Tuple{})
    Base.@_inline_meta
    @inbounds idxs = to_indices(V.parent, _reindex(V.indices, I))
    return Base.unsafe_view(V.parent, idxs...)
end

# reindex was copied from Base because in pre-1.2 it has an incompatible implementation
_reindex(::Tuple{}, ::Tuple{}) = ()

# Skip dropped scalars, so simply peel them off the parent indices and continue
_reindex(idxs::Tuple{Base.ScalarIndex, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) =
    (Base.@_propagate_inbounds_meta; (idxs[1], _reindex(Base.tail(idxs), subidxs)...))

# Slices simply pass their subindices straight through
_reindex(idxs::Tuple{Base.Slice, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (Base.@_propagate_inbounds_meta; (subidxs[1], _reindex(Base.tail(idxs), Base.tail(subidxs))...))

# Re-index into parent vectors with one subindex
_reindex(idxs::Tuple{AbstractVector, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (Base.@_propagate_inbounds_meta; (idxs[1][subidxs[1]], _reindex(Base.tail(idxs), Base.tail(subidxs))...))

# Parent matrices are re-indexed with two sub-indices
_reindex(idxs::Tuple{AbstractMatrix, Vararg{Any}}, subidxs::Tuple{Any, Any, Vararg{Any}}) =
    (Base.@_propagate_inbounds_meta; (idxs[1][subidxs[1], subidxs[2]], _reindex(Base.tail(idxs), Base.tail(Base.tail(subidxs)))...))

# In general, we index N-dimensional parent arrays with N indices
@generated function _reindex(idxs::Tuple{AbstractArray{T,N}, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) where {T,N}
    if length(subidxs.parameters) >= N
        subs = [:(subidxs[$d]) for d in 1:N]
        tail = [:(subidxs[$d]) for d in N+1:length(subidxs.parameters)]
        :(Base.@_propagate_inbounds_meta; (idxs[1][$(subs...)], _reindex(Base.tail(idxs), ($(tail...),))...))
    else
        :(throw(ArgumentError("cannot re-index SubArray with fewer indices than dimensions\nThis should not occur; please submit a bug report.")))
    end
end

# In general, we simply re-index the parent indices by the provided ones
SlowSSubArray{S,T,N,P,I} = SSubArray{S,T,N,P,I,false}
function getindex(V::SSubArray{S,T,N}, I::Vararg{Int,N}) where {S,T,N}
    Base.@_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[_reindex(V.indices, I)...]
    r
end

# overrides default erroring getindex from `StaticArrays`
function getindex(V::SSubArray, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    subind = Base._to_subscript_indices(V, i)
    @inbounds r = V.parent[_reindex(V.indices, subind)...]
    return r
end

# But SSubArrays with fast linear indexing pre-compute a stride and offset
FastSSubArray{S,T,N,P,I} = SSubArray{S,T,N,P,I,true}
function getindex(V::FastSSubArray, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + V.stride1*i]
    r
end

# disambiguation
function getindex(V::FastSSubArray{S,T,1}, i::Int) where {S,T}
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + V.stride1*i]
    r
end

# We can avoid a multiplication if the first parent index is a Colon or AbstractUnitRange,
# or if all the indices are scalars, i.e. the view is for a single value only
FastContiguousSSubArray{S,T,N,P,I<:Union{Tuple{Union{Base.Slice, AbstractUnitRange}, Vararg{Any}},
                                      Tuple{Vararg{Base.ScalarIndex}}}} = SSubArray{S,T,N,P,I,true}
function getindex(V::FastContiguousSSubArray, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + i]
    r
end
# For vector views with linear indexing, we disambiguate to favor the stride/offset
# computation as that'll generally be faster than (or just as fast as) re-indexing into a range.
function getindex(V::FastSSubArray{<:Any, 1}, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + V.stride1*i]
    r
end
function getindex(V::FastContiguousSSubArray{<:Any, 1}, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + i]
    r
end

# Indexed assignment follows the same pattern as `getindex` above
function setindex!(V::SSubArray{T,N}, x, I::Vararg{Int,N}) where {T,N}
    Base.@_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds V.parent[_reindex(V.indices, I)...] = x
    V
end
# overrides the default setindex! from `StaticArrays`
function setindex!(V::SSubArray, x, i::Int)
    Base.@_inline_meta
    subind = Base._to_subscript_indices(V, i)
    @inbounds V.parent[_reindex(V.indices, subind)...] = x
    V
end
function setindex!(V::FastSSubArray, x, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + V.stride1*i] = x
    V
end
function setindex!(V::FastContiguousSSubArray, x, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + i] = x
    V
end
function setindex!(V::FastSSubArray{<:Any, 1}, x, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + V.stride1*i] = x
    V
end
function setindex!(V::FastContiguousSSubArray{<:Any, 1}, x, i::Int)
    Base.@_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + i] = x
    V
end

Base.IndexStyle(::Type{<:FastSSubArray}) = IndexLinear()
Base.IndexStyle(::Type{<:SSubArray}) = IndexCartesian()

if VERSION < v"1.2"
    # substrides copied from base because it's different in Julia < v1.2
    substrides(strds::Tuple{}, ::Tuple{}) = ()
    substrides(strds::NTuple{N,Int}, I::Tuple{Base.ScalarIndex, Vararg{Any}}) where N = (substrides(Base.tail(strds), Base.tail(I))...,)
    substrides(strds::NTuple{N,Int}, I::Tuple{Base.Slice, Vararg{Any}}) where N = (first(strds), substrides(Base.tail(strds), Base.tail(I))...)
    substrides(strds::NTuple{N,Int}, I::Tuple{AbstractRange, Vararg{Any}}) where N = (first(strds)*step(I[1]), substrides(Base.tail(strds), Base.tail(I))...)
    substrides(strds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SSubArrays with indices of type $(typeof(I[1]))"))
else
    using Base: substrides
end

# Strides are the distance in memory between adjacent elements in a given dimension
# which we determine from the strides of the parent
Base.strides(V::SSubArray) = substrides(strides(V.parent), V.indices)

Base.stride(V::SSubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]

Base.elsize(::Type{<:SSubArray{<:Any,<:Any,<:Any,P}}) where {P} = Base.elsize(P)

Base.iscontiguous(A::SSubArray) = Base.iscontiguous(typeof(A))
Base.iscontiguous(::Type{<:SSubArray}) = false
Base.iscontiguous(::Type{<:FastContiguousSSubArray}) = true

Base.first_index(V::FastSSubArray) = V.offset1 + V.stride1 # cached for fast linear SSubArrays
function Base.first_index(V::SSubArray)
    P, I = parent(V), V.indices
    s1 = Base.compute_stride1(P, I)
    s1 + Base.compute_offset1(P, s1, I)
end
Base.unsafe_convert(::Type{Ptr{T}}, V::SSubArray{T,N,P,<:Tuple{Vararg{Base.RangeIndex}}}) where {T,N,P} =
    unsafe_convert(Ptr{T}, V.parent) + (first_index(V)-1)*sizeof(T)

Base.pointer(V::FastSSubArray, i::Int) = pointer(V.parent, V.offset1 + V.stride1*i)
Base.pointer(V::FastContiguousSSubArray, i::Int) = pointer(V.parent, V.offset1 + i)
Base.pointer(V::SSubArray, i::Int) = Base._pointer(V, i)
if isdefined(Base, :_pointer)
    Base._pointer(V::SSubArray{<:Any,1}, i::Int) = pointer(V, (i,))
    Base._pointer(V::SSubArray, i::Int) = pointer(V, Base._ind2sub(axes(V), i))
end

function Base.pointer(V::SSubArray{T,N,<:Array,<:Tuple{Vararg{Base.RangeIndex}}}, is::Tuple{Vararg{Int}}) where {T,N}
    index = Base.first_index(V)
    strds = strides(V)
    for d = 1:length(is)
        index += (is[d]-1)*strds[d]
    end
    return pointer(V.parent, index)
end

if VERSION < v"1.2"
    # copied from Base because in Julia < 1.2 it works in an incompatible way
    _indices_sub(::Real, I...) = (Base.@_inline_meta; _indices_sub(I...))
    _indices_sub() = ()
    function _indices_sub(i1::AbstractArray, I...)
        Base.@_inline_meta
        (Base.unsafe_indices(i1)..., _indices_sub(I...)...)
    end
else
    using Base: _indices_sub
end

axes(S::SSubArray) = (Base.@_inline_meta; _indices_sub(S.indices...))


text:
module MultiLearners

export MultiLearner, transform!, fit!
export multirun

using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils: mergedict
using TSML.System # import LIB_SKL_AVAILABLE

if LIB_SKL_AVAILABLE
    using TSML.SKLearners
end

mutable struct MultiLearner <: TSLearner
    model
    args
    data

    function MultiLearner(args=Dict())
        default_args=Dict(
            :output => :numeric,
            :learner => "RandomForestRegressor",
            :impl_args => Dict()
        )
        new(nothing,mergedict(default_args,args),nothing)
    end
end

function fit!(mlt::MultiLearner, x::T,y::Vector) where {T<:Union{Matrix,Vector}}
    @info "fit: to be implemented"
    @info "matrify x"
    @info "fit(matrified x, y)"
end

function transform!(mlt::MultiLearner,x::T) where {T<:Union{Matrix,Vector}}
    @info "transform: to be implemented"
end

function multirun()
    model = MultiLearner()
    println(model)
end

end


text:
export alloc_opt!, barrier_noop!, gc_invariant_verifier!, lower_exc_handlers!,
       combine_mul_add!, multi_versioning!, propagate_julia_addrsp!, lower_ptls!,
       lower_simdloop!, late_lower_gc_frame!, final_lower_gc!, remove_julia_addrspaces!,
       demote_float16!, remove_ni!, julia_licm!

alloc_opt!(pm::PassManager) =
    API.LLVMAddAllocOptPass(pm)

barrier_noop!(pm::PassManager) =
    API.LLVMAddBarrierNoopPass(pm)

gc_invariant_verifier!(pm::PassManager, strong::Bool=false) =
    API.LLVMAddGCInvariantVerifierPass(pm, convert(Bool, strong))

lower_exc_handlers!(pm::PassManager) =
    API.LLVMAddLowerExcHandlersPass(pm)

combine_mul_add!(pm::PassManager) =
    API.LLVMAddCombineMulAddPass(pm)

multi_versioning!(pm::PassManager) =
    API.LLVMAddMultiVersioningPass(pm)

propagate_julia_addrsp!(pm::PassManager) =
    API.LLVMAddPropagateJuliaAddrspaces(pm)

lower_ptls!(pm::PassManager, imaging_mode::Bool=false) =
    API.LLVMAddLowerPTLSPass(pm, convert(Bool, imaging_mode))

lower_simdloop!(pm::PassManager) =
    API.LLVMAddLowerSimdLoopPass(pm)

late_lower_gc_frame!(pm::PassManager) =
    API.LLVMAddLateLowerGCFramePass(pm)

if VERSION >= v"1.3.0-DEV.95"
    final_lower_gc!(pm::PassManager) =
        API.LLVMAddFinalLowerGCPass(pm)
else
    final_lower_gc!(pm::PassManager) = nothing
end

if VERSION >= v"1.5.0-DEV.802"
    remove_julia_addrspaces!(pm::PassManager) = API.LLVMAddRemoveJuliaAddrspacesPass(pm)
else
    remove_julia_addrspaces!(pm::PassManager) = nothing
end

if VERSION >= v"1.6.0-DEV.1215"
    demote_float16!(pm::PassManager) = API.LLVMAddDemoteFloat16Pass(pm)
else
    demote_float16!(pm::PassManager) = nothing
end

if VERSION >= v"1.6.0-DEV.1476"
    remove_ni!(pm::PassManager) = API.LLVMAddRemoveNIPass(pm)
else
    remove_ni!(pm::PassManager) = nothing # Could implement this in pure Julia if necessary
end

if VERSION >= v"1.6.0-DEV.1477"
    julia_licm!(pm::PassManager) = API.LLVMAddJuliaLICMPass(pm)
else
    julia_licm!(pm::PassManager) = nothing
end


text:
"""
    prefix_pod(name) -> AbstractString

Prepends the prexfix `pod/` to `name` if it does not
already have it.
"""
function prefix_pod(name)
    if !startswith(name, "pod/")
        name = "pod/" * name
    end
    return name
end

"""
    get_status(pod; namespace=nothing) -> String

Get the status of a pod.
"""
function get_status(pod; namespace=nothing)
    pod = prefix_pod(pod)
    ns = isnothing(namespace) ? `` : `--namespace=$namespace`
    return readchomp(`$(kubectl()) get $pod $ns --template=\{\{.status.phase\}\}`)
end

"""
    last_condition(pod; namespace=nothing) -> String

Get the latest "condition" of a pod.
"""
function last_condition(pod; namespace=nothing)
    pod = prefix_pod(pod)
    ns = isnothing(namespace) ? `` : `--namespace=$namespace`
    return JSON3.read(readchomp(`$(kubectl()) get $pod $ns -o jsonpath=\{.status.conditions\[-1:\]\}`))
end

"""
    get_pod_names(labels = ``; namespace=nothing) -> Vector{String}

Get the names of pods in the current namespace (or pass a namespace to specify
a different one). Pass `labels` as e.g. `-l app=myapp` to restrict to specific labels.
"""
function get_pod_names(labels=``; namespace=nothing)
    ns = isnothing(namespace) ? `` : `--namespace=$namespace`
    pods = JSON3.read(readchomp(`$(kubectl()) get pods $ns $labels -o jsonpath="{.items}"`))
    return [pod.metadata.name for pod in pods]
end

"""
    get_current_namespace() -> String

Get the current active Kubernetes namespace.
"""
function get_current_namespace()
    return readchomp(`$(kubectl()) config view --minify --output 'jsonpath={..namespace}'`)
end

"""
    wait_until_pod_ready(pod; namespace=nothing)

Display a progress bar while waiting for `pod` to become `Running`.
"""
function wait_until_pod_ready(pod; namespace=nothing, exit_on_interrupt=false)
    pod = prefix_pod(pod)
    runner(; interrupt_msg="while waiting for $pod to start up", exit_on_interrupt) do
        # We don't want to have to wait to the `kubectl` request to finish in order to display
        # values below the progress bar, or else it will flash annoyingly.
        # So instead we do the `last_condition` calling async from the rest.

        # stores the result
        condition_channel = Channel(1) do c
            while true
                # Put a new value as soon as we can
                put!(c, last_condition(pod; namespace))
            end
        end

        # Show the last condition of the pod dynamically under the ProgressMeter display
        showvalues = () -> begin
            # take the latest, waiting if we need to to get it
            c = take!(condition_channel)
            return pairs(c)
        end

        prog = ProgressUnknown("Starting up $(pod):"; spinner=true)
        while get_status(pod; namespace) == "Pending"
            ProgressMeter.next!(prog; showvalues)
            sleep(2)
        end

        st = get_status(pod; namespace)

        if st == "Succeeded" || st == "Running"
            ProgressMeter.finish!(prog)
        else
            ProgressMeter.finish!(prog; spinner='✗')
            error("""
            Pod did not start up successfully. Got status $(st) and last condition status:

            $(JSON3.pretty(last_condition(pod; namespace)))

            """)
        end
    end
    return nothing
end

"""
    watch_logs(pod; exit_on_interrupt=false, namespace=nothing)

Follows logs from `pod`. If `exit_on_interrupt=true`, exits the Julia session
upon `ctrl-c`. Otherwise throws an error as usual upon interruption.
"""
function watch_logs(pod; exit_on_interrupt=false, namespace=nothing)
    pod = prefix_pod(pod)
    ns = isnothing(namespace) ? `` : `--namespace=$namespace`
    cmd_runner(`$(kubectl()) logs $ns -f $pod`; exit_on_interrupt,
               interrupt_msg="while following the logs of $pod")
    return nothing
end

struct Interrupted <: Exception
    msg::String
end
Base.showerror(io::IO, e::Interrupted) = print(io, "Interrupted ", e.msg)

function cmd_runner(cmd::Cmd; exit_on_interrupt=false, interrupt_msg="Interrupted")
    # We use `wait=false` and `wait_kill_result=true`
    # So that if we ctrl-c we don't get
    # a LoadError (failed process) but rather
    # a usual InterruptException that we can handle.
    runner(() -> run(pipeline(cmd; stdout, stderr); wait=false); exit_on_interrupt,
           interrupt_msg, wait_kill_result=true)
    return nothing
end

function runner(f; exit_on_interrupt=false, interrupt_msg="Interrupted",
                wait_kill_result=false)
    local task
    try
        # Makes interrupts catchable like in the REPL, even if run from a script.
        Base.exit_on_sigint(false)
        task = f()
        wait_kill_result && wait(task)
    catch e
        wait_kill_result && kill(task)
        if e isa InterruptException
            if exit_on_interrupt
                # No need for a stacktrace
                println("Interrupted $(interrupt_msg); exiting...")
                exit(0)
            else
                rethrow(Interrupted(interrupt_msg))
            end
        end
        rethrow()
    finally
        # In non-interative sessions, this defaults to `true`.
        # We should restore that.
        if !Base.isinteractive()
            Base.exit_on_sigint(true)
        end
    end
    return nothing
end

"""
    follow(pod; exit_on_interrupt=false, namespace=nothing)

A simpler helper to follow a `pod` by first waiting until it's
ready ([`wait_until_pod_ready`](@ref)) and then wathing its logs
([`watch_logs`](@ref)).
"""
function follow(pod; exit_on_interrupt=false, namespace=nothing)
    wait_until_pod_ready(pod; exit_on_interrupt, namespace)
    watch_logs(pod; exit_on_interrupt, namespace)
    return nothing
end

"""
    port_forward(pod; remote_port::Int, local_port::Int=remote_port, namespace=nothing) -> Task

Forwards a port from `remote_port` on `pod` to `local_port`. Runs in a detached process,
and so will outlive the current Julia session. Returns the task running this process.
"""
function port_forward(pod; remote_port::Int, local_port::Int=remote_port, namespace=nothing)
    pod = prefix_pod(pod)
    ns = isnothing(namespace) ? `` : `--namespace=$namespace`

    # Launch port forwarding in a separate process group (`detach=true`) so it can outlive this session
    pfd_task = run(Cmd(`$(kubectl()) port-forward $pod $ns $(local_port):$(remote_port)`;
                       detach=true); wait=false)

    println("Fowarding port $(remote_port) of $pod to your local port `$(local_port)`.")
    println("Go to `http://localhost:$(local_port)/` to see the results.")
    return pfd_task
end


text:
# =========================================================================================
"""
    abstract type QuadraturePoints{T <: AbstractFloat} end

Abstract base type for all sets of quadrature points
"""
abstract type QuadraturePoints{T <: AbstractFloat} end

# =========================================================================================
"""
    struct QuadPts2D{T} <: QuadraturePoints{T}
        num   ::Int         # number of points
        x     ::Vector{T}   # x values
        y     ::Vector{T}   # y values
        weight::Vector{T}   # weights
    end

Quadrature points and weights for triangles
"""
struct QuadPts2D{T} <: QuadraturePoints{T}
    """Number of points"""
    num::Int
    """x values"""
    x::Vector{T}
    """y values"""
    y::Vector{T}
    """Weights"""
    weight::Vector{T}
end

# =========================================================================================
"""
    struct QuadPts3D{T} <: QuadraturePoints{T}
        num   ::Int         # number of points
        x     ::Vector{T}   # x values
        y     ::Vector{T}   # y values
        z     ::Vector{T}   # z values
        weight::Vector{T}   # weights
    end

Quadrature points and weights for tetrahedra
"""
struct QuadPts3D{T} <: QuadraturePoints{T}
    """Number of points"""
    num::Int
    """x values"""
    x::Vector{T}
    """y values"""
    y::Vector{T}
    """z values"""
    z::Vector{T}
    """Weights"""
    weight::Vector{T}
end

# =========================================================================================
@doc """
    quadraturepoints(::Type{Tetrahedron{Float64}})
    quadraturepoints(::Type{Tetrahedron{Float32}})
    quadraturepoints(::Type{Triangle{Float64}})
    quadraturepoints(::Type{Triangle{Float32}})

Generator function for quadrature points:
 * *Triangles*: 7 points per element [[Rad48]](@ref Bibliography)
 * *Tetrahedra*: 5 points per element [[Kea86]](@ref Bibliography)

## Return type
`QuadPts2D` or `QuadPts3D`
""" quadraturepoints
for T in [:Float64, :Float32]
    varname = Symbol("triquadpts_", T)
    @eval begin
        const $(varname) = QuadPts2D(7,
            $(T)[1/3, (6+√15)/21, (9-2*√15)/21, (6+√15)/21,
                (6-√15)/21, (9+2*√15)/21, (6-√15)/21],
            $(T)[1/3, (9-2*√15)/21, (6+√15)/21, (6+√15)/21,
                (9+2*√15)/21, (6-√15)/21, (6-√15)/21],
            $(T)[9/80, (155+√15)/2400, (155+√15)/2400, (155+√15)/2400,
                (155-√15)/2400, (155-√15)/2400, (155-√15)/2400]
        )
        quadraturepoints(::Type{Triangle{$(T)}}) = $(varname)
    end
end

for T in [:Float64, :Float32]
    varname = Symbol("tetraquadpts_", T)
    @eval begin
        const $(varname) = QuadPts3D(5,
            $(T)[.25, .5, 1/6, 1/6, 1/6],
            $(T)[.25, 1/6, 1/6, 1/6, .5],
            $(T)[.25, 1/6, 1/6, .5, 1/6],
            $(T)[-.8, .45, .45, .45, .45]/6
        )
        quadraturepoints(::Type{Tetrahedron{$(T)}}) = $(varname)
    end
end


text:
"""
    place_image(p::Point, colmat::VecOrMat; centered = true, alpha = missing)
    -> (upper left point, lower right point)

The color map is centered on p by default, one pixel per value in the matrix or vector.

Width and height is determined by the size of the matrix.

'alpha' mayhave no effect, but is passed on to
   https://www.cairographics.org/manual/cairo-cairo-t.html#cairo-paint-with-alpha

post an issue here if it does in any case.
"""
function place_image(p::Point, colmat::VecOrMat; centered = true, alpha = missing)
    place_image(p, _CairoSurfaceBase(colmat); centered)
end

function place_image(p::Point, colmat::VecOrMat{T}; centered = true, alpha = missing) where T<:Union{Luxor.ARGB32, UInt32}
    # layer to overcome a bug where Luxor would not draw the following strokes
    @layer placeimage(data, p; centered = centered)
    Δp = centered * Point(nx / 2, ny / 2)
    (p - Δp, p - Δp + (nx, ny))
end

"""
    place_image(pos::Point, image::SVGimage,
                  width = missing, height = missing, scalefactor = missing,
                  centered = false)
    place_image(pos::Point, formula::LaTeXString; modifylatex = true, modifysvgcolor = true,
                  width = missing, height = missing, scalefactor = missing,
                  centered = false)
    place_image(pos::Point, data::Luxor.Cairo.CairoSurfaceBase; centered = true, alpha = missing,
                  width = missing, height = missing, scalefactor = missing)
    -> (upper left point, lower right point, scalefactor)

Return the upper left, lower right corners, also the calculated scaling factor from original to placed image.

'image::SVGimage' can be made using MechanicalSketch.readsvg (from Luxor, as usual). It takes files or string input.
'formula::LaTeXString' can be made using MechanicalSketch.@latexify and other methods, see Latexify.
                       Reuse output 'scalefactor' as input for following formula images.
'data::Luxor.Cairo.CairoSurfaceBase' can be made using MechanicalSketch.readpng (from Luxor).

'height', 'width' and 'scalefactor': Uniform scaling should work with quantities or numbers (i.e. pixels or points).

'modifylatex' = true: Calls 'modify_latex', which subscripts bracketed indexes and a number of other adaptions.
'modifysvgcolor' = true: Replace currentColor with the last color and opacity from 'setcolor', 'setopacity' or 'sethue'.
'alpha' = 0.0 - 1.0 : 0.0 is completely translucent, 1.0 is opaque. The effect vary depending on input type.
"""
function place_image(pos::Point, formula::LaTeXString; modifylatex = true, modifysvgcolor = true,
    width = missing, height = missing, scalefactor = missing,
    centered = false)
    svim = readsvg(tex2svg_string(formula; modifylatex, modifysvgcolor))
    place_image(pos, svim ; width, height, scalefactor, centered)
end
function place_image(pos::Point, image::SVGimage;
                      width = missing, height = missing, scalefactor = missing,
                      centered = false)
    scalefac, original_width, original_height = scalingfactor(image, width, height, scalefactor)
    # Destination size
    dest_width_pix, dest_height_pix = scale_to_pt.(scalefac .* (original_width, original_height))
    # Destination upper left corner
    ptupleft = centered ?  pos - 0.5 .* (dest_width_pix, dest_height_pix) : pos
    @layer begin
        translate(ptupleft)
        scale(scalefac)
        placeimage(image, Point(0, 0); centered = false)
    end
    ptupleft, ptupleft + (dest_width_pix, dest_height_pix), scalefac
end

function place_image(pos::Point, data::Luxor.Cairo.CairoSurfaceBase; centered = true, alpha = missing,
                     width = missing, height = missing, scalefactor = missing)
    scalefac, original_width, original_height = scalingfactor(data, width, height, scalefactor)
    # Destination size
    dest_width_pix, dest_height_pix = scale_to_pt.(scalefac .* (original_width, original_height))
    # Destination upper left corner (we'll place the image non-centered afterwards)
    ptupleft = centered ?  pos - 0.5 .* (dest_width_pix, dest_height_pix) : pos
    @layer begin
        translate(ptupleft)
        scale(scalefac)
        ismissing(alpha) ? placeimage(data, Point(0, 0); centered = false) : placeimage(data, Point(0, 0), alpha; centered = false)
    end
    ptupleft, ptupleft + (dest_width_pix, dest_height_pix), scalefac
end
"""
    scalingfactor(image, width, height, scalefactor)
    -> (scalefac, original_width, original_height)
"""
function scalingfactor(image, width, height, scalefactor::Union{Float64, Missing})
    @assert !ismissing(width) + !ismissing(height) + !ismissing(scalefactor) < 2 "Only one of width, height and scalefactor can be specified."
    original_width, original_height = get_width_height(image)
    # Find scaling from input to output
    scalefac = if ismissing(scalefactor)
        if ismissing(width) && ismissing(height)
            1.0
        elseif ismissing(height)
            scale_to_pt(width) / scale_to_pt(original_width)
        elseif ismissing(width)
            scale_to_pt(height) / scale_to_pt(original_height)
        end
    else
        scalefactor
    end
    scalefac, original_width, original_height
end

"""
get_width_height(image)
--> (w, h)::Quantity{Length}

'image' can be of types
- SVGimage (output from 'readsvg'),
- CairoSurfaceBase{UInt32} (output from 'readpng')
- matrix  - it is assumed that matrices follow the image convention: A column represent horizontal pixels.

The width and height of image based on scale_to_pt(m).
Note that the width can be inaccurate, as this is ambiguous for svg input.
"""
get_width_height(image) =  scale_pt_to_unit(m) .* (image.width, image.height)
get_width_height(image::Matrix) =  scale_pt_to_unit(m) .* reverse(size(matrix))

"""
    _CairoSurfaceBase(colmat::VecOrMat{T}) where T

Convert a matrix or vector to a cairo surface suitable for ´placeimage´.
One pixel per element. Relies on ImageIO / ImageMagick because it seems
to work more easily.
"""
function _CairoSurfaceBase(colmat::VecOrMat{T}) where T
    msg = "Can't display matrices of eltype($(eltype(colmat))). \n\tConsider defining a subtypes(AbstractColorLegend) to convert."
    @assert T <: Union{RGB, RGBA, Luxor.Cairo.ARGB32, HSV, HSVA,  HSL, HSLA,LCHab, LCHabA, LCHuv, LCHuvA}  msg
    tempfilename = joinpath(tempname(), "temp_sketch.png")
    # write to a file using ImageIO. Results for transparent colours seem to differ.
    save(File(format"PNG", tempfilename), colmat)
    # Read this image using Cairo's interface to get the right format for placing into other cairo surfaces.
    png = readpng(tempfilename);
    rm(tempfilename)
    png
end


text:
using ArgParse, Core, StatsBase
include("../src/debruijnnets.jl");
include("../src/sampling.jl");
include("../src/motifs.jl")
include("../src/randomwalks.jl")
include("../src/hypa.jl")

function compare_nodes(observed_nodes, sampled_nodes)
    observed = length(intersect(observed_nodes, sampled_nodes))
    unobserved = length(sampled_nodes) - observed
    return observed, unobserved
end

function compare_edges(observed_edges, sampled_edges)
    observed = length(intersect(observed_edges, sampled_edges))
    unobserved = length(sampled_edges) - observed
    return observed, unobserved
end

"""

"""
function sample_properties(input_file, k, frequency, walk_interval, max_intervals, num_samples, output_file)
    print_interval = 100
    println("Reading graph.")
    fo, fo_map, ko, ko_map = from_ngram(input_file, frequency, k)
    rev_fo_map = Dict{Integer, String}(val=>key for (key,val) in fo_map)
    # Observed first-order nodes/edges
    observed_fo_nodes = length(fo_map)
    observed_fo_edges = ne(fo)
    # Observed kth-order nodes
    observed_nodes = Set{Tuple}(keys(ko_map))
    observed_korder_nodes = length(observed_nodes)
    println("Done.")
    println("Computing all possible k-edge walks.")
    all_kedge_walks, all_korder_nodes = get_all_walks_with_nodes(fo, k)
    total_korder_nodes = length(all_korder_nodes)
    total_korder_edges = length(all_kedge_walks)
    println("Done.")
    println("observed_korder_nodes: $(observed_korder_nodes)")
    println("total_korder_nodes: $(total_korder_nodes)")
    observed_edges, weights = walks_from_edges(ko, ko_map, fo_map)
    observed_edgeset = Set{Tuple}(observed_edges)
    observed_ko_edges = length(observed_edges)
    observed_sampled_nodes = zeros(Float64, num_samples, max_intervals)
    unobserved_sampled_nodes = zeros(Float64, num_samples, max_intervals)
    observed_sampled_edges = zeros(Float64, num_samples, max_intervals)
    unobserved_sampled_edges = zeros(Float64, num_samples, max_intervals)
    missing_nodes = zeros(Float64, num_samples, max_intervals)
    missing_edges = zeros(Float64, num_samples, max_intervals)
    Threads.@threads for run in range(1, num_samples)
        println("Run: $run")
        sampled_walks = sample(all_kedge_walks, walk_interval)
        sampled_edgeset = Set{Tuple}(sampled_walks)
        interval_count = 0
        intervals_run = 0
        # Redefining makes this type stable
        non::Integer = 0
        nun::Integer = 0
        noe::Integer = 0
        nue::Integer = 0
        i = 1
        while (nun / total_korder_nodes) < 1.0 && i <= max_intervals
            # Construct a kth-order graph from these walks
            fo_s, ko_s, new_ko_map = from_walks(sampled_walks, k, rev_fo_map)
            sampled_nodes = Set{Tuple}(keys(new_ko_map))

            # Compare nodes/edges from sampled walks to nodes from observed walks
            non, nun = compare_nodes(observed_nodes, sampled_nodes)
            noe, nue = compare_edges(observed_edgeset, sampled_edgeset)

            # Update ko stats
            observed_sampled_nodes[run,i] = non / observed_korder_nodes
            unobserved_sampled_nodes[run,i] = nun / total_korder_nodes
            observed_sampled_edges[run,i] = noe / observed_ko_edges
            unobserved_sampled_edges[run,i] = nue / total_korder_edges

            # Update fo stats
            missing_nodes[run,i] = nv(fo_s) / observed_fo_nodes
            missing_edges[run,i] = ne(fo_s) / observed_fo_edges

            # Compute next set of walks
            nxt_walks = sample(all_kedge_walks, walk_interval)
            append!(sampled_walks, nxt_walks)
            union!(sampled_edgeset, Set(nxt_walks))
            # Print for logging
            interval_count += 1
            if interval_count == print_interval
                println("i: $i, num walks: $(i*walk_interval), proportion unobserved: $(nun/total_korder_nodes)")
                interval_count = 0
            end
            i += 1
        end
    end

    x = [walk_interval*i for i in range(1, max_intervals)]
    observed_sampled_nodes_stats = mean_and_std(observed_sampled_nodes, 1)
    unobserved_sampled_nodes_stats = mean_and_std(unobserved_sampled_nodes, 1)
    observed_sampled_edges_stats = mean_and_std(observed_sampled_edges, 1)
    unobserved_sampled_edges_stats = mean_and_std(unobserved_sampled_edges, 1)

    missing_nodes_stats = mean_and_std(missing_nodes, 1)
    missing_edges_stats = mean_and_std(missing_edges, 1)

    open(output_file, write=true) do file
        out_str = "$(join(x, ','))\n"
        write(file, out_str)
        for arr in [observed_sampled_nodes_stats, unobserved_sampled_nodes_stats, observed_sampled_edges_stats, unobserved_sampled_edges_stats, missing_nodes_stats, missing_edges_stats]
            out_str = "$(join(arr[1], ','))|$(join(arr[2], ','))\n"
            write(file, out_str)
        end
    end
end

function parse_commandline(args)
    s = ArgParseSettings()
    @add_arg_table s begin
	"--file", "-g"
            help = "Full path to input file"
	"--frequency", "-f"
	    help = "Frequencies in ngram?"
	    action = :store_true
        "-k"
	    help = "Order of the motifs."
	    arg_type = Int
    	"-w", "--walk_interval"
	    help = "Number of walks per iteration."
	    arg_type = Int
	"-i", "--iterations"
	    help = "Number of iterations."
	    arg_type = Int
	"-r", "--runs"
	    help = "Number of runs."
	    arg_type = Int
    end
return parse_args(args, s)
end

function main()
    arguments = parse_commandline(ARGS)
    println(arguments)
    input_filepath = arguments["file"]
    println(input_filepath)
    frequency = arguments["frequency"]
    println(frequency)
    k = arguments["k"]
    runs = arguments["runs"]
    walk_interval = arguments["walk_interval"]
    max_intervals = arguments["iterations"]
    num_cpus = Threads.nthreads()
    splitpath = split(input_filepath, '/')
    ngram_filename = splitpath[end]
    println("k: $k")
    println(runs)
    output_file = "../../debruijn-nets/results/motifs/$(ngram_filename)_k-$(k)_samplestats.csv"
    sample_properties(input_filepath, k, frequency, walk_interval, max_intervals, runs, output_file)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
## macro for argument checking

macro check_args(D, cond)
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(D)), ": the condition ", $(string(cond)), " is not satisfied.")))
        end
    end
end

##### Utility functions

isunitvec(v::AbstractVector) = (norm(v) - 1.0) < 1.0e-12

isprobvec(p::AbstractVector{<:Real}) =
    all(x -> x ≥ zero(x), p) && isapprox(sum(p), one(eltype(p)))

pnormalize!(v::AbstractVector{<:Real}) = (v ./= sum(v); v)

add!(x::AbstractArray, y::AbstractVector) = broadcast!(+, x, x, y)
add!(x::AbstractArray, y::Zeros) = x

multiply!(x::AbstractArray, c::Number) = (x .*= c; x)

exp!(x::AbstractArray) = (x .= exp.(x); x)

# get a type wide enough to represent all a distributions's parameters
# (if the distribution is parametric)
# if the distribution is not parametric, we need this to be a float so that
# inplace pdf calculations, etc. allocate storage correctly
@inline partype(::Distribution) = Float64

# for checking the input range of quantile functions
# comparison with NaN is always false, so no explicit check is required
macro checkquantile(p,ex)
    p, ex = esc(p), esc(ex)
    :(zero($p) <= $p <= one($p) ? $ex : NaN)
end

macro checkinvlogcdf(lp,ex)
    lp, ex = esc(lp), esc(ex)
    :($lp <= zero($lp) ? $ex : NaN)
end

# because X == X' keeps failing due to floating point nonsense
function isApproxSymmmetric(a::AbstractMatrix{Float64})
    tmp = true
    for j in 2:size(a, 1)
        for i in 1:(j - 1)
            tmp &= abs(a[i, j] - a[j, i]) < 1e-8
        end
    end
    return tmp
end

# because isposdef keeps giving the wrong answer for samples
# from Wishart and InverseWisharts
hasCholesky(a::Matrix{Float64}) = isa(trycholesky(a), Cholesky)

function trycholesky(a::Matrix{Float64})
    try cholesky(a)
    catch e
        return e
    end
end

"""
    ispossemdef(A, k) -> Bool
Test whether a matrix is positive semi-definite with specified rank `k` by
checking that `k` of its eigenvalues are positive and the rest are zero.
# Examples
```jldoctest; setup = :(using Distributions: ispossemdef)
julia> A = [1 0; 0 0]
2×2 Array{Int64,2}:
 1  0
 0  0
julia> ispossemdef(A, 1)
true
julia> ispossemdef(A, 2)
false
```
"""
function ispossemdef(X::AbstractMatrix, k::Int;
                     atol::Real=0.0,
                     rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
    _check_rank_range(k, minimum(size(X)))
    ishermitian(X) || return false
    dp, dz, dn = eigsigns(Hermitian(X), atol, rtol)
    return dn == 0 && dp == k
end
function ispossemdef(X::AbstractMatrix;
                     atol::Real=0.0,
                     rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
    ishermitian(X) || return false
    dp, dz, dn = eigsigns(Hermitian(X), atol, rtol)
    return dn == 0
end

function _check_rank_range(k::Int, n::Int)
    0 <= k <= n || throw(ArgumentError("rank must be between 0 and $(n) (inclusive)"))
    nothing
end

#  return counts of the number of positive, zero, and negative eigenvalues
function eigsigns(X::AbstractMatrix,
                  atol::Real=0.0,
                  rtol::Real=(minimum(size(X))*eps(real(float(one(eltype(X))))))*iszero(atol))
    eigs = eigvals(X)
    eigsigns(eigs, atol, rtol)
end
function eigsigns(eigs::Vector{<: Real}, atol::Real, rtol::Real)
    tol = max(atol, rtol * eigs[end])
    eigsigns(eigs, tol)
end
function eigsigns(eigs::Vector{<: Real}, tol::Real)
    dp = count(x -> tol < x, eigs)        #  number of positive eigenvalues
    dz = count(x -> -tol < x < tol, eigs) #  number of numerically zero eigenvalues
    dn = count(x -> x < -tol, eigs)       #  number of negative eigenvalues
    return dp, dz, dn
end


text:
#!/usr/bin/env julia
println("helloWorld!")

f(x::Float64, y::Float64) = 2x + y
println(f(2.3,4.5))

max_stars_repo_name: combofish/chips-get

max_issues_repo_name: combofish/chips-get

max_forks_repo_name: combofish/chips-get


text:
# Generate notebooks using literate
using Pkg
Pkg.instantiate()

using SIIPExamples
using Literate

SIIPExamples.literate_all(force = true)


text:
# [src/rkadapt.jl]

"""
Runge-Kutta adaptive-step integration using the specified weights, nodes,
and Runge-Kutta matrix (or Dormand-Prince 5(4) by default).

Implements numerical propagation of an ordinary differential equation
from some initial value over the desired range. This function is similar
to MATLAB's variable-step ODE propagators (e.g., ode45), but is
generalized to any Butcher tableau. Unlike ode45, however, there's no
inherent minimum number of steps that this function will take; if it can
finish integration with a single step, it will, whereas ode45 will always
take at least a certain number of steps, no matter how small the time
interval and how smooth the dynamics. Therefore, this function is
commonly more useful with odehybrid than ode45.

This function is generic for all adaptive-step Runge-Kutta methods. That
is, any adaptive-step Runge-Kutta propagator can be created by passing
the weightes, nodes, and Runge-Kutta matrix (together, the Butcher
tableau) into this function. It will correctly take advantage of the
first-same-as-last property, eliminating one function evaluation per
step when the table has this property. See the examples below.

    t, x = rk4adapt(ode, ts, x0);
    t, x = rk4adapt(ode, ts, x0, options);
    t, x = rk4adapt(ode, ts, x0, a, b, c p);
    t, x = rk4adapt(ode, ts, x0, options, a, b, c, p);

Inputs:

    ode     Ordinary differential equation function (s.t. ẋ = ode(t, x);)
    ts      Time span, [t_start, t_end]
    x0      Initial state (column vector)
    options Options structure from odeset. This function uses the
                InitialStep, MaxStep, OutputFcn, RelTol, and AbsTol fields only.
    a       Runge-Kutta matrix (s-by-s)
    b       Weights (2-by-s), the top row containing the weights for the
                state update
    c       Nodes (s-by-1)
    p       Lowest order method represented by b (e.g., for
                Runge-Kutta-Fehlberg, this is 4).

    (NOTE that if a/b/c/p are not specified, this defaults to Dormand-Prince)

Outputs:

    t       Time history
    x       State history, with each row containing the state corresponding to
                the time in the same row of t.

Example:

    # Simulate an undamped harmonic oscillator for 10s starting from an
    #    initial state of [1; 0] using the default a, b, c, and p
    #    (Dormand-Prince 5(4)):

    t, x = rkadapt((t,x) -> [-x[2]; x[1]], [0 10], [1; 0]);
    plot(t, x);

    # Add options.

    options = ODE_SET(MaxStep = 0.1, InitialStep = 0.05)
    t, x = rkadapt( (t, x) -> [-x[2]; x[1]], [0 10], [1; 0], options)
    plot(t, x);

    # Now simulate with Bogacki-Shampine 3(2).

    a = [  0    0    0    0;
           1/2  0    0    0;
           0    3/4  0    0;
           2/9  1/3  4/9  0];

    b = [ 2/9   1/3  4/9  0;       # For the update
          7/24  1/4  1/3  1/8];    # For the error prediction

    c = [ 0     1/2  3/4  1];

    p = 2;                         # The lower of the two orders

    t, x = rkadapt( (t,x) -> [-x[2]; x[1]], [0 10], [1; 0], a, b, c, p);
    plot(t, x);

See "Runge-Kutta methods" on Wikipedia for discussion of the Butcher
tableau (a, b, and c).

Reference:

Dormand, J. R., and P. J. Prince. "A family of embedded Ruge-Kutta
formulae." _Journal of Computational and Applied Mathematics_ 6.1 (1980):
19-26. Web. 22 Apr. 2014.

See also: odehybrid, ode45, odeset, rkfixed.

Online doc: http://www.anuncommonlab.com/doc/odehybrid/rkadapt.html

Copyright 2014 An Uncommon Lab
"""
function rkadapt(ode, ts, x0, options::ODE_SET, a, b, c, p)

    # Unpack what we care about from the options.
    rel_tol = (options.RelTol == []) ? 1e-3 : options.RelTol
    abs_tol = (options.AbsTol == []) ? 1e-6 : options.AbsTol
    dt_max  = (options.MaxStep == []) ? Inf : options.MaxStep
    dt      = (options.InitialStep == []) ? min(1e-3, dt_max) : min(options.InitialStep, dt_max)

    # See if this table has the first-same-as-last property.
    fsal = all(a[end, :] == b[1, :]);

    ex = 1/(p+1);     # Time step correction exponent
    dt_min = 100 * eps(1.0 * maximum(abs.(ts))); # Minimum time step

    t   = ts[1];       # Time history
    x   = x0[:];       # State history
    tk  = ts[1];       # Current time
    xk  = x0[:];       # Current state
    xkp = xk;          # Proposed state update

    nx = length(x0);    # Number of states
    s  = size(b, 2);    # Length of weights
    d  = zeros(nx, s);  # Matrix of derivatives
    db = b[2, :] - b[1, :];

    # If the user provided an OutputFcn, call it.
    if (options.OutputFcn != [])
        options.OutputFcn(ts, x0, "init");
    end

    # If using a table with first-same-as-last, start the derivative here.
    if fsal
        d[:, 1] = ode(tk, xk);
    end

    # Propagate, updating k to k+1 until we're done until t(k) == ts(2).
    stop = false;
    while tk < (ts[2] - dt_min)

        # Don't step past the end time.
        if tk + dt > ts[2]
            dt = ts[2] - tk;
        end

        # Evaluate the first bit of the derivatives.
        if ~fsal
            d[:, 1] = ode(tk, xk);
        end

        tkp, xkp = [], []  # For scope

        # Try to take a step.
        while true

            # Calculate derivatives.
            for z = 2:s
                dxk = (d[:, 1:z-1])* (dt * a[z, 1:z-1])
                d[:, z] = ode(tk + c[z] * dt, xk + dxk);
            end

            # Update the state and error estimate.
            tkp = tk + dt;
            xkp = xk;
            er  = zeros(nx, 1);
            for z = 1:s
                xkp = xkp + ((dt * b[1, z]) .* d[:, z]);
                er  = er[:, 1] + ((dt * db[z]) .* d[:, z]);
            end

            # Create the scaled error. This single number can now be compared to 1 to see
            # if the worst error satisfies both absolute and relative tolerances. The error is acceptable if
            # it is either less than the relative error or less than the aboslute error.
            scaled_error = maximum(abs.(er) ./ (abs_tol .+ rel_tol .* abs.(xkp)));

            # If the error is small, update the time step and break to the
            # next step.
            if scaled_error <= 1
                dt = min(0.9 * dt * scaled_error^(-ex), min(3.0*dt, dt_max));
                # NOTE - ^This line has a SLIGHTLY different value than the MATLAB equivalent
                break;

            # If the error is too large, update the time step and continue
            # the search for an appropriate time step.
            else

                # If that was already the smallest step we can take, give up.
                if dt == dt_min
                    @warn "Could not meet integration tolerance."
                    break;
                end

                # Cut the time step down faster than we build it up.
                dt = max(0.7 * dt * scaled_error^(-ex), dt_min);
                     # NOTE - Julia and MATLAB slightly differ on last digit
            end

        end

        # Accept the update.
        tk = tkp;
        xk = xkp;

        # For first-same-as-last, we can copy over the derivative at k.
        if fsal
            d[:, 1] = d[:, end];
        end

        # Store.
        t = [t; tk];
        x = [x xk[:]];

        # If the user provided an OutputFcn, call it.
        if (options.OutputFcn != [])
            if (options.OutputFcn(tk, xk, ""))
                break
            end
        end
    end

    # Rectify the ending.
    if ~stop
        t[end] = ts[2];
    end

    # If the user provided an OutputFcn, call it.

    if (options.OutputFcn != [])
        options.OutputFcn([], [], "done")
    end

    return t, x'
end

function dptable()
    """
        Returns default Dormand-Prince Butcher Tableau values
    """

    a = [0           0          0           0        0          0     0;
         1/5         0          0           0        0          0     0;
         3/40        9/40       0           0        0          0     0;
         44/45      -56/15      32/9        0        0          0     0;
         19372/6561 -25360/2187 64448/6561 -212/729  0          0     0;
         9017/3168  -355/33     46732/5247  49/176  -5103/18656 0     0;
         35/384      0          500/1113    125/192 -2187/6784  11/84 0];

    b = [35/384     0  500/1113   125/192  -2187/6784    11/84    0;
         5179/57600 0  7571/16695 393/640  -92097/339200 187/2100 1/40];

    c = [0 1/5 3/10 4/5 8/9 1 1];

    p = 4;

    return a, b, c, p
end

function rkadapt(ode, ts, x0)
    """
        Wrapper function that generates default ODE_SET values and uses the
            Dormand-Prince Butcher Tableau
    """
    options = ODE_SET();
    a, b, c, p = dptable();

    return rkadapt(ode, ts, x0, options, a, b, c, p)
end

function rkadapt(ode, ts, x0, options::ODE_SET)
    """
        Wrapper function that uses the Dormand-Prince Butcher Tableau
    """

    a, b, c, p = dptable();
    return rkadapt(ode, ts, x0, options, a, b, c, p)
end

function rkadapt(ode, ts, x0, a, b, c, p)
    """
        Wrapper function that generates default ODE_SET values
    """
    options = ODE_SET();
    return rkadapt(ode, ts, x0, options, a, b, c, p)
end


text:
######### Stan program example  ###########

using StanSample, Test

eightschools ="
data {
  int<lower=0> J; // number of schools
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
  real mu;
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}
"

schools8data = Dict("J" => 8,
  "y" => [28,  8, -3,  7, -1,  1, 18, 12],
  "sigma" => [15, 10, 16, 11,  9, 11, 10, 18],
  "tau" => 25
)

sm = SampleModel("schools8", eightschools)

rc = stan_sample(sm, data=schools8data)

if success(rc)
  samples = read_samples(sm)

  df = read_summary(sm)
  @test df[df.parameters .== :mu, :mean][1] ≈ 7.64 rtol=0.5
  @test df[df.parameters .== :tau, :mean][1] ≈ 6.50 rtol=0.5
  @test df[df.parameters .== Symbol("theta[1]"), :mean][1] ≈ 11.1 rtol=0.5
  @test df[df.parameters .== Symbol("theta[8]"), :mean][1] ≈ 8.3 rtol=0.5
  @test df[df.parameters .== Symbol("eta[1]"), :mean][1] ≈ 0.4 rtol=0.5
  @test df[df.parameters .== Symbol("eta[8]"), :mean][1] ≈ 0.066 rtol=0.5

end


text:
"""
    MeshUtilModule

Module for mesh utility functions used in other meshing modules.
"""
module MeshUtilModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import Base.BitSet
import LinearAlgebra: norm

mutable struct HyperFaceContainer
    o::BitSet # numbers of the other nodes on the hyperface
    n::Int # new node number generated on the hyperface
end

makecontainer() = Dict{FInt, Array{HyperFaceContainer}}();

function addhyperface!(container,hyperface,newn)
    h=sort([i for i in hyperface])
    anchor=h[1]; other= BitSet(h[2:end]);
    C=get(container,anchor,HyperFaceContainer[]);
    fnd=false;
    for k=1:length(C)
        if C[k].o == other
            fnd=true; break;
        end
    end
    if fnd!=true
        push!(C,HyperFaceContainer(other,newn)); newn=newn+1;#new node added
    end
    container[anchor]= C; # set the new array
    return newn # return the number of the added node, possibly unchanged
end

function findhyperface!(container,hyperface)
    h=sort([i for i in hyperface])
    anchor=h[1]; other= BitSet(h[2:end]);
    C=get(container,anchor,HyperFaceContainer[]);
    for k=1:length(C)
        if C[k].o == other
            return [anchor, C[k].o], C[k].n
        end
    end
    return [], 0
end

function ontosphere(xyz::FFltMat,radius::FFlt)
# % Project nodes onto a sphere of given radius.
# %
# % function fens= onto_sphere(fens,radius,list)
# %
# % fens = finite element node set,
# % radius = radius of the sphere,
# % list = optional argument, if not empty then only
# %	     the nodes in the list are to be moved; otherwise all nodes are moved.

    for j=1:size(xyz,1)
        xyz[j,:] =xyz[j,:]*radius/norm(xyz[j,:]);
    end
    return xyz;
end

"""
    linearspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate linear space.

Generate a linear sequence of numbers between start and top (i. e. sequence
of number with uniform intervals inbetween).
"""
function linearspace(start::T, stop::T, N::Int)  where {T<:Number}
    return range(start, stop = stop, length = N)
end

"""
    gradedspace(start::T, stop::T, N::Int)  where {T<:Number}

Generate quadratic space.

Generate a quadratic sequence of numbers between start and finish.
This sequence corresponds to separation of adjacent numbers that
increases linearly from start to finish.
"""
function gradedspace(start::T, stop::T, N::Int, strength=2)  where {T<:Number}
    x = range(0.0, stop = 1.0, length = N);
    x = x.^strength
    # for i = 1:strength
    #     x = cumsum(x);
    # end
    x = x/maximum(x);
    out = start .* (1.0 .- x) .+ stop .* x;
end

# function inbox(box::FFltVec,sdim::FInt,x::FFltVec)
#     for i=1:sdim
#         if (!inrange(box[2*i-1],box[2*i],x[i]))
#             return false
#         end
#     end
#     return true
# end

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    ModiaMath - Mathematical Utilities for Modia and Modia3D

To define a model use Modia or Modia3D. You can define a model
natively in ModiaMath in the following way:

```julia
  using ModiaMath
  using StaticArrays

  @component Pendulum(;L=1.0, m=1.0, d=0.1, g=9.81) begin
     phi = RealScalar(start=pi/2, unit="rad"    , fixed=true,               numericType=ModiaMath.XD_EXP)
     w   = RealScalar(start=0.0 , unit="rad/s"  , fixed=true, integral=phi, numericType=ModiaMath.XD_EXP)
     a   = RealScalar(            unit="rad/s^2",             integral=w  , numericType=ModiaMath.DER_XD_EXP)
     r   = RealSVector{2}(        unit="m"      ,                           numericType=ModiaMath.WC)
  end;

  function ModiaMath.computeVariables!(p::Pendulum, sim::ModiaMath.SimulationState)
     L = p.L; m = p.m; d = p.d; g = p.g; phi = p.phi.value; w = p.w.value

     p.a.value = (-m*g*L*sin(phi) - d*w) / (m*L^2)

     if ModiaMath.isStoreResult(sim)
        p.r.value = @SVector [L*sin(phi), -L*cos(phi)]
     end
  end;

  simulationModel = ModiaMath.SimulationModel(Pendulum(L=0.8, m=0.5, d=0.2), stopTime=5.0);
```

To simulate a model and plot results:

```julia
  result = ModiaMath.simulate!(simulationModel; log=true);
  ModiaMath.plot(result, [(:phi, :w) :a])
```

[PendulumPlot](PendulumPlot-url)

To run examples:
```julia
  include("$(ModiaMath.path)/examples/Simulate_Pendulum.jl")
  include("$(ModiaMath.path)/examples/Simulate_FreeBodyRotation.jl")
  include("$(ModiaMath.path)/examples/withoutMacros_withoutVariables/Simulate_SimpleStateEvents.jl")
  include("$(ModiaMath.path)/examples/withoutMacros_withoutVariables/Simulate_BouncingBall.jl")
```

To run tests:
```julia
  include("$(ModiaMath.path)/test/runtests.jl")
```

For more information, see (https://github.com/ModiaSim/ModiaMath.jl/blob/master/README.md)
"""
module ModiaMath

"""
    const path

Absolute path of package directory of ModiaMath
"""
const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Time = Float64   # Prepare for later Integer type of time
const Version = "0.5.1"
const Date = "2019-04-07"

println(" \nImporting ModiaMath Version $Version ($Date)")

# Exported symbols
export @component
export RealVariable, RealScalar, RealSVector, RealSVector3
export plot

# Types
"""
    type ModiaMath.AbstractSimulationModel

Struct that is used as simulation model (has field: simulationState)
"""
abstract type AbstractSimulationModel end

"""
    type ModiaMath.AbstractComponentWithVariables

Struct that contains ModiaMath.AbstractVariables as field or as field
in a sub-struct.
"""
abstract type AbstractComponentWithVariables end

" The internal part of a component (has at least fields \"name\" and \"within\") "
abstract type AbstractComponentInternal end

"""
    type ModiaMath.AbstractVariable <: ModiaMath.AbstractComponentWithVariables

A Variable used as element of the DAE model description and is
included in the result (if no residue)
"""
abstract type AbstractVariable <: AbstractComponentWithVariables end

"""
    ModiaMath.AbstractRealVariable <: ModiaMath.AbstractVariable

A real [`ModiaMath.AbstractVariable`](@ref) (either scalar or array)
"""
abstract type AbstractRealVariable <: AbstractVariable end

getVariableAndResidueValues(extraInfo::Any) = nothing    # Return a variable value and a residue value table of nonlinear solver (for error message)

# include sub-modules and make symbols available that have been exported in sub-modules
include("Utilities.jl")
using .Utilities

include("Logging.jl")
using .Logging

include(joinpath("Frames", "_module.jl"))
using .Frames

#include("SparseJacobian.jl")
#using .SparseJacobian

include(joinpath("Result", "_module.jl"))
using .Result

include(joinpath("NonlinearEquations", "_module.jl"))
using .NonlinearEquations
const NonlinearEquationsInfo   = NonlinearEquations.KINSOL.NonlinearEquationsInfo
const solveNonlinearEquations! = NonlinearEquations.KINSOL.solveNonlinearEquations!

include(joinpath("DAE", "_module.jl"))
using .DAE

include(joinpath("SimulationEngine", "_module.jl"))
using .SimulationEngine

include(joinpath("Variables"       , "_module.jl"))
using .Variables

include("ModiaToModiaMath.jl")
using .ModiaToModiaMath

# Import packages that are used in examples and tests
# (in order that there are no requirements on the environment
#  in which the examples and tests are executed).
import DataFrames
import StaticArrays
import Unitful
import LinearAlgebra
import Test

end # module


text:
struct ParticleSwarmOptimization <: QuerySynthesisOptimizer
    eps::Float64
    swarmsize::Int
    maxiter::Int
    minstep::Float64
    minfunc::Float64
    function ParticleSwarmOptimization(; eps=0.1, swarmsize=100, maxiter=100)
        check_epsilon(eps)
        new(eps, swarmsize, maxiter, 1e-8, 1e-8)
    end
end

function query_synthesis_optimize(f::Function, optimizer::ParticleSwarmOptimization, data::Array{T, 2}, labels::Vector{Symbol})::Array{T, 2} where T <: Real
    lb, ub = vec.(data_boundaries(data[:, labels .!= :Lout], optimizer.eps))
    x_opt, _ = pso(x -> vec(-f(x)), lb, ub;
                    swarmsize=optimizer.swarmsize,
                    maxiter=optimizer.maxiter,
                    minstep=optimizer.minstep,
                    minfunc=optimizer.minfunc)
    return reshape(x_opt, size(data, 1), 1)
end

"""
Vectorized version of the partial swarm optimization without constraints.
Original implementation: https://github.com/yuehhua/PSO.jl by Yueh-Hua Tu
"""
function pso(func::Function, lb::Vector, ub::Vector; args=(), kwargs=Dict(),
             swarmsize=100, ω=0.5, ϕp=0.5, ϕg=0.5, maxiter=100, minstep=1e-8, minfunc=1e-8, verbose=false)
    @assert length(ub) == length(lb)
    @assert all(ub .>= lb)

    function update_position!(x, p, fx, fp)
        i_update = (fx .< fp)
        p[:, i_update] .= copy(x[:, i_update])
        fp[i_update] .= fx[i_update]
    end

    obj = x -> func(x, args...; kwargs...)

    # Initialize the particle swarm
    vhigh = abs.(ub .- lb)
    vlow = -vhigh
    S = swarmsize
    D = length(lb)  # the number of dimensions each particle has

    x = lb .+ rand(D, S) .* (ub .- lb)  # particle positions
    v = vlow .+ rand(D, S) .* (vhigh .- vlow)  # particle velocities
    p = zeros(D, S)  # best particle positions

    fx = obj(x)  # current particle function values
    fp = ones(S) * Inf  # best particle function values

    g = copy(x[:, 1])  # best swarm position
    fg = Inf  # best swarm position starting value

    # Store particle's best position (if constraints are satisfied)
    update_position!(x, p, fx, fp)

    # Update swarm's best position
    i_min = argmin(fp)
    if fp[i_min] < fg
        g = copy(p[:, i_min])
        fg = fp[i_min]
    end

    # Iterate until termination criterion met
    it = 1
    while it <= maxiter
        rp = rand(D, S)
        rg = rand(D, S)

        # Update the particles' velocities and positions
        v = ω*v .+ ϕp * rp .*(p .- x) .+ ϕg * rg .*(g .- x)
        x += v
        # Correct for bound violations
        maskl = x .< lb
        masku = x .> ub
        x = x .* (.~(maskl .| masku)) .+ lb .* maskl .+ ub .* masku

        # Update objectives
        fx = obj(x)

        # Store particle's best position (if constraints are satisfied)
        update_position!(x, p, fx, fp)

        # Compare swarm's best position with global best position
        i_min = argmin(fp)
        if fp[i_min] < fg
            verbose && println("New best for swarm at iteration $(it): $(p[i_min, :]) $(fp[i_min])")

            p_min = copy(p[:, i_min])
            stepsize = √(sum((g .- p_min).^2))

            if abs.(fg .- fp[i_min]) <= minfunc
                verbose && println("Stopping search: Swarm best objective change less than $(minfunc)")
                return (g, fg, p, fp)
            end
            if stepsize <= minstep
                verbose && println("Stopping search: Swarm best position change less than $(minstep)")
                return (g, fg, p, fp)
            end

            g = copy(p_min)
            fg = fp[i_min]
        end

        verbose && println("Best after iteration $(it): $(g) $(fg)")
        it += 1
    end

    verbose && println("Stopping search: maximum iterations reached --> $(maxiter)")
    return (g, fg, p, fp)
end


text:
using Pkg

Pkg.activate(dirname(@__FILE__))
Pkg.instantiate()

using JuliaWebAPI

attacks(x1, y1, x2, y2) = (x1 == x2 || y1 == y2 || x1 - y1 == x2 - y2 || x1 + y1 == x2 + y2)

function queens(N, xs...)
    lxs = length(xs)
    c = N - lxs
    c == 0 && return reverse(xs)
    for i=1:N
        any(1:lxs) do j
            attacks(xs[j], N-j+1, i, c)
        end && continue
        v = queens(N, xs..., i)
        v !== nothing && return v
    end
    return nothing
end

webqueens(N)=queens(parse(Int, N))

process(
    JuliaWebAPI.create_responder([
        (webqueens, true)
    ], "tcp://127.0.0.1:9999", true, ""), async=true
)

const apiclnt = APIInvoker("tcp://127.0.0.1:9999")

run_http(apiclnt, 8888)

max_stars_repo_path: Chapter 14/server/server.jl

max_issues_repo_path: Chapter 14/server/server.jl

max_forks_repo_path: Chapter 14/server/server.jl


text:
export AbstractSolver, Solver
export train!

abstract AbstractSolver

# ==========================================================================
# Namespacing for convenient enumerations
module Solver

using ..SupervisedLearning.AbstractSolver

immutable GradientDescent <: AbstractSolver
end

immutable BFGS <: AbstractSolver
end

end


text:
# Verify solutions against matpower results
using KernelAbstractions
using Test
using ExaPF

import ExaPF: PowerSystem

const PS = PowerSystem

@testset "Power flow 9 bus case" begin
    datafile = joinpath(dirname(@__FILE__), "..", "..", "data", "case9.m")
    pf = PS.PowerNetwork(datafile)

    # test impedance matrix entries
    @test isapprox(real(pf.Ybus[1, 1]), 0.0)
    @test isapprox(imag(pf.Ybus[1, 1]), -17.3611111)

    S = [0.0000 + 0.0000im, 1.6300 + 0.0000im, 0.8500 + 0.0000im, 0.0000 + 0.0000im,
        -0.9000 - 0.3000im, 0.0000 + 0.0000im, -1.0000 - 0.3500im, 0.0000 + 0.0000im,
        -1.2500 - 0.5000im]

    @test isapprox(S, pf.sbus)

    polar = PolarForm(pf, CPU())
    x = ExaPF.initial(polar, State())

    cache = ExaPF.get(polar, ExaPF.PhysicalState())
    ExaPF.init_buffer!(polar, cache)
    jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())
    # solve power flow
    convergence = ExaPF.powerflow(polar, jx, cache, NewtonRaphson())
    ExaPF.get!(polar, State(), x, cache)

    x_sol = [0.16875136481876485, 0.0832709533581424,
             -0.04200385129447893, -0.07011446830092488,
             0.033608106889679565, 0.010848015284769322, 0.06630715934781146,
             -0.07592061900861094, 0.9870068781579537, 0.9754722045044448,
             1.003375449839184, 0.9856449067439345,
             0.996185273317625, 0.9576210937650547,
            ]

    @test convergence.n_iterations == 4
    @test isapprox(x_sol, x, atol=1e-7)
end

@testset "Power flow 14 bus case" begin
    datafile = joinpath(dirname(@__FILE__), "..", "..", "data", "case14.m")
    polar = PolarForm(datafile, CPU())

    cache = ExaPF.get(polar, ExaPF.PhysicalState())
    ExaPF.init_buffer!(polar, cache)
    jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())
    # solve power flow
    conv = ExaPF.powerflow(polar, jx, cache, NewtonRaphson())

    @test conv.n_iterations == 2
    @test isapprox(norm(cache.balance, Inf), 1.3158e-10, rtol=1e-4)

    vmag = cache.vmag
    vmag_matpower = [1.060000,1.045000,1.010000,1.017671,1.019514,1.070000,1.061520,1.090000,1.055932,1.050985,1.056907,1.055189,1.050382,1.035530]

    @test isapprox(vmag, vmag_matpower, rtol=1e-6)
end

@testset "Power flow 30 bus case" begin
    datafile = joinpath(dirname(@__FILE__), "..", "..", "data", "case30.m")
    polar = PolarForm(datafile, CPU())

    cache = ExaPF.get(polar, ExaPF.PhysicalState())
    ExaPF.init_buffer!(polar, cache)
    jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())
    # solve power flow
    conv = ExaPF.powerflow(polar, jx, cache, NewtonRaphson())

    @test conv.n_iterations == 3
    @test isapprox(norm(cache.balance, Inf), 9.56998e-10, rtol=1e-4)

    vmag_matpower = [1.000000,1.000000,0.983138,0.980093,0.982406,0.973184,0.967355,0.960624,0.980506,0.984404,0.980506,0.985468,1.000000,0.976677,0.980229,0.977396,0.976865,0.968440,0.965287,0.969166,0.993383,1.000000,1.000000,0.988566,0.990215,0.972194,1.000000,0.974715,0.979597,0.967883]

    @test isapprox(cache.vmag, vmag_matpower, rtol=1e-6)
end

@testset "Power flow 300 bus case" begin
    datafile = joinpath(dirname(@__FILE__), "..", "..", "data", "case300.m")
    polar = PolarForm(datafile, CPU())

    cache = ExaPF.get(polar, ExaPF.PhysicalState())
    ExaPF.init_buffer!(polar, cache)
    jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())
    # solve power flow
    conv = ExaPF.powerflow(polar, jx, cache, NewtonRaphson())

    @test conv.n_iterations == 5
    @test isapprox(norm(cache.balance, Inf), 0.0, atol=1e-10)

    vmag_matpower = [1.028420,1.035340,0.997099,1.030812,1.019109,1.031196,0.993408,1.015300,1.003386,1.020500,1.005657,0.997373,0.997674,0.999170,1.034391,1.031587,1.064906,0.981924,1.001000,0.975168,0.996270,1.050098,1.005656,1.023354,0.998577,0.975035,1.024565,1.041441,0.975688,1.001170,1.020158,1.020312,1.053611,1.021655,1.029283,1.044944,1.000732,1.008735,1.021646,1.034515,0.977877,1.001958,1.047499,1.025388,0.998003,0.996035,1.005135,1.015146,1.033490,0.991822,0.978860,1.024716,0.990654,1.016040,0.958300,0.947956,0.962698,0.951318,0.979391,0.969614,0.977610,0.996488,0.963200,0.983787,0.990023,0.982012,0.987242,1.034127,1.025000,0.987112,0.990818,0.991954,1.015248,1.031724,1.027231,1.052000,1.052000,0.992945,1.018224,1.000000,0.989358,1.005986,1.000708,1.028759,0.995737,1.022267,1.009461,0.990000,0.975245,0.973213,0.974473,0.970155,0.976812,0.960282,1.024861,0.934829,0.929853,1.043500,0.958437,0.987111,0.972796,1.000588,1.023300,1.010300,0.997795,1.000129,1.002406,1.002825,1.019136,0.986142,1.004551,1.001998,1.022076,1.019337,1.047586,1.047088,1.055000,1.011709,1.042991,1.051000,1.015510,1.043500,1.016107,1.008106,1.052800,1.052800,1.057719,1.073500,0.986926,1.004833,1.053500,1.043500,0.966417,1.017724,0.963000,0.984473,0.998709,0.986644,0.999801,1.036082,0.991820,1.041011,0.983914,1.000211,0.997254,0.971492,1.002431,0.987864,0.929000,0.982900,1.024466,0.983654,1.062214,0.973081,1.052200,1.007700,0.939796,0.969910,0.979330,1.051824,1.044628,0.971645,1.038589,1.052200,1.065000,1.065000,1.053282,1.002757,1.055100,1.043500,0.937458,0.998236,1.048984,1.035903,0.973993,0.992473,1.015000,0.954313,0.956174,0.974032,0.990839,1.003359,0.966709,0.985554,1.003768,1.018555,0.999440,1.004774,0.980462,1.001820,1.013262,1.010000,0.991863,0.986632,0.975110,1.021525,1.007547,1.055420,1.008000,1.000000,1.050000,0.996551,1.000254,0.945276,1.018005,1.000000,1.042356,1.049552,1.040000,1.053541,1.041466,1.000000,1.038706,1.009515,1.016500,1.055850,1.010000,1.000000,1.023776,1.050000,0.993000,1.010000,0.992178,0.971140,0.965191,0.969095,0.976999,0.976227,1.020532,1.025125,1.015209,1.014590,1.000433,0.980890,0.974945,0.942873,0.972387,0.960470,1.000921,0.977728,0.958325,1.031028,1.012876,1.024438,1.012197,0.969485,1.050700,1.050700,1.032300,1.014500,1.050700,1.050700,1.050700,1.029000,1.050000,1.014500,1.050700,0.996700,1.021200,1.014500,1.001700,0.989300,1.050700,1.050700,1.014500,1.011774,0.994500,0.983335,0.976825,1.011711,1.002924,0.991387,1.002280,0.988722,0.964884,0.974704,0.970504,0.964756,0.965606,0.931742,0.944074,0.928799,0.997240,0.950422,0.959699,0.957027,0.939160,0.963555,0.950267,0.964683,0.979007,1.000000,0.978627,1.000000,1.000000,1.000000,0.975431,0.980460,0.979888,1.040517]

    @test isapprox(cache.vmag, vmag_matpower, rtol=1e-6)
end


text:
module KhepriRadiance
using KhepriBase

# functions that need specialization
include(khepribase_interface_file())
include("Radiance.jl")

function __init__()
  set_material(radiance, material_basic, radiance_material_neutral)
  set_material(radiance, material_metal, radiance_generic_metal)
  set_material(radiance, material_glass, radiance_generic_glass_80)
  set_material(radiance, material_wood, radiance_light_wood)
  set_material(radiance, material_concrete, radiance_outside_facade_30)
  set_material(radiance, material_plaster, radiance_generic_interior_wall_70)

  add_current_backend(radiance)
end
end


text:
function encrypt(msg::AbstractString, key::AbstractString)
    msg    = uppercase(join(filter(isalpha, collect(msg))))
    key    = uppercase(join(filter(isalpha, collect(key))))
    msglen = length(msg)
    keylen = length(key)

    if keylen < msglen
        key = repeat(key, div(msglen - keylen, keylen) + 2)[1:msglen]
    end

    enc = Vector{Char}(msglen)

    @inbounds for i in 1:length(msg)
        enc[i] = Char((Int(msg[i]) + Int(key[i]) - 130) % 26 + 65)
    end

    return join(enc)
end

function decrypt(enc::AbstractString, key::AbstractString)
    enc    = uppercase(join(filter(isalpha, collect(enc))))
    key    = uppercase(join(filter(isalpha, collect(key))))
    msglen = length(enc)
    keylen = length(key)

    if keylen < msglen
        key = repeat(key, div(msglen - keylen, keylen) + 2)[1:msglen]
    end

    msg = Vector{Char}(msglen)

    @inbounds for i in 1:length(enc)
        msg[i] = Char((Int(enc[i]) - Int(key[i]) + 26) % 26 + 65)
    end

    return join(msg)
end

const messages = ("Attack at dawn.", "Don't attack.", "The war is over.", "Happy new year!")
const key = "LEMON"

for msg in messages
    enc = encrypt(msg, key)
    dec = decrypt(enc, key)
    println("Original: $msg\n -> encrypted: $enc\n -> decrypted: $dec")
end

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia