text:
Base.lock(h::Hist3D) = lock(h.hlock)
Base.unlock(h::Hist3D) = unlock(h.hlock)

"""
    bincounts(h::Hist3D)

Get the bin counts of a histogram.
"""
@inline bincounts(h::Hist3D) = h.hist.weights

"""
    binedges(h::Hist3D)

Get a 3-tuple of the bin edges of a histogram.
"""
@inline binedges(h::Hist3D) = h.hist.edges

"""
    bincenters(h::Hist3D)

Get a 2-tuple of the bin centers of a histogram.
"""
function bincenters(h::Hist3D)
    StatsBase.midpoints.(binedges(h))
end

"""
    binerrors(f::T, h::Hist3D) where T<:Function = f.(h.sumw2)
    binerrors(h::Hist3D) = binerrors(sqrt, h)

Calculate the bin errors from `sumw2` with a Gaussian default.
"""
binerrors(f::T, h::Hist3D) where T<:Function = f.(h.sumw2)
binerrors(h::Hist3D) = binerrors(sqrt, h)

"""
    nbins(h::Hist3D)

Get a 3-tuple of the number of x and y bins of a histogram.
"""
function nbins(h::Hist3D)
    size(bincounts(h))
end

"""
    integral(h::Hist3D)

Get the integral a histogram.
"""
function integral(h::Hist3D)
    sum(bincounts(h))
end

"""
    empty!(h::Hist3D)

Resets a histogram's bin counts and `sumw2`.
"""
function Base.empty!(h::Hist3D{T,E}) where {T,E}
    h.hist.weights .= zero(T)
    h.sumw2 .= 0.0
    return h
end

"""
    push!(h::Hist3D, valx::Real, valy::Real, wgt::Real=1)
    atomic_push!(h::Hist3D, valx::Real, valy::Real, wgt::Real=1)

Adding one value at a time into histogram.
`sumw2` (sum of weights^2) accumulates `wgt^2` with a default weight of 1.
`atomic_push!` is a slower version of `push!` that is thread-safe.

"""
@inline function atomic_push!(h::Hist3D{T,E}, valx::Real, valy::Real, valz::Real, wgt::Real=1) where {T,E}
    lock(h)
    push!(h, valx, valy, valz, wgt)
    unlock(h)
    return nothing
end

@inline function Base.push!(h::Hist3D{T,E}, valx::Real, valy::Real, valz::Real, wgt::Real=1) where {T,E}
    rx, ry, rz = binedges(h)
    Lx, Ly, Lz = nbins(h)
    binidxx = _edge_binindex(rx, valx)
    binidxy = _edge_binindex(ry, valy)
    binidxz = _edge_binindex(rz, valz)
    if h.overflow
        binidxx = clamp(binidxx, 1, Lx)
        binidxy = clamp(binidxy, 1, Ly)
        binidxz = clamp(binidxz, 1, Lz)
        @inbounds h.hist.weights[binidxx,binidxy,binidxz] += wgt
        @inbounds h.sumw2[binidxx,binidxy,binidxz] += wgt^2
    else
        if (unsigned(binidxx - 1) < Lx) && (unsigned(binidxy - 1) < Ly) && (unsigned(binidxz - 1) < Lz)
            @inbounds h.hist.weights[binidxx,binidxy,binidxz] += wgt
            @inbounds h.sumw2[binidxx,binidxy,binidxz] += wgt^2
        end
    end
    return nothing
end

Base.broadcastable(h::Hist3D) = Ref(h)

"""
    Hist3D(elT::Type{T}=Float64; binedges, overflow) where {T}

Initialize an empty histogram with bin content typed as `T` and bin edges.
To be used with [`push!`](@ref). Default overflow behavior (`false`)
will exclude values that are outside of `binedges`.
"""
function Hist3D(elT::Type{T}=Float64; bins, overflow=_default_overflow) where {T}
    counts = zeros(elT, length.(bins) .- 1)
    return Hist3D(Histogram(bins, counts); overflow=overflow)
end

"""
    Hist3D(tuple, edges::NTuple{3,AbstractRange}; overflow)
    Hist3D(tuple, edges::NTuple{3,AbstractVector}; overflow)

Create a `Hist3D` with given bin `edges` and values from
a 2-tuple of arrays of x, y values. Weight for each value is assumed to be 1.
"""
function Hist3D(A::NTuple{3,AbstractVector}, r::NTuple{3,AbstractRange}; overflow=_default_overflow)
    h = Hist3D(Int; bins=r, overflow=overflow)
    push!.(h, A[1], A[2], A[3])
    return h
end
function Hist3D(A::NTuple{3,AbstractVector}, edges::NTuple{3,AbstractVector}; overflow=_default_overflow)
    if all(_is_uniform_bins.(edges))
        r = (range(first(edges[1]), last(edges[1]), length=length(edges[1])),
             range(first(edges[2]), last(edges[2]), length=length(edges[2])),
             range(first(edges[3]), last(edges[3]), length=length(edges[3])))
        return Hist3D(A, r; overflow=overflow)
    else
        h = Hist3D(Int; bins=edges, overflow=overflow)
        push!.(h, A[1], A[2], A[3])
        return h
    end
end

"""
    Hist3D(tuple, wgts::AbstractWeights, edges::NTuple{3,AbstractRange}; overflow)
    Hist3D(tuple, wgts::AbstractWeights, edges::NTuple{3,AbstractVector}; overflow)

Create a `Hist3D` with given bin `edges` and values from
a 2-tuple of arrays of x, y values.
`wgts` should have the same `size` as elements of `tuple`.
"""
function Hist3D(A::NTuple{3,AbstractVector}, wgts::AbstractWeights, r::NTuple{3,AbstractRange}; overflow=_default_overflow)
    @boundscheck @assert size(A[1]) == size(A[2]) == size(A[3]) == size(wgts)
    h = Hist3D(eltype(wgts); bins=r, overflow=overflow)
    push!.(h, A[1], A[2], A[3], wgts)
    return h
end
function Hist3D(A::NTuple{3,AbstractVector}, wgts::AbstractWeights, edges::NTuple{3,AbstractVector}; overflow=_default_overflow)
    if all(_is_uniform_bins.(edges))
        r = (range(first(edges[1]), last(edges[1]), length=length(edges[1])),
             range(first(edges[2]), last(edges[2]), length=length(edges[2])),
             range(first(edges[3]), last(edges[3]), length=length(edges[3])))
        return Hist3D(A, wgts, r; overflow=overflow)
    else
        h = Hist3D(Int; bins=edges, overflow=overflow)
        push!.(h, A[1], A[2], A[3], wgts)
        return h
    end
end

"""
    Hist3D(A::AbstractVector{T}; nbins::NTuple{3,Integer}, overflow) where T
    Hist3D(A::AbstractVector{T}, wgts::AbstractWeights; nbins::NTuple{3,Integer}, overflow) where T

Automatically determine number of bins based on `Sturges` algo.
"""
function Hist3D(A::NTuple{3,AbstractVector{T}};
        nbins::NTuple{3,Integer}=_sturges.(A),
        overflow=_default_overflow,
    ) where {T}
    F = float(T)
    nbinsx, nbinsy, nbinsz = nbins
    lox, hix = minimum(A[1]), maximum(A[1])
    loy, hiy = minimum(A[2]), maximum(A[2])
    loz, hiz = minimum(A[3]), maximum(A[3])
    rx = StatsBase.histrange(F(lox), F(hix), nbinsx)
    ry = StatsBase.histrange(F(loy), F(hiy), nbinsy)
    rz = StatsBase.histrange(F(loz), F(hiz), nbinsz)
    r = (rx, ry, rz)
    return Hist3D(A, r; overflow=overflow)
end

function Hist3D(A::NTuple{3,AbstractVector{T}}, wgts::AbstractWeights;
        nbins::NTuple{3,Integer}=_sturges.(A),
        overflow=_default_overflow,
    ) where {T}
    F = float(T)
    nbinsx, nbinsy, nbinsz = nbins
    lox, hix = minimum(A[1]), maximum(A[1])
    loy, hiy = minimum(A[2]), maximum(A[2])
    loz, hiz = minimum(A[3]), maximum(A[3])
    rx = StatsBase.histrange(F(lox), F(hix), nbinsx)
    ry = StatsBase.histrange(F(loy), F(hiy), nbinsy)
    rz = StatsBase.histrange(F(loz), F(hiz), nbinsz)
    r = (rx, ry, rz)
    return Hist3D(A, wgts, r; overflow=overflow)
end

"""
    function lookup(h::Hist3D, x, y)

For given x/y/z-axis value `x`, `y`, `z`, find the corresponding bin and return the bin content.
If a value is out of the histogram range, return `missing`.
"""
function lookup(h::Hist3D, x, y, z)
    rx, ry, rz = binedges(h)
    !(first(rx) <= x <= last(rx)) && return missing
    !(first(ry) <= y <= last(ry)) && return missing
    !(first(rz) <= z <= last(rz)) && return missing
    return bincounts(h)[_edge_binindex(rx, x), _edge_binindex(ry, y), _edge_binindex(rz, z)]
end

"""
    normalize(h::Hist3D)

Create a normalized histogram via division by `integral(h)`.
"""
function normalize(h::Hist3D)
    return h*(1/integral(h))
end

"""
    project(h::Hist3D, axis::Symbol=:x)
    project(axis::Symbol=:x) = h::Hist3D -> project(h, axis)

Computes the `:x`/`:y`/`:z` axis projection of the 3D histogram by
summing over the specified axis. Returns a `Hist2D`.
"""
function project(h::Hist3D, axis::Symbol=:x)
    @assert axis ∈ (:x, :y, :z)
    dimremove, dimskeep = if axis == :z
        3, [1,2]
    elseif axis == :y
        2, [1,3]
    else
        1, [2,3]
    end
    counts = dropdims(sum(bincounts(h), dims=dimremove), dims=dimremove)
    sumw2 = dropdims(sum(h.sumw2, dims=dimremove), dims=dimremove)
    edges = binedges(h)[dimskeep]
    return Hist2D(Histogram(edges, counts), sumw2; overflow=h.overflow)
end
project(axis::Symbol=:x) = h::Union{Hist2D,Hist3D} -> project(h, axis)


text:
module Recordings

import SignalAnalysis: signal, framerate, nchannels
import Statistics: mean
import DataFrames: DataFrame, DataFrameRow
import TimeZones: localzone, ZonedDateTime, astimezone, @tz_str
import Dates: unix2datetime, datetime2unix, now, UTC

const MAGIC = hton(0x43c04d126f173001)

# defaults from old file format without header
const FRAMERATE = 32000.0
const NCHANNELS = 4

function read(filename; tz=localzone())
  isfile(filename) && return readrec(filename)
  df = DataFrame(time=ZonedDateTime[], filename=String[], duration=Float64[], nchannels=Int[], framerate=Float64[])
  for filename ∈ readdir(filename; join=true)
    m = match(r"/?rec\-(\d+)\.dat", filename)
    if m !== nothing
      t = 0
      filesize = stat(filename).size
      nch = NCHANNELS
      fs = FRAMERATE
      if filesize > 24
        open(filename, "r") do f
          hdr = read(f, RecordingHeader)
          if hdr !== nothing
            t = astimezone(ZonedDateTime(unix2datetime(hdr.millis / 1000), tz"UTC"), tz)
            nch = hdr.nchannels
            fs = hdr.framerate
          end
        end
      end
      d = filesize / (4 * nch * fs)
      t == 0 && (t = astimezone(ZonedDateTime(unix2datetime(parse(Int64, m[1]) / 1000), tz"UTC"), tz))
      push!(df, (t, filename, d, nch, Float64(fs)))
    end
  end
  sort!(df, :time)
end

read(filenames::AbstractVector; tz=localzone()) = sort!(vcat(read.(filenames; tz)...), :time)

read(row::DataFrameRow) = readrec(row.filename)
read(df::DataFrame, i) = readrec(df.filename[i])

Base.@kwdef struct RecordingHeader
  magic::UInt64 = MAGIC
  millis::Int64 = round(Int64, 1000 * time())
  framerate::Int32
  nchannels::Int16 = 1
end

function read(io::IO, ::Type{RecordingHeader})
  p = position(io)
  try
    x = Base.read(io, UInt64)
    if x == MAGIC
      hdr = RecordingHeader(
        millis = Base.read(io, Int64),
        framerate = Base.read(io, Int32),
        nchannels = Base.read(io, Int16))
      [Base.read(io, UInt16) for _ in 1:5]
      return hdr
    end
  catch ex
    ex isa EOFError || @warn "$ex"
  end
  seek(io, p)
  nothing
