text:
# Fractional Brownian motion
immutable FBM <: ContinuousUnivariateStochasticProcess
  t::Vector{Float64}
  n::Int64
  h::Float64 # Hurst index

  function FBM(t::Vector{Float64}, n::Int64, h::Float64)
    t[1] == 0. || error("Starting time point must be equal to 0.0.")
    issorted(t, lt=<=) || error("The time points must be strictly sorted.")
    int64(length(t)) == n || error("Number of time points must be equal to the vector holding the time points.")
    0. < h < 1. || error("Hurst index must be between 0 and 1.")
    new(t, n, h)
  end
end

FBM(t::Vector{Float64}, h::Float64) = FBM(t, int64(length(t)), h)
FBM(t::Ranges, h::Float64) = FBM(collect(t), int64(length(t)), h)
FBM(t::Float64, n::Int64, h::Float64) = FBM(t/n:t/n:t, h)
FBM(t::Float64, h::Float64) = FBM([t], 1, h)

FBM(t::Matrix{Float64}, h::Float64) = FBM[FBM(t[:, i], h) for i = 1:size(t, 2)]
FBM(t::Ranges, np::Int, h::Float64) = FBM[FBM(t, h) for i = 1:np]
FBM(t::Float64, n::Int64, np::Int, h::Float64) = FBM[FBM(t, n, h) for i = 1:np]

# Fractional Gaussian noise
immutable FGN <: ContinuousUnivariateStochasticProcess
  σ::Float64
  h::Float64 # Hurst index

  function FGN(σ::Float64, h::Float64)
    σ > 0. || error("Standard deviation must be positive.")
    0. < h < 1. || error("Hurst index must be between 0 and 1.")
    new(σ, h)
  end
end

FGN(h::Float64) = FGN(1., h)

convert(::Type{FBM}, p::BrownianMotion) = FBM(p.t, p.n, 0.5)
convert(::Type{FGN}, p::BrownianMotion) = FGN((p.t[end]/(p.n-1))^0.5, 0.5)
convert(::Type{FGN}, p::FBM) = FGN((p.t[end]/(p.n-1))^p.h, p.h)

function autocov(p::FGN, i::Int64, j::Int64)
  twoh::Float64 = 2*p.h
  0.5*abs2(p.σ)*(abs(j-i+1)^twoh+abs(j-i-1)^twoh-2*abs(j-i)^twoh)
end

function autocov!(c::Matrix{Float64}, p::FGN)
  n::Int64 = size(c, 1)

  for i = 1:n
    for j = 1:i
      c[i, j] = autocov(p, i, j)
    end
  end

  for i = 1:n
    for j = (i+1):n
      c[i, j] = c[j, i]
    end
  end

  c
end

autocov(p::FGN, maxlag::Int64) = autocov!(Array(Float64, maxlag, maxlag), p)

function autocov!(y::Vector{Float64}, p::FGN, lags::IntegerVector)
  nlags = length(lags)
  twoh::Float64 = 2*p.h

  for i = 1:nlags
    y[i] = 0.5*abs2(p.σ)*(abs(lags[i]+1)^twoh+abs(lags[i]-1)^twoh-2*abs(lags[i])^twoh)
  end

  y
end

autocov(p::FGN, lags::IntegerVector) = autocov!(Array(Float64, length(lags)), p, lags)

function autocov(p::FBM, i::Int64, j::Int64)
  twoh::Float64 = 2*p.h
  0.5*((p.t[i])^twoh+(p.t[j])^twoh-abs(p.t[i]-p.t[j])^twoh)
end

function autocov!(c::Matrix{Float64}, p::FBM)
  n::Int64 = p.n-1

  for i = 1:n
    for j = 1:i
      c[i, j] = autocov(p, i+1, j+1)
    end
  end

  for i = 1:n
    for j = (i+1):n
      c[i, j] = c[j, i]
    end
  end

  c
end

function autocov(p::FBM)
  n::Int64 = p.n-1
  autocov!(Array(Float64, n, n), p)
end

### rand_chol generates FBM using the method based on Cholesky decomposition.
### T. Dieker, Simulation of Fractional Brownian Motion, master thesis, 2004.
### The complexity of the algorithm is O(n^3), where n is the number of FBM samples.
function rand_chol(p::FBM; fbm::Bool=true)
  w = chol(autocov(p), :L)*randn(p.n-1)

  # If fbm is true return FBM, otherwise return FGN
  w = [0., w]
  if !fbm
    w = diff(w)
  end

  w
end

function rand_chol(p::Vector{FBM}; fbm::Bool=true)
  n::Int64 = fbm ? p[1].n : p[1].n-1
  np = length(p)
  x = Array(Float64, n, np)

  for i = 1:np
    x[:, i] = rand_chol(p[i], fbm=fbm)
  end

  x
end

### rand_fft generates FBM using fast Fourier transform (FFT).
### The time interval of FBM is [0, 1] with a stepsize of 2^p, where p is a natural number.
### The algorithm is known as the Davies-Harte method or the method of circular embedding.
### R.B. Davies and D.S. Harte, Tests for Hurst Effect, Biometrika, 74 (1), 1987, pp. 95-102.
### For a more recent publication, see P.F. Craigmile, Simulating a Class of Stationary Gaussian Processes Using the
### Davies–Harte Algorithm, With Application to Long Memory Processes, Journal of Time Series Analysis, 24 (5), 2003,
### pp. 505-511.
### The complexity of the algorithm is O(n*log(n)), where n=2^p is the number of FBM samples.
function rand_fft(p::FBM; fbm::Bool=true)
  # Determine number of points of simulated FBM
  pnmone::Int64 = p.n-1
  n::Int64 = 2^ceil(log2(pnmone))

  # Compute autocovariance sequence of underlying FGN
  c = Array(Float64, n+1)
  autocov!(c, FGN((1/n)^p.h, p.h), 0:n)

  # Compute square root of eigenvalues of circular autocovariance sequence
  l = real(fft([c, c[end-1:-1:2]]))
  all(i->(i>0), l) || error("Non-positive eigenvalues encountered.")
  lsqrt = sqrt(l)

  # Simulate standard random normal variables
  twon::Int64 = 2*n
  z = randn(twon)

  # Generate fractional Gaussian noise (retain only the first p.n-1 values)
  x = sqrt(0.5)*lsqrt[2:n].*complex(z[2*(2:n)-2], z[2*(2:n)-1])
  w = real(bfft([lsqrt[1]*z[1], x, lsqrt[n+1]*z[twon], conj(reverse(x))]))[1:pnmone]/sqrt(twon)

  # If fbm is true return FBM, otherwise return FGN
  if fbm
    w = [0., cumsum(w)]
  end

  w
end

function rand_fft(p::Vector{FBM}; fbm::Bool=true)
  n::Int64 = fbm ? p[1].n : p[1].n-1
  np = length(p)
  x = Array(Float64, n, np)

  for i = 1:np
    x[:, i] = rand_fft(p[i], fbm=fbm)
  end

  x
end

function rand(p::FBM; fbm::Bool=true, rtype::Symbol=:fft)
  if rtype == :chol
    rand_chol(p; fbm=fbm)
  elseif rtype == :fft
    rand_fft(p; fbm=fbm)
  else
  end
end

function rand(p::Vector{FBM}; fbm::Bool=true, rtype::Symbol=:fft)
  if rtype == :chol
    rand_chol(p; fbm=fbm)
  elseif rtype == :fft
    rand_fft(p; fbm=fbm)
  else
  end
end


text:
#-----------------------------------------------------------------------------------------------------
# Tensor Basis
#-----------------------------------------------------------------------------------------------------
export TensorBasis
"""
    TensorBasis

Basis without any symmetries. The index is computed by:

`I(i₁, i₂, ⋯, iₙ) = i₁ B^(n-1) + i₂ B^(n-2) + ⋯ + iₙ`
"""
struct TensorBasis <: AbstractBasis
    dgt::Vector{Int64}
    B::Int64
    TensorBasis(dgt::Vector{Int64}, B::Integer) = new(dgt, Int64(B))
    TensorBasis(dgt::AbstractVector{<:Integer}, B::Integer) = new(collect(Int64, dgt), Int64(B))
    TensorBasis(L::Integer; base::Integer=2) = new(zeros(Int64, L), Int64(base))
end

@inline content(b::TensorBasis) = 1:b.B^length(b.dgt)
@inline content(::TensorBasis, i::Integer) = i
@inline norm(b::TensorBasis) = ones(Int64, B^length(b.dgt))
@inline norm(b::TensorBasis, ::Integer) = 1
@inline eltype(::TensorBasis) = Int64
@inline size(b::TensorBasis, i::Integer) = isone(i) || isequal(i, 2) ? b.B^length(b.dgt) : 1
@inline size(b::TensorBasis) = (l=b.B^length(b.dgt); (l, l))
@inline index(b::TensorBasis) = 1, index(b.dgt, base=b.B)

function schmidt(v::AbstractVector, Ainds::AbstractVector{<:Integer}, b::TensorBasis)
    S = SchmidtMatrix(eltype(v), b, Ainds)
    for i = 1:length(v)
        change!(b, i)
        addto!(S, v[i])
    end
    Matrix(S)
end

#-------------------------------------------------------------------------------------------------------------------------
# Projected Basis
#-------------------------------------------------------------------------------------------------------------------------
export ProjectedBasis
"""
    ProjectedBasis

Basis for subspace that is spanned only by product states.
"""
struct ProjectedBasis <: AbstractBasis
    dgt::Vector{Int64}
    I::Vector{Int64}
    B::Int64
    ProjectedBasis(dgt::Vector{Int64}, I::Vector{Int64}, B::Integer) = new(dgt, I, Int64(B))
end

@inline eltype(::ProjectedBasis) = Int
@inline change!(b::ProjectedBasis, i::Integer) = (change!(b.dgt, b.I[i], base=b.B); 1)
@inline function index(b::ProjectedBasis)
    i = index(b.dgt, base=b.B)
    ind = binary_search(b.I, i)
    ind > 0 ? (1, ind) : error("No such symmetry.")
end

function ProjectedBasis(f, L::Integer; base::Integer=2, alloc::Integer=1000, threaded::Bool=false)
    dgt = zeros(Int64, L)
    I = threaded ? selectindex_threaded(f, L, base=base, alloc=alloc) : selectindex(f, L, 1:base^L, base=base, alloc=alloc)
    ProjectedBasis(dgt, I, base)
end

function schmidt(v::AbstractVector, Ainds::AbstractVector{<:Integer}, b::ProjectedBasis)
    S = SchmidtMatrix(eltype(v), b, Ainds)
    for i = 1:length(v)
        change!(b, i)
        addto!(S, v[i])
    end
    Matrix(S)
end

#-------------------------------------------------------------------------------------------------------------------------
# TranslationalBasis
#-------------------------------------------------------------------------------------------------------------------------
export TranslationalBasis
"""
    TranslationalBasis

Basis for subspace that is spanned by momentum states, which can also incorporate projected restriction.

Properties:
-----------
- `dgt`: Digits.
- `I`  : List of indicies.
- `R`  : List of normalization.
- `C`  : Unit phase factor.
- `B`  : Base.
"""
struct TranslationalBasis{T <: Number} <: AbstractBasis
    dgt::Vector{Int64}
    I::Vector{Int}
    R::Vector{Float64}
    C::Vector{T}
    B::Int64
    TranslationalBasis(dgt::Vector{Int64}, I, R, C::Vector{T}, B::Integer) where T = new{T}(dgt, I, R, C, Int64(B))
end

@inline eltype(::TranslationalBasis{T}) where T = T

"""
    index(b::TranslationalBasis)

To calculate the index of a given digits, we first shift the digits to that of the minimum index, then search for the place in the basis.
Because of the restriction from the momentum, some state with zero normalization (which is not included in the basis) will appear.
To avoid exception in the matrix constructon of `Operation`, we allow the index to not in the basis content.
When this happend, we return index 1, and normalization 0, so it has no effect on the matrix being filled.
"""
@inline function index(b::TranslationalBasis)
    Im, T = translation_index(b.dgt, b.B)
    i = binary_search(b.I, Im)
    N, ind = if iszero(i)
        # Here we allow the case where i is not in the basis.
        # In such case we add 0 to index 1.
        zero(eltype(b)), 1
    else
        b.C[T+1] * b.R[i], i
    end
    N, ind
end

struct TranslationJudge
    F                   # Projective selection
    K::Int              # Momentum
    B::Int              # Base
    C::Vector{Float64}  # Normalization coefficient
end

function (judge::TranslationJudge)(dgt::AbstractVector{<:Integer}, i::Integer)
    Q, N = if judge.F(dgt)
        c, r = translation_check(dgt, i, judge.B)
        if c && iszero( mod(r * judge.K, length(dgt) ) )
            true, judge.C[r]
        else
            false, 0.0
        end
    else
        false, 0.0
    end
    Q, N
end

function TranslationalBasis(f, k::Integer, L::Integer; base::Integer=2, alloc::Integer=1000, threaded::Bool=false)
    dgt = zeros(Int64, L)
    judge = TranslationJudge(f, k, base, [L/sqrt(i) for i = 1:L])
    I, R = if threaded
        selectindexnorm_threaded(judge, L, base=base, alloc=alloc)
    else
        selectindexnorm(judge, L, 1:base^L, base=base, alloc=alloc)
    end
    C = if iszero(k)
        fill(1.0, L)
    elseif isequal(2k, L)
        [iseven(i) ? 1.0 : -1.0 for i=0:L-1]
    else
        [exp(-1im*2π/L*k*i) for i=0:L-1]
    end
    TranslationalBasis(dgt, I, R, C, base)
end

function TranslationalBasis(k::Integer, L::Integer; base::Integer=2, alloc::Integer=1000, threaded::Bool=false)
    TranslationalBasis(x -> true, k, L, base=base, alloc=alloc, threaded=threaded)
end

