text:
model = TrajectoryOptimization.Dynamics.pendulum
n = model.n; m = model.m
model_d = rk3(model)

# cost
Q = Array(1e-3*Diagonal(I,n))
R = Array(1e-3*Diagonal(I,m))
Qf = Q
x0 = zeros(n)
xf = [pi;0.0]

# options
verbose=false
opts_ilqr = iLQRSolverOptions{T}(verbose=false,live_plotting=:off)
opts_al = AugmentedLagrangianSolverOptions{T}(verbose=false,opts_uncon=opts_ilqr,iterations=50,penalty_scaling=10.0)
opts_altro = ALTROSolverOptions{T}(verbose=false,opts_al=opts_al,R_minimum_time=15.0,dt_max=0.15,dt_min=1.0e-3)

# constraints
u_bnd = 3.
bnd = BoundConstraint(n,m,u_min=-u_bnd,u_max=u_bnd)
goal_con = goal_constraint(xf)

# problem
N = 31
U = [ones(m) for k = 1:N-1]
constraints = Constraints(N)
for k = 1:N-1
    constraints[k] += bnd
end
constraints[N] += goal_con
obj = LQRObjective(Q,R,Qf,xf,N)

dt = 0.15
pendulum = Problem(model_d,obj,U,constraints=constraints,dt=dt,x0=x0,xf=xf,N=N)


text:
using SerialPorts
using Base.Test

@show list_serialports()

@linux_only SerialPorts.in_dialout()
SerialPorts.check_serial_access()

ext: jl


text:
include("../src/pwl_approx.jl")
tol = 1e-10
k = 20
r = floor(Int64,(k)/2)
Is = collect(0:k)
# w = min.(Is,4)            # delta-linear

# Playing around with numbers, we see that the asymptotically
# bad behavior of sqrt doesn't set in until extremely large hypereges
w = Is .* (k .- Is)       # clique
w = sqrt.(Is)
# w = log2.(Is .+ 1)
w = Is.^(.8)
w = w[1:r+1]
Is = Is[1:r+1]

pl = scatter(Is,w,legend = false)

# next

epsi = 0.0
Lines, M, X, Y, Ranges = get_pwl_approx(w,epsi,false)
xs = 0:0.01:Is[end]
mlast = M[end]

if mlast < 0
    # negative sloped lines not needed
    # replace with flat line
    M[end] = 0
    X[end] = Is[end]
    Y[end] = w[end]
    pop!(Lines)
elseif Ranges[end][1] == r-1 && Ranges[end-1][2] == r-1 && mlast > 0
    # Covered point r-1 twice. Get rid of it
    M[end] = 0
    X[end] = Is[end]
    Y[end] = w[end]
    pop!(Lines)
elseif mlast > 0
    push!(M,0)
    push!(X,Is[end])
    push!(Y,w[end])
end
push!(Lines, x->w[end]*ones(length(x)))

# construct gadget
b = zeros(length(X)-1)
for t = 1:length(X)-1
    x,y = line_intersection(M[t],M[t+1],X[t],X[t+1],Y[t],Y[t+1])
    b[t] = x
end

a = zeros(length(M)-1)
for i = 1:length(a)
    a[i] = M[i]- M[i+1]
end
J = length(a)

p = x -> sum(a[j].*min.(x,b[j]) for j = 1:J)

wapprox = p(0:r)
@show length(Lines)

# plotting
for l in Lines[1:end]
    plot!(pl,xs,l(xs))
end
plot!(pl,ylimit = [0,w[end]+1])
# xs = 0:0.01:r
# ys = p(xs)
# plot!(xs,ys,linewidth = 2)

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
export create_new_pwq, generate_PiecewiseQuadratic, insert, delete_next
export get_right_endpoint, evalPwq, get_vals, find_minimum, find_minimum_value

"""
    PiecewiseQuadratic{T}

Represents a piecewise-quadratic as a linked list of QuadraticPolynomial{T} and left endpoints.

Note that the first element of the list is sometimes assumed to be a list head with `NaN` polynomial and `left_endpoint=NaN`.

The list/element `x` is the last element in the list when `x.next === x` or `x.next.left_endpoint==Inf`.
"""
mutable struct PiecewiseQuadratic{T}
    p::QuadraticPolynomial{T}
    left_endpoint::T
    next::PiecewiseQuadratic{T}

    PiecewiseQuadratic{T}()  where T = (x=new(); x.left_endpoint=Inf; x.next = x)
    PiecewiseQuadratic{T}(p, left_endpoint, next) where T = new(p, left_endpoint, next)
end

"""
Constructs piece of quadratic polynomial poiting to NULL, i.e.
a rightmost segment, or one that has not been inserted into a piecewise quadratic
"""
function PiecewiseQuadratic(p::QuadraticPolynomial{T}, left_endpoint) where T
    PiecewiseQuadratic{T}(p, left_endpoint, PiecewiseQuadratic{T}())
end
function PiecewiseQuadratic(p::QuadraticPolynomial{T}, left_endpoint, next::PiecewiseQuadratic{T}) where T
    PiecewiseQuadratic{T}(p, left_endpoint, next)
end

islast(λ::PiecewiseQuadratic) = λ.next.left_endpoint == Inf

"""
Creates empty piecewise-quadratic polynomial consisting only of the list head
"""
function create_new_pwq(T)
    return PiecewiseQuadratic(QuadraticPolynomial(T[NaN,NaN,NaN]), T(NaN), PiecewiseQuadratic{T}())
end

"""
Creates piecewise-quadratic polynomial containing one element
"""
function create_new_pwq(p::QuadraticPolynomial{T}) where T
    pwq = PiecewiseQuadratic(p, T(-Inf), PiecewiseQuadratic{T}())
    return PiecewiseQuadratic(QuadraticPolynomial(T[NaN,NaN,NaN]), T(NaN), pwq)
end

function generate_PiecewiseQuadratic(args::Tuple{Vector{T},T}...) where T
    return PiecewiseQuadratic(QuadraticPolynomial(T[NaN,NaN,NaN]), T(NaN), _generate_PiecewiseQuadratic_helper(args...))
end

function _generate_PiecewiseQuadratic_helper(arg)
    return PiecewiseQuadratic(QuadraticPolynomial(arg[1]), arg[2])
end

function _generate_PiecewiseQuadratic_helper(arg, args...)
    return PiecewiseQuadratic(QuadraticPolynomial(arg[1]), arg[2],  _generate_PiecewiseQuadratic_helper(args...))
end

function Base.iterate(pwq::PiecewiseQuadratic{T}, iterstate::PiecewiseQuadratic{T}=pwq.next) where T
    if (iterstate.left_endpoint == Inf)
        return nothing
    end

    return (iterstate, iterstate.next)
end

# For trouble-shooting etc.
function getindex(Λ::PiecewiseQuadratic, n::Integer)
    if n <= 0
        throw(BoundsError(Λ, n))
    end
    # If first element has NaN left_endpoint, assume it to be empty list-head
    λ = isnan(Λ.left_endpoint) ? Λ.next : Λ
    for k=1:n-1
        if islast(λ)
            throw(BoundsError(Λ, n))
        end
        λ = λ.next
    end
    return λ
end

function insert(pwq::PiecewiseQuadratic{T}, p::QuadraticPolynomial, left_endpoint) where T
    pwq.next = PiecewiseQuadratic(p, left_endpoint, pwq.next)
    return pwq.next
end

# Delete the node after pwq and return the node that will follow after
# pwq after the deletion
#OBS This function is unsafe if pwq.next does not exist
function delete_next(pwq::PiecewiseQuadratic{T}) where T
    pwq.next = pwq.next.next
    return pwq.next
end

function get_right_endpoint(λ::PiecewiseQuadratic)
    return λ.next.left_endpoint
end

function length(pwq::PiecewiseQuadratic)
    n = 0
    for x in pwq
        n += 1
    end
    return n
end

function show(io::IO, Λ::PiecewiseQuadratic{T}) where T
    print(io, "PiecewiseQuadratic{$T} with $(length(Λ)) elements:\n")

    for λ in Λ
        if λ.left_endpoint == -Inf
            print(io, "[  -∞ ,")
        else
            @printf(io, "[%3.2f,", λ.left_endpoint)
        end

        if get_right_endpoint(λ) == Inf
            print(io, " ∞  ]")
        else
            @printf(io, " %3.2f]", get_right_endpoint(λ))
        end

        print(io, "\t  :   ")
        show(io, λ.p)
        print(io, "\n")
    end
    return
end

function (Λ::PiecewiseQuadratic)(x::Number)
    for λ in Λ
        if λ.left_endpoint <= x < get_right_endpoint(λ)
            return λ.p(x)
        end
    end
    return NaN
end

function (Λ::PiecewiseQuadratic{T1})(x::AbstractArray) where T1
    y = fill(zero(promote_type(T1,eltype(x))), length(x))
    for λ in Λ
        inds = λ.left_endpoint .<= x .< get_right_endpoint(λ)
        y[inds] .= λ.p.(x[inds])
    end
    return y
end

"""
    x, y, x_all, y_all = get_vals(Λ::PiecewiseQuadratic)
Evaluate `Λ` so that `y .= Λ.(x)` for a set of gridpoints `x`.
Also returns lists x_all, y_all so that `y_all[i] .= Λ[i].(x_all[i])`
where `Λ[i]` is the `i`th segment of the piecewise quadratic.

The gridpoints in x_all[i] cointains points in the interval `[xmin[i]-1,xmax[i]+1]` where
`xmin[i],xmax[i]` are the extrema of the interval for which `Λ[i]` defines `Λ`.

The gridpoints at the ends default to `xmin[1] := xmax[1] - 1` and `xmax[end] := xmin[end] + 1`.
"""
function get_vals(Λ::PiecewiseQuadratic)
    x = Float64[]
    y = Float64[]
    x_all = Array{Float64,1}[]
    y_all = Array{Float64,1}[]
    for λ in Λ
        left_endpoint = λ.left_endpoint
        right_endpoint = get_right_endpoint(λ)

        if left_endpoint != -Inf || right_endpoint != Inf
            if left_endpoint == -Inf
                left_endpoint = right_endpoint - 2.0
            end

            if right_endpoint == Inf
                right_endpoint = left_endpoint + 2.0
            end
        else
            left_endpoint = -10.
            right_endpoint = 10.
        end

        y_grid = range(left_endpoint, stop=right_endpoint, length=50)
        vals =  λ.p.(y_grid)
        append!(x, y_grid)
        append!(y, vals)
        push!(x_all, y_grid)
        push!(y_all, vals)
    end
    return x, y, x_all, y_all
end

#
# function plot_pwq(Λ::PiecewiseQuadratic)
#     x, y, x_all, y_all = get_vals(Λ)
#     p = plot(x,y, l =(3,:red), lab="minimum", size=(1200,800))
#     for i in eachindex(x_all)
#         plot!(p, x_all[i][[1,end]], y_all[i][[1,end]], l = 0, m=(1,:cross,:orange))
#         plot!(p, x_all[i], y_all[i], l=(1,:blue), lab="")
#     end
#     p
#     return p
# end

"""
    p_opt, x_opt, f_opt = find_minimum(Λ::PiecewiseQuadratic)

Find the point `x_opt = argmin Λ(x)` and return `f_opt = Λ(x_opt)` and the quadratic
function `p_opt` which defines `Λ` around the point `x_opt`.
"""
function find_minimum(Λ::PiecewiseQuadratic)
    # May assume that the minimum is at the staionary point of the polynomial
    # TODO: True?

    f_opt = Inf
    p_opt = typeof(Λ.p)()
    x_opt = NaN

    for λ in Λ
        x, f = find_minimum(λ.p)
        if f < f_opt
            f_opt = f
            x_opt = x
            p_opt = λ.p
        end
    end
    return p_opt, x_opt, f_opt
end

function find_minimum_value(Λ::PiecewiseQuadratic)
    f_opt = Inf

    for λ in Λ
        _, f = find_minimum(λ.p)
        if f < f_opt
            f_opt = f
        end
    end
    if isnan(f_opt)
        println("find_minimum_value failed to find a minimum, no minimum exists")
    end
    return f_opt
end


text:
import Base.reverse
reverse(n::Number) = n