end

function readrec(filename)
  open(filename, "r") do f
    nch = NCHANNELS
    fs = FRAMERATE
    hdr = read(f, RecordingHeader)
    if hdr !== nothing
      nch = hdr.nchannels
      fs = hdr.framerate
    end
    bytes = Base.read(f)
    raw = reinterpret(Float32, bytes)
    x = signal(reshape(raw, nch, :)', fs)
    x .- mean(x; dims=1)
  end
end

function write(filename, x; fs=framerate(x))
  open(filename, "w") do io
    hdr = RecordingHeader(
      millis = round(Int64, 1000 * datetime2unix(now(UTC))),
      framerate = round(Int32, fs),
      nchannels = nchannels(x))
    Base.write(io, hdr.magic)
    Base.write(io, hdr.millis)
    Base.write(io, hdr.framerate)
    Base.write(io, hdr.nchannels)
    Base.write(io, UInt16[0, 0, 0, 0, 0])
    Base.write(io, Float32.(x'))
  end
  nothing
end

end # module

ext: jl


text:
__precompile__()

"""
This module includes all the base functions of OpenMendel.
"""
module MendelBase
#
# Required external packages.
#
using CSV
using DataFrames
using Distributions
using GLM
using LinearAlgebra
using Missings
using Printf
using Random
using Statistics
using StatsBase
#
# Required OpenMendel packages and modules.
#
using MendelSearch
using SnpArrays
#
# Define the data structures used by OpenMendel.
#
include("data_structures.jl")
#
# Include functions to process the keywords that specify
# the data files to use and the analysis to perform.
#
include("keywords.jl")
#
# Include useful general utilities.
#
include("general_utilities.jl")
#
# Include useful genetic utilities.
#
include("genetic_utilities.jl")
#
# Include functions to read the genetic data from external files.
#
include("read_data.jl")
#
# Include functions to prepare for and perform a pedigree likelihood evaluation
# via the Elston-Stewart algorithm.
#
include("elston_stewart_preparation.jl")
include("elston_stewart_evaluation.jl")
#
# Method to obtain path to this package's data files
# so they can be used in the documentation and testing routines.
# For example, datadir("Control file.txt") will return
# "/path/to/package/data/Control file.txt"
#
datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)

end # module MendelBase

max_stars_repo_name: OpenMendel/MendelBase

max_issues_repo_name: OpenMendel/MendelBase

max_forks_repo_name: OpenMendel/MendelBase


text:
#####
##### Viscous flux divergences
#####

@inline function ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, clock, closure::AbstractTurbulenceClosure, args...)
    return 1/Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_cᶜᶜᶜ, viscous_flux_ux, clock, closure, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_ζᶠᶠᶜ, viscous_flux_uy, clock, closure, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_ηᶠᶜᵃ, viscous_flux_uz, clock, closure, args...))
end

@inline function ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, clock, closure::AbstractTurbulenceClosure, args...)
    return 1/Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_ζᶠᶠᶜ, viscous_flux_vx, clock, closure, args...) +
                                    δyᵃᶠᵃ(i, j, k, grid, Ay_cᶜᶜᶜ, viscous_flux_vy, clock, closure, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_ξᶜᶠᵃ, viscous_flux_vz, clock, closure, args...))
end

@inline function ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, clock, closure::AbstractTurbulenceClosure, args...)
    return 1/Vᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_ηᶠᶜᶠ, viscous_flux_wx, clock, closure, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_ξᶜᶠᶠ, viscous_flux_wy, clock, closure, args...) +
                                    δzᵃᵃᶠ(i, j, k, grid, Az_cᶜᶜᵃ, viscous_flux_wz, clock, closure, args...))
end

#####
##### Viscosities at different locations
#####

@inline νᶜᶜᶜ(i, j, k, grid, clock, ν::Number) = ν
@inline νᶠᶠᶜ(i, j, k, grid, clock, ν::Number) = ν
@inline νᶠᶜᶠ(i, j, k, grid, clock, ν::Number) = ν
@inline νᶜᶠᶠ(i, j, k, grid, clock, ν::Number) = ν

@inline νᶜᶜᶜ(i, j, k, grid, clock, ν::Function) = ν(xnode(Center, i, grid), ynode(Center, j, grid), znode(Center, k, grid), clock.time)
@inline νᶠᶠᶜ(i, j, k, grid, clock, ν::Function) = ν(xnode(Face,   i, grid), ynode(Face,   j, grid), znode(Center, k, grid), clock.time)
@inline νᶠᶜᶠ(i, j, k, grid, clock, ν::Function) = ν(xnode(Face,   i, grid), ynode(Center, j, grid), znode(Face,   k, grid), clock.time)
@inline νᶜᶠᶠ(i, j, k, grid, clock, ν::Function) = ν(xnode(Center, i, grid), ynode(Face,   j, grid), znode(Face,   k, grid), clock.time)

#####
##### Viscous fluxes
#####

@inline viscous_flux_ux(i, j, k, grid, clock, ν, u) = νᶜᶜᶜ(i, j, k, grid, clock, ν) * ∂xᶜᶜᵃ(i, j, k, grid, u)
@inline viscous_flux_uy(i, j, k, grid, clock, ν, u) = νᶠᶠᶜ(i, j, k, grid, clock, ν) * ∂yᶠᶠᵃ(i, j, k, grid, u)
@inline viscous_flux_uz(i, j, k, grid, clock, ν, u) = νᶠᶜᶠ(i, j, k, grid, clock, ν) * ∂zᵃᵃᶠ(i, j, k, grid, u)

@inline viscous_flux_vx(i, j, k, grid, clock, ν, v) = νᶠᶠᶜ(i, j, k, grid, clock, ν) * ∂xᶠᶠᵃ(i, j, k, grid, v)
@inline viscous_flux_vy(i, j, k, grid, clock, ν, v) = νᶜᶜᶜ(i, j, k, grid, clock, ν) * ∂yᶜᶜᵃ(i, j, k, grid, v)
@inline viscous_flux_vz(i, j, k, grid, clock, ν, v) = νᶜᶠᶠ(i, j, k, grid, clock, ν) * ∂zᵃᵃᶠ(i, j, k, grid, v)

@inline viscous_flux_wx(i, j, k, grid, clock, ν, w) = νᶠᶜᶠ(i, j, k, grid, clock, ν) * ∂xᶠᶜᵃ(i, j, k, grid, w)
@inline viscous_flux_wy(i, j, k, grid, clock, ν, w) = νᶜᶠᶠ(i, j, k, grid, clock, ν) * ∂yᶜᶠᵃ(i, j, k, grid, w)
@inline viscous_flux_wz(i, j, k, grid, clock, ν, w) = νᶜᶜᶜ(i, j, k, grid, clock, ν) * ∂zᵃᵃᶜ(i, j, k, grid, w)

#####
##### Products of viscosity and divergence, vorticity, and vertical momentum gradients
#####

@inline ν_δᶜᶜᶜ(i, j, k, grid, clock, ν, u, v) = @inbounds νᶜᶜᶜ(i, j, k, grid, clock, ν) * div_xyᶜᶜᵃ(i, j, k, grid, u, v)
@inline ν_ζᶠᶠᶜ(i, j, k, grid, clock, ν, u, v) = @inbounds νᶠᶠᶜ(i, j, k, grid, clock, ν) * ζ₃ᶠᶠᵃ(i, j, k, grid, u, v)

@inline ν_uzᶠᶜᶠ(i, j, k, grid, clock, ν, u) = @inbounds νᶠᶜᶠ(i, j, k, grid, clock, ν) * ∂zᵃᵃᶠ(i, j, k, grid, u)
@inline ν_vzᶜᶠᶠ(i, j, k, grid, clock, ν, v) = @inbounds νᶜᶠᶠ(i, j, k, grid, clock, ν) * ∂zᵃᵃᶠ(i, j, k, grid, v)


text:
# Grid Function
@doc """
    CTP_Grid(CTPGrid::Array{Float64,2})

    Returns tmm - minimum measured thickness determined at the time of the inspection
"""
function CTP_Grid(CTPGrid::Array{Float64,2})  # may add plotting feature later - spacings paramter for here
    # Labels for plotting
    #=circumferential_plane_index = Int64.(collect(1:1:s2)) # collect auto rounds - no need to write code to correct
    circ_label = fill("C",length(circumferential_plane_index))
    circ_index_labels = map(string, circ_label, circumferential_plane_index) # join integer and string to string
    circ_spacing_index = collect(0:0.5:(length(circumferential_plane_index)/2)-spacings)
    longitudinal_plane_index = Int64.(collect(1:1:c2)) # collect auto rounds - no need to write code to correct
    long_label = fill("M",length(longitudinal_plane_index))
    long_index_labels = map(string, long_label, longitudinal_plane_index) # join integer and string to string
    long_spacing_index = collect(0:0.5:(length(longitudinal_plane_index)/2)-spacings)
    =#

    # Determine CTPs (Critical Thickness Profiles)
    Circumferential_CTP = zeros(size(CTPGrid,1)) # initialize output
    Longitudinal_CTP = zeros(size(CTPGrid,2)) # initialize output
    col_index = [1:size(CTPGrid,2);]
    row_index = [1:size(CTPGrid,1);]
    @inbounds for i in 1:size(CTPGrid,1)
        Circumferential_CTP[i] = minimum(CTPGrid[row_index[i],1:size(CTPGrid,2)])
    end
    i=2
    @inbounds for i in 1:size(CTPGrid,2)
        Longitudinal_CTP[i] = minimum(CTPGrid[1:size(CTPGrid,1),col_index[i]])
    end

    # find tmm
    minC_CTP = minimum(Circumferential_CTP)
    minL_CTP = minimum(Longitudinal_CTP)
    tmm = minimum([minC_CTP,minL_CTP]) # minimum measured thickness determined at the time of the inspection.
        return tmm, Longitudinal_CTP
end