function schmidt(v::AbstractVector, Ainds::AbstractVector{<:Integer}, b::TranslationalBasis)
    dgt, R, phase = b.dgt, b.R, b.C[2]
    S = SchmidtMatrix(promote_type(eltype(v), eltype(b)), b, Ainds)
    for i = 1:length(v)
        change!(b, i)
        val = v[i] / R[i]
        for j in 1:length(dgt)
            addto!(S, val)
            cyclebits!(dgt)
            val *= phase
        end
    end
    Matrix(S)
end

#-------------------------------------------------------------------------------------------------------------------------
# Translational + Parity Basis
#-------------------------------------------------------------------------------------------------------------------------
export TranslationParityBasis, TranslationFlipBasis
abstract type AbstractTranslationalParityBasis <: AbstractBasis end
"""
    TranslationParityBasis

Basis with translational and reflection symmetries.

Properties:
-----------
- `dgt`: Digits.
- `I`  : Representing states.
- `R`  : Normalization.
- `K`  : {±1}, momentum phase.
- `P`  : {±1}, parity.
- `B`  : Base.
"""
struct TranslationParityBasis <: AbstractTranslationalParityBasis
    dgt::Vector{Int64}    # Digits
    I::Vector{Int}          # Representing states
    R::Vector{Float64}      # Normalization
    C::Vector{Int}          # Momentum phase exp(1im * 0/π) = {±1}
    P::Int                  # {±1}, parity
    B::Int64                # Base
    TranslationParityBasis(dgt, I, R, C, P, B::Integer) = new(dgt, I, R, C, P, Int64(B))
end

"""
    TranslationFlipBasis

Basis with translational and spin-flip symmetries.

Properties:
-----------
- `dgt`: Digits.
- `I`  : Representing states.
- `R`  : Normalization.
- `C`  : Momentum phase.
- `P`  : {±1}, parity.
- `B`  : Base.
"""
struct TranslationFlipBasis{T} <: AbstractTranslationalParityBasis
    dgt::Vector{Int64}    # Digits
    I::Vector{Int}          # Representing states
    R::Vector{Float64}      # Normalization
    C::Vector{T}            # Momentum phase exp(-1im * 2π/L * K)
    P::Int                  # {±1}, parity
    B::Int                  # Base
end

@inline eltype(::TranslationParityBasis) = Float64
@inline eltype(::TranslationFlipBasis{T}) where T = T

function translation_parity_index(parity, b::AbstractTranslationalParityBasis)
    reflect, i, t = translation_parity_index(parity, b.dgt, b.B)
    ind = binary_search(b.I, i)
    N, I = if iszero(ind)
        0.0, 1
    else
        n = reflect ? b.P * b.C[t+1] : b.C[t+1]
        n * b.R[ind], ind
    end
    N, I
end

@inline index(b::TranslationParityBasis) = translation_parity_index(reverse, b)
@inline index(b::TranslationFlipBasis) = translation_parity_index(x -> spinflip(x, b.B), b)

struct TranslationParityJudge
    F                   # Projective selection
    parity              # Parity operation acting on digits
    K::Int64            # Momentum
    P::Int64            # Parity eigenvalue {±1}
    L::Int64            # Length
    B::Int64            # Base
    C::Vector{Float64}  # Normalization coefficients
end

function (judge::TranslationParityJudge)(dgt::AbstractVector{<:Integer}, i::Integer)
    Q, N = if judge.F(dgt)
        isrep, reflect, r, m = translation_parity_check(judge.parity, dgt, i, judge.B)
        if isrep && iszero(mod(judge.K * r, judge.L))
            if reflect
                trans = mod(judge.K * m, judge.L)
                if ( isequal(2 * trans, judge.L) && isequal(judge.P, -1) ) || ( iszero(trans) && isone(judge.P) )
                    true, judge.C[judge.L+r]
                else
                    false, 0.0
                end
            else
                true, judge.C[r]
            end
        else
            false, 0.0
        end
    else
        false, 0.0
    end
    Q, N
end

function TranslationParityBasis(
    f, K::Integer, P::Integer, L::Integer;
    base::Integer=2, alloc::Integer=1000, threaded::Bool=false
)
    @assert iszero(K) || isequal(2K, L) "Momentum incompatible with parity."
    @assert isone(P) || isequal(P, -1) "Invalid parity"
    judge = begin
        N2 = [2L / sqrt(i) for i = 1:L]
        N1 = N2 ./ sqrt(2)
        TranslationParityJudge(f, reverse, K, P, L, base, vcat(N1, N2))
    end

    I, R = if threaded
        selectindexnorm_threaded(judge, L, base=base, alloc=alloc)
    else
        selectindexnorm(judge, L, 1:base^L, base=base, alloc=alloc)
    end

    C = if iszero(K)
        fill(1.0, L)
    else
        [iseven(i) ? 1.0 : -1.0 for i=0:L-1]
    end

    TranslationParityBasis(zeros(Int, L), I, R, C, P, base)
end

TranslationParityBasis(K, P, L; base=2, alloc=1000, threaded=false) = TranslationParityBasis(x -> true, K, P, L, base=base, alloc=alloc, threaded=threaded)

function TranslationFlipBasis(
    f, K::Integer, P::Integer, L::Integer;
    base::Integer=2, alloc::Integer=1000, threaded::Bool=false
)
    @assert isone(P) || isequal(P, -1) "Invalid parity"
    judge = begin
        N2 = [2L / sqrt(i) for i = 1:L]
        N1 = N2 ./ sqrt(2)
        TranslationParityJudge(f, x -> spinflip(x, base), K, P, L, base, vcat(N1, N2))
    end

    I, R = if threaded
        selectindexnorm_threaded(judge, L, base=base, alloc=alloc)
    else
        selectindexnorm(judge, L, 1:base^L, base=base, alloc=alloc)
    end

    C = if iszero(K)
        fill(1.0, L)
    elseif isequal(2K, L)
        [iseven(i) ? 1.0 : -1.0 for i=0:L-1]
    else
        [exp(-1im * 2π/L * K * i) for i=0:L-1]
    end

    TranslationFlipBasis(zeros(Int64, L), I, R, C, P, base)
end

TranslationFlipBasis(K, P, L; base=2, alloc=1000, threaded=false) = TranslationFlipBasis(x -> true, K, P, L, base=base, alloc=alloc, threaded=threaded)

function parity_schmidt(parity, v::AbstractVector, Ainds::AbstractVector{<:Integer}, b::AbstractTranslationalParityBasis)
    dgt, R, phase = b.dgt, b.R, b.C[2]
    S = SchmidtMatrix(promote_type(eltype(v), eltype(b)), b, Ainds)
    for i = 1:length(v)
        change!(b, i)
        val = v[i] / R[i]
        for j in 1:length(dgt)
            addto!(S, val)
            cyclebits!(dgt)
            val *= phase

        end
        dgt .= parity(dgt)
        val *= b.P
        for j in 1:length(dgt)
            addto!(S, val)
            cyclebits!(dgt)
            val *= phase
        end
    end
    Matrix(S)
end

schmidt(v, Ainds, b::TranslationParityBasis) = parity_schmidt(reverse, v, Ainds, b)
schmidt(v, Ainds, b::TranslationFlipBasis) = parity_schmidt(x -> spinflip(x, b.B), v, Ainds, b)

#-----------------------------------------------------------------------------------------------------
# Double Basis
#-----------------------------------------------------------------------------------------------------
export DoubleBasis
"""
    DoubleBasis{Tb1<:AbstractBasis, Tb2<:AbstractBasis}

Basis for constructing transition matrix from one symmetry sector to another.

Properties:
-----------
- `dgt`: Digits.
- `B1` : Basis of the target symmetry sector.
- `B2` : Basis of the starting symmetry sector.
- `B`  : Base.
"""
struct DoubleBasis{Tb1<:AbstractBasis, Tb2<:AbstractBasis} <: AbstractBasis
    dgt::Vector{Int64}
    B1::Tb1
    B2::Tb2
    B::Int64
    DoubleBasis(B1::AbstractBasis, B2::AbstractBasis) = new{typeof(B1), typeof(B2)}(B1.dgt, B1, B2, B2.B)
end

@inline eltype(b::DoubleBasis) = promote_type(eltype(b.B1), eltype(b.B2))
@inline function change!(b::DoubleBasis, i::Integer)
    N = change!(b.B2, i)
    b.dgt .= b.B2.dgt
    N
end
@inline index(b::DoubleBasis) = index(b.B1)
@inline size(b::DoubleBasis) = size(b.B1, 1), size(b.B2, 2)
@inline size(b::DoubleBasis, i::Integer) = isone(i) ? size(b.B1, 1) : isequal(i, 2) ? size(b.B2, 2) : 1

ext: jl


text:
struct Trace{T,S,F} <: Message
    start::T
    X::S
    stop::F
end

function collecttrace(iter, t::Trace)
    start = t.start
    X = t.X
    stop = t.stop
    ϕ = dyniterate(iter, start)
    ϕ === nothing && return X
    x, u = ϕ
    push!(X, x)
    while !stop(x)
        ϕ = dyniterate(iter, u)
        ϕ === nothing && break
        x, u = ϕ
        push!(X, x)
    end
    return X
end

export Trace, collecttrace


text:
using PkgTemplates

t = Template(;
           authors="Robert Rudolph",
           julia=v"1.6.2",
	   dir="C:\\Users\\RR\\julia",
	       plugins=[
               TravisCI(),
               Codecov(),
               AppVeyor(),
               GitHubPages(),
               CirrusCI(),
           ],
       )

print("enter new package name")
pkgName = readline()

generate(t, pkgName)


text:
# run_accuracy.jl

using SPECTrecon

include("project.jl")


text:
using Base.Test
using DataFrames

let
	v = [1, 2, 3, 4]
	dv = DataArray(v, falses(size(v)))

	m = [1 2; 3 4]
	dm = DataArray(m, falses(size(m)))

	t = Array(Int, 2, 2, 2)
	t[1:2, 1:2, 1:2] = 1
	dt = DataArray(t, falses(size(t)))

	dv = DataArray(v)
	dv = DataArray(v, [false, false, false, false])

	dv = DataArray(dv)
	dm = DataArray(dm)
	dt = DataArray(dt)

	dv = DataArray(Int, 2)
	dm = DataArray(Int, 2, 2)
	dt = DataArray(Int, 2, 2, 2)

	similar(dv)
	similar(dm)
	similar(dt)

	similar(dv, 2)
	similar(dm, 2, 2)
	similar(dt, 2, 2, 2)

	# Test vecbind
	a = [1:4]
	d = DataArray(a)
	@assert isequal(vecbind(a,a), [a,a])
	@assert isequal(vecbind(a,1.0 * a), 1.0 * [a,a])
	@assert isequal(vecbind(d,a), [d,d])
	@assert isequal(vecbind(a,d), [d,d])
	@assert isequal(vecbind(1.0 * a,d), 1.0 * [d,d])
	@assert isequal(vecbind(d,d), [d,d])
	@assert isequal(vecbind(a,IndexedVector(a)), [a,a])
	@assert isequal(vecbind(a,IndexedVector(d)), [d,d])
	@assert isequal(vecbind(PooledDataArray(a),IndexedVector(d)), [d,d])
	@assert isequal(vecbind(PooledDataArray(a),IndexedVector(a)), [a,a])
	@assert isequal(vecbind(a,RepeatedVector(a,2)), [a,a,a])
	@assert (vecbind(a,StackedVector({a,1.0*a})) == [a,a,a])
end


text:
@testset "one_component_2d" begin

	for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
		for i = 1:20
			test_image = T.(eval(Symbol("test_image_$(i)")))

			# Call one component 2D algorithm and compare the number of components.
			labels = ImageComponentAnalysis.label_components(OneComponent2D(), test_image)
			num_components = maximum(labels)

			@test num_components == test_image_results[i]
		end
	end

end


text:
# Test for load-order-specific issues with WebIO and IJulia
# when running outside of Jupyter.
#
# See: https://github.com/JuliaGizmos/WebIO.jl/issues/117
using Test
using WebIO
using Conda, IJulia

@testset "IJulia initialization" begin
    @test :ijulia in WebIO.providers_initialised
end

if haskey(ENV, "CI")
    # This takes a long time so only run in CI by default.
    @testset "Jupyter setup" begin
        # We need conda-forge for node.
        Conda.add_channel("conda-forge")
        Conda.add.(["jupyterlab>=2,<3", "nodejs>=10"])

        # Remove path to make sure that we can find Conda if nothing else exists.
        oldpath = ENV["PATH"]
        ENV["PATH"] = ""
        conda_jupyter = WebIO.find_jupyter_cmd()

        # Restore the old PATH because we need to have node.
        pathsep = Sys.iswindows() ? ";" : ":"
        ENV["PATH"] = string(Conda.BINDIR, pathsep, oldpath)
        @assert Sys.which("node") !== nothing

        @testset "JupyterLab setup" begin
            # Uninstall before actual tests begin
            run(Cmd(
                `$conda_jupyter labextension uninstall @webio/jupyter-lab-provider`,
                ignorestatus=true,
            ))
            WebIO.install_jupyter_labextension(conda_jupyter)
            errpipe = Pipe()
            run(pipeline(`$conda_jupyter labextension list`, stderr=errpipe))
            close(errpipe.in)
            installed_extensions = read(errpipe, String)
            @test occursin("@webio/jupyter-lab-provider", installed_extensions)
        end

        @testset "Jupyter Notebook setup" begin
            WebIO.install_jupyter_nbextension()
            installed_extensions = read(`$conda_jupyter nbextension list`, String)
            @test occursin("webio", installed_extensions)
        end

        @testset "Jupyter serverextension setup" begin
            # Note: The server extension should be automatically installed after
            # installing either the lab or notebook plugins.
            installed_extensions = read(`$conda_jupyter serverextension list`, String)
            @test occursin("jlstaticserve", installed_extensions)
        end

        # Restore path for good hygiene
        ENV["PATH"] = oldpath
    end