"""
    ptrace(ρ,dims,sys)

Compute the partial trace of a matrix ρ with subsystems dimensions specified in
the vector sys_dims (e.g.: [dimA,dimB...]) by tracing out the subsystems given
in the vector sys (e.g.: [2,3]).

Function adapted from a MATLAB function by Toby Cubitt available here:
http://www.dr-qubit.org/Matlab_code.html and licensed under GPL2.

# Arguments
* `ρ::Matrix`: matrix
* `dims::Vector{Int}`: vector specifying the dimensions of the subsystems
  that the full matrix contains. A qubit and resonator state (with the Fock
  space truncated to 6) would have sys_dim = [2 6].
* `sys::Union{Vector{Int},Int}`: vector (or single value) specifying the indices
  of the subsystems to trace out.

# Examples
```jldoctest
julia> ρ = [0.5 0 0 0.5; 0 0 0 0; 0 0 0 0; 0.5 0 0 0.5]
julia> dims = [2,2] # two qubits
julia> sys = 1 # trace out the first qubit
julia> ptrace(ρ,dims,sys)
2×2 Array{Float64,2}:
 0.5  0.0
 0.0  0.5
```
"""
function ptrace(ρ::Matrix,dims::Vector{Int},sys::Union{Vector{Int},Int})
    # Make sure the input arguments make sense
    if any(sys.>length(dims))
        error("You have given a system index too high for the dimensions specified in `dims`")
    end
    if prod(dims) != size(ρ,1)
        error("The subsystem dimensions given in `dims` are not consistent with those of the matrix")
    end
    # First, calculate systems, dimensions, etc.
    nsub = length(dims) # Number of subsytems in ρ
    rdims = reverse(dims) # Dimensions of subsystems in reverse order
    keep = setdiff(1:nsub,sys) # Subsystems to keep
    dimstrace = prod(dims[sys]) # Total dim of subsystems to trace out
    dimskeep = prod(dims[keep]) # Total dim of subsystems to keep
    # Reshape density matrix into tensor with one row and one column index
    # for each subsystem, permute traced subsystem indices to the end,
    # reshape again so that first two indices are row and column
    # multi-indices for kept subsystems and third index is a flattened index
    # for traced subsystems, then sum third index over "diagonal" entries
    perm = nsub+1-[reverse(keep);reverse(keep)-nsub;sys;sys-nsub]
    x = reshape(permutedims(reshape(ρ,rdims...,rdims...),perm),dimskeep,dimskeep,dimstrace^2)
    return sum(x[:,:,1:dimstrace+1:dimstrace^2],3)[:,:,1]
end

