text:
"""
    batch_matvec(A::Tensor, X::Matrix)

for A ∈ ℝ^{m × d × n}, X ∈ ℝ^{d × n}. Performs ``n`` matrix-vec multiplications

    A[:,:,i] * X[:,i]

for ``i \\in 1,\\ldots,n``. This is performed via expansion of `X` and can be
efficiently evaluated using BLAS operations, which (I have assumed!) will often
be faster than a loop, especially in relation to AD and GPUs. Cf. `multitask.jl`
"""
batch_matvec(A::AbstractArray{T,3}, X::AbstractArray{T,2}) where T =
    dropdims(sum(A .* unsqueeze(X, 1), dims=2), dims=2)

"""
    batch_matmul(A::MatrixOrArray, X::Array)

for A ∈ ℝ^{n × m × d}, X ∈ ℝ^{m × r × d}. Performs ``n`` matrix multiplications

    A[:,:,i] * X[:,:,i]

for ``i \\in 1,\\ldots,n``. This is performed via expansion of `X` and can be
efficiently evaluated using BLAS operations, which (I have assumed!) will often
be faster than a loop in relation to AD and GPUs.

We also permit A ∈ ℝ^{n × m}, X ∈ ℝ^{m × r × d} => ℝ^{m × r × d}.
"""
batch_matmul(A::AbstractArray{T,3}, X::AbstractArray{T,3}) where T =
    dropdims(sum(unsqueeze(A,3) .* unsqueeze(X, 1), dims=2), dims=2)

batch_matmul(A::AbstractArray{T,2}, X::AbstractArray{T,3}) where T =
    dropdims(sum(unsqueeze(unsqueeze(A,3),4) .* unsqueeze(X, 1), dims=2), dims=2)

ensure_matrix(x::AbstractMatrix) = x
ensure_matrix(x::AbstractVector) = unsqueeze(x, 2)

ensure_tensor(x::AbstractMatrix) = unsqueeze(x, 3)
ensure_tensor(x::AbstractVector) = reshape(x, length(x), 1, 1)
ensure_tensor(x::AbstractArray{T,3}) where T = x

eye(d) = Matrix(I, d, d)

################################################################################
##                                                                            ##
##                         Special Matrix Constructors                        ##
##          *** THIS IS TAKEN DIRECTLY FROM MY JULIA UTILS PKG,               ##
##          BUT COPIED HERE FOR COMPLETENESS / TO SATISFY ANONYMITY ***       ##
##                                                                            ##
################################################################################

#= have a CuArray version in CUDA file. Cannot make more generally available
   because CuArrays.jl will not compile (e.g. on my machine) if no GPU. See
   @requires in [REDACTED].jl and CUDA.jl files in this project. =#

"""
    make_lt(x::AbstractVector{T}, d::Int)

Make lower triangular matrix (``d \\times d``) from vector `x`.
"""
function make_lt(x::AbstractVector{T}, d::Int)::Array{T,2} where T <: Real
    @argcheck (length(x) == Int(d*(d+1)/2))
    make_lt!(zeros(T, d, d), x, d)
end

function make_lt!(out::AbstractMatrix{T}, x::AbstractArray{T,1}, d::Int)::Array{T,2} where T <: Real
    x_i = 1
    for j=1:d, i=j:d
        out[i,j] = x[x_i]
        x_i += 1
    end
    return out
end

# Gradients
function unmake_lt(M::AbstractMatrix{T}, d)::Array{T,1} where T <: Real
    return M[tril!(trues(d,d))]
end

make_lt(x::TrackedArray, d::Int) = Tracker.track(make_lt, x, d)

@grad function make_lt(x, d::Int)
    return make_lt(Tracker.data(x), d), Δ -> (unmake_lt(Δ, d), nothing)
end

"""
    make_lt_strict(x::AbstractVector{T}, d::Int)

Make *strictly* lower triangular matrix (``d \\times d``) from vector `x`.
"""
function make_lt_strict(x::AbstractVector{T}, d::Int)::Array{T,2} where T <: Real
    @argcheck (length(x) == Int(d*(d-1)/2))
    return make_lt_strict!(zeros(T, d, d), x, d)
end

function make_lt_strict!(out::AbstractMatrix{T}, x::AbstractVector, d::Int)::Array{T,2} where T <: Real
    x_i = 1
    for j=1:d-1, i=j+1:d
        out[i,j] = x[x_i]
        x_i += 1
    end
    return out
end

make_lt_strict(x::TrackedArray, d::Int) = Tracker.track(make_lt_strict, x, d)

@grad function make_lt_strict(x, d::Int)
    return make_lt_strict(Tracker.data(x), d), Δ -> (unmake_lt_strict(Δ, d), nothing)
end

function unmake_lt_strict(M::AbstractMatrix{T}, d::Int)::Array{T,1} where T <: Real
    return M[tril!(trues(d,d), -1)]
end

"""
    make_skew(x::AbstractArray{T,1}, d::Int)

Make ``d \\times d`` skew-symmetric matrix from vector `x`.
"""
function make_skew(x::AbstractArray{T,1}, d::Int)::AbstractArray{T,2} where T <: Real
    S = make_lt_strict(x, d)
    return S - S'
end

function cayley_orthog(x::AbstractArray{T,1}, d)::AbstractArray{T,2} where T <: Real
    S = make_skew(x, d)
    I = eye(d)
    return (I - S) / (I + S)  # (I - S)^{-1}(I + S). Always nonsingular.
end

function inverse_cayley_orthog(Q::AbstractMatrix{T}, d::Int=size(Q,1)) where T <: Real
    unmake_lt_strict((I - Q ) / (I + Q), d)
end

"""
    diag0(x::Array{T,1})
Make diagonal matrix from vector `x`. This may have been superseded by `Diag`,
but is here for legacy reasons since Flux's version used to result in Tracked{Tracked}.
"""
function diag0(x::Array{T,1})::Array{T,2} where T <: Real
    d = length(x)
    M = zeros(T, d,d)
    M[diagind(M)] = x
    return M
end

diag0(x::TrackedArray) = Tracker.track(diag0, x)

@grad function diag0(x)
    return diag0(Tracker.data(x)), Δ -> (Δ[diagind(Δ)],)
end

"""
    LDSCell_diag_u{U,V,W}
Simple LDS struct containing three fields: a (transition coeffs), B (input
matrix) and h, the initial state.
"""
mutable struct LDSCell_diag_u{U,V,W}
    a::U
    B::V
    h::W
end

# Operation
function (m::LDSCell_diag_u)(h, x)
    h = m.a .* h + m.B * x
    return h, h
end

Flux.hidden(m::LDSCell_diag_u) = m.h
Flux.@treelike LDSCell_diag_u

"""
    LDSCell_diag_batch_u{U,V,W}
Simple LDS struct containing three fields: a (transition coeffs), B (input
matrix) and h, the initial state. Permits 3D arrays for a, B and performs
batch updates on the inputs.
"""
mutable struct LDSCell_diag_batch_u{U,V,W}
    a::U
    B::V
    h::W
end

# Operation
function (m::LDSCell_diag_batch_u)(h, x)
    h = m.a .* h + batch_matvec(m.B, x)
    return h, h
end

Flux.hidden(m::LDSCell_diag_batch_u) = m.h
Flux.@treelike LDSCell_diag_batch_u

"""
    LDSCell_simple_u{U,V,W}
Simple LDS struct containing three fields: A, B (for transition and input
matrices) and h, the initial state.
"""
mutable struct LDSCell_simple_u{U,V,W}
    A::U
    B::V
    h::W
end

# Operation
function (m::LDSCell_simple_u)(h, x)
    h = m.A * h + m.B * x
    return h, h
end

Flux.hidden(m::LDSCell_simple_u) = m.h
Flux.@treelike LDSCell_simple_u

"""
    LDSCell_batch_u{U,V,W}
Simple LDS struct containing three fields: A, B (for transition and input
matrices) and h, the initial state. Permits 3D arrays for A, B and performs
batch updates on the inputs.
"""
mutable struct LDSCell_batch_u{U,V,W}
    A::U
    B::V
    h::W
end

# Operation
function (m::LDSCell_batch_u)(h, x)
    h = batch_matvec(m.A, h) + batch_matvec(m.B, x)
    return h, h
end

Flux.hidden(m::LDSCell_batch_u) = m.h
Flux.@treelike LDSCell_batch_u

################################################################################
##                                                                            ##
##                         Model specific utilities                           ##
##                                                                            ##
################################################################################

"""
    A(ψ, d)
Construct transition matrix `A` using reduced Cayley form from parameter vector
ψ (extracting relevant elements).
"""
function A(ψ, d)
    n_skew = Int(d*(d-1)/2)
    x_S, x_V = ψ[1:d], ψ[d+1:d+n_skew]
    V = cayley_orthog(x_V/10, d)
    S = diag0(σ.(x_S))
    return S * V
end

"""
    B(ψ, d)
Construct transition matrix `B` using sparse construction from parameter vector
ψ (extracting relevant elements).
"""
function B(ψ::AbstractVector, d, n_u)
    n_skew = Int(d*(d-1)/2)
    i_init = d+n_skew
    ψ_b1 = reshape(ψ[(i_init+1):(i_init+n_u*d)], d, n_u)
    ψ_b2 = reshape(ψ[(i_init + n_u*d +1):(i_init + 2*n_u*d)], d, n_u)
    return σ.(ψ_b1) .* tanh.(ψ_b2)
end

function B(ψ::AbstractMatrix, d, n_u)
    n_skew = Int(d*(d-1)/2)
    n_b = size(ψ, 2)
    i_init = d+n_skew
    ψ_b1 = reshape(ψ[(i_init+1):(i_init+n_u*d), :], d, n_u, n_b)
    ψ_b2 = reshape(ψ[(i_init + n_u*d +1):(i_init + 2*n_u*d), :], d, n_u, n_b)
    return σ.(ψ_b1) .* tanh.(ψ_b2)
end

"""
    C(ψ, d)
Construct emission matrix `C` with no special constructor from parameter vector
ψ (extracting relevant elements).
"""
function C(ψ::AbstractVector, d, n_u, n_out)
    n_skew = Int(d*(d-1)/2)
    i_init = d + n_skew + 2*n_u*d
    ψ_c= ψ[(i_init+1):(i_init+n_out*d)]
    return reshape(ψ_c, n_out, d)
end

function C(ψ::AbstractMatrix, d, n_u, n_out)
    n_skew = Int(d*(d-1)/2)
    n_b = size(ψ, 2)
    i_init = d + n_skew + 2*n_u*d
    ψ_c= ψ[(i_init+1):(i_init+n_out*d), :]
    return reshape(ψ_c, n_out, d, n_b)
end

"""
    dmt(ψ, d)
Construct emission bias `d` with no special constructor from parameter vector
ψ (extracting relevant elements).
"""
function dmt(ψ::AbstractVector, d, n_u, n_out)
    n_skew = Int(d*(d-1)/2)
    i_init = d + n_skew + 2*n_u*d +n_out*d
    return ψ[(i_init+1):(i_init+n_out)]
end

function dmt(ψ::AbstractMatrix, d, n_u, n_out)
    n_skew = Int(d*(d-1)/2)
    i_init = d + n_skew + 2*n_u*d +n_out*d
    return ψ[(i_init+1):(i_init+n_out), :]
end


text:
# averaging.jl
# Small module to test lower bounds on error in Manhattan
# Authored by Arthur Delarue on 11/30/15

manhattan = loadCityGraph()
graph = manhattan.network
travel_times, num_rides, trainTestDf, testing_travel_times, testing_num_rides, testDf = loadInputTravelTimes(manhattan.positions, "average", year=2013, startTime=12, endTime=14, startMonth=1, endMonth=12, radius = 200)
nodeTree, nodePairs = buildNodeKDTree(manhattan.positions)

function average_manhattan_results(fileList)
	travelTimes = zeros(6134, 6134)
	for (i,fileName) in enumerate(fileList)
		newTimes = load(fileName, "times")
		new_graph, new_edge_dists, new_nodes, new_edge_isExpensive = modifyGraphForDijkstra(graph, newTimes, manhattan.positions, turn_cost=0.0)
		old_nodes = getInverseMapping(new_nodes, nv(new_graph))
		new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)
		if i == 1
			travelTimes = new_sp.traveltime
		else
			travelTimes = (travelTimes .* (i - 1) + new_sp.traveltime) ./ i
		end
	end
	return travelTimes
end

newTravelTimes = average_manhattan_results([
	"Outputs/fast_r200_minr1_i50_wd_1214_m112_y2013_relaxed_smooth_ppc3_dataaverage_errorboth_dynConstr_st50000_add0_thr0.7_rem0_thr0.3_every1_gcu_tc0.0/manhattan-times-31.jld",
	"Outputs/fast_r200_minr1_i50_wd_1214_m112_y2013_relaxed_smooth_ppc3_dataaverage_errorboth_dynConstr_st50000_add0_thr0.7_rem0_thr0.3_every2_gcu_tc0.0/manhattan-times-50.jld",
	"Outputs/fast_r200_minr1_i50_wd_1214_m112_y2013_relaxed_smooth_ppc3_dataaverage_errorboth_dynConstr_st50000_add0_thr0.7_rem0_thr0.3_every3_gcu_tc0.0/manhattan-times-50.jld",
	"Outputs/fast_r200_minr1_i50_wd_1214_m112_y2013_relaxed_smooth_ppc3_dataaverage_errorboth_dynConstr_st50000_add0_thr0.7_rem0_thr0.3_every4_gcu_tc0.0/manhattan-times-50.jld"])

avg_sq_err, avg_rel_err, avg_bias = computeTestingError(trainTestDf, nodeTree, nodePairs, newTravelTimes, "average")

ext: jl

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
using LogitNormals
using Documenter

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

makedocs(;
    modules=[LogitNormals],
    authors="Thomas Wutzler <twutz@bgc-jena.mpg.de> and contributors",
    repo="https://github.com/bgctw/LogitNormals.jl/blob/{commit}{path}#{line}",
    sitename="LogitNormals.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/LogitNormals.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bgctw/LogitNormals.jl",
    devbranch="main",
)


text:
# Variates and distribution parameters with physical dimensions.

using Unitful: Unitful, Quantity, Units, NoDims, NoUnits, unit, ustrip