text:
@article{OpenMPIData,
	Author = {{Knopp}, T. and {Szwargulski}, P. and {Griese}, F. and {Gr{\"a}ser}, M.},
	Title = {OpenMPIData: An initiative for freely accessible magnetic particle imaging data},
	Journal = {Data in Brief},
	Year = {2020},
	Volume = {28},
	Pages = {104971},
	Month = {2},
	Doi = {https://doi.org/10.1016/j.dib.2019.104971},
	Url = {https://www.sciencedirect.com/science/article/pii/S2352340919313265?via%3Dihub},
	Keywords = {MPI},
}


text:
# struct to model a element

@tm1def mutable struct Element
  Name::Union{String,Nothing}
  UniqueName::Union{String,Nothing}
  Index::Union{Integer,Nothing}
  # not sure what type to use here
  Type::Union{String,Nothing}
  # not sure what type to use here either so am just going to skip this field
  #Hierarchy::Union{String,Nothing}
  Edges::Union{Vector,Nothing}
  Level::Union{Integer,Nothing}
  Attributes::Union{Dict,Nothing}
  Parents::Union{Vector{Element},Nothing}
  Components::Union{Vector{Element},Nothing}
  LocalizedAttributes::Union{Vector{Dict},Nothing}
end

Element(name::AbstractString) = Element(Dict("Name" => name))

namefield(element::Element) = element.Name

# functions for endpoints

@api_default function get_element(
  api::TM1API,
  dimension_name::AbstractString,
  hierarchy_name::AbstractString,
  element_name::AbstractString;
  options...,
)
  params = Dict("\$expand" => "*")

  result = tm1_get_json(
    api,
    "Dimensions('" *
    dimension_name *
    "')/Hierarchies('" *
    hierarchy_name *
    "')/Elements('" *
    element_name *
    "')";
    params = params,
    options...,
  )

  Element(result)
end

@api_default function delete_element(
  api::TM1API,
  dimension_name::AbstractString,
  hierarchy_name::AbstractString,
  element_name::AbstractString;
  options...,
)

  tm1_delete(
    api,
    "Dimensions('" *
    dimension_name *
    "')/Hierarchies('" *
    hierarchy_name *
    "')/Elements('" *
    element_name *
    "')";
    options...,
  )
end


text:
ssr
##############################################################################
##
## Fit is the only exported function
##
##############################################################################

function fit(m::InteractiveFixedEffectModel,
             f::Formula,
             df::AbstractDataFrame,
             vcov_method::AbstractVcovMethod = VcovSimple();
             method::Symbol = :dogleg,
             lambda::Number = 0.0,
             subset::Union{AbstractVector{Bool}, Void} = nothing,
             weight::Union{Symbol, Void} = nothing,
             maxiter::Integer = 10_000,
             tol::Real = 1e-9,
             save::Union{Bool, Void} = nothing)

    ##############################################################################
    ##
    ## Transform DataFrame -> Matrix
    ##
    ##############################################################################

    ## parse formula
    rf = deepcopy(f)
    (has_absorb, absorb_formula, absorb_terms, has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose!(rf)
    if has_iv
        error("partial_out does not support instrumental variables")
    end
    rt = Terms(rf)
    has_regressors = allvars(rf.rhs) != [] || (rt.intercept == true && !has_absorb)

    # change default if has_regressors
    if save == nothing
        save = !has_regressors
    end
    ## create a dataframe without missing values & negative weights
    vars = allvars(rf)
    vcov_vars = allvars(vcov_method)
    absorb_vars = allvars(absorb_formula)
    factor_vars = [m.id, m.time]
    all_vars = vcat(vars, absorb_vars, factor_vars, vcov_vars)
    all_vars = unique(convert(Vector{Symbol}, all_vars))
    esample = complete_cases(df[all_vars])
    if weight != nothing
        esample &= isnaorneg(df[weight])
        all_vars = unique(vcat(all_vars, weight))
    end
    if subset != nothing
        if length(subset) != size(df, 1)
            error("df has $(size(df, 1)) rows but the subset vector has $(length(subset)) elements")
        end
        esample &= convert(BitArray, subset)
    end
    subdf = df[esample, all_vars]
    main_vars = unique(convert(Vector{Symbol}, vcat(vars, factor_vars)))
    for v in main_vars
        dropUnusedLevels!(subdf[v])
    end

    # Compute data needed for errors
    vcov_method_data = VcovMethodData(vcov_method, subdf)

    # Compute weight
    sqrtw = get_weight(subdf, weight)

    ## Compute factors, an array of AbtractFixedEffects
    if has_absorb
        fes = FixedEffect[FixedEffect(subdf, a, sqrtw) for a in absorb_terms.terms]
        # in case some FixedEffect is aFixedEffectIntercept, remove the intercept
        if any([typeof(f.interaction) <: Ones for f in fes])
            rt.intercept = false
        end
        pfe = FixedEffectProblem(fes, Val{:lsmr})
    else
        pfe = nothing
    end

    iterations = 0
    converged = false
    # get two dimensions
    id = subdf[m.id]
    time = subdf[m.time]

    ##############################################################################
    ##
    ## Construict vector y and matrix X
    ##
    ##############################################################################

    mf = simpleModelFrame(subdf, rt, esample)

    # Compute demeaned X
    if has_regressors
        coef_names = coefnames(mf)
        X = ModelMatrix(mf).m
        broadcast!(*, X, X, sqrtw)
        residualize!(X, pfe, Int[], Bool[])
    end

    # Compute demeaned y
    py = model_response(mf)[:]
    yname = rt.eterms[1]
    if eltype(py) != Float64
        y = convert(py, Float64)
    else
        y = py
    end
    broadcast!(*, y, y, sqrtw)
    oldy = deepcopy(y)
    residualize!(y, pfe, Int[], Bool[])

    ##############################################################################
    ##
    ## Estimate Model on Matrix
    ##
    ##############################################################################

    # initialize factor models at 0.1
    idpool = fill(0.1, length(id.pool), m.rank)
    timepool = fill(0.1, length(time.pool), m.rank)

    if !has_regressors
        fp = FactorModel(y, sqrtw, id.refs, time.refs, m.rank)
        fs = FactorSolution(idpool, timepool)
        # factor model

        (fs, iterations, converged) =
            fit!(Val{method}, fp, fs; maxiter = maxiter, tol = tol, lambda = lambda)
    else
        # interactive fixed effect
        # initialize fs
        coef = X \ y
        fp = FactorModel(y - X * coef, sqrtw, id.refs, time.refs, m.rank)
        fs = FactorSolution(idpool, timepool)
        fit!(Val{:levenberg_marquardt}, fp, fs; maxiter = 100, tol = 1e-3, lambda = lambda)

        fs = InteractiveFixedEffectsSolution(coef, fs.idpool, fs.timepool)
        fp = InteractiveFixedEffectsModel(y, sqrtw, X, id.refs, time.refs, m.rank)
        ym = deepcopy(y)
        Xm = deepcopy(X)

        while true
            # estimate the model
           (fs, iterations, converged) =
                fit!(Val{method}, fp, fs; maxiter = maxiter, tol = tol, lambda = lambda)
            # check that I obtain the same coefficient if I solve
            # y ~ x + γ1 x factors + γ2 x loadings
            # if not, this means fit! ended up on a a local minimum.
            # restart with randomized coefficients, factors, loadings
            newpfe = FixedEffectProblem(getfactors(fp, fs), Val{:lsmr})
            residualize!(ym, newpfe, Int[], Bool[], tol = tol, maxiter = maxiter)
            residualize!(Xm, newpfe, Int[], Bool[], tol = tol, maxiter = maxiter)
            ydiff = Xm * (fs.b - Xm \ ym)
            if iterations >= maxiter || norm(ydiff)  <= 0.01 * norm(y)
                break
            end
            info("Algorithm ended up on a local minimum. Restarting from a new, random, x0.")
            map!(x -> randn() * x, fs, fs)
            copy!(ym, y)
            copy!(Xm, X)
        end
    end

    ##############################################################################
    ##
    ## Compute residuals
    ##
    ##############################################################################

    # compute residuals
    fp = FactorModel(deepcopy(y), sqrtw, id.refs, time.refs, m.rank)
    if has_regressors
        BLAS.gemm!('N', 'N', -1.0, X, fs.b, 1.0, fp.y)
    end
    subtract_factor!(fp, fs)
    broadcast!(/, fp.y, fp.y, sqrtw)
    residuals = fp.y
    ##############################################################################
    ##
    ## Compute errors
    ##
    ##############################################################################

    if !has_regressors
        ess = sumabs2(residuals)
    else
        residualsm = ym .- Xm * fs.b
        crossxm = cholfact!(At_mul_B(Xm, Xm))
        ## compute the right degree of freedom
        df_absorb_fe = 0
        if has_absorb
            df_absorb_fe = 0
            ## poor man adjustement of df for clustedered errors + fe: only if fe name != cluster name
            for fe in fes
                if typeof(vcov_method) == VcovCluster && in(fe.factorname, vcov_vars)
                    df_absorb_fe += 0
                else
                    df_absorb_fe += sum(fe.scale .!= zero(Float64))
                end
            end
        end
        df_absorb_factors = 0
        newfes = getfactors(fp, fs)
        for fe in newfes
            df_absorb_factors +=
                (typeof(vcov_method) == VcovCluster && in(fe.factorname, vcov_vars)) ?
                    0 : sum(fe.scale .!= zero(Float64))
        end
        df_residual = max(size(X, 1) - size(X, 2) - df_absorb_fe - df_absorb_factors, 1)

        ## estimate vcov matrix
        vcov_data = VcovData(Xm, crossxm, residualsm, df_residual)
        matrix_vcov = vcov!(vcov_method_data, vcov_data)

        # compute various r2
        nobs = size(subdf, 1)
        ess = sumabs2(residualsm)
        tss = compute_tss(ym, rt.intercept, sqrtw)
        r2_within = 1 - ess / tss

        ess = sumabs2(residuals)
        tss = compute_tss(oldy, rt.intercept || has_absorb, sqrtw)
        r2 = 1 - ess / tss
        r2_a = 1 - ess / tss * (nobs - rt.intercept) / df_residual
    end

    ##############################################################################
    ##
    ## Save factors and loadings in a dataframe
    ##
    ##############################################################################

    if !save
        augmentdf = DataFrame()
    else
        augmentdf = DataFrame(fp, fs, esample)
        # save residuals in a dataframe
        if all(esample)
            augmentdf[:residuals] = residuals
        else
            augmentdf[:residuals] =  DataArray(Float64, size(augmentdf, 1))
            augmentdf[esample, :residuals] = residuals
        end

        # save fixed effects in a dataframe
        if has_absorb
            # residual before demeaning
            mf = simpleModelFrame(subdf, rt, esample)
            oldresiduals = model_response(mf)[:]
            if has_regressors
                oldX = ModelMatrix(mf).m
                BLAS.gemm!('N', 'N', -1.0, oldX, coef, 1.0, oldresiduals)
            end
            fp = FactorModel(oldresiduals, sqrtw, id.refs, time.refs, m.rank)
            subtract_factor!(fp, fs)
            axpy!(-1.0, residuals, oldresiduals)
            # get fixed effect
            augmentdf = hcat(augmentdf, getfe!(pfe, oldresiduals, esample))
        end
    end

    if !has_regressors
        return FactorResult(esample, augmentdf, ess, iterations, converged)
    else
        return InteractiveFixedEffectsResult(fs.b, matrix_vcov, esample, augmentdf,
            coef_names, yname, f, nobs, df_residual, r2, r2_a, r2_within,
            ess, sum(iterations), all(converged))
    end
end

# Symbol to formula Symbol ~ 0
function fit(m::InteractiveFixedEffectModel,
             variable::Symbol, args...; kwargs...)
    formula = Formula(variable, 0)
    fit(m, formula, args...; kwargs...)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct BlobRestoreRange <: SwaggerModel
    startRange::Any # spec type: Union{ Nothing, String } # spec name: startRange
    endRange::Any # spec type: Union{ Nothing, String } # spec name: endRange

    function BlobRestoreRange(;startRange=nothing, endRange=nothing)
        o = new()
        validate_property(BlobRestoreRange, Symbol("startRange"), startRange)
        setfield!(o, Symbol("startRange"), startRange)
        validate_property(BlobRestoreRange, Symbol("endRange"), endRange)
        setfield!(o, Symbol("endRange"), endRange)
        o
    end
end # type BlobRestoreRange

const _property_map_BlobRestoreRange = Dict{Symbol,Symbol}(Symbol("startRange")=>Symbol("startRange"), Symbol("endRange")=>Symbol("endRange"))
const _property_types_BlobRestoreRange = Dict{Symbol,String}(Symbol("startRange")=>"String", Symbol("endRange")=>"String")
Base.propertynames(::Type{ BlobRestoreRange }) = collect(keys(_property_map_BlobRestoreRange))
Swagger.property_type(::Type{ BlobRestoreRange }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BlobRestoreRange[name]))}
Swagger.field_name(::Type{ BlobRestoreRange }, property_name::Symbol) =  _property_map_BlobRestoreRange[property_name]

function check_required(o::BlobRestoreRange)
    (getproperty(o, Symbol("startRange")) === nothing) && (return false)
    (getproperty(o, Symbol("endRange")) === nothing) && (return false)
    true
end

function validate_property(::Type{ BlobRestoreRange }, name::Symbol, val)
end


text:
@testset "Categorical" begin
    Nclass = 3
    n = 10
    bij_lik = CategoricalLikelihood(
        BijectiveSimplexLink(LogisticSoftMaxLink(zeros(Nclass)))
    )
    @test_skip test_auglik(
        bij_lik;
        n,
        f=[randn(n) for _ in 1:AGPL.nlatent(bij_lik)],
        qf=[Normal.(randn(n), 1.0) for _ in 1:AGPL.nlatent(bij_lik)],
    )
    nonbij_lik = CategoricalLikelihood(LogisticSoftMaxLink(zeros(Nclass)))
    @test_skip test_auglik(
        nonbij_lik,
        n,
        f=[randn(n) for _ in 1:AGPL.nlatent(nonbij_lik)],
        qf=[Normal.(randn(n), 1.0) for _ in 1:AGPL.nlatent(nonbij_lik)],
    )
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://discourse.julialang.org/t/generically-iterating-over-arrays-of-points/65273/2

# %%
using StaticArrays

eachpoint(itr) = itr
eachpoint(ptmatrix::AbstractMatrix) = eachcol(ptmatrix)

function f(pts)
    for pt in eachpoint(pts)
        println(pt)
    end
end

f(SVector(k+1, k+2, k+3) for k in 0:3:12)

# %%
f(reshape(1:15, 3, 5))

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using InvariantMeasures
using ValidatedNumerics

@testset "Ulam assembler" begin

D = Mod1Dynamic(x->2*x)
B = Ulam(8)
P = assemble(B, D)

Ptrue = [
		0.5 0.5 0 0 0 0 0 0;
		0 0 0.5 0.5 0 0 0 0;
		0 0 0 0 0.5 0.5 0 0;
		0 0 0 0 0 0 0.5 0.5;
		0.5 0.5 0 0 0 0 0 0;
		0 0 0.5 0.5 0 0 0 0;
		0 0 0 0 0.5 0.5 0 0;
		0 0 0 0 0 0 0.5 0.5;
		]
Ptrue = Ptrue'

@test all(contains_zero.(P-Ptrue))

# mod1_dynamic with a non-Markov dynamic
D = mod1_dynamic(x->x+0.5)
B = Ulam(8)
P = assemble(B, D)

Ptrue = [
0  0  0  0  1  0  0  0;
0  0  0  0  0  1  0  0;
0  0  0  0  0  0  1  0;
0  0  0  0  0  0  0  1;
1  0  0  0  0  0  0  0;
0  1  0  0  0  0  0  0;
0  0  1  0  0  0  0  0;
0  0  0  1  0  0  0  0;
		]

@test all(contains_zero.(P-Ptrue))

@test opnormbound(L1, DiscretizedOperator(B, D)) == 1
@test opnormbound(Linf, DiscretizedOperator(B, D)) == 1

end


text:
using Documenter, Weave, ONet

for file ∈ readdir(joinpath(dirname(pathof(ONet)), "..", "docs", "jmd"))
      weave(joinpath(dirname(pathof(ONet)), "..", "docs", "jmd", file),
            out_path = joinpath(dirname(pathof(ONet)), "..", "docs", "src"),
            doctype = "github")
end

makedocs(format = Documenter.HTML(),
         modules = [ONet],
         sitename = "ONet.jl",
         pages = ["Introduction" => "index.md",
                  "Manual" => "manual.md",
                  "API" => "api.md"]
    )

deploydocs(repo = "github.com/uva-bi-sdad/ONet.jl.git")


text:
# This file is for functions that are new functions sepecifically for workin with
# NamedDimsArrays, rather than overloads of existing functions

"""
    rename(nda::NamedDimsArray, names)

Returns a new `NameDimsArray` with the given dimension `names`.
`rename` outright replaces the names; while still wrapping the same backing array.
Unlike the constructor, it does not require that new names are compatible
with the old names (though you do still need to match the number of dimensions).
"""
rename(nda::NamedDimsArray, names) = NamedDimsArray(parent(nda), names)

"""
    refine_names(x, names)

Refine the names of the dimensions of `x` to match `names`.
This is like [`rename`](ref), but it only affects unnamed dimensions.
I.e. dimensions of a `NamedDimsArray` called `:_`, or any dimensions of an
`AbstractArray` in general.
"""
@inline function refine_names(orig::NamedDimsArray{L}, names::Tuple) where {L}
    new_names = unify_names(names, L)
    return NamedDimsArray{new_names}(parent(orig))
end
@inline refine_names(orig::AbstractArray, names::Tuple) = NamedDimsArray{names}(orig)

"""
    unname(A::NamedDimsArray) -> AbstractArray
    unname(A::AbstractArray) -> AbstractArray

Return the input array `A` without any dimension names.

For `NamedDimsArray`s this returns the parent array, equivalent to calling `parent`, but for
any other `AbstractArray` simply returns the input.
"""
unname(x::NamedDimsArray) = parent(x)
unname(x::AbstractArray) = x

"""
    dimnames(A) -> Tuple
    dimnames(A, d) -> Symbol

Return the names of all the dimensions of the array `A`,
or just the one for the `d`-th dimension.

Gives wildcards `:_` if this is not a `NamedDimsArray`.
Like `size(A, d)`, it allows `d > ndims(A)`, in this case all the trailing dimension are given the wildcard name (`:_`).
"""
dimnames(::Type{<:NamedDimsArray{L}}) where {L} = L
dimnames(::Type{<:AbstractArray{T, N}}) where {T, N} = ntuple(_->:_, N)
dimnames(x::T) where T<:AbstractArray = dimnames(T)

function dimnames(AT::Type{<:AbstractArray{T,N}}, d::Integer) where {T,N}
    if 1 <= d <= N
        return dimnames(AT)[d]
    elseif d > N
        return :_
    else
        throw(DimensionMismatch("dimnames: dimension out of range"))
    end
end
dimnames(x::T, d::Integer) where T<:AbstractArray = dimnames(T, d)


text:
@test supertype(ImageInTerminal.BigBlocks) <: ImageInTerminal.ImageEncoder
@test supertype(ImageInTerminal.SmallBlocks) <: ImageInTerminal.ImageEncoder

@testset "_charof" begin
    @test @inferred(ImageInTerminal._charof(0.0)) === '⋅'
    @test @inferred(ImageInTerminal._charof(0.2)) === '░'
    @test @inferred(ImageInTerminal._charof(0.5)) === '▒'
    @test @inferred(ImageInTerminal._charof(0.8)) === '▓'
    @test @inferred(ImageInTerminal._charof(1.0)) === '█'
end

@testset "encodeimg 256 small" begin
    @testset "gray square" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), gray_square, 2, 2)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 2
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232;48;5;248m▀\e[38;5;239;48;5;255m▀\e[0m"
    end
    @testset "transparent gray square" begin
        # alpha is ignored for small block encoding.
        # So this yields the exact same results as above.
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), gray_square_alpha, 2, 2)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 2
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232;48;5;248m▀\e[38;5;239;48;5;255m▀\e[0m"
    end
    # the following tests checks the correct use of restrict
    # we compare against a hand checked reference output
    @testset "camera man" begin
        # this checks the correct use of restrict
        # we compare against a hand checked reference output
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), camera_man, 20, 20)
        @test typeof(res) <: Vector{String}
        @test h === 9
        @test w === 17
        @test_reference "reference/camera_small_20x20_256.txt" res
        # too small size
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), camera_man, 1, 1)
        @test typeof(res) <: Vector{String}
        @test h === 3
        @test w === 5
        @test_reference "reference/camera_small_1x1_256.txt" res
        # bigger version
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), camera_man, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 33
        @test_reference "reference/camera_small_60x60_256.txt" res
    end
    @testset "lighthouse" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), lighthouse, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 49
        @test_reference "reference/lighthouse_small_60x60_256.txt" res
    end
    @testset "toucan" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), toucan, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 20
        @test w === 42
        @test_reference "reference/toucan_small_60x60_256.txt" res
    end