else
    @info(
        "Skipping Conda/Jupyter installation tests; set the CI environment "
        * "variable to run these tests as well.",
    )
end

# Test the demo jupyter notebook and make sure all of its cells execute
# without error.
using NBInclude

@testset "Demo notebook" begin
    @nbinclude(joinpath(@__DIR__, "assets", "jupyter-test.ipynb"))
end


text:
tile_size = 50
board_size = 4 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]

chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
img = chessboard

###
using Random, PyPlot
using ScatteredInterpolation
using Dierckx

# params:
point_count = 100

flow = gen_tiled_flow(size(img), 30, 40)

rand_inds = rand(CartesianIndices(flow), point_count)

flow[rand_inds[1]]

## ScatteredInterpolation:
using BenchmarkTools, ScatteredInterpolation

points = reshape([ind[i] for i=1:length(rand_inds[1]) for ind=rand_inds], length(rand_inds), length(rand_inds[1]))'
samples = [flow[I] for I in rand_inds]

idx = rand_inds
@benchmark reduce((x,y)->hcat(x, [y.I...]), idx[2:end]; init=[idx[1].I...])

@benchmark reshape([ind[i] for i=1:length(rand_inds[1]) for ind=rand_inds], length(rand_inds), length(rand_inds[1]))'

# fun
flow_size = size(flow)
x1 = 1:flow_size[1]
y1 = 1:flow_size[2]