function calc_qubittimeevo(ρ,U_array)
    # Time evolve the system density matrix by applying the time evolution
    # operator and measure the qubit state at each step. This is done by
    # tracing out the resonator system and keeping the excited states entry.
    N = div(size(ρ,1),2)
    map(U_array) do U
        real(ptrace(U*ρ*U',[2,N],2)[2,2])
    end
end

function calc_qubittimeevo_rel(ρ,U_rel_array,steps,decimation=10)
    # Time evolve the system density matrix by applying the relative time
    # evolution in succession and measure the qubit state at each step.
    N = div(size(ρ,1),2)
    c = length(U_rel_array)
    ρ_evo = deepcopy(ρ)
    #nvec = collect(0.0:N-1)'

    e_prob = Array{Float64}(ceil(Int,steps/decimation))
    #n_expect = Array{Float64}(ceil(Int,steps/decimation))

    j = 1
    for i = 1:steps
        # Note the use of the modulo operation to apply U in cycles
        ρ_evo = U_rel_array[mod1(i,c)]*ρ_evo*U_rel_array[mod1(i,c)]'
        if (i-1)%decimation==0
            e_prob[j] = real(ptrace(ρ_evo,[2,N],2)[2,2])
            #n_expect[j] = (nvec'*real(diag(ptrace(ρ_evo,[2,N],1))))[1]
            j += 1
        end
    end
    return e_prob
end

function calc_photonnumbers(N,g,times,e_prob)
    M = length(times)
    B = Array{Float64}(M,N)
    @compat B[:,1] .= 1.0
    for n = 2:N
        f = 2*sqrt(n-1)*g
        @compat B[:,n] .= cos.(f.*times)
    end
    return B\(1-2*e_prob)
end

function calc_densitymatrix_resonator(ρ,g,times,U_array)
    # Solving for the density matrix from photon number distributions
    N = div(size(ρ,1),2)
    # We need 2*N displacements to get enough linearly independent
    # equations. We sample the alpha values from a circle of unit radius.
    α = [1.0*complex(cos(θ),sin(θ)) for θ in 0:π/N:2π-π/N]
    # Do the math
    res = pmap(α) do α
        # For each angle/displaced state, we first find the photon number distribution
        # Generating the displacement operator for angle α
        Dα = qeye(2) ⊗ gen_displacementop(α,N)
        # Solving for the displaced state photon number distribution
        e_prob = calc_qubittimeevo(Dα'*ρ*Dα,U_array)
        photon_numbers = calc_photonnumbers(N,g,times,e_prob)
        # We then calculate the displacement operator matrix elements
        M = Array{Complex128}(N,N^2)
        for n=1:N, i=1:N, j=1:N
            # Generating the displacement operator matrix elements
            M[n,(i-1)*N+j] = conj(Dα[i,n]) * Dα[j,n]
        end
        photon_numbers,M
    end
    # Reshaping the result into arrays for solving the system of equations
    photon_numbers = vcat([res[1] for res in res]...)
    M = vcat([res[2] for res in res]...)
    # Solving the system of equations (very slighly overdetermined, there will
    # be some least square fitting). The system to solve is:
    # photon_numbers = M * ρ
    ρ_res = M\photon_numbers # Vector form of the density matrix (ie the elements are 0 0, 0 1, 0 2, ..., N N-1, N N)
    ρ_res = transpose(reshape(ρ_res,N,N)) # Reshaping the vector into an actual matrix
end

function calc_wignerfunction_resonator(ρ,g,times,U_array,w_samples;maxdisp=MAXDISP)
    N = div(size(ρ,1),2)
    # Solving for the Wigner function from photon number distributions
    parity_matrix = [(-1.0)^i for i = 0:N-1]
    # Choose where to plot the function
    disp = linspace(-maxdisp,maxdisp,w_samples)
    # Iterate over all displacements
    W = pmap(product(disp,disp)) do id
        # Generating the displacement operator for that matrix element
        Dα = qeye(2) ⊗ gen_displacementop(complex(id[1],id[2]),N)
        # For each displaced state, we first find the photon number distribution
        e_prob = calc_qubittimeevo(Dα'*ρ*Dα,U_array)
        photon_numbers = calc_photonnumbers(N,g,times,e_prob)
        # We then calculate the entry.
        2/π * sum(parity_matrix .* photon_numbers)
    end
    return convert(Array{Float64,2},reshape(W,w_samples,w_samples))
end

function calc_wignerfunction_resonator(ρ_res,w_samples;maxdisp=MAXDISP)
    N = size(ρ_res,1)
    # Solving for the Wigner function from the density matrix
    parity_matrix = [(-1.0)^i for i = 0:N-1]
    # Choose where to plot the function
    disp = linspace(-maxdisp,maxdisp,w_samples)
    # Iterate over all displacements
    map(product(disp,disp)) do id
        # Generating the displacement operator for that matrix element
        Dα = gen_displacementop(complex(id[1],id[2]),N)
        # We then calculate the entry.
        2/π * sum(parity_matrix .* real(diag(Dα'*ρ_res*Dα)))
    end
end

function cal_densitymatrix_qubit(xup,yup,zup)
    # A qubit in state |ψ> is defined by two values, a & b, where |ψ> =
    # a|↑>+b|↓>. By convention, a is real and positive, and b is in general
    # complex. From the probability of measuring |↑> (that is, measuring "up" in
    # the z basis), we can deduce a, and the magnitude of b
    a = sqrt(zup)
    magb = sqrt(1-zup)
    # The probability of measuring "up" in the x basis is equal to
    # |1/sqrt(2)*(<↑|+<↓|)|ψ>|^2, which is equal to 1/2(a^2+b^2)+a*mag(b)cos(ϕ)
    cosϕ = ( xup - (a^2 + magb^2)/2 ) / (a * magb)
    # The probability of measuring "up" in the y basis is equal to
    # |1/sqrt(2)*(<↑|-i<↓|)|ψ>|^2, which is equal to 1/2(a^2+b^2)+a*mag(b)sin(ϕ)
    sinϕ = ( yup - (a^2 + magb^2)/2 ) / (a * magb)
    b = magb * complex(cosϕ,sinϕ)
    return a, b
end

function ptrace_old(rho,dimOfSubsystemToKeep)
    # First check if the density matrix is Hermitian
    if maximum(abs(rho'-rho))>8E-15
        error("This density matrix is not Hermitian!")
    end
    dimOfSubsystemToTraceOut = div(size(rho)[1] , dimOfSubsystemToKeep)
    matrix = zeros(Complex128,dimOfSubsystemToKeep,dimOfSubsystemToKeep)
    for i=1:dimOfSubsystemToTraceOut
        vector = zeros(Complex128,dimOfSubsystemToTraceOut)
        vector[i] = (1.0+1.0im)/sqrt(2)
        fullVector = qeye(dimOfSubsystemToKeep) ⊗ vector
        matrix = matrix + fullVector' * rho * fullVector
    end
    return matrix
end


text:
#
# xforms.jl --
#
# Implementation of affine transforms which are notably useful for coordinate
# transforms.
#
#------------------------------------------------------------------------------
#
# This file if part of the TwoDimensional Julia package licensed under the MIT
# license (https://github.com/emmt/TwoDimensional.jl).
#
# Copyright (C) 2016-2019, Éric Thiébaut.
#

module AffineTransforms

export
    AffineTransform,
    compose,
    intercept,
    jacobian,
    rotate,
    scale,
    translate

using ..TwoDimensional: AbstractPoint, Point

# Imports for extension.
import Base: +, -, *, ∘, /, \, inv, eltype
import Base: Float16, Float32, Float64
import Base.MPFR: BigFloat
import LinearAlgebra: ⋅, det

"""
# Affine 2D Transforms

An affine 2D transform `A` is defined by 6 real coefficients, `Axx`, `Axy`,
`Ax`, `Ayx`, `Ayy` and `Ay`.  Such a transform maps `(x,y)` as `(xp,yp)` given
by:

```julia
xp = Axx*x + Axy*y + Ax
yp = Ayx*x + Ayy*y + Ay
```

The immutable type `AffineTransform` is used to store an affine 2D transform
`A`, it can be created by:

```julia
I = AffineTransform{T}() # yields the identity with type T
A = AffineTransform{T}(Axx, Axy, Ax, Ayx, Ayy, Ay)
```

The parameter `T` above is used to specify the floating-point type for the
coefficients; if omitted it is guessed from the types of the coefficients.

## Operations with affine 2D transforms

Many operations are available to manage or apply affine transforms:

```julia
(xp, yp) = A(x,y)       # apply affine transform A to coordinates (x,y)
(xp, yp) = A*(x,y)      # idem
(xp, yp) = A(v)         # idem, with v = (x,y)
(xp, yp) = A*v          # idem

A(Point(x,y)) -> Point(xp, yp)
A*Point(x,y)  -> Point(xp, yp)

C = compose(A, B, ...)  # compose 2 (or more) transforms, C = apply B then A
C = A∘B                 # idem
C = A*B                 # idem
C = A⋅B                 # idem

B = translate(x, y, A)  # B = apply A then translate by (x,y)
B = translate(v, A)     # idem with v = (x,y)
B = v + A               # idem

B = translate(A, x, y)  # B = translate by (x,y) then apply A
B = translate(A, v)     # idem with v = (x,y)
B = A + v               # idem

B = rotate(θ, A)   # B = apply A then rotate by angle θ
C = rotate(A, θ)   # C = rotate by angle θ then apply A

B = scale(ρ, A)    # B = apply A then scale by ρ
B = ρ*A            # idem
C = scale(A, ρ)    # C = scale by ρ then apply A
C = A*ρ            # idem

B = inv(A)         # reciprocal coordinate transform
C = A/B            # right division, same as: C = compose(A, inv(B))
C = A\\B            # left division, same as: C = compose(inv(A), B)
```

"`∘`" and "`⋅`" can be typed by `\\circ<tab>` and `\\cdot<tab>`.

## Type conversion

As a general rule, the floating-point type `T` of an `AffineTransform{T}` is
imposed for all operations and for the result.  The floating-point type of the
composition of several coordinate transforms is the promoted type of the
transforms which have been composed.

Calling `eltype(A)` yields floating-point type of the coefficients of the 2D
affine transform `A`.  To convert the floating-point type of the coefficients
of `A` to be `T`, do one of:

```julia
B = T.(A)
B = AffineTransform{T}(A)
B = convert(AffineTransform{T}, A)
```

"""
struct AffineTransform{T<:AbstractFloat} <: Function
    xx::T
    xy::T
    x ::T
    yx::T
    yy::T
    y ::T
    AffineTransform{T}() where T = new{T}(1,0,0, 0,1,0)
    function AffineTransform{T}(Axx::Real, Axy::Real, Ax::Real,
                                Ayx::Real, Ayy::Real, Ay::Real) where T
        new{T}(Axx,Axy,Ax, Ayx,Ayy,Ay)
    end
end

eltype(::AffineTransform{T}) where {T} = T

AffineTransform() = AffineTransform{Float64}()
function AffineTransform(Axx::Txx, Axy::Txy, Ax::Tx,
                         Ayx::Tyx, Ayy::Tyy, Ay::Ty
                         ) where {Txx<:Real,Txy<:Real,Tx<:Real,
                                  Tyx<:Real,Tyy<:Real,Ty<:Real}
    T = float(promote_type(Txx, Txy, Tx, Tyx, Tyy, Ty))
    AffineTransform{T}(Axx,Axy,Ax, Ayx,Ayy,Ay)
end

function AffineTransform(Axx::T, Axy::T, Ax::T,
                         Ayx::T, Ayy::T, Ay::T) where {T<:AbstractFloat}
    AffineTransform{T}(Axx,Axy,Ax, Ayx,Ayy,Ay)
end

AffineTransform(A::AffineTransform) = A
AffineTransform{T}(A::AffineTransform{T}) where {T<:AbstractFloat} = A
AffineTransform{T}(A::AffineTransform) where {T<:AbstractFloat} =
    AffineTransform{T}(A.xx, A.xy, A.x,
                       A.yx, A.yy, A.y)

# Allow for `T.(obj)` to work with `T` a floating-point type.
Broadcast.broadcasted(::Type{T}, A::AffineTransform) where {T<:AbstractFloat} =
    AffineTransform{T}(A)

Base.convert(::Type{T}, A::AffineTransform) where {T<:AffineTransform} = T(A)
Base.promote_type(::Type{AffineTransform{T}}, ::Type{AffineTransform{U}}) where {T,U} =
    AffineTransform{promote_type(T,U)}
Base.promote_type(::Type{AffineTransform{T}}, ::Type{AffineTransform{T}}) where {T} =
    AffineTransform{T}

@deprecate(Float16(A::AffineTransform),  AffineTransform{Float16}(A))
@deprecate(Float32(A::AffineTransform),  AffineTransform{Float32}(A))
@deprecate(Float64(A::AffineTransform),  AffineTransform{Float64}(A))
@deprecate(BigFloat(A::AffineTransform), AffineTransform{BigFloat}(A))

#------------------------------------------------------------------------------
# apply the transform to some coordinates:

(A::AffineTransform{T})(x::T, y::T) where {T<:AbstractFloat} =
    (A.xx*x + A.xy*y + A.x,
     A.yx*x + A.yy*y + A.y)

(A::AffineTransform{T})(x::Real, y::Real) where {T<:AbstractFloat} =
    A(convert(T, x), convert(T, y))

(A::AffineTransform)(v::Tuple{Real,Real}) = A(v[1], v[2])

(A::AffineTransform)(v::Point) = Point(A(v.x, v.y)...)

#------------------------------------------------------------------------------
# Combine a translation with an affine transform.

"""
### Translating an affine transform

Affine transforms can be letf- or right-translated.

```julia
translate(x, y, A)
```
or
```julia
translate((x,y), A)
```

yield an affine transform which translate the output of affine transform `A` by
offsets `x` and `y`.

```julia
translate(A, x, y)
```
or
```julia
translate(A, (x,y))
```

yield an affine transform which translate the input of affine transform `A` by
offsets `x` and `y`.

The same results can be obtained with the `+` operator:

```julia
B = (x,y) + A    # same as: B = translate((x,y), A)
B = A + (x,y)    # same as: B = translate(A, (x,y))
```

See also: [`AffineTransform`](@ref), [`rotate`](@ref), [`scale`](@ref).

""" translate

# Left-translating results in translating the output of the transform.
translate(x::T, y::T, A::AffineTransform{T}) where {T<:AbstractFloat} =
    AffineTransform{T}(A.xx, A.xy, A.x + x,
                       A.yx, A.yy, A.y + y)

translate(x::Real, y::Real, A::AffineTransform{T}) where {T<:AbstractFloat} =
    translate(convert(T, x), convert(T, y), A)

translate(v::Tuple{Real,Real}, A::AffineTransform) =
    translate(v[1], v[2], A)

translate(v::AbstractPoint, A::AffineTransform) =
    translate(v.x, v.y, A)

# Right-translating results in translating the input of the transform.
translate(A::AffineTransform{T}, x::T, y::T) where {T<:AbstractFloat} =
    AffineTransform{T}(A.xx, A.xy, A.xx*x + A.xy*y + A.x,
                       A.yx, A.yy, A.yx*x + A.yy*y + A.y)

translate(A::AffineTransform{T}, x::Real, y::Real) where {T<:AbstractFloat} =
    translate(A, convert(T, x), convert(T, y))

translate(A::AffineTransform, v::Tuple{Real,Real}) =
    translate(A, v[1], v[2])

translate(A::AffineTransform, v::AbstractPoint) =
    translate(A, v.x, v.y)

#------------------------------------------------------------------------------
"""
### Scaling an affine transform

There are two ways to combine a scaling by a factor `ρ` with an affine
transform `A`.  Left-scaling as in:

```julia
B = scale(ρ, A)
```

results in scaling the output of the transform; while right-scaling as in:

```julia
C = scale(A, ρ)
```

results in scaling the input of the transform.  The above examples yield
transforms which behave as:

```julia
B(v) = ρ.*A(v)
C(v) = A(ρ.*v)
```

where `v` is any 2-element tuple.

The same results can be obtained with the `*` operator:

```julia
B = ρ*A    # same as: B = scale(ρ, A)
C = A*ρ    # same as: B = scale(A, ρ)
```

See also: [`AffineTransform`](@ref), [`rotate`](@ref), [`translate`](@ref).

"""
scale(ρ::T, A::AffineTransform{T}) where {T<:AbstractFloat} =
    AffineTransform{T}(ρ*A.xx, ρ*A.xy, ρ*A.x,
                       ρ*A.yx, ρ*A.yy, ρ*A.y)

scale(A::AffineTransform{T}, ρ::T) where {T<:AbstractFloat} =
    AffineTransform{T}(ρ*A.xx, ρ*A.xy, A.x,
                       ρ*A.yx, ρ*A.yy, A.y)

#------------------------------------------------------------------------------
"""
### Rotating an affine transform

There are two ways to combine a rotation by angle `θ` (in radians
counterclockwise) with an affine transform `A`.  Left-rotating as in:

```julia
B = rotate(θ, A)
```

results in rotating the output of the transform; while right-rotating as in:

```julia
C = rotate(A, θ)
```

results in rotating the input of the transform.  The above examples are
similar to:

```julia
B = R∘A
C = A∘R
```

where `R` implements rotation by angle `θ` around `(0,0)`.

See also: [`AffineTransform`](@ref), [`scale`](@ref), [`translate`](@ref).

"""
function rotate(θ::T, A::AffineTransform{T}) where {T<:AbstractFloat}
    cs = cos(θ)
    sn = sin(θ)
    return AffineTransform{T}(cs*A.xx - sn*A.yx,
                              cs*A.xy - sn*A.yy,
                              cs*A.x  - sn*A.y,
                              cs*A.yx + sn*A.xx,
                              cs*A.yy + sn*A.xy,
                              cs*A.y  + sn*A.x)
end

function rotate(A::AffineTransform{T}, θ::T) where {T<:AbstractFloat}
    cs = cos(θ)
    sn = sin(θ)
    return AffineTransform{T}(A.xx*cs + A.xy*sn,
                              A.xy*cs - A.xx*sn,
                              A.x,
                              A.yx*cs + A.yy*sn,
                              A.yy*cs - A.yx*sn,
                              A.y)
end

# Make sure the floating-point type of an affine transform is preserved.
for func in (:scale, :rotate)
    @eval begin
        $func(α::Real, A::AffineTransform{T}) where {T<:AbstractFloat} =
            $func(convert(T, α), A)
        $func(A::AffineTransform{T}, α::Real) where {T<:AbstractFloat} =
            $func(A, convert(T, α))
    end
end

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

"""
`det(A)` returns the determinant of the linear part of the affine
transform `A`.
"""
det(A::AffineTransform) = A.xx*A.yy - A.xy*A.yx

"""
`jacobian(A)` returns the Jacobian of the affine transform `A`, that is the
absolute value of the determinant of its linear part.
"""
jacobian(A::AffineTransform) = abs(det(A))

"""
`inv(A)` returns the inverse of the affine transform `A`.
"""
function inv(A::AffineTransform{T}) where {T}
    Δ = det(A)
    Δ == zero(T) && error("transformation is not invertible")
    α = one(T)/Δ
    Txx =  α*A.yy
    Txy = α*A.xy
    Tyx = α*A.yx
    Tyy =  α*A.xx
    return AffineTransform{T}(+Txx, -Txy, Txy*A.y - Txx*A.x,
                              -Tyx, +Tyy, Tyx*A.x - Tyy*A.y)
end

"""

`compose(A,B)` yields the affine transform which combines the two affine
transforms `A` and `B`, that is the affine transform which applies `B` and then
`A`.  Composition is accessible via: `A∘B`, `A*B` or `A⋅B` ("`∘`" and "`⋅`" can
be typed by `\\circ<tab>` and `\\cdot<tab>`).

It is possible to compose more than two affine transforms.  For instance,
`compose(A,B,C)` yields the affine transform which applies `C` then `B`, then
`A`.

"""
compose(A::AffineTransform) = A

compose(A::AffineTransform{T}, B::AffineTransform{T}) where {T<:AbstractFloat} =
    AffineTransform{T}(A.xx*B.xx + A.xy*B.yx,
                       A.xx*B.xy + A.xy*B.yy,
                       A.xx*B.x  + A.xy*B.y + A.x,
                       A.yx*B.xx + A.yy*B.yx,
                       A.yx*B.xy + A.yy*B.yy,
                       A.yx*B.x  + A.yy*B.y + A.y)

compose(A::AffineTransform, B::AffineTransform) =
    compose(promote(A, B)...)

compose(A::AffineTransform, B::AffineTransform, C::AffineTransform, args::AffineTransform...) =
    compose(compose(A, B), C, args...)

"""

`rightdivide(A,B)` yields `A/B`, the right division of the affine
transform `A` by the affine transform `B`.

"""
rightdivide(A::AffineTransform{T}, B::AffineTransform{T}) where {T<:AbstractFloat} = begin
    Δ = det(B)
    Δ == zero(T) && error("right operand is not invertible")
    α = one(T)/Δ
    Rxx = α*(A.xx*B.yy - A.xy*B.yx)
    Rxy = α*(A.xy*B.xx - A.xx*B.xy)
    Ryx = α*(A.yx*B.yy - A.yy*B.yx)
    Ryy = α*(A.yy*B.xx - A.yx*B.xy)
    AffineTransform{T}(Rxx, Rxy, A.x - (Rxx*B.x + Rxy*B.y),
                       Ryx, Ryy, A.y - (Ryx*B.y + Ryy*B.y))

end

rightdivide(A::AffineTransform, B::AffineTransform) =
    rightdivide(promote(A, B)...)

"""
`leftdivide(A,B)` yields `A\\B`, the left division of the affine
transform `A` by the affine transform `B`.
"""
leftdivide(A::AffineTransform{T}, B::AffineTransform{T}) where {T<:AbstractFloat} = begin
    Δ = det(B)
    Δ == zero(T) && error("left operand is not invertible")
    α = one(T)/Δ
    Txx = α*A.yy
    Txy = α*A.xy
    Tyx = α*A.yx
    Tyy = α*A.xx
    Tx = B.x - A.x
    Ty = B.y - A.y
    AffineTransform{T}(Txx*B.xx - Txy*B.yx,
                       Txx*B.xy - Txy*B.yy,
                       Txx*Tx   - Txy*Ty,
                       Tyy*B.yx - Tyx*B.xx,
                       Tyy*B.yy - Tyx*B.xy,
                       Tyy*Ty   - Tyx*Tx)
end

leftdivide(A::AffineTransform, B::AffineTransform) =
    leftdivide(promote(A, B)...)

"""

`intercept(A)` returns the tuple `(x,y)` such that `A(x,y) = (0,0)`.

"""
function intercept(A::AffineTransform{T}) where {T}
    Δ = det(A)
    Δ == zero(T) && error("transformation is not invertible")
    α = one(T)/Δ
    return (α*(A.xy*A.y - A.yy*A.x), α*(A.yx*A.x - A.xx*A.y))
end

intercept(T::Type{<:Point}, A::AffineTransform) = T(intercept(A)...)

+(v::Union{AbstractPoint,Tuple{Real,Real}}, A::AffineTransform) =
    translate(v, A)

+(A::AffineTransform, v::Union{AbstractPoint,Tuple{Real,Real}}) =
    translate(A, v)

-(A::AffineTransform, v::Tuple{Real,Real}) = A + (-v[1], -v[2])
-(A::AffineTransform, v::AbstractPoint) = A + (-v.x, -v.y)

for op in (:∘, :*, :⋅)
    @eval begin
        $op(A::AffineTransform, B::AffineTransform) = compose(A, B)
    end
end

*(A::AffineTransform, v::Union{Point,Tuple{Real,Real}}) = A(v)

*(ρ::Real, A::AffineTransform) = scale(ρ, A)

*(A::AffineTransform, ρ::Real) = scale(A, ρ)

\(A::AffineTransform, B::AffineTransform) = leftdivide(A, B)

/(A::AffineTransform, B::AffineTransform) = rightdivide(A, B)

Base.show(io::IO, ::MIME"text/plain", A::AffineTransform) =
    print(io, typeof(A),
          "(",   A.xx, ",", A.xy, ",", A.x,
          ",  ", A.yx, ",", A.yy, ",", A.y, ")")

Base.show(io::IO, A::AffineTransform) = show(io, MIME"text/plain"(), A)

Base.print(io::IOBuffer, A::AffineTransform) = show(io, A)

end # module


text:
export GradientBoost

function GradientBoost(prob::Problem{L, B, CS},
                       rate::Real,
                       subsolver::AbstractNextBasisSolver=defaultsolver(prob);
                       ad=false,
                       maxiter::Int=300) where {L, F, B<:AbBasis{F}, CS}
    ens = GBEnsemble(prob)
    enses = BBEnsemble{F,B}[]
    iter = 0
    currloss = trainerr(ens)
    # llog(txt="") = printfmtln("{:3d}) {:.6f} {:.2f}({}) {}", iter, currloss, sum(abs, ens.β), countnz(ens.β), txt)
    llog(t="") = nothing
    llog()

    for iter in 1:maxiter
        prevloss = currloss
        objtv, grad = lossngrad(prob, ens.fit)
        ω, ρ, score = getnextω(prob, subsolver, grad)
        βi = linesearch(prob, ens.fit, ω)
        βi *= rate
        push!(ens, ω, βi)
        currloss = trainerr(ens)

        if ad
            ens, currloss = backpropωs(prob, ens)
            push!(enses, BBEnsemble(ens))
            llog("ad {ω}")
        end

        abs(1-currloss/prevloss) < 1e-5 && (println("loss converged"); break)
        abs(ens.β[end]) < abs(ens.β[1])*1e-4 && (println("β converged"); break)
        llog(format("{:.2f}%", 100*abs(ens.β[end]) / abs(ens.β[1])))
    end

    llog("\nEND")
    ad ? enses : ens
end


text:
using ITensors
using Test
using Compat

@testset "IndexSet" begin
  idim = 2
  jdim = 3
  kdim = 4
  ldim = 5
  i = Index(idim, "i")
  j = Index(jdim, "j")
  k = Index(kdim, "k")
  l = Index(ldim, "l")
  @testset "show" begin
    indset = IndexSet(i, j, k)
    @test length(sprint(show, indset)) > 1
  end
  @testset "Basic constructors" begin
    I = IndexSet(i, j, k)
    @test IndexSet(I) === I
    @test l ∈ IndexSet(I..., l)
    @test l ∈ IndexSet(l, I...)
    @test length(IndexSet(i, j)) == 2
    # construct with function
    ind_list = [i, j, k]
    I = IndexSet(ii -> ind_list[ii], 3)
    @test i ∈ I
    @test j ∈ I
    @test k ∈ I
    I = IndexSet(ii -> ind_list[ii], Order(3))
    @test i ∈ I
    @test j ∈ I
    @test k ∈ I
  end
  @testset "length of IndexSet and friends" begin
    @test length(IndexSet(i, j)) == 2
    @test size(IndexSet(i, j)) == (length(IndexSet(i, j)),)
  end
  @testset "Convert to Index" begin
    @test Index(IndexSet(i)) === i
    @test_throws BoundsError Index(IndexSet(i, j))
  end
  @testset "Index dimensions" begin
    I = IndexSet(i, j, k)
    @test dim(I) == idim * jdim * kdim
    @test dims(I) == [idim, jdim, kdim]
    @test dim(I, 1) == idim
    @test dim(I, 2) == jdim
    @test dim(I, 3) == kdim

    @test maxdim(I) == max(idim, jdim, kdim)
  end

  @testset "Set operations" begin
    I1 = @inferred(IndexSet(i, j, k))
    I2 = @inferred(IndexSet(k, l))
    I3 = @inferred(IndexSet(j, l))
    @test I1 isa Vector{Index{Int}}
    @test @inferred(hassameinds(I1, (k, j, i)))
    @test @inferred(Nothing, getfirst(setdiff(I1, I2, I3))) == i
    @test isnothing(@inferred(Nothing, getfirst(setdiff(I1, IndexSet(k, j, i)))))
    @test @inferred(setdiff(I1, I2)) == [i, j]
    @test hassameinds(@inferred(setdiff(I1, I2)), IndexSet(i, j))
    @test hassameinds(@inferred(setdiff(I1, I2)), (j, i))
    @test I1 ∩ I2 == [k]
    @test hassameinds(I1 ∩ I2, IndexSet(k))
    @test @inferred(Nothing, getfirst(intersect(I1, I2))) == k
    @test isnothing(@inferred(Nothing, getfirst(intersect(I1, IndexSet(l)))))
    @test @inferred(intersect(I1, IndexSet(j, l))) == [j]
    @test hassameinds(@inferred(intersect(I1, IndexSet(j, l))), IndexSet(j))
    @test @inferred(Nothing, getfirst(intersect(I1, IndexSet(j, l)))) == j
    @test @inferred(intersect(I1, IndexSet(j, k))) == [j, k]
    @test hassameinds(@inferred(intersect(I1, (j, k))), IndexSet(j, k))
    @test hassameinds(@inferred(intersect(I1, (j, k, l))), (j, k))
    @test @inferred(filterinds(I1, "i")) == IndexSet(i)
    @test @inferred(filterinds(I1; tags="i")) == IndexSet(i)
    @test @inferred(filterinds(I1; inds=j)) == IndexSet(j)
    @test @inferred(filterinds(I1; tags="i", inds=j)) == IndexSet()
    @test @inferred(filterinds(I1; plev=1, inds=j)) == IndexSet()
    @test @inferred(filterinds(I1; plev=0, inds=k)) == IndexSet(k)
    @test @inferred(filterinds(I1; plev=0)) == IndexSet(i, j, k)
    @test @inferred(filterinds(I1; inds=l)) == IndexSet()
    @test @inferred(hassameinds(filter(I1, "i"), IndexSet(i)))
    @test @inferred(Nothing, getfirst(I1, "j")) == j
    @test isnothing(@inferred(Nothing, getfirst(I1, "l")))
    @test @inferred(Nothing, findfirst(I1, i)) == 1
    @test @inferred(Nothing, findfirst(I1, j)) == 2
    @test @inferred(Nothing, findfirst(I1, k)) == 3
    @test isnothing(@inferred(Nothing, findfirst(I1, Index(2))))
  end

  @testset "Set operations with Order" begin
    i, j, k, l = Index.(2, ("i", "j", "k", "l"))

    Iij = IndexSet(i, j)
    Ijl = IndexSet(j, l)
    Ikl = IndexSet(k, l)
    Iijk = IndexSet(i, j, k)

    #
    # setdiff
    # intersect
    # symdiff
    # union
    # filter
    #

    #
    # setdiff
    #

    @test @inferred(setdiff(Iijk, Ikl)) == [i, j]

    @test @inferred(setdiff(Iij, Iijk)) == Index{Int}[]

    @test @inferred(uniqueinds(Iijk, Ikl; tags="i")) == [i]

    @test @inferred(uniqueinds(Iijk, Ikl; tags=not("i"))) == [j]

    @test @inferred(setdiff(Iijk, Ijl, Ikl)) == [i]

    #
    # intersect
    #

    @test @inferred(intersect(Iijk, Ikl)) == [k]

    @test @inferred(intersect(Iijk, Iij)) == [i, j]

    @test @inferred(commoninds(Iijk, Iij; tags="i")) == [i]

    #
    # symdiff
    #

    @test @inferred(symdiff(Iijk, Ikl)) == [i, j, l]

    @test @inferred(symdiff(Iijk, Iij)) == [k]

    #
    # union
    #

    @test @inferred(union(Iijk, Ikl)) == [i, j, k, l]

    @test @inferred(union(Iijk, Iij)) == [i, j, k]
  end

  @testset "intersect index ordering" begin
    I = IndexSet(i, k, j)
    J = IndexSet(j, l, i)
    # Test that intersect respects the ordering
    # of the indices in the first IndexSet
    @test @inferred(hassameinds(intersect(I, J), IndexSet(i, j)))
    @test @inferred(hassameinds(intersect(J, I), IndexSet(j, i)))
  end
  @testset "adjoint" begin
    I = IndexSet(i, k, j)
    @test adjoint(I) == IndexSet(i', k', j')
  end
  @testset "mapprime" begin
    I = IndexSet(i', k'', j)
    @test mapprime(I, 1, 5) == IndexSet(i^5, k'', j)
    @test mapprime(I, 2, 0) == IndexSet(i', k, j)

    J = IndexSet(i, j, k')
    @test mapprime(J, 0, 2) == IndexSet(i'', j'', k')

    J = mapprime(J, 1, 5)
    @test J == IndexSet(i, j, k^5)
  end
  @testset "strides" begin
    I = IndexSet(i, j)
    @test NDTensors.dim_to_strides(I) == (1, idim)
    @test NDTensors.dim_to_stride(I, 1) == 1
    @test NDTensors.dim_to_stride(I, 2) == idim
  end
  @testset "setprime" begin
    I = IndexSet(i, j)
    J = setprime(I, 2, i)
    @test i'' ∈ J
  end
  @testset "prime" begin
    I = IndexSet(i, j)
    J = prime(I, j)
    @test i ∈ J
    @test j' ∈ J
    J = prime(I; inds=j)
    @test i ∈ J
    @test j' ∈ J
    J = prime(I; inds=not(j))
    @test i' ∈ J
    @test j ∈ J
  end
  @testset "noprime" begin
    I = IndexSet(i'', j')
    J = noprime(I)
    @test i ∈ J
    @test j ∈ J
  end
  @testset "swapprime" begin
    I = IndexSet(i, j)
    @test swapprime(I, 0, 1) == IndexSet(i', j')
    @test swapprime(I, 0, 4) == IndexSet(i^4, j^4)
    I = IndexSet(i, j'')
    @test swapprime(I, 2, 0) == IndexSet(i'', j)
    I = IndexSet(i, j'', k, l)
    @test swapprime(I, 2, 0) == IndexSet(i'', j, k'', l'')
    I = IndexSet(i, k'', j'')
    @test swapprime(I, 2, 1) == IndexSet(i, k', j')
    # In-place version:
    I = IndexSet(i, k'', j''')
    I = swapprime(I, 2, 0)
    @test I == IndexSet(i'', k, j''')
    # With tags specified:
    I = IndexSet(i, k, j)
    @test swapprime(I, 0, 1, "i") == IndexSet(i', k, j)
    @test swapprime(I, 0, 1, "j") == IndexSet(i, k, j')

    I = IndexSet(i, i', j)
    @test swapprime(I, 0, 1, "i") == IndexSet(i', i, j)
    @test swapprime(I, 0, 1, "j") == IndexSet(i, i', j')
  end

  @testset "swaptags" begin
    i1 = Index(2, "Site,A")
    i2 = Index(2, "Site,B")
    is = IndexSet(i1, i2)
    sis = swaptags(is, "Site", "Link")
    for j in sis
      @test !hastags(j, "Site")
      @test hastags(j, "Link")
    end
  end

  @testset "hastags" begin
    i = Index(2, "i, x")
    j = Index(2, "j, x")
    is = IndexSet(i, j)
    @test hastags(is, "i")
    @test anyhastags(is, "i")
    @test !allhastags(is, "i")
    @test allhastags(is, "x")
  end

  @testset "broadcasting" begin
    x = Index([QN(n) => 1 for n in 0:1], "x")
    y = Index([QN(n) => 2 for n in 0:1], "y")
    I = IndexSet(x, y)

    # prime
    J = prime.(I)
    # broken for now
    #@inferred broadcast(prime, I)
    @test J isa IndexSet
    @test x' ∈ J
    @test y' ∈ J

    # prime 2
    J = prime.(I, 2)
    # broken for now
    #@inferred broadcast(prime, I, 2)
    @test J isa IndexSet
    @test x'' ∈ J
    @test y'' ∈ J

    # tag
    J = addtags.(I, "t")
    # broken for now
    #@inferred broadcast(addtags, I, "t")
    @test J isa IndexSet
    @test addtags(x, "t") ∈ J
    @test addtags(y, "t") ∈ J

    # dag
    J = dag.(I)
    # broken for now
    #@inferred broadcast(dag, I)
    @test J isa IndexSet
    @test x ∈ J
    @test y ∈ J
    @test dir(J[1]) == -dir(I[1])
    @test dir(J, x) == -dir(I, x)
    @test dir(J[2]) == -dir(I[2])
    @test dir(J, y) == -dir(I, y)
    @test ITensors.dirs(J, (x, y)) == [-dir(I, x), -dir(I, y)]
    @test ITensors.dirs(J) == [-dir(I, x), -dir(I, y)]

    # dir
    dirsI = dir.(I)
    # broken for now
    #@inferred broadcast(dir, I)
    @test dirsI isa Vector{ITensors.Arrow}
    @test dirsI == [ITensors.Out, ITensors.Out]

    # dims
    dimsI = dim.(I)
    # broken for now
    #@inferred broadcast(dim, I)
    @test dimsI isa Vector{Int}
    @test dimsI == [2, 4]

    # pairs
    J = prime.(I)
    pairsI = I .=> J
    #@inferred broadcast(=>, I, J)
    @test pairsI isa Vector{<:Pair}
    @test pairsI == [x => x', y => y']

    pairsI = I .=> 1
    #@inferred broadcast(=>, I, 1)
    @test pairsI isa Vector{<:Pair}
    @test pairsI == [x => 1, y => 1]

    pairsI = I .=> (1, 2)
    #@inferred broadcast(=>, I, (1, 2))
    @test pairsI isa Vector{<:Pair}
    @test pairsI == [x => 1, y => 2]

    pairsI = I .=> [1, 2]
    #@inferred broadcast(=>, I, [1, 2])
    @test pairsI isa Vector{<:Pair}
    @test pairsI == [x => 1, y => 2]
  end
end

nothing


text:
const PKG_PATH = @__DIR__
const DEFAULT_OUTDIR = joinpath(PKG_PATH, "../output")
mkpath(DEFAULT_OUTDIR)

const R = Float64
const Z = Int64
const G = 32.17 # gravitational acceleration [ft/s^2]

wrap_to_pi(x::R) =  x - (2π * floor((x + π) / (2π)))

"""
Fill array from expression evaluation.
"""
macro array(n, body)
    :([$(esc(body)) for _ in 1:$(esc(n))])
end


text:
function Berry(kx,ky,Ham!,H1,H2,H3)
    Ham!(kx,ky,H1) ; Hx = similar(H1) ; Hy = similar(H1) ;
    val , vec = eigen!(H1)

    curv = similar(val)

    Dx(kx,ky,Ham!,H2,H3,Hx)
    vx = vec'*Hx*vec

    denominator = [ ifelse( n-m ≈ 0, 0, 1/(n-m)^2 ) for n in val , m in val]

    Dy(kx,ky,Ham!,H2,H3,Hy)
    vy = vec'*Hy*vec

    @. vy *= denominator

    @inbounds for c in eachindex(curv)
        curv[c] = -2*imag(sum(vx[c,m]*vy[m,c] for m in 1:length(val)))
    end
    curv
end


text:
module IntegrationAlgorithms

using ....NQCDynamics:
    NQCDynamics,
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    DynamicsMethods,
    DynamicsUtils,
    Calculators,
    natoms, nbeads, ndofs

include("mdef_baoab.jl")
include("bcocb.jl")
include("mint.jl")
include("bcb_electronics.jl")

end # module


text:
"""
    Consecutive(size, step = size)
    Consecutive(size; step = size)

Sliding window of width `size` and interval `step`. Yield tuples.

This transducer is like [`Partition`](@ref) but feeds tuples to the
downstream transducers instead of vectors.

If `step == 1`, this transducer supports parallel reduction for any
collections; i.e., `Consecutive(size, 1)'(op)` is associative if
`op` is associative.

!!! warning
    Currently, in parallel folds, `Consecutive(size, 1)` cannot be
    used with reducing functions that can produce a `Reduced`.

If `step > 1`, this transducer can, in principle, support parallel reduction
if the input colection allows random access (e.g., arrays). However, this
feature is not implemented yet.

# Examples
```jldoctest
julia> using Transducers

julia> 1:8 |> Consecutive(3) |> collect
2-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 2, 3)
 (4, 5, 6)

julia> 1:8 |> Consecutive(3; step = 1) |> collect
6-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 2, 3)
 (2, 3, 4)
 (3, 4, 5)
 (4, 5, 6)
 (5, 6, 7)
 (6, 7, 8)
```
"""
struct Consecutive{Size<:Val,Step<:Val} <: Transducer
    size::Size
    step::Step

    function Consecutive(size::Val, step::Val)
        @argcheck valueof(size) isa Integer
        @argcheck valueof(step) isa Integer
        @argcheck valueof(size) > 0
        @argcheck valueof(step) > 0
        size′ = mapvalue(Int, size)
        step′ = mapvalue(Int, step)
        return new{typeof(size′),typeof(step′)}(size′, step′)
    end
end

Consecutive(size, step) =
    Consecutive(asval(size), asval(step))
Consecutive(size; step = asval(size)) = Consecutive(size, step)

#=
eduction(
    ::Consecutive{Size,Val{step}},
    xs::AbstractArray,
) where {Size<:Val,step} = (
    ntuple(i -> (@inbounds xs[i]), Size) for
    i in firstindex(xs):step:lastindex(xs)-valueof(Size)  # ???
)
=#

mapvalue(f, ::Val{x}) where {x} = Val(f(x))
valueof(::Val{x}) where {x} = x

asval(x::Val) = x
asval(x) = Val(x)

# Special-casing step=1 case for supporting parallel reduction:
const SingleStepConsecutive{Size<:Val} = Consecutive{Size,Val{1}}

function start(rf::R_{SingleStepConsecutive}, acc)
    left = () # first `size-1` items of _this chunk_ saved for `combine`
    buf = () # immediate last `size-1` items for output
    return wrap(rf, (left, buf), start(inner(rf), acc))
end

next(rf::R_{SingleStepConsecutive}, acc, input) =
    wrapping(rf, acc) do (left, buf0), iacc
        output = (buf0..., input)
        if length(output) === valueof(xform(rf).size)
            ((left, tail(output)), next(inner(rf), iacc, output))
        else
            ((output, output), iacc)
        end
    end

function complete(rf::R_{SingleStepConsecutive}, acc)
    (_left, buf), iacc = unwrap(rf, acc)
    return complete(inner(rf), iacc)
end

function combine(rf::R_{SingleStepConsecutive}, a, b)
    (_, bufa), _ = unwrap(rf, a)
    (leftb, bufb), irb = unwrap(rf, b)
    (lefta, _), ira = unwrap(rf, @return_if_reduced foldlargs(rf, a, leftb...))
    irc = combine(inner(rf), ira, irb)
    bufc = lastitems((bufa..., bufb...), Val(valueof(xform(rf).size) - 1))
    return wrap(rf, (lefta, bufc), irc)
end

lastitems(xs::NTuple{N,Any}, ::Val{M}) where {N,M} =
    N >= M ? ntuple(i -> xs[N-M+i], Val(M)) : xs

start(rf::R_{Consecutive}, acc) =
    wrap(rf, ((), 0), start(inner(rf), acc))

next(rf::R_{Consecutive}, acc, input) =
    wrapping(rf, acc) do (buf0, c), iacc
        output = (buf0..., input)
        if length(output) === valueof(xform(rf).size)
            if c == 0
                ((tail(output), valueof(xform(rf).step) - 1), next(inner(rf), iacc, output))
            else
                ((tail(output), c - 1), iacc)
            end
        else
            ((output, 0), iacc)
        end
    end

function complete(rf::R_{Consecutive}, acc)
    _, iacc = unwrap(rf, acc)
    return complete(inner(rf), iacc)
end


text:
mutable struct DataUnit
    id          :: Union{Symbol, String}
    grid        :: Symbol    # could be ('', s) x (T, U, V, W, UV)
    mask        :: Union{Symbol, Nothing}
    data        :: AbstractArray  # original data

    # if data is in spatial structure, then it will oriented as
    # sdata1 = z > x > y
    # sdata2 = x > y > z
    sdata1      :: AbstractArray
    sdata2      :: AbstractArray

    function DataUnit(
        data_table  :: Any,
        id          :: Union{Symbol, String},
        grid        :: Symbol,
        mask        :: Union{Symbol, Nothing},
        data        :: AbstractArray,
    )

        if mask == :mask
            mask = grid
        end

        if ! haskey(data_table.grid_dims, grid)
            throw(ErrorException("Unknown grid: " * string(grid)))
        end

        grid_dim = data_table.grid_dims[grid]
        if length(data) != reduce(*, grid_dim)
            throw(ErrorException("Grid type and data length mismatched: " * string(id)))
        end

        s1data = reshape(data, grid_dim...)
        s2data = PermutedDimsArray(s1data, [2,3,1])

        return new(
            id,
            grid,
            mask,
            data,
            s1data,
            s2data,
        )
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
struct PolicyValueTestResult
    test::PolicyValueTest
    prediction::Vector{Float32}
end

struct PolicyValueTestResultSVG
    result::PolicyValueTestResult
end

svg(x::PolicyValueTestResult) = PolicyValueTestResultSVG(x)
policypredsvg(state::AbstractState, ::PolicyValueTestResult) = @assert false OVERRIDE_MSG

function getdf(T::Type{PolicyValueTestResult}, results, limit=100)
    test_num = "Test #"
    test_pf = "P/F"
    test_result = "Pred (π*)"
    optimal_action = "Actual (π*)"
    d = Dict(
        test_num=>[],
        test_pf=>[],
        test_result=>[],
        optimal_action=>[],
    )
    for result in results
        idx = result.test.idx
        optimal = result.test.actions
        pred = result.prediction
        _, greedy_action = findmax(pred)
        pass = in(greedy_action, optimal)
        from = 1
        to = min(limit, length(pred))
        actions = reverse(sortperm(pred))[from:to]
        prs = pred[actions]
        action_prs = map(x->begin
            action, pr = x
            pr_s = @sprintf "%.2f" pr
            "$action ($(pr_s[2:end]))"
        end, zip(actions, prs))
        result_policy = join(action_prs, ", ")
        push!(d[test_num], idx)
        push!(d[test_pf], pass ? "Pass" : "Fail")
        push!(d[test_result], result_policy)
        push!(d[optimal_action], join(optimal, ","))
    end
    DataFrame(
        test_num=>d[test_num],
        test_pf=>d[test_pf],
        test_result=>d[test_result],
        optimal_action=>d[optimal_action],
    )
end

function Base.show(io::IO, m::MIME"text/plain", x::PolicyValueTestResult)
    println()
    pretty_results(
        io,
        getdf(PolicyValueTestResult, [x]);
        title="Policy Test Result #$(x.test.idx)"
    )
end

function Base.show(io::IO, m::MIME"image/svg+xml", x::PolicyValueTestResultSVG)
    test = x.result.test
    pred = x.result.prediction
    s = test.state
    player_name = s.to_play.name
    actions = test.actions
    state_grid = getgrid(s; header="to play: $player_name", margin_left=0)
    next_states = Rect[]
    from = 1
    len = length(actions)
    to = min(len, 3)
    for (idx,a) in enumerate(actions[from:to])
        push!(
            next_states,
            getgrid(s, a, 0f0, getnextstate(s, a); header="action: $a", margin_left=0)
        )
    end
    if length(actions) > 3
        push!(next_states, Grid(1,1, (r,c)->"..."; line_thickness=0))
    end
    input_grid = Grid(
        [state_grid],
        :,
        1,
        header="INPUT",
        margin_left=0,
        line_thickness=0
    )
    policy_grid_header = "ACTUAL π*"
    policy_grid = Grid(
        next_states,
        :,
        1,
        header=(text=policy_grid_header, height=40),
        margin_left=0,
        line_thickness=0,
    )
    pred_states = Rect[]
    limit = 3
    from = 1
    to = min(limit, length(pred))
    pred_actions = reverse(sortperm(pred))[from:to]
    prs = pred[pred_actions]
    for (a,pr) in zip(pred_actions, prs)
        pr_str = @sprintf "%.2f" pr
        push!(
            pred_states,
            getgrid(
                s,
                a,
                0f0,
                getnextstate(s, a);
                margin_left=0,
                header="action: $a",
                subheader=(text=pr_str, color="#666"),
            )
        )
    end
    pred_grid = Grid(
        pred_states,
        :,
        1;
        header="PRED π*",
        margin_left=0,
        line_thickness=0,
    )
    g = Grid([input_grid,pred_grid,policy_grid]; cell_valign=:top, header="POLICY TEST RESULT #$(test.idx)", line_thickness=0)
    show(io, m, g)
end


text:
"""
Provider module for the `biber` executable.
"""
module Biber

export biber

using Pkg.Artifacts

binary() = joinpath(artifact"biber_bin", "biber$(Sys.iswindows() ? ".exe" : "")")

"""
    biber(f)

Call function `f` with the *biber* binary path as argument.

```julia
biber() do bin
    run(`\$bin file`)
end
```
"""
biber(f) = f(binary())

version() = VersionNumber(lstrip(!isnumeric, biber(bin -> readchomp(`$bin --version`))))

end


text:
type FiveStepFunctionArms <: BanditEnvironment
	shift_time::Int64
end

function arms(env::FiveStepFunctionArms, trial::Int64)
  if trial > env.shift_time
    [
      BernoulliArm(0.1),
      BernoulliArm(0.2),
      BernoulliArm(0.3),
      BernoulliArm(0.4),
      BernoulliArm(0.5)
    ]
  else
    [
      BernoulliArm(0.5),
      BernoulliArm(0.4),
      BernoulliArm(0.3),
      BernoulliArm(0.2),
      BernoulliArm(0.1)
    ]
  end
end

function n_arms(env::FiveStepFunctionArms)
  5
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
using GLAbstraction, Colors, GeometryTypes, GLVisualize, Reactive

function showpath(;follow = true, movie = false, x0 = false, truth = true, obs = true, smooth = true, sample = true, ode = true, nu = true, rotating = false)

    window = glscreen()
    if follow
        ss = 1:2:div(length(tt),2)
    else
        ss = 1:4:5000
    end
    timesignal = loop(ss)

    extra = 10

    Yxyz = collect(Point3f0, Iterators.flatten(Pt[i].Y.yy[1:end-1] for i in 1:m)) # Y
    #Yxyz = collect(Point3f0, Iterators.flatten(Pᵒ[i].V[1:end-1] for i in 1:m)) # V

    XXxyz = collect(Point3f0, Iterators.flatten(XX[i].yy[1:end-1] for i in 1:m))
    XXmeanxyz = collect(Point3f0, Iterators.flatten(XXmean[i].yy[1:end-1] for i in 1:m))
    XXstdr = vcat([Point3f0.( XXstd[i][1:end-1] ) for i in 1:m]...)
    XXscale = vcat([Point3f0.(XXscal[i][1:end-1]) for i in 1:m]...)
    XXrotation = vcat([Vec4f0.(XXrot[i][1:end-1]) for i in 1:m]...)

    Xxyz = Point3f0.(Xtrue.yy)

    Vxyz = Point3f0.(V.yy)

    if follow
        phi = Float32(0); sca = 1.
    else
        sca = 1/15
        phi = Float32(0.2*2pi)
    end

    if rotating
        rotation = map(timesignal) do i
            phi = 25.5tt[i]
            sca_ = sca+ sca*(i>3000)*(2*atan((i-3000)/1000)/pi)
            @SMatrix Float32[ sca_*cos(phi)  -sca_*sin(phi)    0.0     0.0
            sca_*sin(phi)   sca_*cos(phi)    0.0       0.0
            0.0            0.0             sca_*1.0   -sca_/sca*1.25
            0.0            0.0             0.0       1.0]
        end
    else
        rotation = @SMatrix Float32[ sca*cos(phi)  -sca*sin(phi)    0.0       0.0
                                 sca*sin(phi)   sca*cos(phi)    0.0       0.0
                                 0.0            0.0             sca*1.0  -1.25
                                 0.0            0.0             0.0       1.0]

    end

    eyeposition = map(timesignal) do i
        Vec{3,Float32}(XXmeanxyz[i] )
    end

    # create the camera lookat and up vector
    #lookat = Signal(Vec3f0(0))
    lookat = map(timesignal) do i
        Vec{3,Float32}(XXmeanxyz[i+200] - 0XXmeanxyz[i])
    end
    upvector = Signal(Vec3f0(0,-1,0))

    # create a camera from these
    #camera = PerspectiveCamera(window.area, eyeposition, lookat, upvector)
    cam = PerspectiveCamera(
        Signal(Vec3f0(0)), # theta (rotate by x around cam xyz axis)
        Signal(Vec3f0(0)), # translation (translate by translation in the direction of the cam xyz axis)
        lookat, # lookat. We want to look at the middle of the cubes
        eyeposition, # camera position. We want to be on the same height, but further away in y
        upvector, #upvector
        window.area, # window area
        Signal(80f0), # Field of View
        Signal(2f0),  # Min distance (clip distance)
        Signal(100f0), # Max distance (clip distance)
    )
    if movie
        Xvisible = map(timesignal) do i; (500 < i < 950 || i > 3000) end
        Xsmvisible = map(timesignal) do i; (i < 450) end
        Yvisible = map(timesignal) do i; (1500 < i < 1950) end
        XXvisible =  map(timesignal) do i; (2000 < i < 2450) end
        XXmeanvisible = map(timesignal) do i; (i > 2500) end
        XXdistvisible = map(timesignal) do i; (i > 3500) end
        Vvisible = map(timesignal) do i; (1000 < i < 3000) end

    else
        Xvisible =
        Xvisible =
        Xsmvisible =
        Yvisible =
        XXvisible =
        XXmeanvisible =
        XXdistvisible =
        Vvisible = map(timesignal) do i; true end
    end
    X0visible = map(timesignal) do i; true end

    Xcolor = RGBA{Float32}(0.04, 0.15, 0.44, 0.6)
    Xsmcolor = RGBA{Float32}(0.04, 0.15, 0.44, 0.8)
    Ycolor = RGBA{Float32}(0.34, 0.05, 0.14, 0.4)
    XXcolor = RGBA{Float32}(0.04, 0.35, 0.14, 0.4)
    XXmeancolor = RGBA{Float32}(0.04, 0.35, 0.14, 0.15)
    XXdistcolor =  RGBA{Float32}(0.34, 0.05, 0.14, 0.40)
    Vcolor = RGBA{Float32}(0.7, 0.3, 0., 0.5)

    X3d = visualize(
        Xxyz[1:skippoints:end], :lines,
        color = Xcolor,
        model = rotation,
        visible = Xvisible
    )
    Xsm3d = visualize(
        Point3f0.(X2.yy[1:skippoints:end]), :lines,
        color = Xsmcolor,
        model = rotation,
        visible = Xsmvisible
    )

    Y3d = visualize(
        Yxyz[1:skippoints:end], :lines,
        color = Ycolor,
        model = rotation,
        visible = Yvisible
    )

    XX3d = visualize(
        XXxyz[1:skippoints:end], :lines,
        color = XXcolor,
        model = rotation,
        visible = XXvisible
    )

    sphere = Sphere(Point3f0(0,0,0), 1.0f0)
    XXmean3d = visualize(
        (sphere, XXmeanxyz[1:5extra*skippoints:end]),
        #scale = Float32(2*sca)*XXstdr[1:5extra*skippoints:end],
        scale = Float32(1)*XXscale[1:5extra*skippoints:end],
        rotation = XXrotation[1:5extra*skippoints:end],
        color = XXmeancolor,
        model = rotation,
        visible = XXmeanvisible
    )

    XXdist = visualize(
        collect(Iterators.flatten(zip(XXmeanxyz[1:extra*skippoints:end],
         Xxyz[1:extra*skippoints:end-1]))), :linesegment,
        color = XXdistcolor,
        model = rotation,
        visible = XXdistvisible
    )

    sphere2 = Sphere(Point3f0(0), 1f0)
    V3d = visualize(
        (sphere2, Vxyz),
        scale = fill(0.5Point3f0(1,1,1), length(Vxyz)),
        color = Vcolor,
        model = rotation,
        visible = Vvisible
    )

    sphere3 = Sphere(Point3f0(0), 0.1f0)
    X03d = visualize(
        (sphere3, Point3f0.(X0[1:10:end])),
        color = Xcolor,
        model = rotation,
        visible = X0visible
    )

    if follow
        camera = cam
    else
        camera = :perspective
    end
    if smooth
        _view(XXmean3d, window, camera=camera)
    end
    if sample
        _view(XX3d, window, camera=:perspective)
    end
    if x0
        _view(X03d, window, camera=:perspective)
    end
    if nu
        _view(Y3d, window, camera=camera)
    end

    if ode
        _view(Xsm3d, window, camera=camera)
    end
    if truth
        _view(X3d, window, camera=camera)
    end
    if truth && smooth
        _view(XXdist, window, camera=camera)
    end
    if obs
        _view(V3d, window, camera=camera)
    end

    if !movie
        renderloop(window)
    else

        # save video to report dir, or in some tmp dir we'll delete later
        path = "./output"

        isdir(path) || mkdir(path)
        name = path * "/$(randstring()).mkv"

        #@async renderloop(window)

        # create a stream to which we can add frames
        io = GLVisualize.create_video_stream(name, window)
        for i in 1:div(3900,2)
            # do something
               #render current frame
            # if you call @async renderloop(window) you can replace this part with yield
            #yield()
            GLWindow.render_frame(window)
            GLWindow.swapbuffers(window)
            GLWindow.reactive_run_till_now()

            # add the frame from the current window
            GLVisualize.add_frame!(io)
        end
        # closing the stream will trigger writing the video!
        close(io.io)
        GLWindow.destroy!(window)
    end
end

showpath(;follow = false, movie = false, truth = true, obs = false, smooth = true, sample = false, x0 = false, ode = false, nu = true, rotating = false)

#showpath(;follow = false, movie = true, rotating = true)


text:
"Bayesian inference on state space models"
module Baytes

################################################################################
#Import modules
import BaytesCore: BaytesCore, update!, infer, results, init, init!, propose, propose!, propagate!
using BaytesCore:
    BaytesCore,
    AbstractAlgorithm,
    AbstractTune,
    AbstractConfiguration,
    AbstractDiagnostics,
    AbstractConstructor,
    Updater,
    Iterator,
    subset,
    UpdateBool,
    UpdateTrue,
    UpdateFalse,
    DataTune,
    ArrayConfig,
    DataStructure,
    Expanding,
    Rolling,
    Batch,
    SubSampled,
    adjust,
    generate_showvalues,
    TemperingMethod,
    IterationTempering,
    JointTempering,
    ValueHolder,
    PrintDefault,
    ProgressLog,
    SilentLog,
    ConsoleLog,
    ProgressReport,
    SampleDefault

import ModelWrappers:
    ModelWrappers,
    sample,
    sample!

using ModelWrappers:
    ModelWrappers,
    ModelWrapper,
    Tagged,
    Objective,
    DiffObjective,
    AbstractDifferentiableTune,
    ℓObjectiveResult,
    ℓDensityResult,
    ℓGradientResult,
    flatten,
    paramnames,
    FlattenDefault,
    FlattenTypes,
    FlattenAll,
    FlattenContinuous,
    UnflattenTypes,
    UnflattenStrict,
    UnflattenFlexible

using BaytesMCMC, BaytesFilters, BaytesPMCMC, BaytesSMC

#Utility tools
import Base: Base, summary
using DocStringExtensions:
    DocStringExtensions, TYPEDEF, TYPEDFIELDS, FIELDS, SIGNATURES, FUNCTIONNAME
using ArgCheck: ArgCheck, @argcheck
using UnPack: UnPack, @unpack, @pack!
using Random: Random, AbstractRNG, GLOBAL_RNG, randexp
using ProgressMeter: ProgressMeter, @showprogress, Progress, next!
# Post-Processing
using MCMCDiagnosticTools: MCMCDiagnosticTools, ess_rhat, mcse, gelmandiag
using Statistics: Statistics, mean, std, sqrt, quantile, var
using Dates: Dates, today, hour, minute, now
using JLD2: JLD2, jldsave
using PrettyTables: PrettyTables, pretty_table

################################################################################
#Abstract types to be dispatched in Examples section
include("sampling/Sampling.jl")

################################################################################
export
    update!,
    init,
    init!,
    sample,
    sample!,
    propose,
    propose!,
    propagate!,

    ## ModelWrappers ~ For now explicitly import if used with Baytes.jl
#    dynamics,
#    predict,
#    generate,

    #Re-export relevant types/structs from other modules
    ## BaytesCore
    Updater,
    Iterator,
    UpdateBool,
    UpdateTrue,
    UpdateFalse,
    DataTune,
    ArrayConfig,
    DataStructure,
    Expanding,
    Rolling,
    Batch,
    SubSampled,
    TemperingMethod,
    IterationTempering,
    JointTempering,
    ValueHolder,
    PrintDefault,
    ProgressLog,
    SilentLog,
    ConsoleLog,
    ProgressReport,
    SampleDefault,

    ## ModelWrappers
    flatten,
    FlattenDefault,
    FlattenTypes,
    FlattenAll,
    FlattenContinuous,
    UnflattenTypes,
    UnflattenStrict,
    UnflattenFlexible,

    ## BaytesFilters
    ParticleFilter,
    ParticleFilterDefault,
    ParticleFilterConstructor,
    ParticleFilterMemory,

    Markov,
    SemiMarkov,
    SemiMarkovInitiation,
    SemiMarkovTransition,
    Bootstrap,
    Systematic,
    Stratified,
    Residual,
    Multinomial,
    Ancestral,
    Marginal,
    Conditional,

    ## BaytesMCMC
    MCMC,
    MCMCDefault,
    MCMCConstructor,
    NUTS,
    ConfigNUTS,
    HMC,
    ConfigHMC,
    MALA,
    ConfigMALA,
    Metropolis,
    ConfigMetropolis,
    Custom,
    ConfigCustom,

    ConfigStepnumber,
    StepNumberTune,
    KineticEnergy,
    ConfigTuningWindow,
    ConfigStepsize,
    ConfigProposal,
    MatrixMetric,
    MDense,
    MDiagonal,
    MUnit,
    Warmup,
    Adaptionˢˡᵒʷ,
    Adaptionᶠᵃˢᵗ,
    Exploration,

    ## BaytesPMCMC
    PMCMC,
    PMCMCDefault,
    PMCMCConstructor,
    ParticleMetropolis,
    ParticleGibbs,

    ## BaytesSMC
    SMC,
    SMCDefault,
    SMCConstructor,

    SMC2,
    SMC2Constructor,
    SMC2Kernel,
    SMCweight,
    SMCreweight

end


text:
# To see installed packages, do Pkg.status()
using Pkg
Pkg.add("Distributed")
Pkg.add("DataFrames")
Pkg.add("DataFramesMeta")
Pkg.add("DataStructures")
Pkg.add("StatsBase")
Pkg.add("Combinatorics")
Pkg.add("Printf")
Pkg.add("Dates")
Pkg.add("CSV")
Pkg.add("Statistics")
Pkg.add("Distributions")
Pkg.add("HypothesisTests")
Pkg.add("Random")
Pkg.add("QuadGK")
Pkg.add("Measures")
Pkg.add("Test")
Pkg.add("ColorTypes")
Pkg.add("HDF5")
Pkg.add("JLD")
Pkg.add("Colors"      )
Pkg.add("Test")
Pkg.add("Permutations")
Pkg.add("Plots")   # only used in Windows
#Pkg.add("Revise")
using Distributed
using DataFrames
using DataFramesMeta
using StatsBase
using Combinatorics
using Printf
using Dates
using CSV
using Statistics
using Distributions
using HypothesisTests
using Random
using QuadGK
using Measures
using Test
using Permutations
#using Revise


text:
using DelayDiffEq
using Test

# out-of-place problem
function f_notinplace(u,h,p,t)
  [-h(p, t-1/5)[1] + u[1]; -h(p, t-1/3)[2] - h(p, t-1/5)[2]]
end
const prob_notinplace = DDEProblem(f_notinplace, ones(2), (p, t)->zeros(2), (0.0, 100.0),
                                   constant_lags = [1/5, 1/3])

# in-place problem
function f_inplace(du,u,h,p,t)
  du[1] = - h(p, t-1/5)[1] + u[1]
  du[2] = - h(p, t-1/3)[2] - h(p, t-1/5)[2]
end
const prob_inplace = DDEProblem(f_inplace, ones(2), (p, t)->zeros(2), (0.0, 100.0),
                                constant_lags = [1/5, 1/3])

const alg = MethodOfSteps(BS3())

@testset for prob in (prob_notinplace, prob_inplace)
  # reference solution
  dde_int = init(prob, alg)
  sol = solve!(dde_int)

  # save all components
  @testset "all components" begin
    # without keyword argument
    @testset "without keyword" begin
      # solution, solution of DDE integrator, and solution of ODE integrator
      # contain all components
      @test length(sol.u[1]) == 2
      @test length(dde_int.sol.u[1]) == 2
      @test length(dde_int.integrator.sol.u[1]) == 2

      # solution and solution of DDE integrator are equal
      @test sol.t == dde_int.sol.t
      @test sol.u == dde_int.sol.u

      ## interpolation
      @test sol(25:100, idxs=2) ≈ [u[1] for u in sol(25:100, idxs=[2])]
    end

    # with keyword argument
    @testset "with keyword" begin
      dde_int2 = init(prob, alg; save_idxs=[1, 2])
      sol2 = solve!(dde_int2)

      # solution, solution of DDE integrator, and solution of ODE integrator
      # contain all components
      @test length(sol2.u[1]) == 2
      @test length(dde_int2.sol.u[1]) == 2
      @test length(dde_int2.integrator.sol.u[1]) == 2

      # solution and solution of DDE integrator are equal
      @test sol.t == dde_int.sol.t
      @test sol.u == dde_int.sol.u

      # interpolation
      @test sol[2, :] ≈ dde_int.integrator.sol(sol.t, idxs=2)
    end
  end

  # save only second component
  @testset "second component" begin
    # array index
    @testset "array index" begin
      dde_int2 = init(prob, alg; save_idxs=[2])
      sol2 = solve!(dde_int2)

      # solution and solution of DDE integrator contain only second component
      @test length(sol2.u[1]) == 1
      @test length(dde_int2.sol.u[1]) == 1

      # solution of ODE integrator contains both components
      @test length(dde_int2.integrator.sol.u[1]) == 2

      # solution and solution of DDE integrator are equal
      @test sol2.t == dde_int2.sol.t
      @test sol2.u == dde_int2.sol.u

      # interpolation
      @test sol(25:100, idxs=2) ≈ sol2(25:100, idxs=1)
      @test sol(25:100, idxs=[2]) ≈ sol2(25:100, idxs=[1])
    end

    # scalar index
    @testset "scalar index" begin
      dde_int2 = init(prob, alg; save_idxs=2)
      sol2 = solve!(dde_int2)

      # solution and solution of DDE integrator is only vector of floats
      @test typeof(sol2.u) === Vector{Float64}
      @test typeof(dde_int2.sol.u) === Vector{Float64}

      # solution of ODE integrator contains both components
      @test length(dde_int2.integrator.sol.u[1]) == 2

      # solution and solution of DDE integrator are equal
      @test sol2.t == dde_int2.sol.t
      @test sol2.u == dde_int2.sol.u

      # solution equals second component of complete solution
      @test sol.t ≈ sol2.t && sol[2, :] ≈ sol2.u

      # interpolation of solution equals second component of
      # interpolation of complete solution
      @test sol(25:100, idxs=2) ≈ sol2(25:100, idxs=1)
    end
  end
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
#=

  Diet problem in Julia ConstraintSolver.jl

  Standard diet problem.

  Minimize the cost for the products:
   Type of                        Calories   Chocolate    Sugar    Fat
   Food                                      (ounces)     (ounces) (ounces)
   Chocolate Cake (1 slice)       400           3            2      2
   Chocolate ice cream (1 scoop)  200           2            2      4
   Cola (1 bottle)                150           0            4      1
   Pineapple cheesecake (1 piece) 500           0            4      5

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function diet(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    prices = problem[:prices]
    limits = problem[:limits]
    nutrition = problem[:nutrition]
    max_val = problem[:max_val]

    len = length(prices)

    @variable(model, 0 <= x[1:len] <= max_val, Int)
    @variable(model, 1 <= z <= maximum(prices)*max_val, Int)

    for i in 1:len
        scalar_product(model,nutrition[i,:],x,:(>=),limits[i])
    end

    @constraint(model, z == sum(prices.*x) )

    @objective(model,Min,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("z:$z_val")
                println("x:$x_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end

diet_problems = Dict(
:1 => Dict(
     :prices => [ 50, 20, 30, 80.4],  # price in cents for each nutrition
     :limits => [500,  6, 10,  8], #limits, requirements for each nutrition type

     # nutrition for each product
     :nutrition =>
     resize_matrix([[400, 200, 150, 500],  # calories
                  [  3,   2,   0,   0],  # chocolate
                  [  2,   2,   4,   4],  # sugar
                  [  2,   4,   1,   5]] # fat
   ),
     :max_val => 10
   ),
)

@time diet(diet_problems[:1])

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# need one:
# 1. FASTA  with full non redundant sequences
# 2. FASTA with envelop of domain
# return: envelopes to align

function enveloptoalign(filefull; ctype::Symbol=:amino, pos::Bool=true)
    dful = readfull(filefull, ctype=ctype, pos=pos)
    #dful = readenvelop(filefull, ctype=ctype)
    return dful
end

function enveloptoalign(filefull,fileenv,fileins; delta = 0, ctype::Symbol=:amino)
    denv = readenvelop(fileenv, ctype=ctype)
    #println(length(denv))
    #print(denv)
    if ctype == :amino
        dful = readfull(filefull, ctype=ctype, pos=false)
    else
        dful = readfull(filefull, ctype=ctype, pos=true)
    end
    #println(length(dful))
    dins = readenvelop(fileins, ctype= ctype)
    #println(length(dins))
    al = Vector{Tuple{String,String,String,String}}()
    for (name,seq) in dful
        #print("full name ", name)
        l = length(seq)
        if haskey(denv, name)
            envelops = denv[name]
            envinsall = dins[name]
            for env in envelops
                auxins = []
                for ins in envinsall
                    if env[1] == ins[1]
                        auxins = ins
                    end
                end
                if ctype == :amino
                    lb = max(1,env[1]-delta)
                    ub = min(l,env[2]+delta)
                else
                    lb = env[1]
                    ub = env[2]
                end
                if ctype == :amino
                    envseq = seq[lb:ub]
                else
                    envseq = seq
                end
                envhmm = env[3]
                envins = auxins[3]
		if occursin('/', name)
			envname = name
		else
                	envname = name*"/$lb-$ub"
		end
                push!(al,(envname,envseq,envhmm,envins))
            end
        end
    end
    return al
end

function readfull(filefull; ctype::Symbol=:amino, pos::Bool=true)
    ffull = FastaIO.FastaReader(filefull)
    dheader = Dict{String,String}()
    for (_name,seq) in ffull
        name=String(split(_name)[1])
        if occursin('|', name)
            (aux,name) = split(name, "|")
        end
        if occursin('/', name) && !pos
            (name,aux) = split(name,"/")
        end
        if !haskey(dheader,name)
            dheader[name] = seq
        end
    end
    close(ffull)
    dheader
end

function readenvelop(filenv; ctype::Symbol=:amino)
    fenve = FastaIO.FastaReader(filenv)
    dheader = Dict{String, Any}()
    for (name,seq) in fenve
        if occursin('|', name)
            (aux,name) = split(name, "|")
        end
        if occursin('/', name) && ctype == :amino
            name,_envelop = split(name,"/")
            envelop = parse.(Int, split(_envelop,"-"))
        else
            envelop = [0, length(seq)]
        end
        if envelop[2] < envelop[1]
            if ctype == :amino
                error("strange envelop in $name")
            end
        end
        tname = String(name)
        tname = strip(tname)
        if haskey(dheader, tname)
            push!(dheader[tname],[envelop[1],envelop[2],seq])
        else
            dheader[tname] = [[envelop[1],envelop[2],seq],]
        end
    end

    close(fenve)
    dheader
end

"""
function writenewal(fileout, al)
    fp = FastaIO.FastaWriter(fileout)
    for i in eachindex(al)
        header = al[i].name
        seq = al[i].pa
        length(seq) > 0 && writeentry(fp,header,seq)
    end
    close(fp)
end
"""


text:
function modifier(ability)

end

function ability()

end

mutable struct DNDCharacter

end


text:
using BinaryBuilder, SHA
include("../../fancy_toys.jl")

name = "CompilerSupportLibraries"
version = v"0.3.3"

# We are going to need to extract the latest libstdc++ and libgomp from BB
# So let's grab them into tarballs by using preferred_gcc_version:
extraction_script = raw"""
mkdir -p ${libdir}
for d in /opt/${target}/${target}/lib*; do
    # Copy all the libstdc++ and libgomp files:
    cp -av ${d}/libstdc++*.${dlext}* ${libdir} || true
    cp -av ${d}/libgomp*.${dlext}* ${libdir} || true
    # Don't copy `.a` or `.py` files.  >:[
    rm -f ${libdir}/*.a ${libdir}/*.py
done
"""

extraction_platforms = supported_platforms()
extraction_products = [
    LibraryProduct("libstdc++", :libstdcxx),
    LibraryProduct("libgomp", :libgomp),
]

# Don't actually run extraction if we're asking for a JSON, but don't print it either
if any(startswith(a, "--meta-json") for a in ARGS)
    # How delightfully meta, for when we're calculating the meta!  ;D
    self_url = @__FILE__
    self_hash = open(io -> bytes2hex(sha256(io)), self_url)
    build_info = Dict(p => (self_url, self_hash) for p in BinaryBuilder.abi_agnostic.(extraction_platforms))
else
    build_info = autobuild(joinpath(@__DIR__, "build", "extraction"),
        "LatestLibraries",
        version,
        FileSource[],
        extraction_script,
        # Only extract for platforms we're actually going to use
        filter(should_build_platform, extraction_platforms),
        extraction_products,
        Dependency[];
        skip_audit=true,
        preferred_gcc_version=v"100",
        verbose="--verbose" in ARGS,
        debug="--debug" in ARGS,
    )
end

## Now that we've got those tarballs, we're going to use them as sources to overwrite
## the libstdc++ and libgomp that we would otherwise get from our compiler shards:
script = raw"""
# Start by extracting LatestLibraries
tar -zxvf ${WORKSPACE}/srcdir/LatestLibraries*.tar.gz -C ${prefix}

echo ***********************************************************
echo LatestLibraries logs, reproduced here for debuggability:
cat ${prefix}/logs/LatestLibraries.log
echo ***********************************************************
rm -f ${prefix}/logs/LatestLibraries.log

# Make sure expansions aren't empty
shopt -s nullglob

# copy out all the libraries we can find, not clobbering stuff from LL
for d in /opt/${target}/${target}/lib*; do
    cp -uav ${d}/*.${dlext}* ${libdir}/ || true
done

# libwinpthread is a special snowflake and is only within `bin` for some reason
if [[ ${target} == *mingw* ]]; then
	cp -uav /opt/${target}/${target}/sys-root/bin/*.${dlext}* ${libdir}/
fi

# Delete .a and .py files, we don't want those.
rm -f ${libdir}/*.a ${libdir}/*.py

# Delete any `.so` files that are not ELF files, since they're mostly likely linker scripts
for f in ${libdir}/*.so; do
    if [[ "$(file -b "$f")" != ELF* ]]; then
        rm -f "$f"
    fi
done

# change permissions so that rpath succeeds
for l in ${libdir}/*; do
    chmod 0755 "${l}"
done

# libgcc_s.1.dylib receives special treatment for now
if [[ ${target} == *apple* ]]; then
    install_name_tool -id @rpath/libgcc_s.1.dylib ${libdir}/libgcc_s.1.dylib
fi

# Install license (we license these all as GPL3, since they're from GCC)
install_license /usr/share/licenses/GPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgcc_s", "libgcc_s_sjlj", "libgcc_s_seh"], :libgcc_s),
    LibraryProduct("libstdc++", :libstdcxx),
    LibraryProduct("libgfortran", :libgfortran),
    LibraryProduct("libgomp", :libgomp),
]

for platform in platforms
    if should_build_platform(platform)
        # Find the corresponding source for this platform
        tarball_path, tarball_hash = build_info[BinaryBuilder.abi_agnostic(platform)][1:2]
        sources = [
            FileSource(tarball_path, tarball_hash),
        ]
        build_tarballs(ARGS, name, version, sources, script, [platform], products, [])
    end
end


text:
# ---
# title: 1649. Create Sorted Array through Instructions
# id: problem1649
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Indexed Tree, Segment Tree, Ordered Map
# link: <https://leetcode.com/problems/create-sorted-array-through-instructions/description/>
# hidden: true
# ---
#
# Given an integer array `instructions`, you are asked to create a sorted array
# from the elements in `instructions`. You start with an empty container `nums`.
# For each element from **left to right** in `instructions`, insert it into
# `nums`. The **cost** of each insertion is the **minimum** of the following:
#
#   * The number of elements currently in `nums` that are **strictly less than** `instructions[i]`.
#   * The number of elements currently in `nums` that are **strictly greater than** `instructions[i]`.
#
# For example, if inserting element `3` into `nums = [1,2,3,5]`, the **cost** of
# insertion is `min(2, 1)` (elements `1` and `2` are less than `3`, element `5`
# is greater than `3`) and `nums` will become `[1,2,3,3,5]`.
#
# Return _the **total cost** to insert all elements from _`instructions`
# _into_`nums`. Since the answer may be large, return it **modulo** `109 + 7`
#
#
#
# **Example 1:**
#
#
#
#     Input: instructions = [1,5,6,2]
#     Output: 1
#     Explanation: Begin with nums = [].
#     Insert 1 with cost min(0, 0) = 0, now nums = [1].
#     Insert 5 with cost min(1, 0) = 0, now nums = [1,5].
#     Insert 6 with cost min(2, 0) = 0, now nums = [1,5,6].
#     Insert 2 with cost min(1, 2) = 1, now nums = [1,2,5,6].
#     The total cost is 0 + 0 + 0 + 1 = 1.
#
# **Example 2:**
#
#
#
#     Input: instructions = [1,2,3,6,5,4]
#     Output: 3
#     Explanation: Begin with nums = [].
#     Insert 1 with cost min(0, 0) = 0, now nums = [1].
#     Insert 2 with cost min(1, 0) = 0, now nums = [1,2].
#     Insert 3 with cost min(2, 0) = 0, now nums = [1,2,3].
#     Insert 6 with cost min(3, 0) = 0, now nums = [1,2,3,6].
#     Insert 5 with cost min(3, 1) = 1, now nums = [1,2,3,5,6].
#     Insert 4 with cost min(3, 2) = 2, now nums = [1,2,3,4,5,6].
#     The total cost is 0 + 0 + 0 + 0 + 1 + 2 = 3.
#
#
# **Example 3:**
#
#
#
#     Input: instructions = [1,3,3,3,2,4,2,1,2]
#     Output: 4
#     Explanation: Begin with nums = [].
#     Insert 1 with cost min(0, 0) = 0, now nums = [1].
#     Insert 3 with cost min(1, 0) = 0, now nums = [1,3].
#     Insert 3 with cost min(1, 0) = 0, now nums = [1,3,3].
#     Insert 3 with cost min(1, 0) = 0, now nums = [1,3,3,3].
#     Insert 2 with cost min(1, 3) = 1, now nums = [1,2,3,3,3].
#     Insert 4 with cost min(5, 0) = 0, now nums = [1,2,3,3,3,4].
#     Insert 2 with cost min(1, 4) = 1, now nums = [1,2,2,3,3,3,4].
#     Insert 1 with cost min(0, 6) = 0, now nums = [1,1,2,2,3,3,3,4].
#     Insert 2 with cost min(2, 4) = 2, now nums = [1,1,2,2,2,3,3,3,4].
#     The total cost is 0 + 0 + 0 + 0 + 1 + 0 + 1 + 0 + 2 = 4.
#
#
#
#
# **Constraints:**
#
#   * `1 <= instructions.length <= 105`
#   * `1 <= instructions[i] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "postgresql" begin
    dsn = ODBC.DSN("PgSQL-test")

    @testset "basic queries" begin
        dbs = ODBC.query(dsn, "SELECT datname FROM pg_database WHERE datistemplate = false;")
        data = ODBC.query(dsn, "SELECT table_schema, table_name FROM information_schema.tables ORDER BY table_schema,table_name;")
    end

    @testset "test1" begin
        ODBC.execute!(dsn, "drop table if exists test1")
        ODBC.execute!(dsn, "create table test1
                            (test_bigint bigint,
                             test_decimal decimal,
                             test_int integer,
                             test_numeric numeric,
                             test_smallint smallint,
                             test_float real,
                             test_real double precision,
                             test_money money,
                             test_date date,
                             test_timestamp timestamp,
                             test_time time,
                             test_char char(1),
                             test_varchar varchar(16),
                             test_bytea bytea,
                             test_boolean boolean,
                             test_text text,
                             test_array integer[]
                            )")
        data = ODBC.query(dsn, "select * from information_schema.tables where table_name = 'test1'")
        println("Postgres table 'test1' schema:")
        println()
        ODBC.execute!(dsn, "insert into test1 VALUES
                            (1, -- bigint,
                             1.2, -- decimal,
                             2, -- integer,
                             1.4, -- numeric,
                             3, -- smallint,
                             1.6, -- real,
                             1.8, -- double precision,
                             2.0, -- money,
                             '2016-01-01', -- date,
                             '2016-01-01 01:01:01', -- timestamp,
                             '01:01:01', -- time,
                             'A', -- char(1),
                             'hey there sailor', -- varchar(16),
                             NULL, -- bytea,
                             TRUE, -- boolean,
                             'hey there abraham', -- text,
                             ARRAY[1, 2, 3] -- integer array
                            )")
        data = ODBC.query(dsn, "select * from test1")
        @test size(data) == (1,17)
        @test Tables.schema(data).types == (
            Union{Int64, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int32, Missing},
            Union{DecFP.Dec64, Missing},
            Union{Int16, Missing},
            Union{Float32, Missing},
            Union{Float64, Missing},
            Union{Float64, Missing},
            Union{ODBC.API.SQLDate, Missing},
            Union{ODBC.API.SQLTimestamp, Missing},
            Union{ODBC.API.SQLTime, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{Array{UInt8,1}, Missing},
            Union{String, Missing},
            Union{String, Missing},
            Union{String, Missing},
        )
        println("Postgres table 'test1':")
        show(data)
        println()

        @testset "Streaming postgres data to CSV" begin
            # Test exporting test1 to CSV
            temp_filename = "postgres_test1.csv"
            source = ODBC.Query(dsn, "select * from test1")
            CSV.write(temp_filename, source)

            open(temp_filename) do f
                @test readline(f) == (
                    "test_bigint,test_decimal,test_int,test_numeric,test_smallint," *
                    "test_float,test_real,test_money,test_date,test_timestamp,test_time," *
                    "test_char,test_varchar,test_bytea,test_boolean,test_text,test_array"
                )
                @test readline(f) == (
                    "1,1.2,2,1.4,3,1.6,1.8,2.0,2016-01-01,2016-01-01T01:01:01,01:01:01," *
                    "A,hey there sailor,,1,hey there abraham,\"{1,2,3}\""
                )
            end
            rm(temp_filename)
        end

        @testset "Exporting postgres data to SQLite" begin
            # Test exporting test1 to SQLite
            db = SQLite.DB()
            source = ODBC.Query(dsn, "select * from test1")
            SQLite.load!(source, db, "postgres_test1")

            data = SQLite.query(db, "select * from postgres_test1")
            @test size(data) == (1,17)
            @test data[1][1] === 1
            @test data[2][1] === 1.2
            @test data[9][1] === ODBC.API.SQLDate(2016,1,1)
            @test data[17][1] == "{1,2,3}"
        end

        ODBC.execute!(dsn, "drop table if exists test1")
    end

    ODBC.disconnect!(dsn)
end


text:
struct ElasticConstantSolver{T<:CrystalSystem}
    system::T
end

function (::ElasticConstantSolver{Cubic})(
    strains::AbstractVector{<:EngineeringStrain},
    stresses::AbstractVector{<:EngineeringStress},
)
    indices = map(_whichindex, strains)
    strain1 = _select(strains, indices, 1)
    stress1 = _select(stresses, indices, 1)
    c₁₁ = _cij(strain1[1], strain1[2], stress1[1], stress1[2])
    stress2 = _select(stresses, indices, 2)
    c₁₂ = _cij(strain1[1], strain1[2], stress2[1], stress2[2])
    strain4 = _select(strains, indices, 4)
    stress4 = _select(stresses, indices, 4)
    c₄₄ = _cij(strain4[1], strain4[2], stress4[1], stress4[2])
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        [
            c₁₁ c₁₂ c₁₂ 𝟎 𝟎 𝟎
            c₁₂ c₁₁ c₁₂ 𝟎 𝟎 𝟎
            c₁₂ c₁₂ c₁₁ 𝟎 𝟎 𝟎
            𝟎 𝟎 𝟎 c₄₄ 𝟎 𝟎
            𝟎 𝟎 𝟎 𝟎 c₄₄ 𝟎
            𝟎 𝟎 𝟎 𝟎 𝟎 c₄₄
        ],
    )
end

_select(v, indices, index) = v[filter(==(index), indices)]

_cij(ϵᵢ₊, ϵᵢ₋, σⱼ₊, σⱼ₋) = (σⱼ₊ - σⱼ₋) / (ϵᵢ₊ - ϵᵢ₋)

function _whichindex(x)
    indices = findall(!iszero, x)
    @assert length(indices) == 1
    return first(indices)
end


text:
module UW
using Mmap, SeisIO, SeisIO.FastIO, SeisIO.Quake
using Dates: DateTime

include("UW/imports.jl")
include("UW/uwdf.jl")
include("UW/uwpf.jl")
include("UW/uwevt.jl")
include("UW/desc.jl")

# exports
export formats, readuwevt, uwdf, uwdf!, uwpf, uwpf!

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

using .GeoStatsBase
import .GeoStatsBase: preprocess, solve_single

"""
    ImgQuilt(var₁=>param₁, var₂=>param₂, ...)

Image quilting simulation solver as described in Hoffimann et al. 2017.

## Parameters

### Required

* `TI`       - Training image
* `tilesize` - Tile size in x, y and z

### Optional

* `overlap`  - Overlap size in x, y and z (default to (1/6, 1/6, 1/6))
* `path`     - Simulation path (:raster (default), :dilation, or :random)
* `inactive` - Vector of inactive voxels (i.e. tuples (i,j,k)) in the grid
* `soft`     - A vector of `(data,dataTI)` pairs
* `tol`      - Initial relaxation tolerance in (0,1] (default to 0.1)

## Global parameters

### Optional

* `threads`      - Number of threads in FFT (default to number of physical CPU cores)
* `gpu`          - Whether to use the GPU or the CPU (default to false)
* `showprogress` - Whether to show or not the estimated time duration (default to false)
"""
@simsolver ImgQuilt begin
  @param TI
  @param tilesize
  @param overlap       = (1/6, 1/6, 1/6)
  @param path          = :raster
  @param inactive      = nothing
  @param soft          = []
  @param tol           = .1
  @global threads      = cpucores()
  @global gpu          = false
  @global showprogress = false
end

function preprocess(problem::SimulationProblem, solver::ImgQuilt)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  # sanity checks
  @assert pdomain isa RegularGrid "ImgQuilt requires RegularGrid domain"
  @assert ndims(pdomain) ∈ [2,3] "Number of dimensions must be 2 or 3"

  # result of preprocessing
  preproc = Dict{Symbol,Tuple}()

  for (var, V) in variables(problem)
    @assert var ∈ keys(solver.params) "Parameters for variable `$var` not found"

    # get user parameters
    varparams = solver.params[var]

    # add ghost dimension to simulation grid if necessary
    simsize = ndims(pdomain) == 2 ? (size(pdomain)..., 1) : size(pdomain)

    # create hard data object
    hdata = Dict{CartesianIndex{3},Real}()
    for (loc, datloc) in datamap(problem, var)
      push!(hdata, lin2cart(simsize, loc) => pdata[datloc,var])
    end

    # disable inactive voxels
    shape = Dict{CartesianIndex{3},Real}()
    if varparams.inactive ≠ nothing
      for icoords in varparams.inactive
        push!(shape, icoords => NaN)
      end
    end

    preproc[var] = (varparams, simsize, merge(hdata, shape))
  end

  preproc
end

function solve_single(problem::SimulationProblem, var::Symbol,
                      solver::ImgQuilt, preproc)
  # unpack preprocessed parameters
  par, simsize, hard = preproc[var]

  # run image quilting core function
  reals = iqsim(par.TI, par.tilesize, simsize;
                overlap=par.overlap, path=par.path,
                soft=par.soft, hard=hard, tol=par.tol,
                threads=solver.threads, gpu=solver.gpu,
                showprogress=solver.showprogress)

  # flatten result
  reals[1][:]
end


text:
if basename(pwd()) == "aoc"
    cd("2017/25")
end

function transpile(filename::AbstractString)
    """
    function turingmachine()
        cursor = 0
        tape = Set{Int}()
    """ *
    replace(readchomp(filename),
        r"Begin in state (.)\." => s"    state = '\1'",
        r"Perform a diagnostic checksum after (\d+) steps\." => s"    for _ in 1:\1",
        r"In state (.):" => s"if state == '\1'",
        "If the current value is 0:" => "      if !(cursor in tape)",
        "If the current value is 1:" => "      else",
        "- Write the value 1." => "        push!(tape, cursor)",
        "- Write the value 0." => "        delete!(tape, cursor)",
        "- Move one slot to the left." => "        cursor -= 1",
        "- Move one slot to the right." => "        cursor += 1",
        r"- Continue with state (.)\." => s"        state = '\1'",
        r"(?<!steps\.)\r?\n\r?\n" => "\n        end\n    else",
        r"(?!steps\.)\r?\n\r?\n" => "\n\n        "
    ) *
    "\n        end\n    end\n    \nend\n    length(tape)\nend"
end

part1(filename::AbstractString) = eval(Meta.parse(transpile(filename)))()

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
export beads_segment

function beads_segment(img::AbstractArray, threshold::Int64, maxsize = 1000, n_std = 3)
    filtered = zeros(eltype(img), axes(img))

    img_m = img
    m = mean(img)
    s = std(img)
    img_m[img_m.>m+n_std*s] .= 1
    img_m[img_m.!=1] .= 0

    label = label_components(img_m)
    out = component_pixels(label)
    out_fil = filter(out) do item
        length(item) >= threshold && length(item) <= maxsize
    end

    for vec in out_fil
        filtered[vec] = img[vec]
    end

    return filtered
end

function component_pixels(labels)
    n = maximum(labels)
    cp = [CartesianIndex{ndims(labels)}[] for _ = 1:n]
    for i in CartesianIndices(labels)
        c = labels[i]
        if c != 0
            push!(cp[c], i)
        end
    end
    return cp
end