end

# ====================================================================

@testset "encodeimg 256 big" begin
    @testset "gray square" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), gray_square, 4, 4)
        @test typeof(res) <: Vector{String}
        @test h === 2
        @test w === 4
        @test length(res) === 2
        @test res[1] == "\e[0m\e[38;5;232m██\e[38;5;239m██\e[0m"
        @test res[2] == "\e[0m\e[38;5;248m██\e[38;5;255m██\e[0m"
    end
    @testset "transparent gray square" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), gray_square_alpha, 4, 4)
        @test typeof(res) <: Vector{String}
        @test h === 2
        @test w === 4
        @test length(res) === 2
        @test res[1] == "\e[0m\e[38;5;232m██\e[38;5;239m▓▓\e[0m"
        @test res[2] == "\e[0m\e[38;5;248m░░\e[38;5;255m⋅⋅\e[0m"
    end
    # the following tests checks the correct use of restrict
    # we compare against a hand checked reference output
    @testset "camera man" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), camera_man, 40, 40)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 34
        @test_reference "reference/camera_big_20x20_256.txt" res
    end
    @testset "lighthouse" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), lighthouse, 50, 50)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 50
        @test_reference "reference/lighthouse_big_50x50_256.txt" res
    end
    @testset "toucan" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), toucan, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 20
        @test w === 44
        @test_reference "reference/toucan_big_60x60_256.txt" res
    end
end

# ====================================================================

@testset "encodeimg 24bit small" begin
    @testset "gray square" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), gray_square, 2, 2)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 2
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0;48;2;178;178;178m▀\e[38;2;76;76;76;48;2;255;255;255m▀\e[0m"
    end
    @testset "transparent gray square" begin
        # alpha is ignored for small block encoding.
        # So this yields the exact same results as above.
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), gray_square_alpha, 2, 2)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 2
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0;48;2;178;178;178m▀\e[38;2;76;76;76;48;2;255;255;255m▀\e[0m"
    end
    # the following tests checks the correct use of restrict
    # we compare against a hand checked reference output
    @testset "camera man" begin
        # this checks the correct use of restrict
        # we compare against a hand checked reference output
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), camera_man, 20, 20)
        @test typeof(res) <: Vector{String}
        @test h === 9
        @test w === 17
        @test_reference "reference/camera_small_20x20_24bit.txt" res
        # bigger version
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), camera_man, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 33
        @test_reference "reference/camera_small_60x60_24bit.txt" res
    end
    @testset "lighthouse" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), lighthouse, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 49
        @test_reference "reference/lighthouse_small_60x60_24bit.txt" res
    end
    @testset "toucan" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), toucan, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 20
        @test w === 42
        @test_reference "reference/toucan_small_60x60_24bit.txt" res
    end
end

# ====================================================================

@testset "encodeimg 24bit big" begin
    @testset "gray square" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), gray_square, 4, 4)
        @test typeof(res) <: Vector{String}
        @test h === 2
        @test w === 4
        @test length(res) === 2
        @test res[1] == "\e[0m\e[38;2;0;0;0m██\e[38;2;76;76;76m██\e[0m"
        @test res[2] == "\e[0m\e[38;2;178;178;178m██\e[38;2;255;255;255m██\e[0m"
    end
    @testset "transparent gray square" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), gray_square_alpha, 4, 4)
        @test typeof(res) <: Vector{String}
        @test h === 2
        @test w === 4
        @test length(res) === 2
        @test res[1] == "\e[0m\e[38;2;0;0;0m██\e[38;2;76;76;76m▓▓\e[0m"
        @test res[2] == "\e[0m\e[38;2;178;178;178m░░\e[38;2;255;255;255m⋅⋅\e[0m"
    end
    # the following tests checks the correct use of restrict
    # we compare against a hand checked reference output
    @testset "camera man" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), camera_man, 40, 40)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 34
        @test_reference "reference/camera_big_20x20_24bit.txt" res
    end
    @testset "lighthouse" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), lighthouse, 50, 50)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 50
        @test_reference "reference/lighthouse_big_50x50_24bit.txt" res
    end
    @testset "toucan" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), toucan, 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 20
        @test w === 44
        @test_reference "reference/toucan_big_60x60_24bit.txt" res
    end
end

# ====================================================================