X = repeat(x1, flow_size[1])[:]
Y = repeat(y1', flow_size[1])[:]
gridPoints = [X Y]'

itp = interpolate(method, points, samples);
interpolated = ScatteredInterpolation.evaluate(itp, gridPoints)
gridded = reshape(interpolated, flow_size)

# fun end

# methods: RBFs = [Gaussian(2), Multiquadratic(2), InverseQuadratic(2), InverseMultiquadratic(2)]
function test_scatter(points, samples, smoothness; method)
    flow_size = size(flow)
    x1 = 1:flow_size[1]
    y1 = 1:flow_size[2]

    X = repeat(x1, flow_size[1])[:]
    Y = repeat(y1', flow_size[1])[:]
    gridPoints = [X Y]'

    itp = interpolate(method, points, samples);
    interpolated = ScatteredInterpolation.evaluate(itp, gridPoints)
    gridded = reshape(interpolated, flow_size)
    return gridded
end

@btime gridded_scatter = test_scatter($points, $samples, $0)
@benchmark gridded_scatter = test_scatter($points, $samples, $0, method=Shepard())
@benchmark gridded_scatter = test_scatter($points, $samples, $0)
gridded_scatter = test_scatter(points, samples, 0, method=Multiquadratic(2))
gridded_scatter = test_scatter(points, samples, 0, method=Shepard())

showflow(gridded_scatter, figtitle="scatter")

showflow(flow, figtitle="orig")

import LAP_julia: helpers.mse, helpers.angle_rmse, helpers.angle_mae, helpers.vec_len

mse(gridded_scatter, flow)
vec_len(mse(gridded_scatter, flow))
angle_rmse(gridded_scatter, flow)
angle_mae(gridded_scatter, flow)

sum(flow - gridded_scatter)/length(flow)

showflow(flow)
x = [ind[1] for ind in rand_inds]
y = [ind[2] for ind in rand_inds]
scatter(x, y, marker=:x)
gcf()

## Dierckx
using Dierckx

x = float([ind[1] for ind in rand_inds])
y = float([ind[2] for ind in rand_inds])
z1 = [real(flow[ind]) for ind in rand_inds]
z2 = [imag(flow[ind]) for ind in rand_inds]

function test_dier(x, y, z1, z2)
    storage = 1000
    spl1 = Spline2D(x, y, z1; w=ones(length(x)), kx=1, ky=1, s=storage)
    real1 = evalgrid(spl1, 1:size(flow)[1], 1:size(flow)[2])

    spl2 = Spline2D(x, y, z2; w=ones(length(x)), kx=1, ky=1, s=storage)
    imag1 = evalgrid(spl2, 1:size(flow)[1], 1:size(flow)[2])

    gridded = real1 .+ imag1 .* im
    return gridded
end

@benchmark gridded_dier = test_dier($x, $y, $z1, $z2)
gridded_dier = test_dier(x, y, z1, z2)

showflow(gridded_dier, figtitle="dier")
showflow(flow)

### IrregularInterpolate

# using PyPlot
#
# f(x, y) = (x+y-1)^2
#
# x = float([ind[1] for ind in rand_inds]) ./ 200
# y = float([ind[2] for ind in rand_inds]) ./ 200
# a = f.(x, y)
#
# p = Main.IrregularInterpolate.PolyharmonicSpline(2, [x y], a)  # 3 is the order
#
# inds = CartesianIndices(flow)
# x2 = float(vec([I[1] for I in inds])) ./ 200
# y2 = float(vec([I[2] for I in inds])) ./ 200
# data = Main.IrregularInterpolate.interpolate(p, [x2 y2])
#
# gridded = reshape(data, size(flow))
#
# imgshow(gridded)
#
#
# figure()
# plot_surface([I[1] for I in inds], [I[2] for I in inds], gridded)
# gcf()

###

using PyPlot

x = float([ind[1] for ind in rand_inds])
y = float([ind[2] for ind in rand_inds])
z1 = [real(flow[ind]) for ind in rand_inds]
z2 = [imag(flow[ind]) for ind in rand_inds]

function test_code(x, y, z1, z2)
    p1 = Main.IrregularInterpolate.PolyharmonicSpline(2, [x y], z1)  # 3 is the order
    p2 = Main.IrregularInterpolate.PolyharmonicSpline(2, [x y], z2)  # 3 is the order

    inds = CartesianIndices(flow)
    x2 = float(vec([I[1] for I in inds]))
    y2 = float(vec([I[2] for I in inds]))
    data_real = Main.IrregularInterpolate.interpolate(p1, [x2 y2])
    data_imag = Main.IrregularInterpolate.interpolate(p2, [x2 y2])

    gridded = reshape((data_real .+ data_imag .* im), size(flow))
    return gridded
end

@benchmark gridded_code = test_code($x, $y, $z1, $z2)
gridded_code = test_code(x, y, z1, z2)

showflow(gridded_code)

x = [ind[1] for ind in rand_inds]
y = [ind[2] for ind in rand_inds]
scatter(x, y, marker=:x)
gcf()

showflow(flow)

figure()
plot_surface([I[1] for I in inds], [I[2] for I in inds], imag(gridded), alpha=0.5)
gcf()

z1 = [real(point) for point in z]
z2 = [imag(point) for point in z]

### compare results

# prepare data:
flow = gen_tiled_flow(size(img), 30, 60)

# set params:
point_count = 10
rand_inds = rand(CartesianIndices(flow), point_count) # random points

# prepare inputs
points = reshape([ind[i] for i=1:length(rand_inds[1]) for ind=rand_inds], length(rand_inds), length(rand_inds[1]))'
samples = [flow[I] for I in rand_inds]
x = float([ind[1] for ind in rand_inds])
y = float([ind[2] for ind in rand_inds])
z1 = [real(flow[ind]) for ind in rand_inds]
z2 = [imag(flow[ind]) for ind in rand_inds]

# scatter
@benchmark gridded_scatter = test_scatter($points, $samples, $0)
gridded_scatter = test_scatter(points, samples, 0)
# Dierckx
@benchmark gridded_dier = test_dier($x, $y, $z1, $z2)
gridded_dier = test_dier(x, y, z1, z2)
# code
b = @benchmark gridded_code = test_code($x, $y, $z1, $z2)
gridded_code = test_code(x, y, z1, z2)

# show flows:
showflow(flow, figtitle="orig")
showflow(gridded_scatter, figtitle="scatter")
showflow(gridded_dier, figtitle="dier")
showflow(gridded_code, figtitle="code")

# add points
s = scatter(x, y, marker=:x)
gcf()

###
using DataFrames

Methods = ["Dier",
           "Scatter - Thin plate spline",
           "Scatter - Multiquadratic"]

Modes = ["Random_points",
         "Semigridded",
         "Gridded"]

Point_counts = [9,
                16,
                25,
                49,
                81,
                100,
                144]

###
# gridded
test_flow = CartesianIndices(ones(10,10))

flow_size = size(test_flow)
point_count = 10
points_in_dim = round(Int64, sqrt(point_count))
point_spacing = floor.(flow_size ./ points_in_dim)

function gen_rand_points(flow, point_count, mode; max_skew_const=4)
    if mode == "Random_points"
        ret_inds = rand(CartesianIndices(flow), point_count)
    elseif mode == "Gridded"
        points_in_dim = round(Int64, sqrt(point_count))
        point_spacing = floor.(size(flow) ./ points_in_dim)

        ret_inds = Array{CartesianIndex{2}, 1}(undef, 0)
        for ind in CartesianIndices(flow)
            if all((Tuple(ind)) .% point_spacing .== 0)
                add_ind = CartesianIndex(Tuple(ind) .- floor.(Int64, point_spacing ./ 2))
                push!(ret_inds, add_ind)
            end
        end
    elseif mode == "Semigridded"
        points_in_dim = round(Int64, sqrt(point_count))
        point_spacing = floor.(size(flow) ./ points_in_dim)

        ret_inds = Array{CartesianIndex{2}, 1}(undef, 0)
        print(size(ret_inds))
        for ind in CartesianIndices(flow)
            if all((Tuple(ind)) .% point_spacing .== 0)
                skew = (rand((-floor(Int64, point_spacing[1]/max_skew_const):floor(Int64, point_spacing[1]/max_skew_const))),
                        rand((-floor(Int64, point_spacing[1]/max_skew_const):floor(Int64, point_spacing[2]/max_skew_const))))
                add_ind = CartesianIndex(Tuple(ind) .- floor.(Int64, point_spacing ./ 2) .+ skew)
                @assert add_ind in CartesianIndices(flow)

                push!(ret_inds, add_ind)
            end
        end
    end
    return ret_inds
end

points = gen_rand_points(ones(100, 100), 40, "Semigridded")

imgshow(ones(100, 100))
x = [ind[1] for ind in points]
y = [ind[2] for ind in points]
scatter(x, y, marker=:x)
gcf()

# make Semigridded

###

scatter_interp_bench = DataFrame(Methods = )

df = DataFrame(A = 1:2, a = ["ss", "sdfd"], B = [bench, bench], B2 = [showflow(flow, ret="PyObject"), showflow(flow, ret="PyObject")], C =[ 1, 4])

display(df)

showflow(flow, figtitle="orig", ret="PyObject")

import PyCall.PyObject

scatter_interp_bench = DataFrame(
                                Index = Int[],
                                Method = String[],
                                Mode = String[],
                                Point_count = Int[],
                                Median_speed = Float64[],
                                Magnitude_of_MSE = Float64[],
                                Angle_RMS = Float64[],
                                Angle_mean = Float64[],
                                Benchmark = BenchmarkTools.Trial[],
                                Truth_flow = PyObject[],
                                Method_flow = PyObject[],
                                Scatter_points = PyObject[]
                                )

df = scatter_interp_bench
df = DataFrame(A = Int[], B = String[])

push!(df, (1, "shit"))

showtable(df)

insert!(df, 3, undef, :aaArea)

df.aaArea = [s, s]

# test
A = rand(50000, 30)
d = DataFrame(A)

tmp = d[d.x1 .> 0.7, :]
tmp === d[d.x1 .> 0.7, :]

arr = [12.3545235, 12.323, 2.123412]

mat = [12; 12]

round.(arr, digits=2)

vcat([3, 4, 5], arr)

## test try catch
arr = []

for k in 1:12
    try
        tmp = 12
        if k == 4
            @assert true == false
        end
    catch e
        push!(arr, (k, e))
        continue;
    end
    println(tmp)
    println(k)
end

try
    @assert true == false
catch e
    @warn e
    typeof(e)
end

function f(x::Type, y) where Type <: Real
    x .+ y
end

function g2(x::Type, y::mat) where {Type <:Real, mat <: Matrix{<:Complex{<:Real}}}
    x .+ y
end

@code_warntype f(12, fl)

g2(12, 12)

fl = gen_tiled_flow((12,12), 12)

typeof(fl)

point_count_df = DataFrame(Point_c = Int[],
                           Method = String[],
                           Median_speed_mean = Float64[],
                           Magnitude_of_MSE_mean = Float64[],
                           Angle_RMS_mean = Float64[],
                           Angle_Mean_mean = Float64[])

for points in Point_counts
    for method in Methods
        tmp_df = df[((df.Point_count == points) .& (df.Method .== method)), :]
        tmp_df = select(tmp_df, :Median_speed, :Magnitude_of_MSE, :Angle_RMS, :Angle_mean)
        test_count = length(tmp_df.Median_speed)
        tmp_df = describe(tmp_df)
        retrieved_data = select(tmp_df, :mean, copycols=false)
        push!(poin_count_df, [points method retrieved_data])
    end
end

bench = @benchmark 12+1

df = DataFrame(a = Matrix[])

push!(df, Dict(:a => flow))

df[:, :a][1]

typeof(flow)

df = DataFrame()

flow = gen_tiled_flow((100,100), 12)
plot = showflow(flow)

@save "save_df.JLD2" df flow


text:
using App
App.julia_main()


text:
#
# Global state
#

type COMGlobal
    initialized::Bool
    COMGlobal() = new(false)
end

################################################################################
# @vcall macro for virtual function calls
#   usage: @vcall(object, offset, arguments...)
#
#       object::Ptr{Void}
#       offset: vtable offset
#       arguments: arguments in Arg::Type form
################################################################################

macro vcall(this, idx, rtype, args...)
    argnames = {x.args[1] for x in args}
    argtypes = :((Ptr{Void}, $({eval(x.args[2]) for x in args}...)))
    quote
        local fptr = unsafe_load(unsafe_load(pointer(Ptr{Ptr{Void}},$(this))),$(idx)+1)
        ccall(fptr, thiscall, $(esc(rtype)), ($(esc(argtypes))), $(esc(this)),$(argnames...))
    end
end

################################################################################
# Extras
################################################################################

get_default_lcid() = ccall( (:GetSystemDefaultLCID, l_kernel32), LCID, () )

getvtptr(o,idx) = unsafe_load(unsafe_load(pointer(Ptr{Ptr{Void}},o)),idx)

#
# L string helper
#
macro L_str(s)
    quote utf16($s) end
end

#
# Display and conversion
#

show(io::IO, x::GUID) = @printf(io, "%08X-%04hX-%04hX-%02hhX%02hhX%02hhX%02hhX%02hhX%02hhX%02hhX%02hhX",
                     x.data1, x.data2, x.data3,
                     x.data4.d1,x.data4.d2,x.data4.d3,x.data4.d4,x.data4.d5,x.data4.d6,
                     x.data4.d7,x.data4.d8)

#show(io::IO, x::CLSID) = print(io, "CLSID{", x.guid, "}")

convert(T::Type{Ptr{CLSID}}, x::CLSID) = reinterpret(Ptr{CLSID}, pointer(x.guid))
convert(::Type{Ptr{GUID}}, x::CLSID) = pointer(x.guid)
convert(::Type{Ptr{IID}}, x::IID) = pointer(x.guid)
convert(::Type{IID}, x::CLSID) = IID(x.guid)
convert(::Type{GUID}, x::CLSID) = x.guid

convert(::Type{Ptr{Void}}, x::IUnknown) = x.ptr


text:
# Generate center coordinates for specific index
"""
    st_coords(ga; mid)

`st_coords` is for regular AbstractGeoArray.

# Return
- `X,Y`: A matrix with the dimension of [NX, NY]
"""
function st_coords(ga::AbstractGeoArray; mid::Vector{Int} = [1, 1])
    x, y = st_dim(ga; mid)
    meshgrid(x, y)
end

# This script is modified from:
#   `https://github.com/evetion/AbstractGeoArrays.jl/blob/master/src/AbstractGeoArray.jl`
# Copyright (c) 2018 Maarten Pronk, MIT license

"""
    st_coordinates(ga)

`st_coordinates` is similar as `st_dim`, but for irregular AbstractGeoArray.
"""
function st_coordinates(ga::AbstractGeoArray, p::SVector{2, Int})
    ga.f(p .- 0.5)
end
st_coordinates(ga::AbstractGeoArray, p::Vector{Int}) = st_coordinates(ga, SVector{2}(p))

function st_coordinates(ga::AbstractGeoArray)
    # (ui, uj) = size(ga)[1:2]
    # ci = [st_coordinates(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj]
    nrow, ncol = size(ga)[1:2]
    X = zeros(nrow, ncol)
    Y = zeros(nrow, ncol)
    @inbounds for i = 1:nrow, j = 1:ncol
        ci = st_coordinates(ga, SVector{2}(i,j))
        X[i, j] = ci[1]
        Y[i, j] = ci[2]
    end
    X, Y
end

function st_coordinates(ga::AbstractGeoArray, dim::Symbol)
    if is_rotated(ga)
        error("This method cannot be used for a rotated AbstractGeoArray")
    end
    if dim==:x
        ui = size(ga,1)
        ci = [st_coordinates(ga, SVector{2}(i,1))[1] for i in 1:ui]
    elseif dim==:y
        uj = size(ga,2)
        ci = [st_coordinates(ga, SVector{2}(1,j))[2] for j in 1:uj]
    else
        error("Use :x or :y as second argument")
    end
    return ci
end

"""
    st_coordinates!(ga, x::AbstractUnitRange, y::AbstractUnitRange)

Set AffineMap of `AbstractGeoArray` by specifying the *center coordinates* for each x, y
dimension by a `UnitRange`.
"""
function st_coordinates!(ga, x::AbstractUnitRange, y::AbstractUnitRange)
    size(ga)[1:2] != (length(x), length(y)) && error("Size of `AbstractGeoArray` $(size(ga)) does not match size of (x,y): $((length(x),length(y))). Note that this function takes *center coordinates*.")
    ga.f = unitrange_to_affine(x, y)
    ga
end

export st_coords, st_coordinates, st_coordinates!


text:
# Test that the macro call `m` throws an error exception during pre-compilation
macro test_macro_throws(errortype, m)
    # See https://discourse.julialang.org/t/test-throws-with-macros-after-pr-23533/5878
    :(@test_throws $errortype try @eval $m catch err; throw(err.error) end)
end


text:
"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct AbstractParametersVPRK{IT, DT, TT, D, S, ET <: NamedTuple, PT <: NamedTuple} <: Parameters{DT,TT}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    pparams::PT

    t̄::TT
    q̄::Vector{DT}
    p̄::Vector{DT}

    function AbstractParametersVPRK{IT,DT,D}(equs::ET, tab::TableauVPRK{TT}, Δt::TT, pparams::PT=NamedTuple()) where {IT, DT, TT, D, ET <: NamedTuple, PT <: NamedTuple}
        new{IT, DT, TT, D, tab.s, ET, PT}(equs, tab, Δt, pparams, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

function update_params!(params::AbstractParametersVPRK, sol::Union{AtomicSolutionPODE, AtomicSolutionPDAE})
    # set time for nonlinear solver and copy previous solution
    params.t̄  = sol.t
    params.q̄ .= sol.q
    params.p̄ .= sol.p
end

@inline equation(int::AbstractIntegratorVPRK, i::Symbol) = parameters(int).equ[i]
@inline equations(int::AbstractIntegratorVPRK) = parameters(int).equ
@inline timestep(int::AbstractIntegratorVPRK) = parameters(int).Δt
@inline tableau(int::AbstractIntegratorVPRK) = parameters(int).tab

function Integrators.IntegratorCache{ST}(params::AbstractParametersVPRK{IT,DT,TT,D,S}; kwargs...) where {IT,ST,DT,TT,D,S}
    IntegratorCacheVPRK{ST,D,S}(true; kwargs...)
end

@inline Integrators.CacheType(ST, params::AbstractParametersVPRK{IT,DT,TT,D,S}) where {IT,DT,TT,D,S} = IntegratorCacheVPRK{ST,D,S}


text:
using Distributions, Plots; pyplot()

alpha, beta = 8, 2
prior(lam) = pdf(Gamma(alpha, 1/beta), lam)
data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

like(lam) = *([pdf(Poisson(lam),x) for x in data]...)
posteriorUpToK(lam) = like(lam)*prior(lam)

sig = 0.5
foldedNormalPDF(x,mu) = (1/sqrt(2*pi*sig^2))*(exp(-(x-mu)^2/2sig^2)
                                                + exp(-(x+mu)^2/2sig^2))
foldedNormalRV(mu) = abs(rand(Normal(mu,sig)))

function sampler(piProb,qProp,rvProp)
    lam = 1
    warmN, N = 10^5, 10^6
    samples = zeros(N-warmN)

    for t in 1:N
        while true
            lamTry = rvProp(lam)
            L = piProb(lamTry)/piProb(lam)
            H = min(1,L*qProp(lam,lamTry)/qProp(lamTry,lam))
            if rand() < H
                lam = lamTry
                if t > warmN
                    samples[t-warmN] = lam
                end
                break
            end
        end
    end
    return samples
end

mcmcSamples = sampler(posteriorUpToK,foldedNormalPDF,foldedNormalRV)
println("MCMC Bayes Estimate: ",mean(mcmcSamples))

stephist(mcmcSamples, bins=100,
	c=:black, normed=true, label="Histogram of MCMC samples")

lamRange = 0:0.01:10
plot!(lamRange, prior.(lamRange),
	c=:blue, label="Prior distribution")

closedFormPosterior(lam)=pdf(Gamma(alpha + sum(data),1/(beta+length(data))),lam)
plot!(lamRange, closedFormPosterior.(lamRange),
	c=:red, label="Posterior distribution",
	xlims=(0, 10), ylims=(0, 1.2),
    xlabel=L"\lambda",ylabel="Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#Prezentacja wyników do sprawozdania: Zadanie P1.2
#Kacper Kulczak 279079

#Pomocnicza funkcja pomgająca. Potegowanie
function pow(n,w)
    result = one(n)
    for i = one(w):w
        result *= n
    end
    return result
end

# Funkcja obliczajaca stala e za pomoca przyblizenia granicy (1+1/n)^n
function definition_e(n)
    result = (one(n) + one(n) / n)
    for i = one(n):log2(n)
        result *= result
    end
    return result
end
# Funkcja obliczajaca stala e za pomoca szeregu taylora sum(1/ n!)
function taylor_series(n)
    result = zero(n)
    for i = zero(n):n
        result += one(n) / factorial(n-i)
    end
    return result
end
# Funkcja obliczajaca stala e za pomoca ulamka lancuchowego
function continued_fraction(n)
    i = n + n
    result = i
    result = one(n) + one(n) / result
    result = one(n) + one(n) / result
    i -= one(n) + one(n)
    while i > one(n) + one(n)
        result = i + one(n) / result
        result = one(n) + one(n) / result
        result = one(n) + one(n) / result
        i -= one(n) + one(n)
    end
    result = one(n) + one(n) + one(n) / result
    result = one(n) + one(n) / result
    result = one(n) + one(n) + one(n) / result
    return result

end
#Drugi szereg SUM (3k+1)/(3k)!
function second_series(n)
    result = zero(n)
    for i = zero(n):n
        result += (one(n) + (convert(typeof(n),3) * (n-i) ) * ( convert(typeof(n),3) * (n-i) ) )/ factorial(convert(typeof(n),3) *(n-i))
    end
    return result
end

# Trzeci szreg SUM (4k + 3) / 2^(2k+1) * (2k+1)!
function third_series(n)
    result = zero(n)
    for i = zero(n):n
        temp = convert(typeof(n),4) * (n-i) + one(n) + one(n) + one(n)
        temp /= factorial(convert(typeof(n),2) *(n-i) + one(n)) * pow( convert(typeof(n),2), convert(typeof(n),2)*(n-i) + one(n))
        result += temp
    end
    return result * result
end
#Test zapisujace wyniki na zadanym zasiegu 1->range
function test(fun,range)
    container= [[],[],[],[]]
    for i = one(range) : range
        temp = fun(i)
        err = convert(typeof(range),eu)-temp
        push!(container[1], i)
        push!(container[2], temp)
        push!(container[3],err)
        push!(container[4], -log10(abs( err )) )
    end
    return container
end

#Test dla definicji, przedstawiajacy tylko niektore iteracje
function test_def(fun,range)
    container= [[],[],[],[]]
    x = one(range)
    for i = one(range) : range
        x *= convert(typeof(range),8.)
        temp = fun(x)
        err = convert(typeof(range),eu)-temp
        push!(container[1], x)
        push!(container[2], temp)
        push!(container[3],err)
        push!(container[4], -log10(abs( err )) )
    end
    return container
end

#Funkcja formatujaca wyjscie
function myshow(list)
    for i = 1:length(list[1])
        @printf "n=%d val=%.f error=%.4e #cyfr=%.0f\n" list[1][i] list[2][i] list[3][i] list[4][i]
    end
end

setprecision(256)
@show precision(BigFloat)
range = BigFloat(60.)
data1 = test_def(definition_e, BigFloat(30.) )
data2 = test(taylor_series, range)
data3 = test(continued_fraction, range )
data4 = test(second_series, range)
data5 = test(third_series, range)
@printf "\n\n\n"
println("Definition of e")
myshow(data1)
@printf "\n\n\n"
println("Taylor series")
myshow(data2)
@printf "\n\n\n"
println("Contiuned fraction")
myshow(data3)
@printf "\n\n\n"
println("Second Series")
myshow(data4)
@printf "\n\n\n"
println("Third Series")
myshow(data5)

#Funckje do obliczenia ilości mnożeń i dzieleń w zależnosci od iteracji w wykorzytanych
# przeze mnie metodach
function Tdef(n)
    return one(n)+ one(n)+log2(n)
end
function TCF(n)
    return n * convert(typeof(n),3)
end
function Tser1(n)
    return n * convert(typeof(n),2)
end
function Tser2(n)
    return n * convert(typeof(n),4)
end
function Tser3(n)
    return n * convert(typeof(n),5) + one(n)
end


text:
#utilities for QUtilities
#some of these are candidates for Finometrics
const CRSP_MISSINGS = [
  "A", "B", "C", "D", "E", "S", "T", "P",
  -55.0, -66.0, -77.0, -88.0, -99.0,
  -55, -66, -77, -88, -99,
  "-55.0", "-66.0", "-77.0", "-88.0", "-99.0"]

const CRSP_PARSED_MISSINGS = (-55.0, -66.0, -77.0, -88.0, -99.0)

const FRET = :ret
const CQUART_FVALS = [:mknegcash, :bknegcash, :mkflev, :bkflev, :mkliab, :bkliab]
const CQUART_FCONTROLS = :LLret12m
const CQUART_FWEIGHTONS = nothing
const CQUART_FNAMES = (s::Symbol->Symbol("P_",s)).(CQUART_FVALS)
const CQUART_FGROUPS = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES)

const CQUART_FCONTROLS_ALT = :Lmkat
const CQUART_FNAMES_ALT = (s::Symbol->Symbol(s,"_ac")).(CQUART_FNAMES)
const CQUART_FGROUPS_ALT = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_ALT)

const CQUART_FVALS_D = [:Dmknegcash, :Dbknegcash, :Dmkflev, :Dbkflev, :Dmkliab, :Dbkliab]
const CQUART_FCONTROLS_D = :LLret12m
const CQUART_FWEIGHTONS_D = nothing
const CQUART_FNAMES_D = (s::Symbol->Symbol("P_",s)).(CQUART_FVALS_D)
const CQUART_FGROUPS_D = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_D)

const CQUART_FCONTROLS_ALT_D = :Lmkat
const CQUART_FNAMES_ALT_D = (s::Symbol->Symbol(s,"_ac")).(CQUART_FNAMES_D)
const CQUART_FGROUPS_ALT_D = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_ALT_D)