"""
Type of the variate of a distribution that supports both pure `Real`s and real Unitful
Quantities. If a parameter of the distribution is dimensioned, it is also a `RealQuantity`.
```julia
const RealQuantity = Union{Real, Unitful.Quantity{<:Real}}
```
_i.e._ a `RealQuantity` is either a pure `Real`, or a Unitful `Quantity` with a `Real` as
numeric backing type.
"""
const RealQuantity = Union{Real, Unitful.Quantity{<:Real}}

"""
    unit(d::Distribution)

Retrieve the physical units of the distribution's variate, as a
[`Unitful.Units`](https://painterqubits.github.io/Unitful.jl/stable/types/).

 - If the distribution does not have units attached (*i.e.* it is over real values), this
   returns `NoUnits`, which is a 0-size singleton that gets pruned away in mathematical
   expressions at compile-time.

 - If the distribution does have units attached, this is `unit(eltype(d))`.
"""
Unitful.unit(::D) where {D <: Distribution} = unit(eltype(D))


text:
module CUDArt_gen  # the generated module

using Compat

# In the runtime API, these are all used only inside Ptrs,
# so these typealiases are safe (if you don't need access to
# struct elements)
typealias cudaUUID_t Void

function checkerror(code::Cuint)
    if code == cudaSuccess
        return nothing
    end
    # Because try/finally may disguise the source of the problem,
    # let's show a backtrace here
    warn("CUDA error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    throw(bytestring(cudaGetErrorString(code)))
end

include("../gen-6.5/gen_libcudart_h.jl")

typealias cudaError_t cudaError

if is_windows()
    # location of cudart64_xx.dll or cudart32_xx.dll have to be in PATH env var
    # ex: C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v6.5\bin
    # (by default, it is done by CUDA toolkit installer)

    const dllnames = (Base.Sys.WORD_SIZE==64) ?
        ["cudart64_70", "cudart64_65", "cudart64_60", "cudart64_55", "cudart64_50", "cudart64_50_35"] :
        ["cudart32_70", "cudart32_65", "cudart32_60", "cudart32_55", "cudart32_50", "cudart32_50_35"]
    const libcudart = Libdl.find_library(dllnames, [""])
else # linux or mac
    const libdir = (Base.Sys.WORD_SIZE==64) ? "lib64" : "lib"
    const libcudart = Libdl.find_library(["libcudart", "cudart"], ["/usr/local/cuda/$libdir", "/usr/local/cuda-6.5/$libdir", "/usr/local/cuda-6.0/$libdir", "/usr/local/cuda-5.5/$libdir", "/usr/local/cuda-5.0/$libdir", "/usr/local/cuda-7.0/$libdir"])
end

if isempty(libcudart)
    error("CUDA runtime API library cannot be found")
end

include("../gen-6.5/gen_libcudart.jl")

# # Make cudaSetupArgument more convenient
# function cudaSetupArgument(arg, size, offset)
#   checkerror(ccall( (:cudaSetupArgument, libcudart), cudaError_t, (Ptr{Void}, Csize_t, Csize_t), &arg, size, offset))
# end

# Fix issues stemming from the inability to pass structs as args
const libwrapcuda = Libdl.find_library(["libwrapcuda"],[joinpath(dirname(@__FILE__), "..", "deps")])
isempty(libwrapcuda) && error("Cannot find libwrapcuda")

function cudaMalloc3D(p::Array{cudaPitchedPtr,1}, ext::cudaExtent)
    checkerror(ccall((:wrapcudaMalloc3D, libwrapcuda), cudaError_t, (Ptr{cudaPitchedPtr}, Ptr{cudaExtent}), p, &ext))
end
function cudaMemset3D(pitchedDevPtr::cudaPitchedPtr, value::UInt8, extent::cudaExtent)
  checkerror(ccall( (:wrapcudaMemset3D, libwrapcuda), cudaError_t, (Ptr{cudaPitchedPtr}, Cint, Ptr{cudaExtent}), &pitchedDevPtr, value, &extent))
end

# # Issue warning, not error, on bad cudaFree
# function cudaFree(devPtr)
#     code = ccall( (:cudaFree, libcudart), cudaError_t, (Ptr{Void},), devPtr)
#     if code != cudaSuccess
#         warn(string(bytestring(cudaGetErrorString(code)), " on ", devPtr))
#         Base.show_backtrace(STDOUT, backtrace())
#     end
# end

# Convenience constructors
cudaPitchedPtr() = cudaPitchedPtr(C_NULL, 0, 0, 0)
cudaPitchedPtr(pp::cudaPitchedPtr, xsize::Integer) = cudaPitchedPtr(pp.ptr, pp.pitch, xsize, pp.ysize)

end


text:
# SigmoidFit
function SigmoidFit(u,t,p=zeros(2))
  m = (x, p) -> @. p[1]/(1+exp(x-p[2]))
  c_f = curve_fit(m,t,u,p)
  Curvefit{true}(u,t,m,c_f)
end

# HillFit
function HillFit(u,t,p=zeros(2),n=nothing)
  if n == nothing
    push!(p,rand(1)[1])
    m = (x, p) -> @. p[1] * inv(1.0 + p[2]/(x^p[3]))
  else
    m = (x, p) -> @. p[1] * inv(1.0 + p[2]/(x^n))
  end
  c_f = curve_fit(m,t,u,p)
  Curvefit{true}(u,t,m,c_f)
end

# WeibullFit
function WeibullFit(u,t,p=zeros(3))
  m = weibull_fun
  c_f = curve_fit(m,t,u,p)
  Curvefit{true}(u,t,m,c_f)
end

# OneCompartmentPKFit
function OneCompartmentPKFit(u,t,p=zeros(5))
  m = (x,p) -> @. p[1] * p[2] * p[3] * (exp(-p[4]*x) - exp(-p[1]*x)) / (p[5] * (p[1] - p[4]))
  c_f = curve_fit(m,t,u,p)
  Curvefit{true}(u,t,m,c_f)
end


text:
import Base: size, length, eltype, show

using Oceananigans: AbstractGrid

"""
    VerticallyStretchedCartesianGrid{FT<:AbstractFloat, R<:AbstractRange, A<:AbstractArray} <: AbstractGrid{FT}

A Cartesian grid with with constant horizontal grid spacings `Δx` and `Δy`, and
non-uniform or stretched vertical grid spacing `Δz` between cell centers and cell faces.
"""
struct VerticallyStretchedCartesianGrid{FT<:AbstractFloat, R<:AbstractRange, A<:AbstractArray} <: AbstractGrid{FT}
    # Number of grid points in (x,y,z).
    Nx::Int
    Ny::Int
    Nz::Int
    # Halo size in (x,y,z).
    Hx::Int
    Hy::Int
    Hz::Int
    # Total number of grid points (including halo regions).
    Tx::Int
    Ty::Int
    Tz::Int
    # Domain size [m].
    Lx::FT
    Ly::FT
    Lz::FT
    # Grid spacing [m].
    Δx::FT
    Δy::FT
    ΔzF::A
    ΔzC::A
    # Range of coordinates at the centers of the cells.
    xC::R
    yC::R
    zC::A
    # Range of grid coordinates at the faces of the cells.
    # Note: there are Nx+1 faces in the x-dimension, Ny+1 in the y, and Nz+1 in the z.
    xF::R
    yF::R
    zF::A
end

function VerticallyStretchedCartesianGrid(FT=Float64, arch=CPU();
        size, halo=(1, 1, 1), length=nothing, x=nothing, y=nothing, z=nothing, zF=nothing)

    # Hack that allows us to use `size` and `length` as keyword arguments but then also
    # use the `size` and `length` functions.
    sz, len = size, length
    length = Base.length

    Lx, Ly, Lz, x, y, z = validate_grid_size_and_length(sz, len, halo, x, y, z)

    Nx, Ny, Nz = sz
    Hx, Hy, Hz = halo

    Tx = Nx + 2*Hx
    Ty = Ny + 2*Hy
    Tz = Nz + 2*Hz

    Δx = convert(FT, Lx / Nx)
    Δy = convert(FT, Ly / Ny)

    x₁, x₂ = convert.(FT, [x[1], x[2]])
    y₁, y₂ = convert.(FT, [y[1], y[2]])
    z₁, z₂ = convert.(FT, [z[1], z[2]])

    xF = range(x₁, x₂; length=Nx+1)
    yF = range(y₁, y₂; length=Ny+1)

    xC = range(x₁ + Δx/2, x₂ - Δx/2; length=Nx)
    yC = range(y₁ + Δy/2, y₂ - Δy/2; length=Ny)

    zF, zC, ΔzF, ΔzC = validate_and_generate_variable_grid_spacing(FT, zF, Nz, z₁, z₂)

    VerticallyStretchedCartesianGrid{FT, typeof(xF), typeof(zF)}(Nx, Ny, Nz, Hx, Hy, Hz, Tx, Ty, Tz, Lx, Ly, Lz,
                                                                 Δx, Δy, ΔzF, ΔzC, xC, yC, zC, xF, yF, zF)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Exiv2_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Exiv2")
JLLWrappers.@generate_main_file("Exiv2", UUID("d2dc9c37-b624-5257-85b1-2eb67ce3a11b"))
end  # module Exiv2_jll


text:
using MaximumLikelihoodPower
const MLE = MaximumLikelihoodPower
import MaximumLikelihoodPower.Example

using Test

include("testmle.jl")


text:
using Test
using TestSetExtensions

using SparseAlgebra

@info "Testing started"

@testset "All the tests" begin

    @includetests [:sparsik_tests, :wiedemannchik_tests,
                    :dok_sparsik_tests,
                    :find_basis_tests]

end

@info "All (the) tests OK"


text:
#На прямоугольном поле произвольных размеров расставить маркеры в виде "шахматных" клеток, начиная с юго-западного угла поля, когда каждая отдельная "шахматная" клетка имеет размер n x n клеток поля (n - это параметр функции).
#Начальное положение Робота - произвольное, конечное - совпадает с начальным. Клетки на севере и востоке могут получаться "обрезанными" - зависит от соотношения размеров поля и "шахматных" клеток.
function chess(r, n)
    z = 0
    m = 0
    z, m = corner(r)
    side = Ost;
    x = 0;
    y = 0;
    if ((x ÷ n) + (y ÷ n)) % 2 == 0
        putmarker!(r)
    end
    while isborder(r,Nord) == false
        x, y = chess_line!(side, n, x, y, r)
        move!(r,Nord);
        y = y + 1;
        side = revers(side);
        put!(r, x, y, n)
    end
    x, y = chess_line!(side, n, x, y, r)
    go_home!(r, z, m)
end
function movements_mark!(r,side)
    while (isborder(r,side) ==  false) && (ismarker(r) == true)
        move!(r,side)
    end
end
function revers(side)
    side = HorizonSide((Int(side)+2)%4)
    return side
end
function count_move!(r,side,count)
    for i = 1:count
        move!(r,side)
    end
end
function movements_to_side!(r,side)
    count = 0
    while isborder(r,side) ==  false
        move!(r,side)
        count = count + 1
    end
    return count
end
function corner(r)
    x = 0
    y = 0
    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        y = y + movements_to_side!(r, Sud)
        x = x + movements_to_side!(r, West)
    end
    return x, y
end
function go_home!(r, x, y)
    corner(r)
    count_move!(r, Nord, y)
    count_move!(r, Ost, x)
end
function chess_line!(side, n, x, y, r)
    while isborder(r, side) == false
        move!(r, side)
        if side == Ost
            x = x + 1;
        else
            x = x - 1;
        end

        put!(r, x, y, n)
    end
    return x, y
end
function put!(r, x, y, n)
    if ((x ÷ n) + (y ÷ n)) % 2 == 0
        putmarker!(r)
    end
end
n = 3
chess(r,n)


text:
import Distributions
import MLJBase
import MLJModelInterface

struct MixedVariate <: Distributions.VariateForm end
struct MixedSupport <: Distributions.ValueSupport end

abstract type AbstractSossMLJPredictor{M} <: Distributions.Distribution{MixedVariate, MixedSupport}
end

struct SossMLJPredictor{M,PO,PR,A} <: AbstractSossMLJPredictor{M}
    model::M
    post::PO
    pred::PR
    args::A
end

abstract type AbstractSossMLJModel{P, M} <: MLJModelInterface.JointProbabilistic
end

mutable struct SossMLJModel{P, M, H, I, R, T} <: AbstractSossMLJModel{P, M}
    # the first type parameter P is for the predictor_type
    # then we do the model, which is of type M
    model::M
    # after the model, we put the rest of the fields in alphabetical order for simplicity
    hyperparams::H
    infer::I
    response::R
    transform::T
end

@inline function default_transform(tbl)
    return (X = MLJModelInterface.matrix(tbl),)
end

function SossMLJModel(;
                      predictor::Type{P} = SossMLJPredictor,
                      model::M,
                      hyperparams::H = NamedTuple(),
                      infer::I = Soss.dynamicHMC,
                      response::R = :y,
                      transform::T = default_transform) where P where M where H where I where R where T
    result = SossMLJModel{P, M, H, I, R, T}(
        model,
        hyperparams,
        infer,
        response,
        transform,
    )
    return result
end

struct RMSDistribution <: MLJBase.Measure end
struct RMSExpected <: MLJBase.Measure end
struct RMSMedian <: MLJBase.Measure end


text:
module BacenSGS

using HTTP, Dates, DataFrames,JSON3, Gumbo, StringEncodings, ProgressLogging

const SEARCH_URL = Dict("pt"=>"https://www3.bcb.gov.br/sgspub/index.jsp?idIdioma=P", "en"=>"https://www3.bcb.gov.br/sgspub/index.jsp")

function timeserie_value(indicador::Integer;dtInicial::Dates.Date=Date(2018),dtFinal=missing)
    linkSgs = "http://api.bcb.gov.br/dados/serie/bcdata.sgs.$(indicador)/dados?formato=json&dataInicial=$(Dates.format(dtInicial,dateformat"dd/mm/yyyy"))&dataFinal=$(Dates.format(coalesce(dtFinal,Dates.today()),dateformat"dd/mm/yyyy"))"
    resp = HTTP.get(linkSgs,require_ssl_verification = false)
    if resp.status != 200
        return DataFrame()
    end
    r = resp.body |> String |> JSON3.read |> DataFrame
    r[!,:date] = r.data .|> x-> Dates.Date(x,dateformat"dd/mm/yyyy")
    r[!,:indicador] .= indicador
    r[!,:value] = r.valor .|> x->0.01parse(Float64,x)
    select(r,Not([:data,:valor]))
end

get_tipo_busca(codigo::Integer) = ("localizarSeriesPorCodigo",4,codigo,nothing)
get_tipo_busca(str::AbstractString) = ("localizarSeriesPorTexto",6,nothing,str)
function _get_params(searchMethod,tipo,codigo,texto)
    Dict(
        "method"=>searchMethod,
        "periodicidade"=>4,
        "codigo"=>codigo,
        "fonte"=>341,
        "texto"=>texto,
        "hdFiltro"=>nothing,
        "hdOidGrupoSelecionado"=>nothing,
        "hdSeqGrupoSelecionado"=>nothing,
        "hdNomeGrupoSelecionado"=>nothing,
        "hdTipoPesquisa"=>tipo,
        "hdTipoOrdenacao"=>0,
        "hdNumPagina"=>nothing,
        "hdPeriodicidade"=>"Todas",
        "hdSeriesMarcadas"=>nothing,
        "hdMarcarTodos"=>nothing,
        "hdFonte"=>nothing,
        "hdOidSerieMetadados"=>nothing,
        "hdNumeracao"=>nothing,
        "hdOidSeriesLocalizadas"=>nothing,
        "linkRetorno"=>"/sgspub/consultarvalores/telaCvsSelecionarSeries.paint",
        "linkCriarFiltros"=>"/sgspub/manterfiltros/telaMfsCriarFiltro.paint",
    )
end

function init_search_session(language="pt")
    session = HTTP.get(SEARCH_URL[language]; cookies = true)
end

function serialize_form(params)
    r = []
    for k in keys(params)
        push!(r,"$(k)=$(params[k]==nothing ? "" : params[k])")
    end
    join(r,"&")
end

"""
    function search_timeseries(query)

Localiza séries por `query` - código(Integer) ou texto(String)
"""
function search_timeseries(query)
    init_search_session()
    (searchMethod, tipo,codigo,texto) = get_tipo_busca(query)
    url = "https://www3.bcb.gov.br/sgspub/localizarseries/localizarSeries.do?"
    params = _get_params(searchMethod,tipo,codigo,texto)
    sf = serialize_form(params)
    r = HTTP.post(url*"&$sf",[],cookies=true)
    decode(r.body,enc"Latin1") |> parse_resp_timeseries
end

"""
    function parse_resp_timeseries(resp)

Localiza campos a partir da resposta da query de times series
"""
function parse_resp_timeseries(resp)
    html = replace(replace(replace(resp,"\n"=>""),"\t"=>""),"\r"=>"") |> x-> replace(x, r".+(<table[^>]+id=.tabelaSeries.+</table[^>]*>).*"=>s"\1") |> parsehtml
    reduce(html.root.children[2].children[1].children[2].children,init=DataFrame()) do a,b
        push!(a,(;
            id=parse(Int64,b.children[2].children[1].text),
            descrição=b.children[3].children[1].text,
            periodicidade=b.children[5].children[1].children[1].text
        ))
    end

end

"""
    function timeseries_values(timeseries::Vector;dtInicial::Dates.Date=Date(2018),dtFinal=missing, progress=nothing)

Retorna os valores de um conjunto de `timeseries` a partir de uma data inicial `dtInicial` até a data final `dtFinal`
"""
function timeseries_values(timeseries::Vector;dtInicial::Dates.Date=Date(2018),dtFinal=missing)
    ProgressLogging.progress() do id
        progresso = 0.0
        reduce(timeseries,init=DataFrame()) do a,b
            progresso += 1/length(timeseries)
            r = append!(a,timeserie_value(b;dtInicial=dtInicial,dtFinal=dtFinal))
            @info "Downloading series: " progress=progresso _id=id
            r
        end
    end
end

export timeserie_value, timeseries_values, search_timeseries

end


text:
using Statistics
using PyCall
pygui(:tk)
import PyPlot; const plt = PyPlot
import InteratomicPotentials as Potentials
using LinearAlgebra
n=50
########################################################################################
########################################################################################

n = 50
pe_true = zeros(5000, n)
pe_snap = zeros(5000, n)
pe_lj = zeros(5000, n)
for j = 1:n
    t, pe = Potentials.get_pe("examples/Argon13/DATA_training/true_ensemble/seed_$j/tmp.pe", 0.005)
    pe_true[:, j] = pe

    t, pe = Potentials.get_pe("examples/Argon13/fit_lj/samples/sample_1/seed_$j/tmp.pe", 0.005)
    pe_lj[:, j] = pe

   t, pe = Potentials.get_pe("examples/Argon13/fit_snap/samples/sample_1/seed_$j/tmp.pe", 0.005)
   pe_snap[:, j] = pe

end

fig = plt.figure(figsize=(15, 8))
ax = fig.add_subplot(1, 2, 1)

ax.hist(vec(pe_true), bins = 100, density = true, stacked = true, color = "red", alpha = 0.75, label = "True LJ")
ax.hist(vec(pe_lj), bins = 100, density = true, stacked = true, color = "blue", alpha = 0.5, label = "Fitted LJ")
ax.hist(vec(pe_snap), bins = 100, density = true, stacked = true, color = "black", alpha = 0.25, label = "Fitted SNAP")
ax.legend()
ax.set_ylabel("Probability")
ax.set_xlabel("Energy (LJ Units)")
ax.set_title("Distribution of Energy over seeds")
plt.grid("on")

# pe_snap = zeros(5000, 4*12)
# pe_lj = zeros(5000, 4*12)
# count = 1
# for j = 1:4
#     for i = 1:12
#         t, pe = Potentials.get_pe("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.pe", 0.005)
#         pe_lj[:, count] = pe

#         t, pe = Potentials.get_pe("examples/Argon13/fit_snap/samples/sample_$i/seed_$j/tmp.pe", 0.005)
#         pe_snap[:, count] = pe
#         count += 1
#     end
# end

# ax = fig.add_subplot(1, 3, 2)
# # ax.hist(vec(pe_true), bins = 100, density = true, stacked = true, color = "red", alpha = 0.25, label = "True LJ")
# ax.hist(vec(pe_lj), bins = 100, density = true, stacked = true, color = "blue", alpha = 0.5, label = "Fitted LJ")
# ax.hist(vec(pe_snap), bins = 100, density = true, stacked = true, color = "black", alpha = 0.25, label = "Fitted SNAP")
# ax.set_xlabel("Energy (LJ Units)")
# ax.set_title("Distribution of Energy over samples")
# plt.grid("on")

pe_snap = zeros(5000, 40*12)
pe_lj = zeros(5000, 40*12)
count = 1
for j = 1:50
    # println(j)
    for i = 1:12
        # println(i)
        try Potentials.get_rdf("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.rdf")
            # println("SNAP")
            t, pe = Potentials.get_pe("examples/Argon13/fit_snap/samples/sample_$i/seed_$j/tmp.pe", 0.005)
            pe_snap[:, count] = pe

            # println("LJ")
            t, pe = Potentials.get_pe("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.pe", 0.005)
            pe_lj[:, count] = pe
            count += 1
        catch
            continue
        end

        if count > 40*12
            break
        end
    end

end

ax = fig.add_subplot(1, 2, 2)
# ax.hist(vec(pe_true), bins = 100, density = true, stacked = true, color = "red", alpha = 0.25, label = "True LJ")
ax.hist(vec(pe_lj), bins = 100, density = true, stacked = true, color = "blue", alpha = 0.5, label = "Fitted LJ")
ax.hist(vec(pe_snap), bins = 100, density = true, stacked = true, color = "black", alpha = 0.25, label = "Fitted SNAP")
ax.set_xlabel("Energy (LJ Units)")
ax.set_title("Variability of Energy over seeds and samples")
plt.grid("on")

plt.gcf()

plt.savefig("variability_of_pe.png", dpi = 600, bbox_inches = "tight")

########################################################################################
########################################################################################
lj = Potentials.LennardJones(1.0, 1.0)

c = Potentials.load_lammps("examples/Argon13/DATA_start", atom_names = [:Ar], radii = [3.1], weights = [1.0], boundary_type = ["f", "f", "f"])
Nt = 500000
dT = 1000
dt = 0.01
T = dT:dT:Nt

n = 50
pe_true = zeros(500, n)
pe_snap = zeros(500, n)
pe_lj = zeros(500, n)
for j = 1:n
    println(j)
    r = Potentials.get_positions(c, "examples/Argon13/DATA_training/true_ensemble/seed_$j/", T)
    pe = Potentials.force(r, lj)
    pe_true[:, j] = norm.(pe)

    r = Potentials.get_positions(c, "examples/Argon13/fit_lj/samples/sample_1/seed_$j/", T)
    pe = Potentials.force(r, lj)
    pe_lj[:, j] = norm.(pe)

   r = Potentials.get_positions(c, "examples/Argon13/fit_snap/samples/sample_1/seed_$j/", T)
   pe = Potentials.force(r, lj)
   pe_snap[:, j] = norm.(pe)

end

fig = plt.figure(figsize=(15, 8))
ax = fig.add_subplot(1, 2, 1)

ax.hist(vec(pe_true), bins = 100, density = true, stacked = true, color = "red", alpha = 0.75, label = "True LJ")
ax.hist(vec(pe_lj), bins = 100, density = true, stacked = true, color = "blue", alpha = 0.5, label = "Fitted LJ")
ax.hist(vec(pe_snap), bins = 100, density = true, stacked = true, color = "black", alpha = 0.25, label = "Fitted SNAP")
ax.legend()
ax.set_ylabel("Probability")
ax.set_xlabel("Force (LJ Units)")
ax.set_title("Distribution of Force over seeds")
plt.grid("on")

# pe_snap = zeros(5000, 4*12)
# pe_lj = zeros(5000, 4*12)
# count = 1
# for j = 1:4
#     for i = 1:12
#         t, pe = Potentials.get_pe("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.pe", 0.005)
#         pe_lj[:, count] = pe

#         t, pe = Potentials.get_pe("examples/Argon13/fit_snap/samples/sample_$i/seed_$j/tmp.pe", 0.005)
#         pe_snap[:, count] = pe
#         count += 1
#     end
# end

# ax = fig.add_subplot(1, 3, 2)
# # ax.hist(vec(pe_true), bins = 100, density = true, stacked = true, color = "red", alpha = 0.25, label = "True LJ")
# ax.hist(vec(pe_lj), bins = 100, density = true, stacked = true, color = "blue", alpha = 0.5, label = "Fitted LJ")
# ax.hist(vec(pe_snap), bins = 100, density = true, stacked = true, color = "black", alpha = 0.25, label = "Fitted SNAP")
# ax.set_xlabel("Energy (LJ Units)")
# ax.set_title("Distribution of Energy over samples")
# plt.grid("on")

pe_snap = zeros(500, 20*12)
pe_lj = zeros(500, 20*12)
count = 1
for j = 1:30
    # println(j)
    for i = 1:12
        # println(i)
        try Potentials.get_rdf("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.rdf")
            # println("SNAP")
            r = Potentials.get_positions(c, "examples/Argon13/fit_lj/samples/sample_$i/seed_$j/", T)
            pe = Potentials.force(r, lj)
            pe_lj[:, count] = norm.(pe)

            # println("LJ")
            r = Potentials.get_positions(c, "examples/Argon13/fit_snap/samples/sample_$i/seed_$j/", T)
            pe = Potentials.force(r, lj)
            pe_snap[:, count] = norm.(pe)
            count += 1
        catch
            continue
        end

        if count > 20*12
            break
        end
    end

end

ax = fig.add_subplot(1, 2, 2)
ax.hist(vec(pe_true), bins = 100, density = true, stacked = true, color = "red", alpha = 0.25, label = "True LJ")
ax.hist(vec(pe_lj), bins = 100, density = true, stacked = true, color = "blue", alpha = 0.5, label = "Fitted LJ")
ax.hist(vec(pe_snap), bins = 100, density = true, stacked = true, color = "black", alpha = 0.25, label = "Fitted SNAP")
ax.set_xlabel("Force (LJ Units)")
ax.set_title("Distribution of Force over seeds and samples")
plt.grid("on")

plt.gcf()

plt.savefig("variability_of_lj_force.png", dpi = 600, bbox_inches = "tight")

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


text:
using Maker
using Base.Test

Maker.phony("phony")

Maker.variable("a", "phony") do   # phony limits this to one run
    global COUNT += 1
    pi
end

Maker.variable("b", "in1.csv") do
    global COUNT += 1
    readcsv("in1.csv", skipstart=1)
end

Maker.file("in1.csv")
Maker.file("in2.csv")

Maker.file("c.csv", ["a", "b"]) do
    global COUNT += 1
    writecsv("c.csv", a * b)
end

Maker.variable("c", "c.csv") do
    global COUNT += 1
    readcsv("c.csv")
end

Maker.variable("d", "in2.csv") do
    global COUNT += 1
    readcsv("in2.csv", skipstart=1)
end

Maker.variable("e", ["c", "d"]) do
    global COUNT += 1
    c .* d
end

Maker.file("e.csv", "e") do
    global COUNT += 1
    writecsv("e.csv", e)
end


text:
include("../utils.jl")

@testset "FunctionCurve manifold" begin
    M = FunctionCurveSpace(Euclidean(2))
    @testset "FunctionCurve Basics" begin
        @test manifold_dimension(M) === Inf
    end

    x = t -> SA[t, t - 1]
    y = t -> SA[t^2, sin(t)]
    z = t -> SA[2 - 2, cos(t)]

    test_manifold(
        M,
        [x, y, z],
        test_forward_diff=false,
        test_reverse_diff=false,
        test_vector_spaces=true,
        test_project_tangent=false,
        test_musical_isomorphisms=false,
        test_default_vector_transport=false,
        test_tangent_vector_broadcasting=false,
        test_representation_size=false,
        is_mutating=false,
        expected_dimension_type=Float64,
    )
end


text:
function fib100()
#
#  test polynomial suggested by Goedecker
#
    n = 100;
    p = fib(n)
    z = [-.9891099736226092;
         -.9871753104718534-.06190245872746059*im;  -.9871753104718534+.06190245872746059*im;
         -.9813788056446001-.1235640063982463*im;  -.9813788056446001+.1235640063982463*im;
         -.9717428840956687-.1847446602565987*im; -.9717428840956687+.1847446602565987*im;
         -.9583048244621342-.2452062905225305*im; -.9583048244621342+.2452062905225305*im;
         -.9411166150274026-.3047135378280153*im; -.9411166150274026+.3047135378280153*im;
         -.9202447528605484-.3630347198287647*im; -.9202447528605484+.3630347198287647*im;
         -.8957699869173070-.4199427234459649*im; -.8957699869173070+.4199427234459649*im;
         -.8677870061092108-.4752158792455410*im; -.8677870061092108+.4752158792455410*im;
         -.8364040735640849-.5286388145285358*im; -.8364040735640849+.5286388145285358*im;
         -.8017426085137990-.5800032817847029*im; -.8017426085137990+.5800032817847029*im;
         -.7639367174532022-.6291089592520593*im; -.7639367174532022+.6291089592520593*im;
         -.7231326764169818-.6757642204274601*im; -.7231326764169818+.6757642204274601*im;
         -.6794883664183336-.7197868694867672*im; -.6794883664183336+.7197868694867672*im;
         -.6331726642844369-.7610048396973263*im; -.6331726642844369+.7610048396973263*im;
         -.5843647913085710-.7992568520396451*im; -.5843647913085710+.7992568520396451*im;
         -.5332536223172150-.8343930313987525*im; -.5332536223172150+.8343930313987525*im;
         -.4800369579227939-.8662754778380540*im; -.4800369579227939+.8662754778380540*im;
         -.4249207628992477-.8947787906289414*im; -.4249207628992477+.8947787906289414*im;
         -.3681183737790589-.9197905428773492*im; -.3681183737790589+.9197905428773492*im;
         -.3098496789278928-.9412117047633791*im; -.3098496789278928+.9412117047633791*im;
         -.2503402745082667-.9589570135917125*im; -.2503402745082667+.9589570135917125*im;
         -.1898205998989927-.9729552890388070*im; -.1898205998989927+.9729552890388070*im;
         -.1285250562957008-.9831496921783526*im; -.1285250562957008+.9831496921783526*im;
         -.06669111238377778-.9894979270705116*im; -.06669111238377778+.9894979270705116*im;
         -.004558401154059704-.9919723839157337*im; -.004558401154059704+.9919723839157337*im;
         .05763218786929409-.9905602230050737*im; .05763218786929409+.9905602230050737*im;
         .1196394164912155-.9852633989536015*im; .1196394164912155+.9852633989536015*im;
         .1812226010328281-.9760986249939283*im; .1812226010328281+.9760986249939283*im;
         .2421425057566853-.9630972774521361*im; .2421425057566853+.9630972774521361*im;
         .3021622210555556-.9463052409577104*im; .3021622210555556+.9463052409577104*im;
         .3610480204611693-.9257826954965046*im; .3610480204611693+.9257826954965046*im;
         .4185701899698170-.9016038471668328*im; .4185701899698170+.9016038471668328*im;
         .4745038225339292-.8738566055392182*im; .4745038225339292+.8738566055392182*im;
         .5286295698615102-.8426422119868863*im; .5286295698615102+.8426422119868863*im;
         .5807343429795305-.8080748254361926*im; .5807343429795305+.8080748254361926*im;
         .6306119525307706-.7702810749337576*im; .6306119525307706+.7702810749337576*im;
         .6780636798327496-.7293995925445015*im; .6780636798327496+.7293995925445015*im;
         .7228987709578736-.6855805456977909*im; .7228987709578736+.6855805456977909*im;
         .7649348495515867-.6389851953915993*im; .7649348495515867+.6389851953915993*im;
         .8039982514449404-.5897855154613771*im; .8039982514449404+.5897855154613771*im;
         .8399242976881979-.5381639173058062*im; .8399242976881979+.5381639173058062*im;
         .8725575452576575-.4843131310830981*im; .8725575452576575+.4843131310830981*im;
         .9017520886609776-.4284362924194566*im; .9017520886609776+.4284362924194566*im;
         .9273720302113955-.3707472629394290*im; .9273720302113955+.3707472629394290*im;
         .9492922833419482-.3114711598054638*im; .9492922833419482+.3114711598054638*im;
         .9673998997347004-.2508449722447981*im; .9673998997347004+.2508449722447981*im;
         .9815960786437149-.1891180053506043*im; .9815960786437149+.1891180053506043*im;
         .9917988816037989-.1265517518780266*im; .9917988816037989+.1265517518780266*im;
         .9979464229745987-.06341873580445710*im; .9979464229745987+.06341873580445710*im;
         2.000000000000000];
    z = [z ones(n)];
    p, PolyZeros(z)
end


text:
module Logs

using Statistics, DelimitedFiles, BeastUtils.MatrixUtils, JLD, CSV, DataFrames

export get_log,
       get_log_match,
       import_log,
    list_cols,
    make_meanmatrix,
    make_meanvector,
    make_diagonal_vector,
    make_vector,
    make_correlationMatrix,
    confidence_intervals,
    HDP_intervals,
    log2corr,
    find_cols,
    get_cols,
    combine_logs,
    condense_logs,
    condense_trees,
    make_log,
    compress_log,
    compress_log!,
    read_time

const DELIM_CHAR = '\t'

function compress_log(log_path::String, compressed_path::String)
    cols, data = get_log(log_path, burnin = 0.0)
    save(compressed_path, "cols", cols, "data", data)
end

function compress_log!(log_path::String, compressed_path::String)
    compress_log(log_path, compressed_path)
    rm(log_path)
end

function uncompress_log(compressed_path::String, log_path::String)
    d = load(compressed_path)
    cols = d["cols"]
    data = d["data"]
    make_log(log_path, data, cols, includes_states = true)
end

function uncompress_log!(compressed_path::String, log_path::String)
    uncompress_log(compressed_path, log_path)
    rm(compressed_path)
end

function uncompress_log!(path::String)
    split_path = split(path, '.')
    ext = last(split_path)
    if ext != "jld"
        error("Compressed file must have `.jld` file extension.")
    end

    split_path[end] = "log"
    uncompress_log!(path, join(split_path, '.'))
end

function make_log(path::String, data::AbstractMatrix{Float64},
            col_labels::Vector{String}; includes_states::Bool = false)
    n, p = size(data)
    if length(col_labels) != p
        error("The number of column labels must be the same as the number of columns.")
    end
    open(path, "w") do f
        if !includes_states
            write(f, "state$DELIM_CHAR")
        end
        write(f, join(col_labels, DELIM_CHAR))
        write(f, '\n')
        for i = 1:n
            if !includes_states
                write(f, "$(i - 1)$DELIM_CHAR")
            end
            write(f, join(data[i, :], DELIM_CHAR))
            write(f, '\n')
        end
    end
end

# function below should eventually deprecate get_log()
function import_log(path::String; burnin::Float64 = 0.0)
    header_row = find_header_row(path)
    all_lines = countlines(path)
    n = all_lines - header_row
    skipto = header_row + Int(round(n * burnin)) + 1

    df = CSV.read(path,
                  DataFrame,
                  header=header_row,
                  skipto = skipto,
                  delim='\t')

    return df
end

function get_log(inpath::String; burnin::Float64 = 0.1)
    if endswith(inpath, ".jld")
        d = load(inpath)
        cols = d["cols"]
        data = d["data"]
        n = size(data, 1)
        n_burn = Int(round(burnin * n))

        return cols, data[(n_burn + 1):end, :]
    else

        col_labels, labels_ind = get_cols_and_ind(inpath)
        p = length(col_labels)
        f = open(inpath)
        m = countlines(f) - labels_ind
        close(f)
        n_burn = Int(round(burnin * m))
        n = m - n_burn
        data = readdlm(inpath, skipstart = labels_ind + n_burn, dims = (n, p))
        return string.(col_labels), data
    end
end

function read_last_line(inpath::String)
    col_labels, labels_ind = get_cols_and_ind(inpath)
    p = length(col_labels)
    f = open(inpath)
    m = countlines(f)

    data = readdlm(inpath, skipstart = m-1, dims = (1, p))
    return string.(col_labels), vec(data)
end

function get_log(inpath::String, cols::Array{Int}; burnin::Float64 = 0.1)
    col_labels, data = get_log(inpath, burnin = burnin)
    return col_labels[cols], data[:, cols]
end

function get_log(inpath::String, header::Union{Regex, String};
                    burnin::Float64 = 0.1)
    Base.depwarn("get_log(path::String, header::$(typeof(header))) is deprecated. " *
                "Use get_log_match(path, header) instead.", :get_log)
    all_cols = get_cols(inpath)
    f, l = find_cols(all_cols, header)
    cols, data = get_log(inpath, collect(f:l), burnin = burnin)
end

function get_log_match(inpath::String, header::Union{Regex, String};
                        burnin::Float64 = 0.1)
    all_cols = get_cols(inpath)
    inds = findall(x -> matches(x, header), all_cols)
    cols, data = get_log(inpath, inds, burnin = burnin)
end

function get_cols(inpath::String)
    if endswith(inpath, ".jld")
        cols = load(inpath, "cols")
        return cols
    else
        cols, ind = get_cols_and_ind(inpath)
        return cols
    end
end

function get_cols_and_ind(inpath::String)
    f = open(inpath)
    col_labels = ""
    ind = 0
    for line in eachline(f)
        ind += 1
        if startswith(line, "state")
            col_labels = line
            break
        end
    end
    close(f)
    cols = split(col_labels)
    return cols, ind
end

function find_header_row(inpath::String)
    counter = 0
    open(inpath, "r") do f
        for line in eachline(inpath)
            counter += 1
            if startswith(line, "state")
                break
            end
        end
    end
    return counter
end

function list_cols(cols::Vector{String})
    p = length(cols)
    for i = 1:p
        println("$i: $(cols[i])")
    end
end

function find_cols(cols::Vector{S}, col::Union{String, Regex}) where S <: AbstractString
    first = 0
    last = 0
    p = length(cols)
    for i = 1:p
        if matches(cols[i], col)
            first = i
            break
        end
    end
    if first == 0
        error("$col not found in column names.")
    end
    for i = p:-1:first
        if matches(cols[i], col)
            last = i
            break
        end
    end
    return first, last
end

function matches(col::S, header::String) where S <: AbstractString
    if startswith(col, header)
        return true
    end
    return false
end

function matches(col::S, header::Regex) where S <: AbstractString
    if occursin(header, col)
        return true
    end
    return false
end

function make_meanmatrix(data::Matrix{T}, cols::Vector{String},
        col::String) where T <: Real

    first, last = find_cols(cols, col)
    return make_meanmatrix(data, first, last)
end

function make_meanmatrix(data::Matrix{T}, first::Int, last::Int) where T <: Real
    t_data = data[:, (first:last)]
    means = mean(t_data, dims = 1)
    q = length(means)
    q = Int(sqrt(q))
    mat = reshape(means, q, q)
    return mat
end

function make_meanmatrix(data::Matrix{T}) where T <: Real
    means = mean(data, dims = 1)
    q = length(means)
    q = Int(sqrt(q))
    mat = reshape(means, q, q)
    return mat
end

function make_meanvector(data::Matrix{T}, first::Int, last::Int) where T <: Real

    means = [mean(data[:, i]) for i in first:last]

    return means
end

function make_diagonal_vector(data::Matrix{T}, first::Int, last::Int) where T <: Real
    q = last - first + 1
    p = Int(sqrt(q))
    n = size(data, 1)
    new_mat = zeros(n, p)
    for i = 1:p
        new_mat[:, i] .= data[:, first + (i - 1) * (p + 1)]
    end
    return new_mat
end

function make_vector(data::Matrix{T}, first::Int, last::Int) where T <: Real
    t_data = data[:, (first:last)]
    return t_data
end

function make_correlationMatrix(X::Matrix{T}) where T <: Real
    n, p = size(X)
    @assert n == p
    C = zeros(n, n)
    for i = 1:n
        for j = 1:n
            C[i, j] = X[i, j] / sqrt(X[i, i] * X[j, j])
        end
    end
    return C
end

function confidence_intervals(d::Matrix{Float64}, p::Int)
    x = d[:, p]
    n = length(x)
    sort!(x)
    bottom = Int(floor(0.025 * n))
    top = Int(ceil(0.975 * n))
    return x[bottom], x[top]
end

function HPD_intervals(d::Matrix{Float64}, p::Int; conf::Float64 = 0.95)
    x = d[:, p]
    return HDP_intervals(x, conf = conf)
end

function HPD_intervals(x::Vector{Float64}; conf::Float64 = 0.95)
    n = length(x)
    sort!(x)
    s = Int(ceil(conf * n))
    smallest = Inf
    smallest_inds = [0, 0]
    for i = 1:(n - s + 1)
        diff = x[i + s - 1] - x[i]
        if diff < smallest
            smallest = diff
            smallest_inds[1] = i
            smallest_inds[2] = i + s - 1
        end
    end
    return x[smallest_inds]
end

function HPD_intervals(d::Matrix{Float64}, ps::Vector{Int}; conf::Float64 = 0.95)
    hdps = zeros(length(ps), 2)
    for i = 1:length(ps)
        hdpi = HDP_intervals(d, ps[i], conf = conf)
        hdps[i, 1] = hdpi[1]
        hdps[i, 2] = hdpi[2]
    end
    return hdps
end

function log2corr(X::Matrix{Float64}; invert::Bool = false)
    n, q = size(X)
    p = Int(sqrt(q))
    M = zeros(p, p)
    new_X = zeros(n, q)
    for i = 1:n
        for j = 1:q
            M[j] = X[i, j]
        end
        if invert
            M = inv(M)
        end
        corrs = cov2corr(M)
        for j = 1:q
            new_X[i, j] = corrs[j]
        end
    end
    return new_X
end

function get_header(path::String)
    f = open(path)
    lines = Vector{String}(undef, 0)
    for line in eachline(f)
        if startswith(line, "state")
            break
        end
        push!(lines, line)
    end
    close(f)
    return lines
end

function combine_logs(paths::Vector{String}, outpath::String;
        force::Bool = false)
    n = length(paths)
    cols = get_cols(paths[1])
    for i = 2:n
        test_cols = get_cols(paths[i])
        @assert test_cols == cols
    end
    lines = get_header(paths[1])
    push!(lines, join(cols, "\t"))
    state = 0
    all_data = Matrix{Float64}(undef, 0, length(cols))
    for i = 1:n
        cols, data = get_log(paths[i], burnin = 0.0)
        m = size(data, 1)
        for i = 1:m
            data[i, 1] += state
        end
        if i != 1
            data = data[2:end, :]
        end
        all_data = [all_data; data]
        state = data[end, 1]
    end
    m = size(all_data, 1)
    string_data = Vector{String}(undef, m)
    for i = 1:m
        string_data[i] = join(all_data[i, :], "\t")
    end
    lines = [lines; string_data]
    lines = join(lines, "\n")
    write(outpath, lines)
end

function writeline(io::IOStream, s::AbstractString)
    write(io, s)
    write(io, "\n")
end

function condense_logfiles(in_path::String, out_path::String;
                       every::Int = 10,
                       is_valid = is_valid_logfile
    )

    counter = 0

    open(in_path, "r") do input
        open(out_path, "w") do output

            for line in eachline(input)
                if is_valid(line)
                    if rem(counter, every) == 0
                        writeline(output, line)
                    end
                    counter += 1
                else
                    writeline(output, line)
                end
            end
        end
    end
end

function is_valid_logfile(line::AbstractString)::Bool
    first_tab = findfirst('\t', line)
    if isnothing(first_tab)
        return false
    end

    sub_line = @view line[1:(first_tab - 1)]

    try parse(Int, sub_line)
        return true
    catch
        return false
    end
end

function is_valid_treefile(line::AbstractString)::Bool
    return startswith(line, "tree STATE_")
end

function condense_logs(in_path::String, out_path::String; every::Int = 10)
    return condense_logfiles(in_path, out_path, every = every, is_valid = is_valid_logfile)
end

function condense_trees(in_path::String, out_path::String; every::Int = 10)
    return condense_logfiles(in_path, out_path, every = every, is_valid = is_valid_treefile)
end

const TIME_DICT = Dict("seconds" => 1.0,
                       "minutes" => 60.0,
                       "hours" => 3600.0,
                       "days" => 86400.0)

function read_time(path::String; unit::String = "seconds")
    t, u = split(read(path, String))
    t = parse(Float64, t)
    return t * TIME_DICT[u] / TIME_DICT[unit]
end

end


text:
@testset "$(rpad("Tableau coefficients",80))" begin

    # test Gauss-Legendre Runge-Kutta coefficients

    glrk2_tab1 = getCoefficientsGLRK(2)
    glrk2_tab2 = getCoefficientsGLRK(2, high_precision=false)
    @test glrk2_tab1.a ≈ glrk2_tab2.a atol=4eps()
    @test glrk2_tab1.b == glrk2_tab2.b
    @test glrk2_tab1.c ≈ glrk2_tab2.c atol=4eps()

    glrk3_tab1 = getCoefficientsGLRK(3)
    glrk3_tab2 = getCoefficientsGLRK(3, high_precision=false)
    @test glrk3_tab1.a ≈ glrk3_tab2.a atol=4eps()
    @test glrk3_tab1.b == glrk3_tab2.b
    @test glrk3_tab1.c ≈ glrk3_tab2.c atol=4eps()

    @test getCoefficientsGLRK(1) == getCoefficientsGLRK1()
    @test getCoefficientsGLRK(2) == getCoefficientsGLRK2()
    @test getCoefficientsGLRK(3) == getCoefficientsGLRK3()
    @test getCoefficientsGLRK(4) == getCoefficientsGLRK4()
    @test getCoefficientsGLRK(5) == getCoefficientsGLRK5()

    # test Gauss-Lobatto Runge-Kutta coefficients

    @test getCoefficientsLobIII(2) == getCoefficientsLobIII2()
    @test getCoefficientsLobIII(3) == getCoefficientsLobIII3()
    @test getCoefficientsLobIII(4) == getCoefficientsLobIII4()

    @test getCoefficientsLobIIIA(2) == getCoefficientsLobIIIA2()
    @test getCoefficientsLobIIIA(3) == getCoefficientsLobIIIA3()
    @test getCoefficientsLobIIIA(4) == getCoefficientsLobIIIA4()

    @test getCoefficientsLobIIIB(2) == getCoefficientsLobIIIB2()
    @test getCoefficientsLobIIIB(3) == getCoefficientsLobIIIB3()
    @test getCoefficientsLobIIIB(4) == getCoefficientsLobIIIB4()

    @test getCoefficientsLobIIIC(2) == getCoefficientsLobIIIC2()
    @test getCoefficientsLobIIIC(3) == getCoefficientsLobIIIC3()
    @test getCoefficientsLobIIIC(4) == getCoefficientsLobIIIC4()

    @test getCoefficientsLobIIID(2) == getCoefficientsLobIIID2()
    @test getCoefficientsLobIIID(3) == getCoefficientsLobIIID3()
    @test getCoefficientsLobIIID(4) == getCoefficientsLobIIID4()

    @test getCoefficientsLobIIIE(2) == getCoefficientsLobIIIE2()
    @test getCoefficientsLobIIIE(3) == getCoefficientsLobIIIE3()
    @test getCoefficientsLobIIIE(4) == getCoefficientsLobIIIE4()

    # test PGLRK coefficients

    pglrk2 = getCoefficientsPGLRK(2)
    @test typeof(pglrk2) <: CoefficientsPGLRK

end


text:
using ITensorGaussianMPS
using ITensors
using LinearAlgebra
using Test

@testset "Basic" begin
  # Test Givens rotations
  v = randn(6)
  g, r = ITensorGaussianMPS.givens_rotations(v)
  @test g * v ≈ r * [n == 1 ? 1 : 0 for n in 1:length(v)]
end

@testset "Fermion" begin
  N = 10
  Nf = N ÷ 2

  # Hopping
  t = 1.0

  # Hopping Hamiltonian
  h = Hermitian(diagm(1 => fill(-t, N - 1), -1 => fill(-t, N - 1)))
  e, u = eigen(h)

  @test h * u ≈ u * Diagonal(e)

  E = sum(e[1:Nf])

  # Get the Slater determinant
  Φ = u[:, 1:Nf]
  @test h * Φ ≈ Φ * Diagonal(e[1:Nf])

  # Diagonalize the correlation matrix as a
  # Gaussian MPS (GMPS)
  n, gmps = slater_determinant_to_gmps(Φ; maxblocksize=4)

  ns = round.(Int, n)
  @test sum(ns) == Nf

  Λ = conj(Φ) * transpose(Φ)
  @test gmps * Λ * gmps' ≈ Diagonal(ns) rtol = 1e-2
  @test gmps' * Diagonal(ns) * gmps ≈ Λ rtol = 1e-2

  # Form the MPS
  s = siteinds("Fermion", N; conserve_qns=true)
  ψ = slater_determinant_to_mps(s, Φ; blocksize=4)

  ampo = AutoMPO()
  for i in 1:N, j in 1:N
    if h[i, j] ≠ 0
      ampo .+= h[i, j], "Cdag", i, "C", j
    end
  end
  H = MPO(ampo, s)

  @test inner(ψ, H, ψ) ≈ E rtol = 1e-5

  # Compare to DMRG
  sweeps = Sweeps(10)
  maxdim!(sweeps, 10, 20, 40, 60)
  cutoff!(sweeps, 1E-12)
  energy, ψ̃ = dmrg(H, productMPS(s, n -> n ≤ Nf ? "1" : "0"), sweeps; outputlevel=0)

  # Create an mps
  @test abs(inner(ψ, ψ̃)) ≈ 1 rtol = 1e-5
  @test inner(ψ̃, H, ψ̃) ≈ inner(ψ, H, ψ) rtol = 1e-5
  @test E ≈ energy
end

@testset "Fermion (complex)" begin
  N = 10
  Nf = N ÷ 2

  # Hopping
  θ = π / 8
  t = exp(im * θ)

  # Hopping Hamiltonian
  h = Hermitian(diagm(1 => fill(-t, N - 1), -1 => fill(-conj(t), N - 1)))
  e, u = eigen(h)

  @test h * u ≈ u * Diagonal(e)

  E = sum(e[1:Nf])

  # Get the Slater determinant
  Φ = u[:, 1:Nf]
  @test h * Φ ≈ Φ * Diagonal(e[1:Nf])

  # Diagonalize the correlation matrix as a
  # Gaussian MPS (GMPS)
  n, gmps = slater_determinant_to_gmps(Φ; maxblocksize=4)

  ns = round.(Int, n)
  @test sum(ns) == Nf

  Λ = conj(Φ) * transpose(Φ)
  @test gmps * Λ * gmps' ≈ Diagonal(ns) rtol = 1e-2
  @test gmps' * Diagonal(ns) * gmps ≈ Λ rtol = 1e-2

  # Form the MPS
  s = siteinds("Fermion", N; conserve_qns=true)
  ψ = slater_determinant_to_mps(s, Φ; blocksize=4)

  ampo = AutoMPO()
  for i in 1:N, j in 1:N
    if h[i, j] ≠ 0
      ampo .+= h[i, j], "Cdag", i, "C", j
    end
  end
  H = MPO(ampo, s)

  @test inner(ψ, H, ψ) ≈ E rtol = 1e-5
  @test inner(ψ, H, ψ) / norm(ψ) ≈ E rtol = 1e-5

  # Compare to DMRG
  sweeps = Sweeps(10)
  maxdim!(sweeps, 10, 20, 40, 60)
  cutoff!(sweeps, 1E-12)
  energy, ψ̃ = dmrg(H, productMPS(s, n -> n ≤ Nf ? "1" : "0"), sweeps; outputlevel=0)

  # Create an mps
  @test abs(inner(ψ, ψ̃)) ≈ 1 rtol = 1e-5
  @test inner(ψ̃, H, ψ̃) ≈ inner(ψ, H, ψ) rtol = 1e-5
  @test E ≈ energy
end


text:
# Metatheory.options[:printiter] = true

or_alg = @theory begin
    ((p ∨ q) ∨ r)       ==  (p ∨ (q ∨ r))
    (p ∨ q)             ==  (q ∨ p)
    (p ∨ p)             =>  p
    (p ∨ true)          =>  true
    (p ∨ false)         =>  p
end

and_alg = @theory begin
    ((p ∧ q) ∧ r)       ==  (p ∧ (q ∧ r))
    (p ∧ q)             ==  (q ∧ p)
    (p ∧ p)             =>  p
    (p ∧ true)          =>  p
    (p ∧ false)         =>  false
end

comb = @theory begin
    # DeMorgan
    ¬(p ∨ q)            ==  (¬p ∧ ¬q)
    ¬(p ∧ q)            ==  (¬p ∨ ¬q)
    # distrib
    (p ∧ (q ∨ r))       ==  ((p ∧ q) ∨ (p ∧ r))
    (p ∨ (q ∧ r))       ==  ((p ∨ q) ∧ (p ∨ r))
    # absorb
    (p ∧ (p ∨ q))       =>  p
    (p ∨ (p ∧ q))       =>  p
    # complement
    (p ∧ (¬p ∨ q))      =>  p ∧ q
    (p ∨ (¬p ∧ q))      =>  p ∨ q
end

negt = @theory begin
    (p ∧ ¬p)            =>  false
    (p ∨ ¬(p))          =>  true
    ¬(¬p)               ==  p
end

impl = @theory begin
    (p == ¬p)           =>  false
    (p == p)            =>  true
    (p == q)            =>  (¬p ∨ q) ∧ (¬q ∨ p)
    (p => q)            =>  (¬p ∨ q)
end

fold = @theory begin
    (p::Bool == q::Bool)    |>     (p == q)
    (p::Bool ∨ q::Bool)     |>     (p || q)
    (p::Bool => q::Bool)    |>     ((p || q) == q)
    (p::Bool ∧ q::Bool)     |>     (p && q)
    ¬(p::Bool)              |>     (!p)
end

# t = or_alg ∪ and_alg ∪ neg_alg ∪ demorgan ∪ and_or_distrib ∪
#     absorption ∪ calc

t = or_alg ∪ and_alg ∪ comb ∪ negt ∪ impl ∪ fold

@test @areequal t true ((¬p == p) == false)
@test @areequal t true ((¬p == ¬p) == true)
@test @areequal t true ((¬p ∨ ¬p) == ¬p) (¬p ∨ p) ¬(¬p ∧ p)
@test @areequal t true ((p => (p ∨ p)) == true)
@test @areequal t true ((p => (p ∨ p)) == ((¬(p) ∧ q) => q)) == true

# Frege's theorem
@test @areequal t true (p => (q => r)) => ((p => q) => (p => r))

# Demorgan's
@test @areequal t true (¬(p ∨ q) == (¬p ∧ ¬q))

# Consensus theorem
@test @areequal t ((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z))   ((x ∧ y) ∨ (¬x ∧ z))

# @timev areequal(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)), :((x ∧ y) ∨ (¬x ∧ z)))
# @timev areequal(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)), :((x ∧ y) ∨ (¬babo ∧ z)))
# @timev areequalmagic(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)),   :((x ∧ y) ∨ (¬x ∧ z)))
# @timev areequalmagic(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)),   :((babo ∧ y) ∨ (¬x ∧ z)))

# TODO proof strategies?
# Constructive Dilemma

# @test @areequal (t ∪ [@rule :p => true]) true (((p => q) ∧ (r => s)) ∧ (p ∨ r)) => (q ∨ s)

# @test areequal(t, true, :(¬(((¬p ∨ q) ∧ (¬r ∨ s)) ∧ (p ∨ r)) ∨ (q ∨ s)))

function prove(t, ex, steps)
    params = SaturationParams(timeout=8, sizeout=5300, schedulerparams=(8,2))
    hist = UInt64[]
    push!(hist, hash(ex))
    for i ∈ 1:steps
        g = EGraph(ex)
        saturate!(g, t, params)
        ex = extract!(g, astsize)
        println(ex)
        if !TermInterface.istree(ex)
            return ex
        end
        if hash(ex) ∈ hist
            println("loop detected")
            return ex
        end
        push!(hist, hash(ex))
    end
    return ex
end

ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
@test prove(t, ex, 3)

# using Metatheory.EGraphs.Schedulers
# ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
# g = EGraph(ex)
# params = SaturationParams(timeout=10, sizeout=2^15, scheduler=ScoredScheduler)
# @timev saturate!(g, t, params)
#
# extran = addanalysis!(g, ExtractionAnalysis, astsize)
# ex = extract!(g, extran)
# println(ex)
# #
# #
# ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
# g = EGraph(ex)
# @timev areequal(g, t, ex, true; timeout=10, sizeout=2^15, scheduler=ScoredScheduler, schedulerparams=(32, 1))
#

# @profiler saturate!(g, t; timeout=8, sizeout=2^15)
# exit(0)

# ex = rewrite(:(((p => p) ∧ (r => z) ∧ (p ∨ r)) => (q ∨ s)), impl)
# @test false == prove(t, ex, 4)

# g = EGraph(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)))
# @time saturate!(g, t; timeout=30, sizeout=Inf)
#
# in_same_set(g.U, g.root, addexpr!(g, true).id) |> println
#
# struct LogicAnalysis <: AbstractAnalysis
#     egraph::EGraph
#     logic_egraph::EGraph
# end


text:
# Includes alternate versions which print more information
using Combinatorics

# test_degen_complexity actually tests that the different versions of complexity agree.
function rand_X( numinputs::Int64, numints::Int64 )
  MyInt_bits = MyIntBits( MyInt )
  X = rand(collect(MyInt(0):MyInt(2^numinputs-1)),numints)
end

function rand_tc( numinputs::Int64, numints::Int64 )
  test_complexity( rand_X( numinputs, numints ), numinputs )
end

function test_complexity( X::Vector{MyInt},  numinputs::Int64 )
  gbX = get_bits(X,numinputs)
  println(to_binary(X,2^numinputs))
  #println("X: ",X,"  gbX: ",gbX,"  cmplx: ",complexity0(X,numinputs))
  #println("===================")
  println("X: ",X,"  gbX: ",gbX,"  cmplx: ",complexity(X,numinputs))
end

# Version of complexity5 which prints more information
function complexity( X::Vector{MyInt}, numinputs::Int64; base=Float64=2.0 )
  n = length(X)
  gbX = get_bits( X, numinputs )
  ent_X = entropy(gbX,base=base)
  println("X: ",X,"  gbX: ",gbX," ent_X: ",ent_X)
  for k = 1:(length(X))
    #println("k: ",k,"  ")
    #for s in combinations(X,k)
      #println("  s: ",s,"  gbs: ",get_bits(s,numinputs),"  ents: ",entropy(get_bits(s,numinputs),base=base))
    #end
  end
  ents = [map(x->entropy(x,base=base),map(x->get_bits(x,numinputs),[s for s in combinations(X,k)])) for k = 1:(length(X))]
  #println("ents: ",ents)
  ents_avg = map(x->sum(x)/length(x), ents )
  println("ents_avg: ",ents_avg)
  summand = [ents_avg[k] - k/n*ent_X for k = 1:n]
  println("summand: ",summand)
  sum(summand)
end

# Version of complexity6 which prints more information
function complexity1( X::Vector{MyInt}, numinputs::Int64; base::Float64=2.0 )
  n = length(X)
  Xinds = collect(1:length(X))
  ssum = 0.0
  for k = 1:Int(floor(n/2))
    subset_pairs = [(s,setdiff(Xinds,s)) for s in combinations(Xinds,k)]
    #println("k: ",k,"  subset_pairs: ",subset_pairs)
    X_pairs = map( i->( X[subset_pairs[i][1]], X[subset_pairs[i][2]] ), collect(1:length(subset_pairs)))
    #println("X_pairs: ",X_pairs)
    gbX_pairs = [ (get_bits(Xp[1],numinputs), get_bits(Xp[2],numinputs)) for Xp in X_pairs ]
    #println("gbX_pairs: ",gbX_pairs)
    mutints = [ mutual_information(get_bits(Xp[1],numinputs), get_bits(Xp[2],numinputs)) for Xp in X_pairs ]
    #println("mutints: ",mutints)
    summand = sum(mutints)/length(mutints)
    #println("k: ",k,"  sum mutual informaiton: ",summand)
    println(summand)
    if k == Int(ceil(n/2))
      summand /= 2
    end
    ssum += summand
  end
  ssum
end

# Version of complexity7 which prints more information
function complexity0( X::Vector{MyInt}, numinputs::Int64; base::Float64=2.0 )
  n = length(X)
  IX = integration( X, numinputs, base=base )
  println("X: ",X,"  IX: ",IX)
  Xinds = collect(1:length(X))
  ssum = 0.0
  for k = 1:(n-1)
    #println("k: ",k)
    for s in combinations(Xinds,k)
      compl = X[setdiff(Xinds,s)]
      IXc = integration0(compl,numinputs,base=base)
      #println("s: ",s,"  compl: ",compl,"  IXc ",IXc)
    end
    subsets = [X[setdiff(Xinds,s)] for s in combinations(Xinds,k)]
    IXk = map( x->integration(x,numinputs,base=base), subsets )
    #println("IXk: ",IXk)
    avgIXk = sum(IXk)/length(IXk)
    println("k: ",k,"  avgIXk: ",avgIXk)
    #println("k: ",k,"  avgIXk: ",avgIXk,"  summand: ",k/n*IX - avgIXk)
    ssum += k/n*IX - avgIXk
  end
  ssum
end

# Version of integration() with get_bits() that prints more information
function integration0( X::Vector{MyInt}, numinputs::Int64; base::Float64=2.0 )
  gb_list = [ get_bits([x],numinputs) for x in X]
  ent_list =  [ entropy(get_bits([x],numinputs)) for x in X ]
  #println("X: ",X,"  entX: ",entropy(get_bits(X,numinputs)))
  #println("gb_list: ",gb_list,"  ent_list: ",ent_list,"  sum: ",sum(ent_list))
  sum(ent_list) - entropy(get_bits(X,numinputs))
end

# Version of integration() without get_bits() that prints more information
function integration0( X::Vector{Int64}; base::Float64=2.0 )
  ent_list =  [ entropy([x]) for x in X ]
  println("X: ",X,"  entX: ",entropy(X))
  println("ent_list: ",ent_list,"  sum: ",sum(ent_list))
  sum(ent_list) - entropy(X)
end

# Integration of XI without get_bits
function integration0( XI::Vector{Vector{Int64}} ; base::Float64=2.0 )
  X = vcat( XI... )   # Combines all of the lists in XI into a long list
  println("X: ",X,"  entX: ",entropy(X))
  ent_list = [entropy(x) for x in XI]
  println("ent_list: ",ent_list,"  sum: ",sum(ent_list))
  sum( entropy(x) for x in XI ) - entropy( X)
end

# get_bits for a single MyInt
function get_bits1( x::Main.CGP.MyInt, numinputs::Int64 )
  numbits = 2^numinputs
  result = zeros(MyInt,2^numinputs)
  shift = numbits-1
  mask = MyInt(1) << shift
  for j = collect(numbits:-1:1)
    result[j] |= (mask & x) >> shift
    shift -= 1
    mask >>= 1
  end
  result
end

# Agrees with the results of get_bits
function get_bits1( v::Vector{MyInt}, numinputs::Int64 )
  numbits = 2^numinputs
  result = zeros(MyInt,numbits)
  mask_shift = length(v)-1
  for i = 1:length(v)
    shift = numbits-1
    mask = MyInt(1) << shift
    j = numbits-1
    while j > 1
    #for j = collect(numbits:-1:1)
      rr = (mask & v[i]) >> shift
      result[j] |= rr << mask_shift
      shift -= 1
      mask >>= 1
      j -= 1
    end
    mask_shift -= 1
  end
  result
end

# [([z,x,y], get_bits([z,x,y],2),[entropy(get_bits([p],2)) for p in [z,y,x]],(fract_part(integration([z,x,y],2)+entropy(get_bits([z,x,y],2))))) for y = 0x0:0x7  for x = 0x0:0x7  for z = 0x0:0x7]


text:
# This file includes an example to complare FIF and HiddenFIF

using FractalTools
using Plots

# Generate interpolation data
xi, dx, xf = 0, 0.01, 1
x = collect(xi : dx : xf)
y = sin.(2π * x)
z = cos.(2π * x)
z2 = tan.(2π * x)
u = 0.1
d = u * ones(length(x) - 1)
h = 0.1 * ones(length(x) - 1)
l = 0.1 * ones(length(x) - 1)
m = 0.1 * ones(length(x) - 1)

# Interplate data
interp = interpolate(x, y, d)
hiddeninterp = hiddenfif(x, y, z, d, h, l, m)
hiddeninterp2 = hiddenfif(x, y, z2, d, h, l, m)

# Calculate interpolation data
xd = collect(xi : 0.01dx : xf)
yd = interp.(xd)
hidden_yd = hiddeninterp.(xd)
hidden_yd2 = hiddeninterp2.(xd)

# Plots
plt = plot(layout=(2,1))
plot!(xd, getindex.(yd, 1), label="interpolate", subplot=1, lw=5)
scatter!(x, y, marker=(:circle, 5), subplot=1)
plot!(xd, getindex.(hidden_yd, 1), label="hidden", subplot=1)
scatter!(x, y, marker=(:circle, 5), subplot=1)
plot!(xd, getindex.(hidden_yd2, 1), label="hidden2", subplot=1)
scatter!(x, y, marker=(:circle, 5), subplot=1)
plot!(xd, getindex.(hidden_yd, 2), label="hidden_extra", subplot=2)
scatter!(x, z, marker=(:circle, 5), subplot=2)
plot!(xd, getindex.(hidden_yd2, 2), label="hidden_extra", subplot=2)
scatter!(x, z2, marker=(:circle, 5), subplot=2)


text:
using CSV, DataFrames, Gadfly, StatsBase
using Cairo, Fontconfig

INDIR = "/home/valentin/Development/fuzzing_exps/1112/tiff_out2_long"
OUTDIR = "figures"

base_dir = split(INDIR, "/")[end]
println("Analyzing: $(base_dir)")

################################################################################
################################# 1. Make Map ##################################

function calcNorm(x1, x2, y1, y2)
    xd = x2 - x1
    yd = y2 - y1
    return sqrt(xd*xd + yd*yd)
end
function prepareData(df::DataFrame, xSym::Symbol, ySym::Symbol, sigX::Float64 ,sigY::Float64)
    seeds = filter(e->e[:kind]=="seed", df)
    otherDims = filter(e-> !(e in [:index, xSym, ySym]), names(seeds))
    deletecols!(seeds, otherDims)
    rename!(seeds, Dict(xSym => :x, ySym => :y))
    #
    centers = filter(e->e[:kind]=="center", df)
    deletecols!(centers, otherDims)
    rename!(centers, Dict(xSym => :xend, ySym => :yend))

    data = join(seeds, centers, on = :index)
    data[:norm] = map((x1, x2, y1, y2) -> calcNorm(x1, x2, y1, y2),
        data[:x], data[:xend], data[:y], data[:yend])
    return data
end
function doPlotFlow(data::DataFrame, xl::String, yl::String)
    l = layer(data,
        x = :x, y = :y, xend = :xend, yend = :yend, color = :norm,
        Geom.segment(arrow=true), Geom.point,
    )

    xmin = min(min(data[:xend]...), min(data[:x]...))
    xmax = max(max(data[:xend]...), max(data[:x]...))
    ymin = min(min(data[:yend]...), min(data[:y]...))
    ymax = max(max(data[:yend]...), max(data[:y]...))

    coord = Coord.cartesian(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)
    xsc  = Scale.x_continuous(minvalue = xmin, maxvalue = xmax)
    ysc  = Scale.y_continuous(minvalue = ymin, maxvalue = ymax)
    colsc = Scale.color_continuous(minvalue = min(data[:norm]...), maxvalue = max(data[:norm]...))
    return plot(l, xsc, ysc, colsc, coord, Guide.xlabel(xl), Guide.ylabel(yl))
end
function plotFlow(df::DataFrame, xSym::Symbol, ySym::Symbol, varX::Float64 ,varY::Float64)
    data = prepareData(df, xSym, ySym, sqrt(varX), sqrt(varY))
    return doPlotFlow(data, string(xSym), string(ySym))
end

df = CSV.read("$(INDIR)/coords.csv")
vars = filter(e->e[:kind]=="variance", df)[1, 3:end]

pcs = names(df)[3:end]
pcN = 6 #length(pcs)
plotN = convert(Int, pcN * (pcN - 1) / 2)
println("#Plot: $(plotN)")
plots = Array{Plot}(undef, plotN)

plotI = 1
for i in 1:pcN
    for j in (i+1):pcN
        global plotI
        plots[plotI] = plotFlow(df, pcs[i], pcs[j], vars[i], vars[j])
        plotI += 1
    end
end

set_default_plot_size(20cm, plotN*20cm)
p = vstack(plots);

p |> PDF("$(OUTDIR)/$(base_dir)_multidim_map.pdf")

################################################################################
############################ 2. Divergence Analysis ############################

function lincor(x, y)
    n = length(x)
    prod = sum(x .* y)
    ux, sx = mean_and_std(x)
    uy, sy = mean_and_std(y)
    cor = prod / n - ux * uy
    cor = cor / (sx * sy)
    return cor
end
function getKind(df::DataFrame, k::String)
    tmp = filter(e->e[:kind]==k, df)
    deletecols!(tmp, :kind)
    rename!(tmp, Dict(:value => Symbol(k)))
    sort!(tmp, Symbol(k))
    return tmp
end

df = CSV.read("$(INDIR)/distances.csv")
df = filter(e -> !isinf(e[:value]) , df)

df_divs = join(getKind(df, "divergence"), getKind(df, "mle_divergence"),
                   on = [:index1, :index2])
df_divs = join(df_divs, getKind(df, "hist_divergence"), on = [:index1, :index2])
#
s_cor = corspearman(df_divs[:divergence], df_divs[:mle_divergence])
p_cor = lincor(df_divs[:divergence], df_divs[:mle_divergence])
println("[Normal-MLE divergences correlation] Spearman: $(s_cor)\tPearson: $(p_cor) .")
#
s_cor = corspearman(df_divs[:divergence], df_divs[:hist_divergence])
p_cor = lincor(df_divs[:divergence], df_divs[:hist_divergence])
println("[Normal-Hist divergences correlation]\tSpearman: $(s_cor)\tPearson: $(p_cor) .")
#
s_cor = corspearman(df_divs[:mle_divergence], df_divs[:hist_divergence])
p_cor = lincor(df_divs[:mle_divergence], df_divs[:hist_divergence])
println("[MLE-Hist divergences correlation]\tSpearman: $(s_cor)\tPearson: $(p_cor) .")
println("Size of df_divs: $(size(df_divs))")

symbolsStr = ["c2c_full_eucli", "c2c_proj_eucli", "c2c_maha",
    "s2s_full_eucli", "s2s_proj_eucli", "s2s_maha"]
println("\nDistance-MLE_Div correlations:")
for str in symbolsStr
    sym = Symbol(str)
    div_df = join(getKind(df, str), getKind(df, "mle_divergence"), on = [:index1, :index2])
    div_df[:id] = 1:(size(div_df)[1])
    s_cor = corspearman(div_df[sym], div_df[:mle_divergence])
    p_cor = lincor(div_df[sym], div_df[:mle_divergence])
    println("[$(str)]\tspearman: $(s_cor)\tp_cor: $(p_cor) .")
end
#
println("\nDistance-Normal_Div correlations:")
for str in symbolsStr
    sym = Symbol(str)
    div_df = join(getKind(df, str), getKind(df, "divergence"), on = [:index1, :index2])
    div_df[:id] = 1:(size(div_df)[1])
    s_cor = corspearman(div_df[sym], div_df[:divergence])
    p_cor = lincor(div_df[sym], div_df[:divergence])
    println("[$(str)]\tspearman: $(s_cor)\tp_cor: $(p_cor) .")
end
#
println("\nDistance-Hist_Div correlations:")
for str in symbolsStr
    sym = Symbol(str)
    div_df = join(getKind(df, str), getKind(df, "hist_divergence"), on = [:index1, :index2])
    div_df[:id] = 1:(size(div_df)[1])
    s_cor = corspearman(div_df[sym], div_df[:hist_divergence])
    p_cor = lincor(div_df[sym], div_df[:hist_divergence])
    println("[$(str)]\tspearman: $(s_cor)\tp_cor: $(p_cor) .")
end

eucli_df = join(getKind(df, "c2c_full_eucli"), getKind(df, "c2c_proj_eucli"), on = [:index1, :index2])
pca_loss = mean(eucli_df[:c2c_proj_eucli] ./ eucli_df[:c2c_full_eucli])
println("\nPCA loss: $(100*(1-pca_loss))")

max_stars_repo_name: Jiliac/hemipt

max_issues_repo_name: Jiliac/hemipt

max_forks_repo_name: Jiliac/hemipt


text:
using Documenter: Selectors, Expanders, Markdown
using Documenter.Markdown: Link, Paragraph
struct DatabaseLookup <: Expanders.ExpanderPipeline end

Selectors.order(::Type{DatabaseLookup}) = 0.5
Selectors.matcher(::Type{DatabaseLookup}, node, page, doc) = false

const regex_pattern = r"\@example_database\(([\"a-zA-Z_0-9.+ ]+)(?:, )?(plot|code|stepperplot)?\)"

const atomics = (
    heatmap,
    image,
    lines,
    linesegments,
    mesh,
    meshscatter,
    scatter,
    surface,
    text,
    Makie.volume
)

match_kw(x::String) = occursin(regex_pattern, x)
match_kw(x::Paragraph) = any(match_kw, x.content)
match_kw(x::Any) = false

Selectors.matcher(::Type{DatabaseLookup}, node, page, doc) = match_kw(node)

# ============================================= Simon's implementation
function look_up_source(database_key)
    entries = findall(x-> x.title == database_key, database)
    # current implementation finds titles, but we can also search for tags too
    isempty(entries) && error("No entry found for database reference $database_key")
    length(entries) > 1 && error("Multiple entries found for database reference $database_key")
    sprint() do io
        idx = entries[1]
        print_code(
            io, database[entries[1]],
            scope_start = "",
            scope_end = "",
            indent = "",
            resolution = (entry)-> "resolution = (500, 500)",
            outputfile = (entry, ending)-> Pkg.dir("Makie", "docs", "media", string(entry.unique_name, ending))
        )
    end
end

function Selectors.runner(::Type{DatabaseLookup}, x, page, doc)
    matched = nothing
    for elem in x.content
        if isa(elem, AbstractString)
            matched = match(regex_pattern, elem)
            matched != nothing && break
        end
    end
    matched == nothing && error("No match: $x")
    capture = matched[1]
    embed = matched[2]

    # The sandboxed module -- either a new one or a cached one from this page.
    database_keys = filter(x-> !(x in ("", " ")), split(capture, '"'))
    # info("$database_keys with embed option $embed")
    map(database_keys) do database_key
        is_stepper = false

        # buffer for content
        content = []

        # bibliographic stuff
        idx = findall(x-> x.title == database_key, database)
        entry = database[idx[1]]
        uname = string(entry.unique_name)
        lines = entry.file_range
        tags = entry.tags

        # check whether the example is a stepper or not
        is_stepper = any(x-> occursin("stepper", x), tags)

        if is_stepper && isequal(embed, "plot")
            @warn("you tried to output an @example_database(\"entry\", plot) but the entry type is a stepper!")
        end

        # print source code for the example
        if embed == nothing || isequal(embed, "code")
            source = Markdown.Code("julia", look_up_source(database_key))
            src_code = Markdown.MD(source)
            push!(content, src_code)
        end

        # embed plot for the example
        if (embed == nothing && !is_stepper) || isequal(embed, "plot")
            # print to buffer
            io = IOBuffer()
            embed_plot(io, uname, mediapath, buildpath; src_lines = lines, pure_html = true)
            str = String(take!(io))

            # print code for embedding plot
            src_plot = Documenter.Documents.RawHTML(str)
            push!(content, src_plot)
        end

        # use stepper for the example
        if (embed == nothing && is_stepper) || isequal(embed, "stepperplot")

            steps = enumerate_stepper_examples(mediapath, uname; filter = "thumb")
            steppermediapath = joinpath(mediapath, uname)

            # print to buffer
            io = IOBuffer()

            for step = 1:steps
				divblock = """

"""
	            caption = "Step $step
"
	            imgpath = "media/$uname/thumb-$uname-$step.jpg"
	            imgsrc = """

"""
	            println(io, divblock * caption * imgsrc)
                str = String(take!(io))

                # print code for embedding plot
                src_plot = Documenter.Documents.RawHTML(str)
                push!(content, src_plot)
            end
        end

        # finally, map the content back to the page
        page.mapping[x] = content
    end
end

"""
    enumerate_stepper_examples(filter!(x -> !startswith(x, "thumb"), plots_list))

Enumerates the stepper plots from the `mediapath` from the example with the name `uname`.
Accepts an optional `filter` argument to filter for specific files starting with this text.
"""
function enumerate_stepper_examples(mediapath::AbstractString, uname::String; filter = nothing)
    # get all of the plots from the stepper
    steppermediapath = joinpath(mediapath, uname)
    plots_list = readdir(steppermediapath)
    if filter != nothing
        filter!(x -> !startswith(x, "$filter"), plots_list)
    end

    # the plot step indices always start with 1
    steps = length(plots_list)
end

function enumerate_stepper_examples(mediapath::AbstractString, uname::Symbol; filter = nothing)
    enumerate_stepper_examples(mediapath, String(uname); filter = filter)
end

"""
    embed_video(relpath::AbstractString[; pure_html::Bool = false])

Generates a MD-formatted string for embedding video into Documenter Markdown files
(since `Documenter.jl` doesn't support directly embedding mp4's using ![]() syntax).
`raw_mode` specifies whether to print the @raw block for use in Documenter.
"""
function embed_video(relpath::AbstractString; pure_html::Bool = false)
    embed_code = """
        <video controls autoplay loop muted>
          <source src="$(relpath)" type="video/mp4">
          Your browser does not support mp4. Please use a modern browser like Chrome or Firefox.
        </video>
    """
    if !pure_html
        return str = "```@raw html\n" * embed_code * "```"
    else
        return embed_code
    end
end

"""
    embed_thumbnail(io::IO, func::Function, currpath::AbstractString)

Insert thumbnails matching a search tag.
"""
function embed_thumbnail(io::IO, func::Function, currpath::AbstractString)
    indices = find_indices(func)
    !ispath(currpath) && @warn("currepath does not exist!")
    for idx in indices
        uname = database[idx].unique_name
        title = database[idx].title
        # TODO: currently exporting video thumbnails as .jpg because of ImageMagick issue#120
        testpath1 = joinpath(srcmediapath, "thumb-$uname.png")
        testpath2 = joinpath(srcmediapath, "thumb-$uname.jpg")
        if isfile(testpath1)
            embedpath = relpath(testpath1, currpath)
            println(io, "![]($(embedpath))")
            # [![Alt text](/path/to/img.jpg)](http://example.net/)
            # println(io, "[![$title]($(embedpath))](@ref)")
        elseif isfile(testpath2)
            embedpath = relpath(testpath2, currpath)
            println(io, "![]($(embedpath))")
            # println(io, "[![$title]($(embedpath))](@ref)")
        else
            @warn("thumbnail for index $idx with uname $uname not found")
            embedpath = "not_found"
        end
        embedpath = []
    end
end

embed_thumbnail(io::IO, func::Function) = embed_thumbnail(io::IO, func::Function, atomicspath)

"""
    embed_thumbnail_link(io::IO, func::Function, currpath::AbstractString, tarpath::AbstractString)

Insert thumbnails matching a search tag.
"""
function embed_thumbnail_link(io::IO, func::Function, currpath::AbstractString, tarpath::AbstractString)
    indices = find_indices(func)
    !ispath(currpath) && @warn("currepath does not exist!")
    !ispath(tarpath) && @warn("$(tarpath) does not exist! Note that on your first run of docs generation and before you `makedocs`, you will likely get this error.")
    for idx in indices
        is_stepper = false

        entry = database[idx]
        uname = entry.unique_name
        title = entry.title
        src_lines = entry.file_range
        tags = entry.tags

        # check whether the example is a stepper or not
        is_stepper = any(x-> occursin("stepper", x), tags)

        if is_stepper
            steps = enumerate_stepper_examples(mediapath, uname; filter = "thumb")
            println(io, "```@raw html\n")
            for i = 1:steps
                divblock = """<div style="display:inline-block"><p style="display:inline-block; text-align: center">"""
                caption = "Step $i<br>"
                imgpath = "media/$uname/thumb-$uname-$i.jpg"
                imgsrc = """<img src="$imgpath" alt="$caption" /></p></div>"""
                println(io, divblock * caption * imgsrc)
            end
            println(io, "```")
            break
        else
            # TODO: currently exporting video thumbnails as .jpg because of ImageMagick issue#120
            testpath1 = joinpath(mediapath, "thumb-$uname.png")
            testpath2 = joinpath(mediapath, "thumb-$uname.jpg")
            link = relpath(tarpath, currpath)
            if isfile(testpath1)
                embedpath = relpath(testpath1, currpath)
                println(io, "[![library lines $(src_lines)]($(embedpath))]($(link))")
            elseif isfile(testpath2)
                embedpath = relpath(testpath2, currpath)
                println(io, "[![library lines $(src_lines)]($(embedpath))]($(link))")
            else
                @warn("thumbnail for index $idx with uname $uname not found")
                embedpath = "not_found"
            end
        end
        embedpath = []
    end
end

# embed_thumbnail_link(io::IO, func::Function) = embed_thumbnail_link(io::IO, func::Function, atomicspath)

"""
    embed_plot(io::IO, uname::AbstractString, mediapath::AbstractString, buildpath::AbstractString[;
    pure_html::Bool = false])

Outputs markdown code for embedding plots in `Documenter.jl`.
"""
function embed_plot(
        io::IO,
        uname::AbstractString,
        mediapath::AbstractString,
        buildpath::AbstractString;
        src_lines::AbstractRange = nothing,
        pure_html::Bool = false
    )
    isa(uname, AbstractString) ? nothing : error("uname must be a string!")
    isa(mediapath, AbstractString) ? nothing : error("mediapath must be a string!")
    isa(buildpath, AbstractString) ? nothing : error("buildpath must be a string!")
    medialist = readdir(mediapath)

    if ("$(uname)") in medialist
        steps = enumerate_stepper_examples(mediapath, uname; filter = "thumb")
        println(io, "```@raw html\n")
        for i = 1:steps
            divblock = """<div style="display:inline-block"><p style="display:inline-block; text-align: center">"""
            caption = "Step $i<br>"
            imgpath = "media/$uname/thumb-$uname-$i.jpg"
            imgsrc = """<img src="$imgpath" alt="$caption" /></p></div>"""
            println(io, divblock * caption * imgsrc)
        end
        println(io, "```")
    end
    extensions = [".jpg", ".gif"]
    for i in extensions
        if ("$(uname)" * i) in medialist
            println("$(uname)" * i)
            embedpath = joinpath(relpath(mediapath, buildpath), "$(uname)" * i)
            if pure_html
                embed_code = """

                """
                println(io, embed_code)
            else
                println(io, "![library lines $(src_lines)]($(embedpath))")
            end
            break
        elseif "$(uname).mp4" in medialist
            println("$(uname).mp4")
            embedcode = embed_video(joinpath(relpath(mediapath, buildpath), "$(uname).mp4"); pure_html = pure_html)
            println(io, embedcode)
            break
        else
            @warn("file $(uname)$i with unknown extension in mediapath, or file nonexistent")
        end
    end
    print(io, "\n")
end

"""
    print_table(io::IO, dict::Dict)

Print a Markdown-formatted table with the entries from `dict` to specified `io`.
"""
function print_table(io::IO, dict::Dict)
    # get max length of the keys
    k = string.("`", collect(keys(dict)), "`")
    maxlen_k = max(length.(k)...)

    # get max length of the values
    v = string.(collect(values(dict)))
    maxlen_v = max(length.(v)...)

    j = sort(collect(dict), by = x -> x[1])

    # column labels
    labels = ["Symbol", "Description"]

    # print top header
    print(io, "|")
    print(io, "$(labels[1])")
    print(io, " "^(maxlen_k - length(labels[1])))
    print(io, "|")
    print(io, "$(labels[2])")
    print(io, " "^(maxlen_v - length(labels[2])))
    print(io, "|")
    print(io, "\n")

    # print second line (toprule)
    print(io, "|")
    print(io, "-"^maxlen_k)
    print(io, "|")
    print(io, "-"^maxlen_v)
    print(io, "|")
    print(io, "\n")

    for (idx, entry) in enumerate(j)
        print(io, "|")
        print(io, "`$(entry[1])`")
        print(io, " "^(maxlen_k - length(string(entry[1])) - 2))
        print(io, "|")
        print(io, "$(entry[2])")
        print(io, " "^(maxlen_v - length(entry[2])))
        print(io, "|")
        print(io, "\n")
    end
end

using Makie, ImageTransformations, FileIO
using ImageFiltering  # needed for Gaussian-filtering images during resize

function rescale_image(path::AbstractString, target_path::AbstractString, sz::Int = 200)
    !isfile(path) && @warn("Input argument must be a file!")
    img = FileIO.load(path)

    # calculate new image size `newsz`
    (height, width) = size(img)
    (scale_height, scale_width) = sz ./ (height, width)
    scale = min(scale_height, scale_width)
    newsz = round.(Int, (height, width) .* scale)

    # filter image + resize image
    gaussfactor = 0.4
    σ = map((o,n) -> gaussfactor*o/n, size(img), newsz)
    kern = KernelFactors.gaussian(σ)   # from ImageFiltering
    imgf = ImageFiltering.imfilter(img, kern, NA())
    newimg = ImageTransformations.imresize(imgf, newsz)
    # save image
    FileIO.save(target_path, newimg)
end

"""
    generate_thumbnail(path::AbstractString, target_path, thumb_size::Int = 200)

Generates a (proportionally-scaled) thumbnail with maximum side dimension `sz`.
`sz` must be an integer, and the default value is 200 pixels.
"""
function generate_thumbnail(path, thumb_path, thumb_size = 128)
    if any(ext-> endswith(path, ext), (".png", ".jpeg", ".jpg"))
        rescale_image(path, thumb_path, thumb_size)
    elseif any(ext-> endswith(path, ext), (".gif", ".mp4", ".webm"))
        seektime = get_video_duration(path) / 2
        println(thumb_path)
        run(`ffmpeg -loglevel quiet -ss $seektime -i $path -vframes 1 -vf "scale=$(thumb_size):-2" -y -f image2 $thumb_path`)
    else
        @warn("Unsupported return file format in $path")
    end
end

"""
    get_video_duration(path::AbstractString)

Returns the duration of the video in seconds (Float32).
Accepted file types: mp4, mkv, and gif.

Requires `ffprobe` (usually comes installed with `ffmpeg`).

Note that while this accepts gif, it will not work to get duration of the gif
(`ffprobe` doesn't support that), so it will just fallback to return 0.5 sec.
"""
function get_video_duration(path::AbstractString)
    !isfile(path) && error("input is not a file!")
    accepted_exts = ("mp4", "gif", "mkv")
    filename = basename(path)
    !(split(filename, ".")[2] in accepted_exts) && error("accepted file types are mp4 and gif!")
    try
        dur = read(`ffprobe -loglevel quiet -print_format compact=print_section=0:nokey=1:escape=csv -show_entries format=duration -i "$(path)"`, String)
        dur = parse(Float32, dur)
    catch e
        @warn("`get_video_duration` on $filename did not work, using fallback video duration of 0.5 seconds")
        dur = 0.5
    end
end


text:
import FFTW
import GPUifyLoops: @launch, @loop, @synchronize

function init_poisson_solver(::CPU, g::Grid, tmp_rhs)
    tmp_rhs.data .= rand(Float64, g.Nx, g.Ny, g.Nz)
    PoissonSolver(g, tmp_rhs, FFTW.MEASURE)
end

function init_poisson_solver(::GPU, g::Grid, tmp_rhs)
    tmp_rhs.data .= CuArray{Complex{Float64}}(rand(Float64, g.Nx, g.Ny, g.Nz))
    PoissonSolverGPU(g, tmp_rhs)
end

# Translations to print FFT timing.
let pf2s = Dict(FFTW.ESTIMATE   => "FFTW.ESTIMATE",
                FFTW.MEASURE    => "FFTW.MEASURE",
                FFTW.PATIENT    => "FFTW.PATIENT",
                FFTW.EXHAUSTIVE => "FFTW.EXHAUSTIVE")
    global plannerflag2string
    plannerflag2string(k::Integer) = pf2s[Int(k)]
end

"""
    PoissonSolver(grid, example_field, planner_flag; verbose=false)

Return a `PoissonSolver` on `grid`, using `example_field` and `planner_flag`
to plan fast transforms.
"""
struct PoissonSolver{T<:AbstractArray} <: AbstractPoissonSolver
    kx²::T
    ky²::T
    kz²::T
    FFT!
    DCT!
    IFFT!
    IDCT!
end

function PoissonSolver(g::Grid, exfield::CellField, planner_flag=FFTW.PATIENT; verbose=false)
    kx² = zeros(eltype(g), g.Nx)
    ky² = zeros(eltype(g), g.Ny)
    kz² = zeros(eltype(g), g.Nz)

    for i in 1:g.Nx; kx²[i] = (2sin((i-1)*π/g.Nx)    / (g.Lx/g.Nx))^2; end
    for j in 1:g.Ny; ky²[j] = (2sin((j-1)*π/g.Ny)    / (g.Ly/g.Ny))^2; end
    for k in 1:g.Nz; kz²[k] = (2sin((k-1)*π/(2g.Nz)) / (g.Lz/g.Nz))^2; end

    if verbose
        print("Planning Fourier transforms... (planner_flag=$(plannerflag2string(planner_flag)))\n")
        print("FFT!:  "); @time FFT!  = FFTW.plan_fft!(exfield.data, [1, 2]; flags=planner_flag)
        print("IFFT!: "); @time IFFT! = FFTW.plan_ifft!(exfield.data, [1, 2]; flags=planner_flag)
        print("DCT!:  "); @time DCT!  = FFTW.plan_r2r!(exfield.data, FFTW.REDFT10, 3; flags=planner_flag)
        print("IDCT!: "); @time IDCT! = FFTW.plan_r2r!(exfield.data, FFTW.REDFT01, 3; flags=planner_flag)
    else
        FFT!  = FFTW.plan_fft!(exfield.data, [1, 2]; flags=planner_flag)
        IFFT! = FFTW.plan_ifft!(exfield.data, [1, 2]; flags=planner_flag)
        DCT!  = FFTW.plan_r2r!(exfield.data, FFTW.REDFT10, 3; flags=planner_flag)
        IDCT! = FFTW.plan_r2r!(exfield.data, FFTW.REDFT01, 3; flags=planner_flag)
    end

    PoissonSolver{Array{eltype(g),1}}(kx², ky², kz², FFT!, DCT!, IFFT!, IDCT!)
end

"""
    solve_poisson_3d_ppn_planned!(args...)

Solve Poisson equation with Periodic, Periodic, Neumann boundary conditions in x, y, z using planned
FFTs and DCTs.

  Args
  ----

      solver : PoissonSolver
           g : solver grid
           f : RHS to Poisson equation
           ϕ : Solution to Poisson equation
"""
function solve_poisson_3d_ppn_planned!(solver::PoissonSolver, g::RegularCartesianGrid, f::CellField, ϕ::CellField)
    solver.DCT!*f.data  # Calculate DCTᶻ(f) in place.
    solver.FFT!*f.data  # Calculate FFTˣʸ(f) in place.

    for k in 1:g.Nz, j in 1:g.Ny, i in 1:g.Nx
        @inbounds ϕ.data[i, j, k] = -f.data[i, j, k] / (solver.kx²[i] + solver.ky²[j] + solver.kz²[k])
    end
    ϕ.data[1, 1, 1] = 0

    solver.IFFT!*ϕ.data  # Calculate IFFTˣʸ(ϕ) in place.
    solver.IDCT!*ϕ.data  # Calculate IDCTᶻ(ϕ) in place.
    @. ϕ.data = ϕ.data / (2*g.Nz)
    nothing
end

"""
    PoissonSolverGPU(grid, example_field)

Return a `PoissonSolverGPU` on `grid`, using `example_field` to plan
CuFFTs on a GPU.
"""
struct PoissonSolverGPU{T<:AbstractArray} <: AbstractPoissonSolver
    kx²
    ky²
    kz²
    dct_factors
    idct_bfactors
    FFT_xy!
    FFT_z!
    IFFT_xy!
    IFFT_z!
end

function PoissonSolverGPU(g::Grid, exfield::CellField; verbose=false)
    kx² = CuArray{Float64}(undef, g.Nx)
    ky² = CuArray{Float64}(undef, g.Ny)
    kz² = CuArray{Float64}(undef, g.Nz)

    for i in 1:g.Nx; kx²[i] = (2sin((i-1)*π/g.Nx)    / (g.Lx/g.Nx))^2; end
    for j in 1:g.Ny; ky²[j] = (2sin((j-1)*π/g.Ny)    / (g.Ly/g.Ny))^2; end
    for k in 1:g.Nz; kz²[k] = (2sin((k-1)*π/(2g.Nz)) / (g.Lz/g.Nz))^2; end

    # Exponential factors required to calculate the DCT on the GPU.
    factors = 2 * exp.(collect(-1im*π*(0:g.Nz-1) / (2*g.Nz)))
    dct_factors = CuArray{Complex{Float64}}(repeat(reshape(factors, 1, 1, g.Nz), g.Nx, g.Ny, 1))

    # "Backward" exponential factors required to calculate the IDCT on the GPU.
    bfactors = exp.(collect(1im*π*(0:g.Nz-1) / (2*g.Nz)))
    bfactors[1] *= 0.5
    idct_bfactors = CuArray{Complex{Float64}}(repeat(reshape(bfactors, 1, 1, g.Nz), g.Nx, g.Ny, 1))

    if verbose
        print("Creating CuFFT plans...\n")
        print("FFT_xy!:  "); @time FFT_xy!  = plan_fft!(exfield.data, [1, 2])
        print("FFT_z!:   "); @time FFT_z!   = plan_fft!(exfield.data, 3)
        print("IFFT_xy!: "); @time IFFT_xy! = plan_ifft!(exfield.data, [1, 2])
        print("IFFT_z!:  "); @time IFFT_z!  = plan_ifft!(exfield.data, 3)
    else
        FFT_xy!  = plan_fft!(exfield.data, [1, 2])
        FFT_z!   = plan_fft!(exfield.data, 3)
        IFFT_xy! = plan_ifft!(exfield.data, [1, 2])
        IFFT_z!  = plan_ifft!(exfield.data, 3)
    end

    PoissonSolverGPU{CuArray{Float64}}(kx², ky², kz², dct_factors, idct_bfactors, FFT_xy!, FFT_z!, IFFT_xy!, IFFT_z!)
end

"""
    solve_poisson_3d_ppn_gpu_planned!(args...)

Solve Poisson equation with Periodic, Periodic, Neumann boundary conditions in x, y, z using planned
CuFFTs on a GPU.

  Args
  ----

      Tx, Ty : Thread size in x, y
  Bx, By, Bz : Block size in x, y, z
      solver : PoissonSolverGPU
           g : solver grid
           f : RHS to Poisson equation
           ϕ : Solution to Poisson equation
"""
function solve_poisson_3d_ppn_gpu_planned!(Tx, Ty, Bx, By, Bz, solver::PoissonSolverGPU, g::RegularCartesianGrid, f::CellField, ϕ::CellField)
    # Calculate DCTᶻ(f) in place using the FFT.
    solver.FFT_z! * f.data
    f.data .*= solver.dct_factors
    @. f.data = real(f.data)

    solver.FFT_xy! * f.data  # Calculate FFTˣʸ(f) in place.

    @launch device(GPU()) f2ϕ!(g, f.data, ϕ.data, solver.kx², solver.ky², solver.kz², threads=(Tx, Ty), blocks=(Bx, By, Bz))
    ϕ.data[1, 1, 1] = 0

    solver.IFFT_xy! * ϕ.data  # Calculate IFFTˣʸ(ϕ̂) in place.

    # Calculate IDCTᶻ(ϕ̂) in place using the FFT.
    ϕ.data .*= solver.idct_bfactors
    solver.IFFT_z! * ϕ.data
    nothing
end

"Kernel for computing the solution `ϕ` to Poisson equation for source term `f` on a GPU."
function f2ϕ!(grid::Grid, f, ϕ, kx², ky², kz²)
    @loop for k in (1:grid.Nz; blockIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds ϕ[i, j, k] = -f[i, j, k] / (kx²[i] + ky²[j] + kz²[k])
            end
        end
    end
    @synchronize
end


text:
using Unitful, ImagineInterface, Imagine
using Base.Test

#modified from ImagineInterface test
srate = 10000*Unitful.s^-1
pmin = 0.0*Unitful.μm
pmax = 200.0*Unitful.μm
stack_time = 3.0*Unitful.s
reset_time = 2.0*Unitful.s
exp_time = 0.011*Unitful.s
flash_frac = 0.1 #fraction of time to keep laser on during exposure
z_spacing = 3.1*Unitful.μm
z_pad = 5.0*Unitful.μm
nstacks = 10
d = gen_unidirectional_stack(pmin, pmax, z_spacing, stack_time, stack_time, exp_time, srate, flash_frac; z_pad = z_pad)

ocpi2 = rigtemplate("ocpi-2"; sample_rate = srate)
pos = getpositioners(ocpi2)[1]
pos_mon = getname(ocpi2, monitor_name(pos))
append!(pos, "uni_stack_pos", d["positioner"])
replicate!(pos, nstacks-1)
#write_commands("test.json", ocpi2, nstacks, nframes, exp_time; isbidi = false)

aos = [pos;]
#ais = [pos_mon;]
#ais = ImagineSignal[]

#prepare_ao(aos, bufsz::Int, trigger_terminal::String)
ImagineWorker._set_device("Dev2/")
#tsk = ImagineWorker.prepare_ao(aos, 100000, "PFI1")
#print("starting...\n")
#start(tsk)

ready_chan = RemoteChannel(()->Channel{Int}(1))
ImagineWorker._output_analog_signals(aos, 50000, "PFI1", ready_chan)

#and execute this from another julia instance:
#using NIDAQ, Imagine; Imagine.ttl_pulse(; line_name = "Port2/Line0")


text:
using Plotly

# Signin with your userid and password
# Plotly.signin("myuserid", "abc32def7g")

Plotly.signin("sherrinm", "b8vc9bz7en")

trace1 = [
  "x" => [0, 1, 2, 3, 4, 5, 6, 7, 8],
  "y" => [8, 7, 6, 5, 4, 3, 2, 1, 0],
  "type" => "scatter"
];
trace2 = [
  "x" => [0, 1, 2, 3, 4, 5, 6, 7, 8],
  "y" => [0, 1, 2, 3, 4, 5, 6, 7, 8],
  "type" => "scatter"
];
data = [trace1, trace2];

layout = [
  "xaxis" => ["type" => "log", "autorange" => true],
  "yaxis" => ["type" => "log", "autorange" => true]
];

response = Plotly.plot(data,
		["layout" => layout,
		 "filename" => "plotly-log-axes",
		 "fileopt" => "overwrite"]);

plot_url = response["url"]