@testset "encodeimg 256 small" begin
    @testset "gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), gray_line, 10)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 4
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232m█\e[38;5;239m█\e[38;5;248m█\e[38;5;255m█\e[0m"
    end
    @testset "transparent gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), gray_line_alpha, 10)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 4
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232m█\e[38;5;239m▓\e[38;5;248m░\e[38;5;255m⋅\e[0m"
    end
    @testset "rgb line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), rgb_line, 8)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 6
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;18m█\e[38;5;56m█\e[38;5;91m█\e[38;5;126m█\e[38;5;161m█\e[38;5;88m█\e[0m"
    end
end

# ====================================================================

@testset "encodeimg 256 big" begin
    @testset "gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), gray_line, 9)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232m██ \e[0m … \e[38;5;255m██ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), gray_line, 12)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 12
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232m██ \e[38;5;239m██ \e[38;5;248m██ \e[38;5;255m██ \e[0m"
    end
    @testset "transparent gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), gray_line_alpha, 10)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232m██ \e[0m … \e[38;5;255m⋅⋅ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), gray_line_alpha, 12)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 12
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;232m██ \e[38;5;239m▓▓ \e[38;5;248m░░ \e[38;5;255m⋅⋅ \e[0m"
    end
    @testset "rgb line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), rgb_line, 9)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;21m██ \e[0m … \e[38;5;196m██ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), rgb_line, 22)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 21
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;21m██ \e[38;5;21m██ \e[38;5;56m██ \e[0m … \e[38;5;161m██ \e[38;5;196m██ \e[38;5;196m██ \e[0m"
    end
end

# ====================================================================

@testset "encodeimg 24bit small" begin
    @testset "gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), gray_line, 10)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 4
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0m█\e[38;2;76;76;76m█\e[38;2;178;178;178m█\e[38;2;255;255;255m█\e[0m"
    end
    @testset "transparent gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), gray_line_alpha, 10)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 4
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0m█\e[38;2;76;76;76m▓\e[38;2;178;178;178m░\e[38;2;255;255;255m⋅\e[0m"
    end
    @testset "rgb line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor24bit(), rgb_line, 8)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 6
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;6;0;122m█\e[38;2;47;0;208m█\e[38;2;101;0;154m█\e[38;2;154;0;101m█\e[38;2;208;0;47m█\e[38;2;122;0;6m█\e[0m"
    end
end

# ====================================================================

@testset "encodeimg 24bit big" begin
    @testset "gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), gray_line, 9)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0m██ \e[0m … \e[38;2;255;255;255m██ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), gray_line, 12)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 12
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0m██ \e[38;2;76;76;76m██ \e[38;2;178;178;178m██ \e[38;2;255;255;255m██ \e[0m"
    end
    @testset "transparent gray line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), gray_line_alpha, 10)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0m██ \e[0m … \e[38;2;255;255;255m⋅⋅ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), gray_line_alpha, 12)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 12
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;0m██ \e[38;2;76;76;76m▓▓ \e[38;2;178;178;178m░░ \e[38;2;255;255;255m⋅⋅ \e[0m"
    end
    @testset "rgb line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), rgb_line, 9)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;255m██ \e[0m … \e[38;2;255;0;0m██ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor24bit(), rgb_line, 22)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 21
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;2;0;0;255m██ \e[38;2;13;0;242m██ \e[38;2;27;0;228m██ \e[0m … \e[38;2;228;0;27m██ \e[38;2;242;0;13m██ \e[38;2;255;0;0m██ \e[0m"
    end
end

@testset "non-1 indexing" begin
    @testset "rgb line" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), OffsetArray(rgb_line, (-1,)), 8)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 6
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;18m█\e[38;5;56m█\e[38;5;91m█\e[38;5;126m█\e[38;5;161m█\e[38;5;88m█\e[0m"
    end
    @testset "rgb line2" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), OffsetArray(rgb_line, (2,)), 9)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 9
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;21m██ \e[0m … \e[38;5;196m██ \e[0m"
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.BigBlocks(), ImageInTerminal.TermColor256(), OffsetArray(rgb_line, (-2,)), 22)
        @test typeof(res) <: Vector{String}
        @test h === 1
        @test w === 21
        @test length(res) === 1
        @test res[1] == "\e[0m\e[38;5;21m██ \e[38;5;21m██ \e[38;5;56m██ \e[0m … \e[38;5;161m██ \e[38;5;196m██ \e[38;5;196m██ \e[0m"
    end
    @testset "lighthouse" begin
        res, h, w = ensurecolor(ImageInTerminal.encodeimg, ImageInTerminal.SmallBlocks(), ImageInTerminal.TermColor256(), OffsetArray(lighthouse, (2,-10)), 60, 60)
        @test typeof(res) <: Vector{String}
        @test h === 17
        @test w === 49
        @test_reference "reference/lighthouse_small_60x60_256.txt" res
    end
end


text:
n_values  = [200, 300, 400]
k_values  = [200, 300, 400]
m_values  = [200, 300, 400]

testset_name_suffix = "(main)"

for T ∈ (Float64,Float32,Int64,Int32)
  @time test_real(T, m_values, k_values, n_values, testset_name_suffix)
  @time test_complex(T, m_values, k_values, n_values, testset_name_suffix)
end

A = rand(2,2); B = rand(2,2); AB = A*B; C = fill(NaN, 2, 2);
@test Octavian.matmul!(C, A, B, true, false) ≈ AB
@test Octavian.matmul!(C, A, B, true, true) ≈ 2AB


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

#import HDF5

@testset "hdf5" begin
end


text:
##COMPARE TIME TAKEN BY ADAPTIVE AND NON-ADAPTIVE RE-ABC IN ABAKALIKI EXAMPLE
using RareEventABC
using Distributions
using DataFrames
using ProgressMeter
using Plots; pgfplots();
using StatPlots;

"
A pseudo-marginal Metropolis-Hastings algorithm which uses RE-ABC likelihood estimates.

## Arguments

- `a` Details of the application.
- `schedule` Sequence of epsilon values. Use this option for non-adaptive RE-ABC.
- `targetε` Target ε value. Use this option for adaptive RE-ABC.
- `N` Number of iterations to do.
- `nparticles` Number of particles to use.
- `prior` Prior distribution.
- `θinit` Initial state of chain.
- `inc` Distribution of proposal increments
- `k` Number of particles to accept at each iteration of adaptive RE-ABC.
"
function timedPMMH(a::RareEventABC.ApplicationDetails, schedule::Array{Float64, 1}, N::Int, nparticles::Int, prior::MultivariateDistribution, θinit::Array{Float64, 1}, inc::MultivariateDistribution)
    p = length(θinit)
    states = zeros(Float64, (N,p))
    loglikelihoods = zeros(Float64, N)
    states[1,:] = θ = copy(θinit)
    times = DataFrame(time=Float64[], stop=String[])
    tic()
    loglikelihoods[1] = ll = doSMC(a, θ, schedule, nparticles=nparticles, quiet=true).llest
    t = toq()
    push!(times, (t, "complete"))
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            tic()
            p = doSMC(a, θprop, schedule, nparticles=nparticles, llBound=llBound, quiet=true)
            t = toq()
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
                push!(times, (t, "complete"))
            else
                push!(times, (t, "early"))
            end
        end
        states[i,:], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods, times
end

function timedPMMH(a::RareEventABC.ApplicationDetails, targetε::Float64, N::Int, nparticles::Int, prior::MultivariateDistribution, θinit::Array{Float64, 1}, inc::MultivariateDistribution; k=ceil(Int, nparticles/2))
    p = length(θinit)
    states = zeros(Float64, (N,p))
    loglikelihoods = zeros(Float64, N)
    states[1,:] = θ = copy(θinit)
    times = DataFrame(time=Float64[], stop=String[])
    tic()
    loglikelihoods[1] = ll = doSMC(a, θ, targetε, nparticles=nparticles, quiet=true, k=k).llest
    t = toq()
    push!(times, (t, "complete"))
    logprior = logpdf(prior, θ)
    @showprogress for i in 2:N
        θprop = θ + rand(inc)
        logpriorProp = logpdf(prior, θprop)
        llBound = log(rand()) + ll + logprior - logpriorProp
        if llBound < Inf
            tic()
            p = doSMC(a, θprop, targetε, nparticles=nparticles, llBound=llBound, quiet=true, k=k)
            t = toq()
            if !p.reachedBound ##Acceptance
                θ, ll = θprop, p.llest
                push!(times, (t, "complete"))
            else
                push!(times, (t, "early"))
            end
        end
        states[i,:], loglikelihoods[i] = θ, ll
    end
    states, loglikelihoods, times
end

##ABAKALIKI DATA
size_obs = 30;
popsize = 120;
robs = cumsum([0,13,7,2,3,0,0,1,4,5,3,2,0,2,0,5,3,1,4,0,1,1,1,2,0,1,5,0,5,5]);
duration_obs = maximum(robs);
robs = [robs; [Inf for i in 1:90]];

##Model (SIR with gamma infectious period) and prior
sir_gamma_I = RareEventABC.SEIRDetails(robs, popsize, θ->RareEventABC.Exponential(1.), θ->RareEventABC.DiracDelta(0.0), θ->Distributions.Gamma(θ[2], θ[3]), x->x);
prior = MvIndependent([Exponential(10.), Exponential(10.), Exponential(10.)]);

##Tuning choices (taken from original analysis)
start = [0.103064, 4.11688, 4.28655];
inc = MvNormal(zeros(3), [0.00139816   -0.113794    0.0590885;
                         -0.113794     19.9754    -10.3812;
                          0.0590885   -10.3812      7.47516]);

##NON-ADAPTIVE RE-ABC
##Tuning
srand(1);
schedule = doSMC(sir_gamma_I, start, 15.0, nparticles=300, d=size_obs, quiet=false).schedule;

##Main run
thetas, ll, times = timedPMMH(sir_gamma_I, schedule, 2000, 300, prior, start, inc);

##ADAPTIVE RE-ABC
thetas2, ll2, times2 = timedPMMH(sir_gamma_I, 15.0, 2000, 300, prior, start, inc);

alltimes = vcat(hcat(times, fill("non-adaptive", size(times)[1])), hcat(times2, fill("adaptive", size(times2)[1])))
rename!(alltimes, :x1, :method)

writetable("pmmh_Abakaliki_timed.csv", alltimes) ##To plot using R


text:
########################################################################
# Atoms of the nuclear-norm ball.
########################################################################
"""
    NucBallAtom(m, n, u, v)

Atom in the NucBall of dimension m×n.
"""
mutable struct NucBallAtom{T1<:Int64, T2<:Vector{Float64}} <: AbstractAtom
    m::T1
    n::T1
    u::T2
    v::T2
    function NucBallAtom(u::Vector{Float64}, v::Vector{Float64})
        m = length(u)
        n = length(v)
        new{Int64, Vector{Float64}}(m, n, u, v)
    end
end

"""
    *(M::AbstractLinearOp, a::NucBallAtom) -> Y

Multiply a nuclear-norm atom by a linear map, and return the matrix `Y`.
"""
Base.:(*)(M::LinearOp, a::NucBallAtom{Int64, Vector{Float64}}) = M*vec(a.u*a.v')
LinearAlgebra.mul!(Ma::Vector{Float64}, M::LinearOp, a::NucBallAtom{Int64, Vector{Float64}}) = mul!(Ma, M, vec(a.u*a.v'))
Base.vec(a::NucBallAtom{Int64, Vector{Float64}}) = vec(a.u*a.v')
Base.length(a::NucBallAtom{Int64, Vector{Float64}}) = a.m * a.n
Base.size(a::NucBallAtom{Int64, Vector{Float64}}) = a.m, a.n

########################################################################
# Nuclear-norm ball.
########################################################################
"""
    NucBall(m, n, maxrank=min(m,n))

Atomic set defined by nuclear-norm ball for the space of `m`-by-`n` matrices.
The atomic set takes two optional parameters:

`maxrank` is the maximum dimension of the face exposed by a vector.
"""
struct NucBall{T<:Int64} <: AbstractAtomicSet
    m::T
    n::T
    maxrank::T
    function NucBall(m::Int64, n::Int64, maxrank::Int64)
        m ≥ 1 || throw(DomainError(m,"m must be a positive integer"))
        n ≥ 1 || throw(DomainError(n,"n must be a positive integer"))
        min(m,n) ≥ maxrank ≥ 0 || throw(DomainError(maxrank,"maxrank must be ≥ 0"))
        new{Int64}(m, n, maxrank)
    end
end
NucBall(m::Int64, n::Int64; maxrank=min(m, n)) = NucBall(m, n, maxrank)

"""
    gauge(A::NucBall, x::Matrix)

Gives the nuclear-norm of `x`.
"""
gauge(::NucBall{Int64}, x::Matrix{Float64}) = sum(svdvals(x))
gauge(A::NucBall{Int64}, x::Vector{Float64}) = gauge(A, reshape(x, A.m, A.n))

"""
    support(A::NucBall, z::AbstractMatrix)

Gives the largest singular value of `z`.
"""
function support(A::NucBall{Int64}, z::AbstractMatrix{Float64})
    r = A.maxrank
    if r == 0
        s = [0]
    else
        ~, s, ~ = svds(z, nsv=1)[1]
    end
    return s[1]
end

support(A::NucBall{Int64}, z::AbstractVector{Float64}) = support(A, reshape(z, A.m, A.n))

"""
    expose!(A::NucBall, z::Matrix, a::NucBallAtom)

Gives the top singular vectors of `z`.
"""
function expose!(::NucBall{Int64}, z::AbstractMatrix{Float64}, a::NucBallAtom{Int64, Vector{Float64}})
    u, ~, v = svds(z, nsv=1)[1]
    a.u .= u[:]
    a.v .= v[:]
    return nothing
end
expose!(A::NucBall{Int64}, z::AbstractVector{Float64}, a::NucBallAtom{Int64, Vector{Float64}}) = expose!(A, reshape(z, A.m, A.n), a)

function expose(A::NucBall{Int64}, z::AbstractMatrix{Float64})
    m, n = size(A)
    a = NucBallAtom(zeros(m), zeros(n))
    expose!(A, z, a)
    return a
end
expose(A::NucBall{Int64}, z::AbstractVector{Float64}) = expose(A, reshape(z, A.m, A.n))

Base.length(A::NucBall{Int64}) = A.m*A.n
Base.size(A::NucBall{Int64}) = A.m, A.n
rank(A::NucBall{Int64}) = A.maxrank
atom_name(A::NucBall{Int64}) = "nuclear-norm ball"
atom_description(A::NucBall{Int64}) = "{ uv^T | u ∈ ℝᵐ, v ∈ ℝⁿ, ||u||₂ = ||v||₂ ≤ 1 }"
atom_parameters(A::NucBall{Int64}) = "m = $(A.m), n = $(A.n), maxrank = $(A.maxrank)"

########################################################################
# Face of the nuclear-norm ball.
########################################################################

mutable struct NucBallFace{T1<:Int64, T2<:Matrix{Float64}, T3<:Adjoint{Float64, Matrix{Float64}}} <: AbstractFace
    m::T1
    n::T1
    r::T1
    U::T2
    V::T3
    function NucBallFace(U::Matrix{Float64}, V::Adjoint{Float64, Matrix{Float64}})
        m = size(U, 1)
        n = size(V, 1)
        r = size(U, 2)
        new{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}(m, n, r, U, V)
    end
end

"""
    face(A, z) -> F

Return a face `F` of the atomic set `A` exposed by the vector `z`. The dimension of
the exposed face is limited by `k=maxrank(A)`. If the dimension is being limited
by `k`, then the routine returns a set of rank `k` atoms that define the subset
of the exposed face.
"""
function face!(A::NucBall{Int64}, z::AbstractMatrix{Float64}, F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}})
    r = A.maxrank
    if r == 0
        U = zeros(A.m, 1)
        V = zeros(1, A.n)'
    elseif r < min(A.m, A.n) - 1
        # Arpack won't accept nsv too large
        U, _, V = svds(z, nsv=r)[1]
        F.U .= U
        F.V .= V
    else
        Fa = svd(z)
        F.U .= Fa[1].U[:, 1:r]
        F.V .= Fa[1].V[:, 1:r]
    end
    return nothing
end
face!(A::NucBall{Int64}, z::AbstractVector{Float64}, F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = face!(A, reshape(z, A.m, A.n), F)

function face(A::NucBall{Int64}, z::AbstractMatrix{Float64})
    m, n, r = A.m, A.n, max(A.maxrank, 1)
    U = zeros(m, r)
    V = zeros(r, n)'
    F = NucBallFace(U, V)
    face!(A, z, F)
    return F
end
face(A::NucBall{Int64}, z::AbstractVector{Float64}) = face(A, reshape(z, A.m, A.n))

"""
Given a face and a set of weights, reveal the corresponding
point on the face.
"""
Base.:(*)(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}, S::Matrix{Float64}) = vec(F.U*S*F.V')
Base.:(*)(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}, c::Vector{Float64}) = F*reshape(c, F.r, F.r)

"""
    *(M::LinearOp, F::NucBallFace)

Return a LinearMap `L` whose forward and adjoing operators
are defined as

    L*p = M*vec(U*matrix(p)*V')
    L'*q = vec(U'*matrix(M'*q)*V)

where the `n`-by-`k` matrices `U` and `V` are the left and right
singular vectors that define the face `F` of the nuclear-norm ball.
"""
function Base.:(*)(M::LinearOp, F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}})
    m, n, r = F.m, F.n, F.r
    f = p->M*vec(F.U*reshape(p, r, r)*F.V')
    fc = q->vec(F.U'*reshape(M'*q, m, n)*F.V)
    return LinearMap(f, fc, size(M, 1), r*r)
end

Base.:(*)(λ::Real, F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = λ*vec(F)
Base.length(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = F.m*F.n
Base.size(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = F.m, F.n
rank(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = F.r^2
vec(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = I(F.m*F.n)*F
face_name(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = "face of nuclear-norm ball"
face_parameters(F::NucBallFace{Int64, Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}}) = "m = $(F.m), n = $(F.n), rank = $(rank(F))"


text:
#This program fits a regression model with one
#unknown change point using generated data
using DataFrames
using Gadfly
using Distributions
dataFaked = false

function generateChangeData(ChangePointLengths,vecOfPoissonLambda)
	if length(ChangePointLengths) != length(vecOfPoissonLambda)
		error("length(ChangePointLengths) != length(vecOfPoissonLambda)")
	else
		asdf = [rand(Poisson(vecOfPoissonLambda[ii] ),ChangePointLengths[ii]) for ii = [1:length(vecOfPoissonLambda)]]
		reduce(vcat,asdf)
	end
end

if dataFaked
	ChangePointLengths = [20, 70]
	vecOfPoissonLambda = [1.2, 3.2]
	nobs = sum(ChangePointLengths)
	df = DataFrame(deaths = generateChangeData(ChangePointLengths,vecOfPoissonLambda))

else
	df = readtable("coaldata.txt")
	nobs = size(df)[1]
end
#Set conditions for experiment and generate the data
#Set up the prior values;
#NOTE: In this exercise, we are using the routine "gam_rnd",
#which is parameterized differently than the routine
#in the book. In particular, the second argument
#of the density is the reciprocal of how it is
#parameterized in our book. Thus, we change the
#hyperparameters below to accomodate for this
#fact, and alter the conditional distributions
#accordingly.
a1 = 1; a2 = 1;
d1 = 1; d2 = 1;

#Begin the Gibbs sampler
iter =1000000;
burn = 2000;
lambda_grid = [1:1:nobs-1]';
lambda_date = [1851:1:1961]';

gamma_final = zeros(iter-burn,1);	#Initialize the chains
delta_final = zeros(iter-burn,1);
lambda_final = zeros(iter-burn,1);

lambda = 50;	#First guess is that it is at time 50

for ii = 1:iter
    #----------------------------------------------------
    #Sample gamma, the parameter of the "first" Poisson density
    #----------------------------------------------------
    y_gamma = df[:deaths][1:lambda];
    n_lambda = length(y_gamma);
    gammas = rand(Gamma((a1 + sum(y_gamma)),1/(a2 + n_lambda)))
    #------------------------------------------
    #Sample delta, the parameter of the "second" Poisson density
    #----------------------------------------
    y_delta = df[:deaths][lambda+1:nobs];
    deltas = rand(Gamma((d1 + sum(y_delta)),1/(d2 + nobs-n_lambda)));

    #-----------------------------------
    #Sample the changepoint lambda
    #-----------------------------------

    log_dens_unnorm = zeros(nobs-1,1);
    for jj = 1:(nobs-1); #loop over the number of discrete values for Lambda
       ypart1 = df[:deaths][1:jj];
       ypart2 = df[:deaths][(jj+1):nobs];
       n_temp = length(ypart1);
       log_dens_unnorm[jj] = (sum(ypart1)*log(gammas) - n_temp*gammas + sum(ypart2)*log(deltas) - (nobs-n_temp)*deltas)[1];
       #dens_unnorm(jj,1) = (gammas^(sum(ypart1)))*exp(-n_temp*gammas)*(deltas^(sum(ypart2)))*exp(-(nobs-n_temp)*deltas);
    end;
    d = maximum(log_dens_unnorm);
    log_dens_unnorm = log_dens_unnorm -d;
    dens_unnorm = exp(log_dens_unnorm);
    dens_norm = dens_unnorm/sum(dens_unnorm);
    lambda = rand(Categorical(vec(dens_norm)));
    lambda_keep = lambda_date[lambda];

    if ii > burn;
        gamma_final[ii-burn] = gammas;
        delta_final[ii-burn] = deltas;
        lambda_final[ii-burn] = lambda_keep;
    end;
end;

plot(x=gamma_final,Geom.histogram(bincount = 100))
plot(x=delta_final,Geom.histogram(bincount = 100))
plot(x=lambda_final,Geom.histogram())

#disp('Post means and std deviatiations');
#disp('For betas, thetas, sig, tau and lambda');
#[mean(gamma_final')' std(gamma_final')']
#[mean(delta_final')' std(delta_final')']
#[mean(lambda_final) std(lambda_final)]

max_stars_repo_name: Darf-Ferrara/bayes-scripts

max_issues_repo_name: Darf-Ferrara/bayes-scripts

max_forks_repo_name: Darf-Ferrara/bayes-scripts


text:
module svgPath
using Graphics, Cairo
export elliptic_arc, parsePathData, drawPath, ScalePath

# ======================================================================================
# Adapted from: https://gitlab.gnome.org/GNOME/librsvg/blob/master/rsvg_internals/src/path_builder.rs
# ======================================================================================
approx_eq(a,b) = (abs(a-b) < 0.00390625) #

function elliptic_arc( ctx, x1::Float64, y1::Float64, rx::Float64, ry::Float64,
                   x_axis_rotation::Float64, is_large_arc::Float64, sweep::Float64,
                   x2::Float64, y2::Float64)

    approx_eq(x1, x2) && approx_eq(y1, y2) &&  return

    is_positive_sweep = (sweep == 1.0)

    # X-axis
    f = x_axis_rotation * pi / 180.0
    sinf, cosf = sin(f), cos(f)
    rx, ry = abs(rx), abs(ry)

    # A bit further down we divide by the square of the radii. Check that we won't divide by zero.
    if ((rx * rx) < 2) || ((ry * ry) < 2)
        line_to(ctx, x2, y2)
        return
    end

    k1 = (x1 - x2) / 2.0
    k2 = (y1 - y2) / 2.0

    x1_ =  cosf * k1 + sinf * k2
    y1_ = -sinf * k1 + cosf * k2

    gamma = (x1_ * x1_) / (rx * rx) + (y1_ * y1_) / (ry * ry)
    if gamma > 1.0
        rx *= sqrt(gamma)
        ry *= sqrt(gamma)
    end

    # Compute the center
    k1 = rx * rx * y1_ * y1_ + ry * ry * x1_ * x1_
    approx_eq(k1, 0.0) && return
    k1 = abs(sqrt(((rx * rx * ry * ry) / k1 - 1.0)))
    #println(k1)
    # println(is_positive_sweep == is_large_arc)
    # println(is_positive_sweep)
    # println(is_large_arc)
    is_positive_sweep == is_large_arc && (k1 = -k1)

    cx_ =  k1 * rx * y1_ / ry
    cy_ = -k1 * ry * x1_ / rx

    cx = cosf * cx_ - sinf * cy_ + (x1 + x2) / 2.0
    cy = sinf * cx_ + cosf * cy_ + (y1 + y2) / 2.0

    # Compute start angle
    k1 = (x1_ - cx_) / rx
    k2 = (y1_ - cy_) / ry
    k3 = (-x1_ - cx_) / rx
    k4 = (-y1_ - cy_) / ry

    k5 = abs(sqrt((k1 * k1 + k2 * k2)))
    approx_eq(k5, 0.0) && return

    k5 = k1 / k5
    k5 = clamp(k5, -1.0, 1.0)
    theta1 = acos(k5)
    k2 < 0.0 && (theta1 = -theta1)

    # Compute delta_theta
    k5 = abs(sqrt(((k1 * k1 + k2 * k2) * (k3 * k3 + k4 * k4))))
    approx_eq(k5, 0.0) && return

    k5 = (k1 * k3 + k2 * k4) / k5
    k5 = clamp(k5, -1.0, 1.0)
    delta_theta = acos(k5)
    if k1 * k4 - k3 * k2 < 0.0
        delta_theta = -delta_theta
    end

    if is_positive_sweep && delta_theta < 0.0
        delta_theta += pi * 2.0
    elseif !is_positive_sweep && delta_theta > 0.0
        delta_theta -= pi * 2.0
    end

    # Now draw the arc
    n_segs::Int32 = abs(ceil((delta_theta / (pi * 0.5 + 0.001)))) # as i32
    n_segs_dbl::Float64 = n_segs

    for i in 0:n_segs-1
        f::Float64 = i
        arc_segment(ctx, cx, cy,
            theta1 + f * delta_theta / n_segs_dbl,
            theta1 + (f + 1.0) * delta_theta / n_segs_dbl,
            rx, ry,
            x_axis_rotation )
    end
end
# ======================================================================================
# .
# ======================================================================================
function arc_segment( ctx, xc::Float64, yc::Float64, th0::Float64, th1::Float64,
                        rx::Float64, ry::Float64, x_axis_rotation::Float64)

    f = x_axis_rotation * pi / 180.0
    sinf = sin(f)
    cosf = cos(f)

    th_half = 0.5 * (th1 - th0)
    t = (8.0 / 3.0) * sin(th_half * 0.5) * sin(th_half * 0.5) / sin(th_half)
    x1 = rx * (cos(th0) - t * sin(th0))
    y1 = ry * (sin(th0) + t * cos(th0))
    x3 = rx * cos(th1)
    y3 = ry * sin(th1)
    x2 = x3 + rx * (t * sin(th1))
    y2 = y3 + ry * (-t * cos(th1))

    curve_to(ctx,   xc + cosf * x1 - sinf * y1,
                    yc + sinf * x1 + cosf * y1,
                    xc + cosf * x2 - sinf * y2,
                    yc + sinf * x2 + cosf * y2,
                    xc + cosf * x3 - sinf * y3,
                    yc + sinf * x3 + cosf * y3 )
end
#///////////////////////////////////////////////////////////////////////////////
function parsePathData(str)
           steps = matchall(r"([a-zA-Z][,\s-\.\d]*)", str)
           s = []
    for step in steps
        values = matchall(r"(-?\d*\.*\d+)", step)
        fvals = map(x->parse(Float64,x), values)
        l = step[1]
        push!(s, (l, fvals))
    end
    return s
end
#///////////////////////////////////////////////////////////////////////////////
# "M50,50   # Start x/y
#  A30,50   # arc (rx/ry)
#  0        # x-axis-rotation
#  0,1      # large-arc/sweep flags
#  100,100" # End
#///////////////////////////////////////////////////////////////////////////////
function drawPath(ctx, path)
    here = nothing

    for step in path
        l = step[1]
        if length(step[2]) > 0
            data = step[2]
        else
            data = [0,0]
        end
        if l == 'M'
            move_to(ctx, data[1], data[2])
            here = data[end-1:end]
        elseif l == 'm'
            rel_move_to(ctx, data[1], data[2])
            here = [here[1]+data[end-1] here[2]+data[end]]

        elseif l == 'L'
            for i in 1:2:length(data)-2
                line_to(ctx, data[i], data[i+1])
                here = data[end-1:end]
            end
        elseif l == 'l'
            for i in 1:2:length(data)-2
                rel_line_to(ctx, data[i], data[i+1])
                here = [here[1]+data[end-1] here[2]+data[end]]
            end

        elseif l == 'H'
                line_to(ctx, data[1], here[2])
                here[1] = data[1]
        elseif l == 'h'
                rel_line_to(ctx, data[1], 0)
                here[1] = here[1]+data[1]

        elseif l == 'V'
                line_to(ctx, here[1], data[1])
                here[2] = data[1]
        elseif l == 'v'
                rel_line_to(ctx, 0, data[1])
                here[2] = here[2]+data[1]

        elseif l == 'C'
            curve_to(ctx, data... );
            here = data[end-1:end]
        elseif l == 'c'
            rel_curve_to(ctx, data... );
            here = [here[1]+data[end-1] here[2]+data[end]]

        # elseif step[1] == 'q' curve_to(ctx, )
    elseif l == 'A'
            for i in 1:7:length(data)
                elliptic_arc( ctx, here[1], here[1],
                    data[i  ], data[i+1],
                    data[i+2], data[i+3], data[i+4],
                    data[i+5], data[i+6])
            end
            here = data[end-1:end]
        elseif l == 'a'
            for i in 1:7:length(data)
                elliptic_arc( ctx, here[1], here[1],
                    here[1]+data[i  ], here[2]+data[i+1],
                    data[i+2], data[i+3], data[i+4],
                    here[1]+data[i+5], here[2]+data[i+6])
            end
            here = [here[1]+data[end-1] here[2]+data[end]]

        elseif l == 'z' ||  l == 'Z'
            close_path(ctx) # stroke(ctx)
        end

    end

end
#///////////////////////////////////////////////////////////////////////////////
function ScalePath(path, scale)
    newpath = copy(path)

    for step in newpath
        l, data = step[1], step[2]
        if (l == 'M' || l == 'm' || l == 'L' || l == 'l' || l == 'C' || l == 'c' ||
           l == 'V' || l == 'v' || l == 'Q' || l == 'q' || l == 'S' || l == 's')
            data .*= scale
        elseif l == 'a' ||  l == 'A'
            for i in 1:7:length(data) # print(b[i:i+6])
                data[i  ] *= scale
                data[i+1] *= scale
                data[i+5] *= scale
                data[i+6] *= scale
            end
        end

    end
    return newpath
end

#///////////////////////////////////////////////////////////////////////////////
#
#///////////////////////////////////////////////////////////////////////////////

# function get_current_point(ctx::CairoContext)
# #    ccall((:cairo_get_current_point, Cairo._jl_libcairo), nothing, (Ptr{nothing},Ptr{Cdouble},Ptr{Cdouble}), ctx.ptr, x, y)
#     w::Vector{Cdouble,2} = []
#     w = ccall((:cairo_get_current_point, Cairo._jl_libcairo), Vector{Cdouble,2}, (Ptr{nothing},), ctx.ptr)
#     return w
# end

end # module


text:
function large01()
#
#  test polynomial suggested by Jenkins and Traub
#
    p = [
         1.0;
         0;
        -0.28000000000000;
        -0.59200000000000;
        -0.96890000000000;
        -0.62020000000000;
        0.54382200000000;
        0.08842680000000;
        0.43204191000000;
        -0.15503004000000;
        0.44119337520000;
        -0.65443195750000;
        -0.35708097260000;
        -0.55597791620000;
        0.63916879370000;
        0.14637576020000;
        -0.06709193010000;
        -0.03193754010000;
        0.50857405010000;
        -0.14968258190000;
        0.63273024920000
        ];
    z = [
        -1.0 + 0.3im;
        -1.0 - 0.30im;
        -0.9 + 0.4im;
        -0.9 - 0.4im;
        -0.7 + 0.7im;
        -0.7 - 0.7im;
        -0.4 + 0.9im;
        -0.4 - 0.9im;
        -0.0 + 1.10000000000im;
        -0.000000000000 - 1.10000000000im;
        1.2;
        1.00000000000000;
        0.9 + 0.4im;
        0.9 - 0.4im;
        0.6 + 0.600000000000im;
        0.6 - 0.600000000000im;
        0.4 + 0.900000000000im;
        0.4 - 0.900000000000im;
        0.00000000000 + 0.8im
        0.00000000000 - 0.8im
        ];
    z = [z ones(20)]
    p, PolyZeros(z)
end


text:
@testset "filters.jl" begin
    @testset "FullSF" begin
        μ = [0.1, 0.3, 0.5, -0.3]
        Σ = [   1.3  0.3  0.2  0.4;
                0.3  1.0 -0.4  0.8;
                0.2 -0.4  0.6  0.1;
                0.4  0.8  0.1  1.8;  ]

        fstate = FilterState(μ, Σ)

        smodel = OUModel(1., 1.)
        omodel = ExpGainModel(0.05, 0.1)

        x = [1, 1, 0, 1]
        y = 1
        obs = NeuronObs(x, y)
        dt = 1e-3

        filter = FullSF(4, smodel, omodel)

        update!(fstate, filter, obs, dt)

        @test fstate.μ ≈ [0.2998895344899539, 0.5096890112144515, 0.4895005232755023, 0.00028430173493088073]
        @test fstate.Σ ≈ [1.2993979068979908 0.2993978022428903 0.19960010465510047 0.3991968603469862; 0.2993978022428903 0.9999976923550348 -0.3991998901121445 0.7983967033643355; 0.19960010465510047 -0.3991998901121445 0.600799994767245 0.0998001569826507; 0.3991968603469862 0.7983967033643355 0.0998001569826507 1.7983952905204792]
    end#FullSF

    @testset "DiagSF" begin
        using LinearAlgebra: diag, Diagonal

        for dim in 10:10:1000
            μ = rand(dim)
            μ2 = copy(μ)
            Σ = rand(dim)
            Σ2 = Matrix(Diagonal(Σ))

            fstate = FilterState(μ, Σ)
            fstate2 = FilterState(μ2, Σ2)

            smodel = OUModel(1., 1.)
            omodel = ExpGainModel(0.05, 0.1)

            x = rand(dim)
            y = 1
            obs = NeuronObs(x, y)
            dt = 1e-3

            filter = DiagSF(dim, smodel, omodel)
            filter2 = FullSF(dim, smodel, omodel)

            update!(fstate, filter, obs, dt)
            update!(fstate2, filter2, obs, dt)

            @test fstate.μ ≈ fstate2.μ
            @test maximum(abs.(fstate.Σ .- diag(fstate2.Σ))) < 1e-6
        end
    end#DiagSF

    @testset "BlockSF" begin
        for blocksize in [2, 4, 8, 16, 32]
            numblocks = 5
            dim = blocksize * numblocks

            μ = rand(dim)
            μ2 = copy(μ)

            Σ = rand(blocksize, blocksize, numblocks)
            Σ2 = zeros(dim, dim)

            for i in 1:numblocks, j in 1:blocksize, k in 1:blocksize
                Σ2[(i-1)*blocksize + j, (i-1)*blocksize + k] = Σ[j, k, i]
            end

            fstate = FilterState(μ, Σ)
            fstate2 = FilterState(μ2, Σ2)

            smodel = OUModel(1., 1.)
            omodel = ExpGainModel(0.05, 0.1)

            x = rand(dim)
            y = 1
            obs = NeuronObs(x, y)
            dt = 1e-3

            filter = BlockSF(numblocks, blocksize, smodel, omodel)
            filter2 = FullSF(dim, smodel, omodel)

            update!(fstate, filter, obs, dt)
            update!(fstate2, filter2, obs, dt)

            @test fstate.μ ≈ fstate2.μ

            for i in 1:numblocks, j in 1:blocksize, k in 1:blocksize
                @test fstate2.Σ[(i-1)*blocksize+j, (i-1)*blocksize+k] ≈ fstate.Σ[j, k, i]
            end
        end
    end#BlockSF

    if RUN_BENCHMARKS
    @testset "Benchmark FullSF" begin
        using SynapticFilter: _filter_update!

        dim = 1024
        μ  = rand(dim)
        Σ  = rand(dim, dim)
        x  = rand(dim)
        y  = 1
        g0 = 0.1
        β  = 0.01
        τ  = 1.
        σs = 1.
        dt = 1e-3

        println("")
        println("Benchmarking one filter update step for FullSF")
        display(@benchmark _filter_update!($μ, $Σ, $τ, $σs, $g0, $β, $x, $y, $dt))
        println("")
        println("")
    end#Benchmark FullSF

    @testset "Benchmark BlockSF" begin
        using SynapticFilter: _filter_update!

        dim = 1024
        numblocks = 128
        blocksize = 8
        μ  = rand(dim)
        Σ  = rand(blocksize, blocksize, numblocks)
        x  = rand(dim)
        y  = 1
        g0 = 0.1
        β  = 0.01
        τ  = 1.
        σs = 1.
        dt = 1e-3

        println("")
        println("Benchmarking one filter update step for BlockSF")
        display(@benchmark _filter_update!($μ, $Σ, $τ, $σs, $g0, $β, $x, $y, $dt))
        println("")
        println("")
    end#Benchmark BlockSF

    @testset "Benchmark DiagSF" begin
        using SynapticFilter: _filter_update!

        dim = 1024
        μ  = rand(dim)
        Σ  = rand(dim)
        x  = rand(dim)
        y  = 1
        g0 = 0.1
        β  = 0.01
        τ  = 1.
        σs = 1.
        dt = 1e-3

        println("")
        println("Benchmarking one filter update step for DiagSF")
        display(@benchmark _filter_update!($μ, $Σ, $τ, $σs, $g0, $β, $x, $y, $dt))
        println("")
        println("")
    end#Benchmark DiagSF

    @testset "Benchmark FullSF on GPU" begin
        using SynapticFilter: _filter_update!
        using CUDA

        if CUDA.functional()
            CUDA.allowscalar(false)
            dim = 1024
            μ  = CUDA.rand(dim)
            Σ  = CUDA.rand(dim, dim)
            x  = CUDA.rand(dim)
            y  = 1
            g0 = 0.1f0
            β  = 0.01f0
            τ  = 1f0
            σs = 1f0
            dt = 1f-3

            println("")
            println("Benchmarking one filter update step for FullSF on GPU")
            display(
                @benchmark CUDA.@sync _filter_update!($μ, $Σ, $τ, $σs, $g0, $β, $x, $y, $dt)
            )
            println("")
            println("")
        end
    end#Benchmark FullSF on GPU
    end#end if
end#filters.jl


text:
using Polynomials
using StatsBase
using Random

export dcca, log_space

"""
    log_space(start, stop, num)

Returns a linear spacing in loglog scale. The parameter `start` and `stop` define the range and `num` corresponds to the number of desired points.
"""
function log_space(start::Int,stop::Int,num::Int)
    spacing = map(x -> round(Int,x), exp10.(range(log10(start), stop=log10(stop), length=num)))
    return sort(unique(spacing))
end

"""
    partitioning(x, box_size; overlap = div(box_size,2))

Partitions `x` into several segment of length `box_size`. The default behavior is that the windowed data have a 50% overlap.
To change this, change the `overlap` parameter.

returns an array of size (n,box_size), n being the total number of segments.
"""
function partitioning(x::Array{Float64,1},box_size::Int64; overlap::Int = div(box_size,2))
    nb_windows = div(length(x) - box_size, box_size - overlap) + 1
    partitionned_data = zeros(nb_windows, box_size)
    compteur = 1
    for i in 0:box_size - overlap:(length(x) - box_size)
        for j in 1:box_size
            partitionned_data[compteur,j] = x[i+j]
        end
        compteur += 1
    end
    return partitionned_data
end

"""
    detrending(x; order = 1)

Performs a linear detrending of `x`. You can change the order of the polynomials to a higher order for a non-linear detrending.
"""
function detrending1(values::Array{Float64,1}; order::Int = 1)
    position = collect(1:length(values))
    fit = polyfit(position,values,order)
    return values .-  polyval(fit,position)
end

function detrending(values::Array{Float64,1}; order::Int = 1)
    x = collect(1:length(values))
    X = hcat(ones(length(values)),x)
    A = (X'*X)\X'*values
    return values .- (x.*A[2] .+ A[1])
end

function detrending_optimized(values::Array{Float64,1}, x::Array{Float64}, X::Array{Float64}; order::Int = 1)
    A = (X'*X)\X'*values
    return values .- (x.*A[2] .+ A[1])
end

function integrate(x::Array{Float64,1})
    return cumsum(x)
end

"""
    dcca(x,y; box_start = 3, box_stop = div(length(x),10), nb_pts = 30)

Performs the DCCA analysis of `x` and `y`. The default analysis starts with a window size of 3 up to one tenth of the total length of `x` for statistical reasons.

returns the dcca coefficients.
"""
function dcca(x::Array{Float64,1},y::Array{Float64,1};
     box_start::Int = 3, box_stop::Int = div(length(x),10), nb_pts::Int = 30)
    if box_start < 3
        print("ERROR : size of windows must be greater than 3")
        return
    end
    window_sizes = log_space(box_start,box_stop,nb_pts)
    rho_DCCA = zeros(length(window_sizes))
    ffi, ff1i, ff2i = 0, 0, 0
    for (index,i) in enumerate(window_sizes)
        xi = partitioning(integrate(x),i)
        yi = partitioning(integrate(y),i)
        n = size(xi,1)
        m = size(xi,2)
        x = Array{Float64}(collect(1:m))
        X = hcat(ones(m),x)

        for j in 1:n
            #testing another version of detrending
            seg_x = detrending_optimized(xi[j,:], x, X)
            seg_y = detrending1(yi[j,:])
            ffi += (1/n)*((1/i)*seg_x'seg_y)
            ff1i += (1/n)*((1/i)*seg_x'seg_x)
            ff2i += (1/n)*((1/i)*seg_y'seg_y)
        end
        rho_DCCA[index] = ffi/(sqrt(ff1i)*sqrt(ff2i))
    end
    return  rho_DCCA
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
# Fractal Surface Interpolation demo

using Makie
using Random
using FractalTools
using GeometryBasics
using LinearAlgebra

# Get interpolation data
trig  = Triangle(
    Point(BigFloat(-1.), BigFloat(-1)),
    Point(BigFloat(1.), BigFloat(-1.)),
    Point(BigFloat(0.), BigFloat(2.))
)
f(x, y) = x^2 + y^2 + 1
pnts3d = getdata(f, trig, 100)

# Construct test points
tpnts3d = getdata(f, trig, 100)
tpnts = project(tpnts3d)

α = Node(0.001)
niter = Node(1)
err = lift(niter, α) do n, α
    itp = interpolate(pnts3d, α=α, maxiters=n)
    fval = map(pnt -> f(pnt...), tpnts)
    ival = map(pnt -> itp(pnt...), tpnts)
    val = abs.(fval - ival) ./ abs.(fval) * 100
    val
end

fig = Figure(title="Fractal Interpolation")
ax = Axis(fig[1,1], xlabel="npoints", ylabel="Relative Error (%)", title="Fractal Surface Interpolation Error Demo")
plt = lines!(ax, err)

gl = GridLayout()
gl2 = GridLayout()
gl[1, 1] = slider = Slider(fig, range=0.001 : 0.0001 : 0.01, startvalue=0.001, tellwidth=false)
gl2[1, 1] = button1 = Button(fig, label="Increment", tellwidth=false)
gl2[1, 2] = button2 = Button(fig, label="Decrement", tellwidth=false)
gl[2, 1] = gl2
gl[1, 2] = label1 = Label(fig, "alpha = $(slider.value[])")
gl[2, 2] = label2 = Label(fig, "niter = $(niter[])")
fig[2, 1] = gl

on(button1.clicks) do n
    niter[] += 1
    label2.text[] = "niter = " * string(niter[])
    autolimits!(ax)
end

on(button2.clicks) do n
    niter[] -= 1
    label2.text[] = "niter = " * string(niter[])
    autolimits!(ax)
end

on(slider.value) do val
    α[] = val
    label1.text[] = "alpha = " * string(val)
    autolimits!(ax)
end

display(fig)

# # For video recording.
# fps = 60
# record(fig.scene, "demo.mp4"; framerate = fps) do io
#     for i = 1 : 50 * fps
#         sleep(1/fps)
#         recordframe!(io)
#     end
# end


text:
using Scallop.Underlying

@assert 12 == xunderlying(6)


text:
# this is a annotation
#definiton a function
#only use in console
#v(x) = (print(x);x)
#v(x)
function f(x,y)
    x + y
end
print(f(3,2))
print("\n")
print(typeof(f))
#type limit
function func(x::Int,explaintion::String)
    print(x)
    print(explaintion)
    # this is return value
    x
end
func(3,"fdsf")
#throw error
#func("fsdf",3)

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
using PolyChaos, Test, StaticArrays

deg = 5
op = genLaguerreOrthoPoly(deg,1.23)

n = 10
X = [ exp(3 + rand()), exp.(3 .+ rand(n)) ]

α, β, αs, βs = op.α, op.β, SVector(op.α...), SVector(op.β...)
combinations = Iterators.product([α, αs], [β, βs])

@testset "evaluate univariate basis at point(s)" begin
    for x in X
        [ @test evaluate(collect(0:deg),x,op) == evaluate(collect(0:deg), x, a, b)  for (a, b) in combinations ]
        [ @test evaluate(SVector(0:deg...),x,op) == evaluate(SVector(0:deg...), x, a, b)  for (a, b) in combinations ]
        [ @test evaluate(0:deg,x,op) == evaluate(0:deg, x, a, b)  for (a, b) in combinations ]
        [ @test evaluate(0:2:deg,x,op) == evaluate(0:2:deg, x, a, b)  for (a, b) in combinations ]
        for d in 0:deg
            if typeof(x) <: Real
                [ @test evaluate(d,x,a,b) == evaluate(d,x,op) == evaluate(d,[x],a,b) for (a, b) in combinations ]
            else
                [ @test evaluate(d,x,a,b) == evaluate(d,x,op) for (a, b) in combinations ]
            end
        end

        res = hcat(map(i -> evaluate(i,x,op.α,op.β),collect(0:op.deg))...)
        if typeof(x) <: Real
            [ @test evaluate(x,op) == evaluate([x], op) == evaluate(collect(0:op.deg),[x],op) == evaluate(collect(0:op.deg),[x],a,b) == res for (a, b) in  combinations ]
            [ @test evaluate(x,op) == evaluate([x], op) == evaluate(0:op.deg,[x],op) == evaluate(0:op.deg,[x],a,b) == res for (a, b) in  combinations ]
        else
            [ @test evaluate(x, op) == evaluate(collect(0:op.deg),x,op) == evaluate(collect(0:op.deg),x,a,b) == res for (a, b) in combinations ]
            [ @test evaluate(x, op) == evaluate(0:op.deg,x,op) == evaluate(0:op.deg,x,a,b) == res for (a, b) in combinations ]
        end
    end
end

ops = [  genLaguerreOrthoPoly(deg,1.23),
        genHermiteOrthoPoly(deg,3.),
        LegendreOrthoPoly(deg),
        JacobiOrthoPoly(deg,4.3,10.)
     ]

mop = MultiOrthoPoly(ops,deg)

nop = length(ops)
X = [ exp.(3 .+ rand(nop)) , exp.(3 .+ rand(n,nop)) ]

α, β = coeffs(mop)
αs, βs = SVector(α...), SVector(β...)
combinations = Iterators.product([α, αs], [β, βs])

@testset "evaluate multivariate basis at point(s)" begin
    for x in X
        [ @test evaluate(mop.ind,x,a,b) == evaluate(mop.ind,x,mop) == evaluate(x,mop) for (a, b) in combinations ]
        if typeof(x) <: AbstractVector{<:Real}
            [ @test evaluate(mop.ind,x,a,b) == evaluate(mop.ind,x,mop) == evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...) for (a, b) in combinations]
            [ @test evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),mop) == evaluate(mop.ind,reshape(x,1,length(x)),a,b) for (a, b) in combinations ]
        else
            [ @test evaluate(x,mop) == evaluate(mop.ind,x,mop) == evaluate(mop.ind,x,a,b) for (a, b) in combinations ]
        end
        if VERSION >= VersionNumber("1.1.1")
            for d in eachrow(mop.ind[1:3,:])
                [ @test evaluate(collect(d),x,a,b) == evaluate(collect(d),x,mop) for (a, b) in combinations ]
                [ @test evaluate(d,x,a,b) == evaluate(collect(d),x,mop) for (a, b) in combinations ]
            end
        end
    end
end


text:
# Tests concerning cross-entropy

arr21 = rand(100)
arr22 = rand(100)
arr41 = rand(10000)
arr42 = rand(10000)

# The cross-entropy between two copies of an array
@test cross_entropy(arr41, arr41) ≈ log2(100) atol = 0.01
println("Cross-entropy between array and itself passed.")

# The cross-entropy between an array and itself is equal to the regular entropy of that array
@test cross_entropy(arr21, arr21) == get_entropy(arr21)
println("Equality of cross-entropy between array and itself and entropy passed.")

# The cross-entropy between two small samples from the same distribution will be somewhat larger than the regular entropy
@test cross_entropy(arr21, arr22) ≈ log2(10) atol = 0.4
println("Cross-entropy between two small sample arrays passed.")
@test cross_entropy(arr41, arr42) ≈ log2(100) atol = 0.02
println("Cross-entropy between two large sample arrays passed.")


text:
module  TestCross

using PopGenCore
using PopGenSims
using Test

cats = @nancycats ;

@testset "Crosses" begin
    @testset "sample genotypes" begin
        @test length(PopGenSims.sample_genotype(cats.genodata.genotype[30], 2)) == 2
        @test eltype(PopGenSims.sample_genotype(cats.genodata.genotype[30],1)) == eltype(cats.genodata.genotype[30])
        @test length(PopGenSims.sample_genotype(missing, 2)) == 1
        @test first(PopGenSims.sample_genotype(missing, 2)) === missing
    end

    @testset "crosses" begin
        f1 = cross(cats, "N111", "N107", n = 10)
        @test f1 isa PopData
        @test f1.metadata.samples == 10
        @test length(f1.genodata.name) == 90
        f1 = cross(cats, "N111", "N107", n = 10, generation = "firstgen")
        @test f1 isa PopData
        @test f1.metadata.samples == 10
        @test length(f1.genodata.name) == 90
        f2 = cross(cats => "N111", f1 => "firstgen_10", n = 10)
        @test f2 isa PopData
        @test f2.metadata.samples == 10
        @test length(f2.genodata.name) == 90
        f2 = cross(cats => "N111", f1 => "firstgen_10", n = 10, generation = "F2")
        @test f2 isa PopData
        @test f2.metadata.samples == 10
        @test length(f2.genodata.name) == 90
    end
end
end # module TestCross


text:
@testset "dynamic_distribution_model" begin
    f(s, a) = (nextstate = 1, reward = 1.0, prob = 1.0)
    m = DynamicDistributionModel(f, 2, 2)

    @test length(observation_space(m)) == 2
    @test length(action_space(m)) == 2

    for s = 1:2
        for a = 1:2
            @test m(s, a) == f(s, a)
        end
    end
end