const CQUART_FVALS_RMD = [:RDmkflev, :MDmkflev, :RDmkliab, :MDmkliab]
const CQUART_FCONTROLS_RMD = :LLret12m
const CQUART_FWEIGHTONS_RMD = nothing
const CQUART_FNAMES_RMD = (s::Symbol->Symbol("P_",s)).(CQUART_FVALS_RMD)
const CQUART_FGROUPS_RMD = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_RMD)

const CQUART_FCONTROLS_ALT_RMD = :Lmkat
const CQUART_FNAMES_ALT_RMD = (s::Symbol->Symbol(s,"_ac")).(CQUART_FNAMES_RMD)
const CQUART_FGROUPS_ALT_RMD = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_ALT_RMD)

const CQUART_FNAMES_FD = (s::Symbol->Symbol("FP_",s)).(CQUART_FVALS_D)
const CQUART_FGROUPS_FD = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_FD)

const USED_COMPA_VALUE_FIELDS = [:che, :ch, :at, :dltt, :dlc, :ceq, :txditc,
  :txdb, :pstk, :lt]
const USED_COMPQ_VALUE_FIELDS = [:rdq]
#const USED_COMP_VALUE_FIELDS = [USED_COMPA_VALUE_FIELDS; USED_COMPQ_VALUE_FIELDS;]
const USED_CRSP_FIELDS = [:ret, :price, :shrout]

#final columns for the decompressed file
const RETAINED_COLS_CCM = [:gvkey, :lpermno, :linkeffdate, :linkenddate, :cusip]

const RETAINED_COLS_COMP = [:gvkey; :fyear; :fyrmonth;
  :linkeffdate; :linkenddate; :Nadate;
  :ddate; :adate; :lpermno; :begindate; :enddate; USED_COMPA_VALUE_FIELDS;
  :bknegcash; :cash; :netdebt; :bkequity; :bkflev; :bkliab; :mkequity; :nop;]

#not sure if we really need much of these
const RETAINED_COLS_UNIV = [RETAINED_COLS_COMP;
  :compid; :permno; :date; :ret; :mktcap; :net]

#Only keep what we need, also serves as a check to make sure we ahve what we think
const RETAINED_COLS_DATA_SERIES = unique([:gvkey; :fyear; :fyrmonth;
  :linkeffdate; :linkenddate; :Nadate;
  :ddate; :adate; :lpermno; :begindate; :enddate; :mkequity; :Lmkequity; :mkat;:lmkequity; :Llmkequity;
  CQUART_FCONTROLS; CQUART_FCONTROLS_ALT;
  CQUART_FVALS; (s->Symbol(:G_,s)).(CQUART_FNAMES); (s->Symbol(:G_,s)).(CQUART_FNAMES_ALT);
  CQUART_FCONTROLS_D; CQUART_FCONTROLS_ALT_D;
  CQUART_FVALS_D; (s->Symbol(:G_,s)).(CQUART_FNAMES_D); (s->Symbol(:G_,s)).(CQUART_FNAMES_ALT_D);
  CQUART_FCONTROLS_RMD; CQUART_FCONTROLS_ALT_RMD;
  CQUART_FVALS_RMD; (s->Symbol(:G_,s)).(CQUART_FNAMES_RMD); (s->Symbol(:G_,s)).(CQUART_FNAMES_ALT_RMD);
  CQUART_FGROUPS_FD;
  :vol252d; :Lvol252d; :LLvol252d; :LDvol252d;
  :Dlat; :LDlat; :ret12m; :Lret12m; :LLret12m; :bm; :Lbm; :nop; :Lnop;
  (s->Symbol(:LD, s)).(CQUART_FVALS)])

#Table controls
const CONTROL0 = Vector{Symbol}()
const T1_CONTROL4 = [:Lbm; :Llmkequity; :LDlat; :Lnop]
const T1_CONTROL10 = [T1_CONTROL4; :L12ret12m; :L24ret12m;
  :fyret12m; :L12vol252d; :L24vol252d; :fyvol252d]

@inline function finiteormissing!(df::AbstractDataFrame, s::Symbol)::Nothing
  df[!,s] .= (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(df[!,s])

  return nothing
end

#@inline Base.isfinite(::Missing)=missing

#NOTE: The below adjusts for logged vs unlogged returns
#const RETAINED_FACT_COLS = (( RETURN_FIELD == :lret) ?
#  [RETURN_FIELD; FACTOR_FIELDS; :lrfr, :lmkt] :
#  [RETURN_FIELD; FACTOR_FIELDS; :rfr, :mkt]

#consider promiting to the below to finometrics
@inline obliterate!(df::DataFrame)::DataFrame = select!(df, Not(names(df)))

#checks if a column or multiple columns are a key
@inline function checkdfkey(df::AbstractDataFrame, keys::Union{Symbol, Vector{Symbol}})
  N::Int = size(df,1)
  Nunique::Int = length(groupby(df, keys))

  (N==Nunique) || error("$keys is not a key!!")

  return nothing
end

#converts gzip to lz4
function gz2lz4(file::String; filetype::String = "csv", location::String = "",
  gzextension::String = "gz", lz4extension::String = "lz4",
  deleteold::Bool = false)

  infile::String = "$location\\$file.$filetype.$gzextension"
  outfile::String = "$location\\$file.$filetype.$lz4extension"

  isfile(infile) || error("file not found. Full path: $location\\$file.$filetype.$gzextension")

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4CompressorStream(outdata))
  write(s, indata)
  close(s)
end

function archivefile(file::String; inpath::String=DATA_PATH, inext::String = "",
  archivepath::String = ARCHIVE_PATH)

  (length(inext) > 0) && (file = "$file.$inext")
  inp::String = "$inpath\\$file"
  outp::String = "$archivepath\\$file.gz"

  open(GzipCompressorStream, outp, "w") do sout
    open(inp) do sin
      write(sout, sin)
    end
  end

  return nothing
end

function unarchivefile(file::String; outpath::String=DATA_PATH, outext::String = "",
  archivepath::String = ARCHIVE_PATH)

  (length(outext) > 0) && (file = "$file.$outext")
  outp::String = "$outpath\\$file"
  inp::String = "$archivepath\\$file.gz"

  open(GzipDecompressorStream, inp) do sin
    open(outp, "w") do sout
      write(sout, sin)
    end
  end

end

function recompressfiles(files = FILES_TO_RECOMPRESS; datapath=DATA_PATH)::Nothing

  for f ∈ files
    gz2lz4(f, location=datapath)
  end

  return nothing
end

recompressfile(file; datapath=DATA_PATH)::Nothing = recompressfiles(
  [file]; datapath=datapath)

function uniquebyyear(df::DataFrame; Frequired::Union{Nothing, Vector{Symbol}} = nothing,
  Fidentifier::Symbol = :permno, Fgroup::Symbol = :fyear)

  for sdf::SubDataFrame ∈ groupby(df, Fgroup)
    local ssdf::SubDataFrame

    if !isnothing(Frequired)
      ssdf = view(sdf, completecases(sdf[!, Frequired]), :)
    else
      ssdf = sdf
    end

    N::Int = length(unique(ssdf[!, Fidentifier]))
    N==0 && continue
    println("$(ssdf[1, Fgroup]): $N unique records by $Fidentifier ",
      isnothing(Frequired) ? "" : "containing $Frequired")
  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
include("genericDistribution.jl")
module uniform
<:genericDistribution #genericModel

    """
    ```Inputs
    Range: assumes Probability Density function is defined, & Continuous on range

    ```
    """
    function pdf(x, range)
        return pdf(x, range)
    end

    function cdf(x, range) #TODO: approx Integral of pdf(x)
        return sum(pdf(x, range)) #sum of all points #assumes all po
    end

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using Test
using MPI
using ClimateMachine.Mesh.Topologies

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    @assert csize == 5

    topology = BrickTopology(
        comm,
        (0:10,);
        boundary = ((1, 2),),
        periodicity = (true,),
    )

    elems = topology.elems
    realelems = topology.realelems
    ghostelems = topology.ghostelems
    sendelems = topology.sendelems
    elemtocoord = topology.elemtocoord
    elemtoelem = topology.elemtoelem
    elemtoface = topology.elemtoface
    elemtoordr = topology.elemtoordr
    elemtobndy = topology.elemtobndy
    nabrtorank = topology.nabrtorank
    nabrtorecv = topology.nabrtorecv
    nabrtosend = topology.nabrtosend

    globalelemtoelem = [
        10 1 2 3 4 5 6 7 8 9
        2 3 4 5 6 7 8 9 10 1
    ]

    globalelemtoface = [
        2 2 2 2 2 2 2 2 2 2
        1 1 1 1 1 1 1 1 1 1
    ]

    globalelemtoordr = ones(Int, size(globalelemtoface))
    globalelemtobndy = zeros(Int, size(globalelemtoface))

    globalelemtocoord = Array{Int}(undef, 1, 2, 10)
    globalelemtocoord[:, :, 1] = [0 1]
    globalelemtocoord[:, :, 2] = [1 2]
    globalelemtocoord[:, :, 3] = [2 3]
    globalelemtocoord[:, :, 4] = [3 4]
    globalelemtocoord[:, :, 5] = [4 5]
    globalelemtocoord[:, :, 6] = [5 6]
    globalelemtocoord[:, :, 7] = [6 7]
    globalelemtocoord[:, :, 8] = [7 8]
    globalelemtocoord[:, :, 9] = [8 9]
    globalelemtocoord[:, :, 10] = [9 10]

    @assert csize == 5
    nrealelem = 2
    if crank == 0
        globalelems = [1, 2, 3, 10]
        elemtoelem_expect = [4 1 3 4; 2 3 3 4]
        nabrtorank_expect = [1, 4]
    elseif crank == 1
        globalelems = [3, 4, 2, 5]
        elemtoelem_expect = [3 1 3 4; 2 4 3 4]
        nabrtorank_expect = [0, 2]
    elseif crank == 2
        globalelems = [5, 6, 4, 7]
        elemtoelem_expect = [3 1 3 4; 2 4 3 4]
        nabrtorank_expect = [1, 3]
    elseif crank == 3
        globalelems = [7, 8, 6, 9]
        elemtoelem_expect = [3 1 3 4; 2 4 3 4]
        nabrtorank_expect = [2, 4]
    elseif crank == 4
        globalelems = [9, 10, 1, 8]
        elemtoelem_expect = [4 1 3 4; 2 3 3 4]
        nabrtorank_expect = [0, 3]
    end
    nabrtorecv_expect = UnitRange{Int}[1:1, 2:2]
    nabrtosend_expect = UnitRange{Int}[1:1, 2:2]

    @test elems == 1:length(globalelems)
    @test realelems == 1:nrealelem
    @test ghostelems == (nrealelem + 1):length(globalelems)
    @test elemtocoord == globalelemtocoord[:, :, globalelems]
    @test elemtoface[:, realelems] ==
          globalelemtoface[:, globalelems[realelems]]
    @test elemtoelem == elemtoelem_expect
    @test elemtobndy == globalelemtobndy[:, globalelems]
    @test elemtoordr == globalelemtoordr[:, globalelems]
    @test nabrtorank == nabrtorank_expect
    @test nabrtorecv == nabrtorecv_expect
    @test nabrtosend == nabrtosend_expect

    @test collect(realelems) ==
          sort(union(topology.exteriorelems, topology.interiorelems))
    @test unique(sort(sendelems)) == topology.exteriorelems
    @test length(intersect(topology.exteriorelems, topology.interiorelems)) == 0
end

main()


text:
function makedoubleindex(files)
    idx = Dict{String, Dict}()
    for file in files
        str = lowercase(read(file, String))
        words = split(str, r"\s+")
        for word in words
            if !haskey(idx, word)
                idx[word] = Dict{String, Int}()
            elseif !haskey(idx[word], file)
                (idx[word])[file] = 1
            else
                (idx[word])[file] += 1
            end
        end
    end
    idx
end

function wordsearch(dict, words::Vector{String})
    for word in words
        if haskey(dict, word)
            for (f, n) in dict[word]
                println("File $f contains $n instances of <$word>.")
            end
        else
            println("No instances of \"$word\" were found.")
        end
    end
end
wordsearch(dict, word::String) = wordsearch(dict, [word])

const filenames = ["file1.txt", "file2.txt", "file3.txt"]
const didx = makedoubleindex(filenames)
const searchterms = ["forehead", "of", "hand", "a", "foot"]
wordsearch(didx, searchterms)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# unary
# TODO: deprecate
@i @inline function Base.:-(a!::GVar)
    -(a!.x)
    -(a!.g)
end

@i @inline function NEG(a!::GVar)
    NEG(a!.x)
    NEG(a!.g)
end

function INV(x!::GVar{T}) where T
    x2 = x!.x ^ 2
    GVar(INV(x!.x), -x!.g * x2)
end

@i @inline function DEC(a!::GVar)
    DEC(a!.x)
end

# +-
@i @inline function ⊖(identity)(a!::GVar, b::GVar)
    a!.x -= b.x
    b.g += a!.g
end
@nograd ⊖(identity)(a!::Real, b::GVar)
@nograd ⊖(identity)(a!::GVar, b::Real)

# inv
@eval @i @inline function ⊖(inv)(out!::GVar{T}, y::GVar) where T
    out!.x -= inv(y.x)
    @routine @invcheckoff begin
        @zeros T a1
        a1 += y.x ^ 2
    end
    y.g -= out!.g / a1
    ~@routine
end
@nograd ⊖(inv)(a!::Real, b::GVar)
@nograd ⊖(inv)(a!::GVar, b::Real)

# +- (triple)
@i @inline function ⊖(+)(out!::GVar, x::GVar, y::GVar)
    out!.x -= x.x + y.x
    x.g += out! |> grad
    y.g += out! |> grad
end

@i @inline function ⊖(+)(out!::GVar, x::GVar, y::Real)
    out!.x -= (x |> value) + (y |> value)
    x.g += out! |> grad
end

@i @inline function ⊖(+)(out!::GVar, x::Real, y::GVar)
    out!.x -= (x |> value) + (y |> value)
    y.g += out! |> grad
end

@i @inline function ⊖(-)(out!::GVar, x::GVar, y::GVar)
    out!.x -= x.x - y.x
    x.g += out! |> grad
    y.g -= out! |> grad
end

@i @inline function ⊖(-)(out!::GVar, x::Real, y::GVar)
    out!.x -= (x |> value) - y.x
    y.g -= out!.g
end

@i @inline function ⊖(-)(out!::GVar, x::GVar, y::Real)
    out!.x -= x.x - (y |> value)
    x.g += out! |> grad
end

# NOTE: it will error on `SWAP(a!::GVar, b)` or `SWAP(a!, b:GVar)`
@i @inline function SWAP(a!::GVar, b!::GVar)
    SWAP(a! |> value,  b! |> value)
    SWAP(a!.g, b!.g)
end

# */
@i @inline function ⊖(*)(out!::GVar, x::GVar, y::GVar)
    out!.x -= x.x * y.x
    x.g += out!.g * y.x
    y.g += x.x * out!.g
end

@i @inline function ⊖(*)(out!::GVar, x::Real, y::GVar)
    out!.x -= (x |> value) * y.x
    y.g += (x |> value) * out!.g
end

@i @inline function ⊖(*)(out!::GVar, x::GVar, y::Real)
    out!.x -= x.x * (y |> value)
    x.g += out!.g * (y |> value)
end

for DIV in [:/, :÷]
@eval @i @inline function ⊖($DIV)(out!::GVar{T}, x::GVar, y::GVar) where T
    out!.x -= $DIV(x.x, y.x)
    @routine @invcheckoff begin
        a1 ← zero(out! |> grad)
        a2 ← zero(out! |> grad)
        a1 += x.x * out!.g
        a2 += $DIV(a1, y.x)
    end
    x.g += $DIV(out!.g, y.x)
    y.g -= $DIV(a2, y.x)
    ~@routine
end

@eval @i @inline function ⊖($DIV)(out!::GVar{T}, x::Real, y::GVar) where T
    out!.x -= $DIV(x, y.x)
    @routine @invcheckoff begin
        a1 ← zero(out!.g)
        a2 ← zero(out!.g)
        a1 += x * out!.g
        a2 += $DIV(a1, y.x)
    end
    y.g -= $DIV(a2, y.x)
    ~@routine
end

@eval @i @inline function ⊖($DIV)(out!::GVar, x::GVar, y::Real)
    out!.x -= $DIV(x.x, y)
    x.g += $DIV(out!.g, y)
end
end

@i @inline function ⊖(^)(out!::GVar{T}, x::GVar, n::GVar) where T
    ⊖(^)(out!.x, x.x, n.x)

    # grad x
    @routine @invcheckoff begin
        @zeros T anc1 anc2 anc3 jac1 jac2
        DEC(n.x)
        anc1 += x.x^n.x
        INC(n.x)
        jac1 += anc1 * n.x

        # get grad of n
        anc2 += log(x.x)
        anc3 += x.x ^ n.x
        jac2 += anc3*anc2
    end
    x.g += out!.g * jac1
    n.g += out!.g * jac2
    ~@routine
end

@i @inline function ⊖(^)(out!::GVar{T}, x::GVar, n::Real) where T
    ⊖(^)(out!.x, x.x, n)
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        jac ← zero(x.x)

        DEC(n |> value)
        anc1 += x.x ^ n
        INC(n |> value)
        jac += anc1 * n
    end
    x.g += out!.g * jac
    ~@routine
end

@i @inline function ⊖(^)(out!::GVar{T}, x::Real, n::GVar) where T
    ⊖(^)(out!.x, x, n.x)
    # get jac of n
    @routine @invcheckoff begin
        anc1 ← zero(x)
        anc2 ← zero(x)
        jac ← zero(x)

        anc1 += log(x)
        anc2 += x ^ n.x
        jac += anc1*anc2
    end
    n.g += out!.g * jac
    ~@routine
end

@i @inline function ⊖(atan)(out!::GVar{T}, y::GVar, x::GVar) where T
    ⊖(atan)(out!.x, y.x, x.x)
    @routine @invcheckoff begin
        @zeros T xy2 jac_x jac_y
        xy2 += abs2(x.x)
        xy2 += abs2(y.x)
        jac_y += x.x / xy2
        jac_x += (-y.x) / xy2
    end
    y.g += out!.g * jac_y
    x.g += out!.g * jac_x
    ~@routine
end

@i @inline function ⊖(atan)(out!::GVar{T}, x::GVar) where T
    ⊖(atan)(out!.x, x.x)
    @routine @invcheckoff begin
        xy2 ← one(T)
        xy2 += abs2(x.x)
    end
    x.g += out!.g / xy2
    ~@routine
end

@i @inline function ⊖(abs)(out!::GVar, x::GVar{T}) where T
    out!.x -= abs(x.x)
    if (x > 0, ~)
        x.g += out!.g
    else
        x.g -= out!.g
    end
end

@i @inline function ⊖(abs2)(out!::GVar, x::GVar{T}) where T
    out!.x -= abs2(x.x)
    x.g += out!.g * x.x
    x.g += out!.g * x.x
end
@nograd ⊖(abs2)(a!::GVar, b::Real)
@nograd ⊖(abs2)(a!::Real, b::GVar)

for op in [:*, :/, :^, :+, :-]
    @eval @nograd ⊖($op)(out!::GVar, x::Real, y::Real)
    @eval @nograd ⊖($op)(out!::Real, x::Real, y::GVar)
    @eval @nograd ⊖($op)(out!::Real, x::GVar, y::GVar)
    @eval @nograd ⊖($op)(out!::Real, x::GVar, y::Real)
end

@i @inline function ⊖(sqrt)(out!::GVar, x::GVar{T}) where T
    if x.x != 0
        @routine @invcheckoff begin
            @zeros T anc1 anc2
            anc1 += sqrt(x.x)
            anc2 += 2 * anc1
        end
        out!.x -= anc1
        x.g += out!.g / anc2
        ~@routine
    end
end

@i @inline function ⊖(exp)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        anc1 ← zero(T)
        anc1 += exp(x.x)
    end
    out!.x -= anc1
    x.g += out!.g * anc1
    ~@routine
end

@i @inline function ⊖(log)(out!::GVar, x::GVar{T}) where T
    out!.x -= log(x.x)
    x.g += out!.g / x.x
end

@i @inline function ⊖(sin)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        @zeros T s c
        (s, c) += sincos(x.x)
    end
    out!.x -= s
    x.g += out!.g * c
    ~@routine
end

@i @inline function ⊖(sinh)(out!::GVar, x::GVar{T}) where T
    out!.x -= sinh(x.x)
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc1 += cosh(x.x)
    end
    x.g += out!.g * anc1
    ~@routine
end

@i @inline function (:-=)(asin)(out!::GVar, x::GVar{T}) where T
    out!.x -= asin(x.x)
    @routine @invcheckoff begin
        @zeros T sqrt_1_x2 x2
        x2 += x.x ^ 2
        sqrt_1_x2 += sqrt(x2 |> NEG |> AddConst(1))
    end
    x.g += out!.g / sqrt_1_x2
    ~@routine
end

@i @inline function (:-=)(cos)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        @zeros T s c
        (s, c) += sincos(x.x)
    end
    out!.x -= c
    x.g -= out!.g * s
    ~@routine
end

@i @inline function ⊖(cosh)(out!::GVar, x::GVar{T}) where T
    out!.x -= cosh(x.x)
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc1 += sinh(x.x)
    end
    x.g += out!.g * anc1
    ~@routine
end

@i @inline function ⊖(acos)(out!::GVar, x::GVar{T}) where T
    out!.x -= acos(x.x)
    @routine @invcheckoff begin
        @zeros T sqrt_1_x2 x2
        x2 += x.x ^ 2
        sqrt_1_x2 += sqrt(x2 |> NEG |> AddConst(1))
    end
    x.g -= out!.g / sqrt_1_x2
    ~@routine
end

@i @inline function ⊖(tan)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc2 ← one(x.x)
        anc1 += tan(x.x)
        anc2 += anc1^2
    end
    out!.x -= anc1
    x.g += out!.g * anc2
    ~@routine
end

@i @inline function ⊖(tanh)(out!::GVar, x::GVar{T}) where T
    @routine @invcheckoff begin
        anc1 ← zero(x.x)
        anc2 ← one(x.x)
        anc1 += tanh(x.x)
        anc2 -= anc1^2
    end
    out!.x -= anc1
    x.g += out!.g * anc2
    ~@routine
end

@i @inline function ⊖(sincos)(out!::Tuple{T1,T1}, x::GVar{T}) where {T1<:GVar, T}
    @routine @invcheckoff begin
        s ← zero(T)
        c ← zero(T)
        (s, c) += sincos(x.x)
    end
    (out! .|> value) -= (s, c)
    x.g += (out!.:1 |> grad) * c
    x.g -= (out!.:2 |> grad) * s
    ~@routine
end

for op in [:sqrt, :exp, :log, :sin, :cos, :tanh]
    @eval @nograd ⊖($op)(out!::Real, x::GVar)
    @eval @nograd ⊖($op)(out!::GVar, x::Real)
end

@nograd ⊖(sincos)(out!::Tuple{<:Real,<:Real}, x::GVar)
@nograd ⊖(sincos)(out!::Tuple{<:GVar,<:GVar}, x::Real)

@i @inline function IROT(a!::GVar, b!::GVar, θ::GVar)
    IROT(a!.x, b!.x, θ.x)
    NEG(θ |> value)
    θ.x -= π/2
    ROT(a!.g, b!.g, θ.x)
    θ.g += a!.x * a!.g
    θ.g += b!.x * b!.g
    θ.x += π/2
    NEG(θ |> value)
    ROT(a!.g, b!.g, π/2)
end

@i @inline function IROT(a!::GVar, b!::GVar, θ::Real)
    IROT(a!.x, b!.x, θ)
    NEG(θ)
    θ -= π/2
    ROT(a!.g, b!.g, θ)
    θ += π/2
    NEG(θ)
    ROT(a!.g, b!.g, π/2)
end

@nograd IROT(a!::Real, b!::Real, θ::GVar)

export primitive_grad
function primitive_grad end

@i function (mf::MinusEq)(out!::GVar, args...; kwargs...)
    out!.x -= mf.f((args .|> value)...; kwargs...)
    (args .|> grad) .+= (@skip! out!.g) .* (@skip! primitive_grad(mf.f, (args .|> value)...; kwargs...))  # unsafe statement, error on recursive gradient
end

@i function (mf::MinusEq)(out!::GVar, x::GVar; kwargs...)
    out!.x -= mf.f(x .|> value; kwargs...)
    x.g += (@skip! out!.g) * (@skip! primitive_grad(mf.f, x.x; kwargs...))  # unsafe statement
end

@i function :(-=)(convert)(out!::GVar{Tx, Tg}, y::GVar) where {Tx, Tg}
    out!.x -= convert(y.x)
    y.g += convert(out!.g)
end

@i function HADAMARD(x::GVar, y::GVar)
    HADAMARD(x.x, y.x)
    HADAMARD(x.g, y.g)
end

@i function (f::AddConst)(y::GVar)
    y.x += f.x
end

# more data views
for (DT, OP, NOP) in [(:AddConst, :+, :-, :add), (:SubConst, :-, :+)]
    @eval chfield(x::GVar, ac::$DT, xval::GVar) = GVar($NOP(xval.x, ac.x), xval.g)
end

#chfield(x::T, ::typeof(INV), xval::T) where T<:GVar = GVar(INV(xval.x), -xval.g*(xval.x^2))
#chfield(x::T, ::typeof(NEG), xval::T) where T<:GVar = GVar(-xval.x, -xval.g)

for F in [:INV, :NEG, :FLIP, :INC, :DEC]
    @eval NiLangCore.chfield(x::T, ::typeof($F), xval::T) where T<:GVar = (~$F)(xval)
end

# accumulation on arrays: initially for Bennett algorithm
# TODO: also define it for composite types. or maybe a macro for it.
@i function :(+=)(identity)(target::AbstractArray, source::AbstractArray)
    @safe @assert length(target) == length(source)
    @inbounds for i=1:length(target)
        target[i] += source[i]
    end
end


text:
"""
    read_nnet(fname::String; last_layer_activation = Id())

Read in neural net from a `.nnet` file and return Network struct.
The `.nnet` format is borrowed from [NNet](https://github.com/sisl/NNet).
The format assumes all hidden layers have ReLU activation.
Keyword argument `last_layer_activation` sets the activation of the last
layer, and defaults to `Id()`, (i.e. a linear output layer).
"""
function read_nnet(fname::String; last_layer_activation = Id())
    f = open(fname)
    line = readline(f)
    while occursin("//", line) #skip comments
        line = readline(f)
    end
    # number of layers
    nlayers = parse(Int64, split(line, ",")[1])
    # read in layer sizes
    layer_sizes = parse.(Int64, split(readline(f), ",")[1:nlayers+1])
    # read past additonal information
    for i in 1:5
        line = readline(f)
    end
    # i=1 corresponds to the input dimension, so it's ignored
    layers = Layer[read_layer(dim, f) for dim in layer_sizes[2:end-1]]
    push!(layers, read_layer(last(layer_sizes), f, last_layer_activation))

    return Network(layers)
end

"""
    read_layer(output_dim::Int, f::IOStream, [act = ReLU()])

Read in layer from nnet file and return a `Layer` containing its weights/biases.
Optional argument `act` sets the activation function for the layer.
"""
function read_layer(output_dim::Int64, f::IOStream, act = ReLU())

    rowparse(splitrow) = parse.(Float64, splitrow[findall(!isempty, splitrow)])
     # first read in weights
     W_str_vec = [rowparse(split(readline(f), ",")) for i in 1:output_dim]
     weights = vcat(W_str_vec'...)
     # now read in bias
     bias_string = [split(readline(f), ",")[1] for j in 1:output_dim]
     bias = rowparse(bias_string)
     # activation function is set to ReLU as default
     return Layer(weights, bias, act)
end

"""
Prepend `//` to each line of a string.
"""
to_comment(txt) = "//"*replace(txt, "\n"=>"\n//")

"""
    print_layer(file::IOStream, layer)

Print to `file` an object implementing `weights(layer)` and `bias(layer)`
"""
function print_layer(file::IOStream, layer)
   print_row(W, i) = println(file, join(W[i,:], ", "), ",")
   W = layer.weights
   b = layer.bias
   [print_row(W, row) for row in axes(W, 1)]
   [println(file, b[row], ",") for row in axes(W, 1)]
end

"""
    print_header(file::IOStream, network[; header_text])

The NNet format has a particular header containing information about the network size and training data.
`print_header` does not take training-related information into account (subject to change).
"""
function print_header(file::IOStream, network; header_text="")
   println(file, to_comment(header_text))
   layer_sizes = [size(layer.weights, 1) for layer in network.layers] # doesn't include the input layer
   pushfirst!(layer_sizes, size(network.layers[1].weights, 1)) # add the input layer

   # num layers, num inputs, num outputs, max layer size
   num_layers = length(network.layers)
   num_inputs = layer_sizes[1]
   num_outputs = layer_sizes[end]
   max_layer = maximum(layer_sizes[1:end-1]) # chop off the output layer for the maximum,
   println(file, join([num_layers, num_inputs, num_outputs, max_layer], ", "), ",")
   #layer sizes input, ..., output
   println(file, join(layer_sizes, ", "), ",")
   # empty
   println(file, "This line extraneous")
   # minimum vals of inputs
   println(file, string(join(fill(-floatmax(Float16), num_inputs), ","), ","))
   # maximum vals of inputs
   println(file, string(join(fill(floatmax(Float16), num_inputs), ","), ","))
   # mean vals of inputs + 1 for output
   println(file, string(join(fill(0.0, num_inputs+1), ","), ","))
   # range vals of inputs + 1 for output
   println(file, string(join(fill(1.0, num_inputs+1), ","), ","))
   return nothing
end

"""
    write_nnet(filename, network[; header_text])

Write `network` to \$filename.nnet.
Note: Does not perform safety checks on inputs, so use with caution.

Based on python code at https://github.com/sisl/NNet/blob/master/utils/writeNNet.py
and follows .nnet format given here: https://github.com/sisl/NNet.
"""
function write_nnet(outfile, network; header_text="Default header text.\nShould replace with the real deal.")
    name, ext = splitext(outfile)
    outfile = name*".nnet"
    open(outfile, "w") do f
        print_header(f, network, header_text=header_text)
        for layer in network.layers
            print_layer(f, layer)
        end
    end
    nothing
end
"""
    compute_output(nnet::Network, input::Vector{Float64})

Propagate a given vector through a nnet and compute the output.
"""
function compute_output(nnet::Network, input)
    curr_value = input
    for layer in nnet.layers # layers does not include input layer (which has no weights/biases)
        curr_value = layer.activation(affine_map(layer, curr_value))
    end
    return curr_value # would another name be better?
end

"""
    get_activation(L, x::Vector)
Finds the activation pattern of a vector `x` subject to the activation function given by the layer `L`.
Returns a Vector{Bool} where `true` denotes the node is "active". In the sense of ReLU, this would be `x[i] >= 0`.
"""
get_activation(L::Layer{ReLU}, x::Vector) = x .>= 0.0
get_activation(L::Layer{Id}, args...) = trues(n_nodes(L))

"""
    get_activation(nnet::Network, x::Vector)

Given a network, find the activation pattern of all neurons at a given point x.
Returns Vector{Vector{Bool}}. Each Vector{Bool} refers to the activation pattern of a particular layer.
"""
function get_activation(nnet::Network, x::Vector{Float64})
    act_pattern = Vector{Vector{Bool}}(undef, length(nnet.layers))
    curr_value = x
    for (i, layer) in enumerate(nnet.layers)
        curr_value = affine_map(layer, curr_value)
        act_pattern[i] = get_activation(layer, curr_value)
        curr_value = layer.activation(curr_value)
    end
    return act_pattern
end

"""
    get_activation(nnet::Network, input::Hyperrectangle)

Given a network, find the activation pattern of all neurons for a given input set.
Assume ReLU.
return Vector{Vector{Int64}}.
- 1: activated
- 0: undetermined
- -1: not activated
"""
function get_activation(nnet::Network, input::Hyperrectangle)
    bounds = get_bounds(nnet, input, before_act = true)
    return get_activation(nnet, bounds)
end

"""
    get_activation(nnet::Network, bounds::Vector{Hyperrectangle})

Given a network, find the activation pattern of all neurons given the node-wise bounds.
Assume ReLU. Assume pre-activation bounds where the bounds on the input are given by the first
hyperrectangle, the first hidden layer by the second hyperrectangle, and so on.
return Vector{Vector{Int64}}.
- 1: activated
- 0: undetermined
- -1: not activated
"""
function get_activation(nnet::Network, bounds::Vector{Hyperrectangle})
    act_pattern = Vector{Vector{Int}}(undef, length(nnet.layers))
    for (i, layer) in enumerate(nnet.layers)
        act_pattern[i] = get_activation(layer, bounds[i+1])
    end
    return act_pattern
end

"""
    get_activation(L::Layer{ReLU}, bounds::Hyperrectangle)

Given a layer, find the activation pattern of all neurons in the layer given the node-wise bounds.
Assume ReLU. Assume bounds is the pre-activation bounds for each ReLU in the layer.
return Vector{Vector{Int64}}.
- 1: activated
- 0: undetermined
- -1: not activated
"""
function get_activation(L::Layer{ReLU}, bounds::Hyperrectangle)
    lower = low(bounds)
    upper = high(bounds)
    act_pattern = zeros(n_nodes(L))
    for j in 1:n_nodes(L) # For evey node
        if lower[j] > 0.0
            act_pattern[j] = 1
        elseif upper[j] < 0.0
            act_pattern[j] = -1
        end
    end
    return act_pattern
end

"""
    get_gradient(nnet::Network, x::Vector)

Given a network, find the gradient at the input x
"""
function get_gradient(nnet::Network, x::Vector)
    z = x
    gradient = Matrix(1.0I, length(x), length(x))
    for (i, layer) in enumerate(nnet.layers)
        z_hat = affine_map(layer, z)
        σ_gradient = act_gradient(layer.activation, z_hat)
        gradient = Diagonal(σ_gradient) * layer.weights * gradient
        z = layer.activation(z_hat)
    end
    return gradient
end

"""
    act_gradient(act, z_hat::Vector{N}) where N

Compute the gradient of an activation function at point z_hat.
Currently only supports ReLU and Id.
"""
act_gradient(act::ReLU, z_hat::Vector) = z_hat .>= 0.0
act_gradient(act::Id,   z_hat::Vector) = trues(length(z_hat))

"""
    relaxed_relu_gradient(l::Real, u::Real)

Return the slope of a ReLU activation based on its lower and upper bounds

Returns `0` if u<0, `1` if l>0, `u/(u-l)` otherwise
"""
function relaxed_relu_gradient(l::Real, u::Real)
    u <= 0.0 && return 0.0
    l >= 0.0 && return 1.0
    return u / (u - l)
end

"""
    act_gradient_bounds(nnet::Network, input::AbstractPolytope)

Compute the bounds on the gradient of all activation functions given an input set.
Currently only support ReLU.
Return:
- `LΛ, UΛ::NTuple{2, Vector{BitVector}}`: lower and upper bounds on activation
"""
function act_gradient_bounds(nnet::Network, input::AbstractPolytope)
    # get the pre-activation bounds, and get rid of the input set
    bounds = get_bounds(nnet, input, before_act=true)
    popfirst!(bounds)

    LΛ = Vector{BitVector}(undef, 0)
    UΛ = Vector{BitVector}(undef, 0)
    for (i, layer) in enumerate(nnet.layers)
        l = act_gradient(layer.activation, low(bounds[i]))
        u = act_gradient(layer.activation, high(bounds[i]))
        push!(LΛ, l)
        push!(UΛ, u)
    end
    return (LΛ, UΛ)
end

"""
    get_gradient_bounds(nnet::Network, LΛ::Vector{AbstractVector}, UΛ::Vector{AbstractVector})
    get_gradient_bounds(nnet::Network, input::AbstractPolytope)

Get lower and upper bounds on network gradient for given gradient bounds on activations, or given an input set.
Return:
- `(LG, UG)::NTuple{2, Matrix{Float64}` lower and upper bounds.
"""
function get_gradient_bounds(nnet::Network, input::AbstractPolytope)
    LΛ, UΛ = act_gradient_bounds(nnet, input)
    return get_gradient_bounds(nnet, LΛ, UΛ)
end
function get_gradient_bounds(nnet::Network, LΛ::Vector{<:AbstractVector}, UΛ::Vector{<:AbstractVector})
    n_input = size(nnet.layers[1].weights, 2)
    LG = Matrix(1.0I, n_input, n_input)
    UG = Matrix(1.0I, n_input, n_input)
    for (i, layer) in enumerate(nnet.layers)
        LG_hat, UG_hat = interval_map(layer.weights, LG, UG)
        LG = Diagonal(LΛ[i]) * max.(LG_hat, 0) + Diagonal(UΛ[i]) * min.(LG_hat, 0)
        UG = Diagonal(LΛ[i]) * min.(UG_hat, 0) + Diagonal(UΛ[i]) * max.(UG_hat, 0)
    end
    return (LG, UG)
end

"""
    interval_map(W::Matrix, l::AbstractVecOrMat, u::AbstractVecOrMat)

Simple linear mapping on intervals.
`L, U := ([W]₊*l + [W]₋*u), ([W]₊*u + [W]₋*l)`

Outputs:
- `(lbound, ubound)` (after the mapping)
"""
function interval_map(W::AbstractMatrix{N}, l::AbstractVecOrMat, u::AbstractVecOrMat) where N
    l_new = max.(W, zero(N)) * l + min.(W, zero(N)) * u
    u_new = max.(W, zero(N)) * u + min.(W, zero(N)) * l
    return (l_new, u_new)
end

"""
    get_bounds(problem::Problem)
    get_bounds(nnet::Network, input::Hyperrectangle, [true])

Computes node-wise bounds given a input set. The optional last
argument determines whether the bounds are pre- or post-activation.

Return:
- `Vector{Hyperrectangle}`: bounds for all nodes. `bounds[1]` is the input set.
"""
function get_bounds(nnet::Network, input; before_act::Bool = false) # NOTE there is another function by the same name in convDual. Should reconsider dispatch
    input = overapproximate(input, Hyperrectangle)
    bounds = Vector{Hyperrectangle}(undef, length(nnet.layers) + 1)
    bounds[1] = input
    b = input
    for (i, layer) in enumerate(nnet.layers)
        if before_act
            bounds[i+1] = approximate_affine_map(layer, b)
            b = approximate_act_map(layer, bounds[i+1])
        else
            b = approximate_affine_map(layer, bounds[i])
            bounds[i+1] = approximate_act_map(layer, b)
        end
    end
    return bounds
end
get_bounds(problem::Problem; kwargs...) = get_bounds(problem.network, problem.input; kwargs...)

"""
    affine_map(layer, x)

Compute W*x ⊕ b for a vector or LazySet `x`
"""
affine_map(layer::Layer, x) = layer.weights*x + layer.bias
function affine_map(layer::Layer, x::LazySet)
    LazySets.affine_map(layer.weights, x, layer.bias)
end

"""
   approximate_affine_map(layer, input::Hyperrectangle)

Returns a Hyperrectangle overapproximation of the affine map of the input.
"""
function approximate_affine_map(layer::Layer, input::Hyperrectangle)
    c = affine_map(layer, input.center)
    r = abs.(layer.weights) * input.radius
    return Hyperrectangle(c, r)
end

"""
   approximate_act_map(layer, input::Hyperrectangle)

Returns a Hyperrectangle overapproximation of the activation map of the input.
`act`must be monotonic.
"""
function approximate_act_map(act::ActivationFunction, input::Hyperrectangle)
    β    = act.(input.center)
    βmax = act.(high(input))
    βmin = act.(low(input))
    c    = (βmax + βmin)/2
    r    = (βmax - βmin)/2
    return Hyperrectangle(c, r)
end

approximate_act_map(layer::Layer, input::Hyperrectangle) = approximate_act_map(layer.activation, input)

"""
    split_interval(dom, i)

Split a set into two at the given index.

Inputs:
- `dom::Hyperrectangle`: the set to be split
- `i`: the index to split at
Return:
- `(left, right)::Tuple{Hyperrectangle, Hyperrectangle}`: two sets after split
"""
function split_interval(dom::Hyperrectangle, i::Int64)
    input_lower, input_upper = low(dom), high(dom)

    input_upper[i] = dom.center[i]
    input_split_left = Hyperrectangle(low = input_lower, high = input_upper)

    input_lower[i] = dom.center[i]
    input_upper[i] = dom.center[i] + dom.radius[i]
    input_split_right = Hyperrectangle(low = input_lower, high = input_upper)
    return (input_split_left, input_split_right)
end

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

function isbounded(input)
    if input isa HPolytope
        return LazySets.isbounded(input, false)
    else
        return LazySets.isbounded(input)
    end
end

is_hypercube(set::Hyperrectangle) = all(iszero.(set.radius .- set.radius[1]))
is_halfspace_equivalent(set) = length(constraints_list(set)) == 1


text:
#=
day13:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-13
=#

using AoC, Test

parseCoordinate(line) = split(line,",") .|> c -> parse(Int, c)
@test parseCoordinate("12,56")==[12,56]

function parseFold(line)
    matches = match(r"fold along (.)=(.*)",line)
    return (matches.captures[1] == "x" ? 1 : 2, parse(Int, matches.captures[2]))
end
@test parseFold("fold along y=7") == (2,7)
@test parseFold("fold along x=5") == (1,5)

function parseFile(lines)
    dots = Set()
    folds = []

    mode = :dots
    for line in lines
        if line == ""
            mode = :folds
        elseif mode == :dots
            push!(dots, parseCoordinate(line))
        else
            push!(folds, parseFold(line))
        end
    end
    return (dots, folds)
end
@test parseFile(["3,4","5,6","","fold along x=7"]) == (Set([[3,4],[5,6]]), [(1,7)])

fold(dots, f) = fold(dots, f...)
function fold(dots, direction, location)
    result = Set()
    for dot in dots
        resultDot = [dot...]
        if dot[direction] > location
            resultDot[direction] = 2*location - dot[direction]
        end
        push!(result, resultDot)
    end
    return result
end
@test fold(Set([[3,5]]), (1,2)) == Set([[1,5],])
@test fold(Set([[3,5],(1,5)]), (1,2)) == Set([[1,5],])
@test fold(Set([[3,5]]), (1,4)) == Set([[3,5],])

function part1(lines)
    (dots, folds) = parseFile(lines)

    dots = fold(dots, folds[1])

    return length(dots)
end
@test part1(exampleLines(13,1)) == 17

using FileIO, ColorTypes, FixedPointNumbers

# do not transpose RGB
Base.transpose(x::RGB) = x

function part2(lines)
    (dots, folds) = parseFile(lines)

    for f in folds
        dots = fold(dots, f)
    end

    origin = CartesianIndex(1,1)

    image = zeros(RGB{N0f8}, max(get.(dots,1,0)...)+1, max(get.(dots,2,0)...)+1)
    for dot in dots
        image[CartesianIndex(dot...) + origin] = RGB(1,0,0)
    end

    image = transpose(image)

    for x in 1:size(image,1)
        line = ""
        for y in 1:size(image,2)
            line *= image[x,y] == RGB{N0f8}(0) ? "  " : "██"
        end
        println(line)
    end

    scale = 10
    image2 = similar(image, size(image).*scale)
    for I in CartesianIndices(image2)
        image2[I] = image[CartesianIndex(Tuple(I-origin) .÷ scale) + origin]
    end

    save("target/day13.gif", image2)
end
@test part1(exampleLines(13,1)) == 17

@show lines(13) |> ll -> @time part1(ll)

@time part2(lines(13))


text:
# playTime.jl

#global flog=open("log.playTime","a+")

##文法
include("grammer.jl")

qvar=[誰 何 何処]

#魔法世界
const QB=:QB
const ほむら=:ほむら
const 暁美=:暁美
const まどか=:まどか
const さやか=:さやか
const 上条=:上条
const 杏子=:杏子
const まみ=:まみ
const ひとみ=:ひとみ
const 魔法少女=:魔法少女
const 魔女=:魔女
const ソウルジェム=:ソウルジェム
const ワルプルギスの夜=:ワルプルギスの夜

魔法会議_0=[
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [暁美 の ソウルジェム は きれい],
 [さやか の ソウルジェム は きれい],
 [まみ の ソウルジェム は きれい]
]
魔法会議_1=[
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [暁美 の ソウルジェム は きれい],
 [さやか の ソウルジェム は 濁っている],
 [まみ の ソウルジェム は きれい]
]
魔法会議_2=[
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [暁美 の ソウルジェム は きれい],
 [さやか の ソウルジェム は 濁っている],
 [まみ の ソウルジェム は きれい],
 [さやか は 魔女 です],
]

魔法会議_1n=[
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [暁美 の ソウルジェム は きれい],
 [さやか の ソウルジェム は きれい でない],
 [まみ の ソウルジェム は きれい]
]
魔法会議_2n=[
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [暁美 の ソウルジェム は きれい],
 [さやか の ソウルジェム は きれい でない],
 [まみ の ソウルジェム は きれい],
 [さやか は 魔女 です],
]

魔法会議_01=[魔法会議_0;魔法会議_1]
魔法会議_012=[魔法会議_01;魔法会議_2]

m0=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は きれい]],魔法会議_0)
m1=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は 濁っている]],魔法会議_1)
m2=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は 濁っている]],魔法会議_2)
m3=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は 濁っている],[:x は 魔女 です]],魔法会議_2)

mn0=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は きれい でない]],魔法会議_0)
m1n=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は きれい でない]],魔法会議_1n)
m2n=solveQ([:x,:y,:z],[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は きれい でない],[:x は 魔女 です]],魔法会議_2n)

### 推論
#==
1. [:x は きれい でない] <=> [:y の ソウルジェム は きれい　でない] => :x <- [:y の ソウルジェム]
2. [:x は 濁っている] <= [:x は きれい でない]
3. [:x は きれい でない] => [:y の ソウルジェム は きれい でない]

==#
N=[:x,:y,:z]
Q0=[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は 濁っている],[:x は 魔女 です]]
Q1=[[:x は ソウルジェム を 持っている],[:x の ソウルジェム は きれい でない],[:x は 魔女 です]]

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using MacroTools: @q
using GeneralizedGenerated
using MLStyle

using SymbolicUtils
using SymbolicUtils: Sym, Symbolic, symtype, term

using MappedArrays
import SpecialFunctions
using SpecialFunctions: logfactorial

using NestedTuples: schema
import NestedTuples

using FillArrays

export schema

export symlogdensity

symlogdensity(d, x::Symbolic) = logdensity(d,x)

function symlogdensity(d::ProductMeasure{<:AbstractArray}, x::Symbolic{A}) where {A <: AbstractArray}
    dims = size(d)

    iters = Sym{Int}.(gensym.(Symbol.(:i, 1:length(dims))))

    marginals = d.data

    # To begin, the result is just the summand
    result = getsummand(marginals, x, iters)

    # Then we wrap in a summation index for each dimension
    for i in 1:length(dims)
        result = Sum(result, iters[i], 1, dims[i])
    end

    return result
end

function getsummand(marginals::AbstractMappedArray, x::Symbolic, iters)
    T = eltype(marginals.data)
    symlogdensity(marginals.f(term(getindex, marginals.data, iters...; type=T)), x[iters...])
end

function getsummand(marginals::Fill, x::Symbolic, iters)
    symlogdensity(marginals.value, x[iters...])
end

function NestedTuples.schema(cm::ConditionalModel)
    trace = simulate(cm; trace_assignments=true).trace
    types = schema(merge(trace, argvals(cm)))
    return types
end

function symlogdensity(cm::ConditionalModel{A,B,M}; noinline=()) where {A,B,M}
    trace = simulate(cm; trace_assignments=true).trace
    types = schema(merge(trace, argvals(cm)))
    # m = symify(cm.model)
    m = cm.model
    dict = symdict(cm; noinline=noinline)
    symlogdensity(m, types, dict, trace)
end

function symlogdensity(m::Model{A,B,M}, types, dict, trace) where {A,B,M}
    s = SymCall(convert(Dict,trace))(_symlogdensity,M, m, to_type(types))
    s = rewrite(s)
    return SymbolicCodegen.foldconstants(s, dict)
end

# Convert a named tuple to a dictionary for symbolic substitution
symdict(nt::NamedTuple; noinline=()) = Dict((k => v for (k,v) in pairs(nt) if k ∉ noinline))
symdict(cm::ConditionalModel; noinline=()) = symdict(merge(cm.argvals, cm.obs); noinline=noinline)

function sourceSymlogdensity(cm::ConditionalModel{A,B,M}) where {A,B,M}
    types = schema(cm)
    return sourceSymlogdensity(types)(Model(cm))
end

# Convert a type into the SymbolicUtils type we'll use to represent it
# for example,
#     julia> SymbolicCodegen.sym(Int)
#     :(SymbolicCodegen.Sym{Int64})
#
#     julia> SymbolicCodegen.sym(Int, :n)
#     :(SymbolicCodegen.Sym{Int64}(:n))
#

sym(T::Type) = :(Soss.Sym{$T})

sym(T::Type, s::Symbol) = :($(sym(T))($(QuoteNode(s))))

export sourceSymlogdensity
function sourceSymlogdensity(types)
    sym(s::Symbol) = Soss.sym(getproperty(types, s), s)

    function(_m::Model)
        function proc(_m, st :: Assign)
            x = st.x
            xname = QuoteNode(x)
            xsym = sym(x)
            return :($x = $xsym)
            # rhs = st.rhs
            # return :($x = $rhs)
        end

        function proc(_m, st :: Sample)
            x = st.x
            xname = QuoteNode(x)
            xsym = sym(x)
            s = @q begin
                $x = $xsym
                _ℓ += symlogdensity($(st.rhs), $x)
            end
        end
        proc(_m, st :: Return)     = nothing
        proc(_m, st :: LineNumber) = nothing

        function wrap(kernel)
            q = @q begin
                _ℓ = 0.0
            end

            for x in arguments(_m)
                xname = QuoteNode(x)
                xsym = sym(x)
                push!(q.args, :($x = $xsym))
            end

            @q begin
                $q
                $kernel
                return _ℓ
            end
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

@gg function _symlogdensity(M::Type{<:TypeLevel}, _m::Model, ::Type{T}) where {T}
    types = GeneralizedGenerated.from_type(T)
    Sym = SymbolicUtils.Sym
    body = type2model(_m) |> sourceSymlogdensity(types)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
using Documenter, PoreMatMod, PlutoSliderServer

# run the Pluto example notebooks and export them to HTML for inclusion in the docs
cd("examples") # next line fails if not actually in examples/
PlutoSliderServer.export_directory() # runs each notebook in examples/ and exports to HTML
rm("index.html") # previous line makes additional table-of-contents page (not needed)
export_path = "../docs/src/examples/"
if !isdir(export_path)
    mkdir(export_path)
end
for file in readdir() # loop over files in examples/
    if length(split(file, ".html")) > 1 # select only the HTML files
        @info "Staging file $file"
        mv(file, export_path * file, force=true) # move HTML files to docs source folder for build
    end
end
cd("..") # return to root for running deploydocs

# build the docs
makedocs( # to test docs, run this call to `makedocs` in the REPL
    root=joinpath(dirname(pathof(PoreMatMod)), "..", "docs"),
    modules=[PoreMatMod, Xtals],
    sitename="PoreMatMod.jl",
    clean=true,
    pages=[
        "PoreMatMod" => "index.md",
        "Manual" => [
            "Getting Started" => "manual/start.md",
            "Loading Data" => "manual/inputs.md",
            "Substructure Search" => "manual/find.md",
            "Substructure Find/Replace" => "manual/replace.md"],
        "Examples" => "examples.md",
        "PoreMatModGO" => "PoreMatModGO.md",
        "Contribute/Report Issues" => "collab.md"
    ],
    format=Documenter.HTML(assets = ["assets/flux.css"]),
    push_preview=true,
    doctest=false # doctests are run in testing; running them here is redundant and slow
)

# deploy the docs
deploydocs(repo = "github.com/SimonEnsemble/PoreMatMod.jl.git", versions = ["latest" => "v^", "v#.#.#", "dev" => "master"])


text:
__precompile__()

module FixedPolynomials
    import MultivariatePolynomials
    const MP = MultivariatePolynomials
    import Base: ==, convert, promote_rule

    abstract type AbstractPolySystem{T} end
    export AbstractPolySystem

    #import Base: gradient

    include("poly.jl")
    include("show.jl")
    include("convert_promote.jl")
    include("tables.jl")
    include("config.jl")

end


text:
"""
    eom_logistic4(u, p, t)

Equations of motion for a 4D transitive causal chain of unidirectionally
coupled logistic maps, where X1 -> X2 -> X3 -> X4 [1].

# References
Ye, Hao, et al. "Distinguishing time-delayed causal interactions using
convergent cross mapping." Scientific reports 5 (2015): 14750
"""
function eom_logistic4(u, p, t)
    r₁, r₂, r₃, r₄, c₂, c₃, c₄  = (p...,)
    y₁, y₂, y₃, y₄ = (u...,)

    dy₁ = y₁*(r₁ - r₁*y₁)
    dy₂ = y₂*(r₂ - c₂*y₁ - r₂*y₂)
    dy₃ = y₃*(r₃ - c₃*y₂ - r₃*y₃)
    dy₄ = y₄*(r₄ - c₄*y₃ - r₄*y₄)
    return SVector{4}(dy₁, dy₂, dy₃, dy₄)
end

function logistic4(u₀, r₁, r₂, r₃, r₄, c₂, c₃, c₄)
    p = [r₁, r₂, r₃, r₄, c₂, c₃, c₄]
    DiscreteDynamicalSystem(eom_logistic4, u₀, p)
end

"""
    logistic4(;u₀ = rand(4), r₁ = 3.9, r₂ = 3.6, r₃ = 3.6, r₄ = 3.8,
        c₂ = 0.4, c₃ = 0.4, c₄ = 0.35) -> DiscreteDynamicalSystem

Initialise a system of a a transitive chain of four unidirectionally coupled
logistic maps, where X1 -> X2 -> X3 -> X4 [1]. Default parameters are as in
Ye et al. (2015).

# References
Ye, Hao, et al. "Distinguishing time-delayed causal interactions using convergent cross mapping." Scientific reports 5 (2015): 14750
"""
logistic4(;u₀ = rand(4),
            r₁ = 3.9, r₂ = 3.6, r₃ = 3.6, r₄ = 3.8,
            c₂ = 0.4, c₃ = 0.4, c₄ = 0.35) =
    logistic4(u₀, r₁, r₂, r₃, r₄, c₂, c₃, c₄)