text:
module TestView
using Compat
using Compat.Test
using CategoricalArrays

const ≅ = isequal

@testset "view($(CategoricalArray{Union{T, eltype(a)}}), $inds), ordered=$order construction" for
    T in (Union{}, Missing), order in (true, false),
    a in (1:10, 10:-1:1, ["a", "c", "b", "b", "a"]),
    inds in [1:2, :, 1, []]

    x = CategoricalArray{Union{T, eltype(a)}}(a, ordered=order)
    v = view(x, inds)
    @test levels(v) === levels(x)
    @test unique(v) == (ndims(v) > 0 ? unique(a[inds]) : [a[inds]])
    @test isordered(v) === isordered(x)
end

@testset "views comparison" begin
    ca1 = CategoricalArray([1, 2, 3])
    ca2 = CategoricalArray{Union{Int, Missing}}([1, 2, 3])
    ca3 = CategoricalArray([1, 2, missing])
    ca4 = CategoricalArray([4, 3, 2])
    ca5 = CategoricalArray([1 2; 3 4])

    @test view(ca1, 1:2) == view(ca1, 1:2)
    @test view(ca2, 1:2) == view(ca2, 1:2)
    @test view(ca1, 1:2) == view(ca2, 1:2)
    @test view(ca1, 1:2) == view(ca3, 1:2)
    @test view(ca1, 1:2) != view(ca2, 1:3)
    @test view(ca1, 1:2) != view(ca5, 1:2, 1:1)
end

@testset "fill! on view" for
    a in (categorical(["c", "a", "b"]), categorical(["c", "a", missing]))
    v = view(a, 1:2)
    @test fill!(v, a[1]) == ["c", "c"]
    @test fill!(v, "a") == ["a", "a"]
    @test fill!(v, "d") == ["d", "d"]
    if eltype(a) >: Missing
        @test fill!(v, missing) ≅ [missing, missing]
    else
        @test_throws MethodError fill!(v, missing)
    end
end

end


text:
# batched
"""
    sum_autobatch(v)

it's equivalent to sum(v, dims=1) except for the case
of a vector, where it will return a scalar instead of
a 1-element vector.
"""
sum_autobatch(v) = sum(v, dims=1)
sum_autobatch(v::Vector) = sum(v)

"""
    _batched_outer_prod!(R, vb, wb)

Efficiently performs the outer product R[i] .= vb[i] .* wb[i]'
along the batch dimension i, assuming that the batch dimension
is the last.
Internally uses the fact that R is a StridedView
"""
@inline function _batched_outer_prod!(R::StridedView, vb, wb)
    #@unsafe_strided R begin
        @inbounds @simd for i=1:size(R, 3)
            for j=1:size(wb, 1)
                for k=1:size(vb, 1)
                    R[k,j,i] = vb[k,i]*conj(wb[j,i])
                end
            end
        end
    #end

    #=@unsafe_strided R vb wb begin
        for i=1:size(R, 3)
            BLAS.ger!(1.0, vb[:,i], wb[:,i], R[:,:,i])
        end
    end=#
    return R
end

ext: jl


text:
module paqueteFourier

import Base: +, -, *, /

import Gadfly

# Notar que los métodos de +, -, *, / y diff se exportan solos
export trigPoly, innerProduct, nodes, vals, tpfplot

# El objeto básico de este paquete es un trigPoly, cuyo campo coefs contiene
# los coeficientes de la representación de un polinomio trigonométrico como
# combinación lineal de 1, exp(-im θ), exp(im θ), exp(-2im θ), exp(2im θ), …,
# en ese orden.
# Junto a la definición de este 'type' Julia automáticamente deja disponible un
# constructor, que es una función que recibe un vector de números complejos y
# devuelve un trigPoly con ese vector en el campo coefs.
# A continuación justificaré mi elección de ordenamiento para los coeficientes
# con la ayuda de una tabla.  Denoto al orden elegido (0, -1, 1, -2, 2, …) por
# A, al orden visto en clase (-N/2, -N/2+1, …, -1, 0, 1, …, N/2-2, N/2) por B y
# al orden por defecto de la FFT (0, 1, …, N/2-1, -N/2, -N/2+1, …, -1) por C (B
# y C son ligeramente distintos si N es impar, pero la tabla sigue valiendo).
#                                 ┌───┬───┬───┐
#                                 │ A │ B │ C │
# ┌───────────────────────────────┼───┼───┼───┤
# │                          Suma │ ✓ │ ✗ │ ✗ │
# │Interpretación de coeficientes │ ✓ │ ✗ │ ~ │
# │     Vector de números de onda │ ~ │ ✓ │ ✗ │
# │           Interpolación (FFT) │ ✗ │ ~ │ ✓ │
# ├───────────────────────────────┴───┴───┴───┤
# │   ✓ = Fácil, ~ = Moderado, ✗ = Difícil    │
# └───────────────────────────────────────────┘
"""
    trigPoly(x)

Polinomio trigonométrico que es combinación lineal de 1, exp(-im θ), exp(im θ), exp(-2im θ), … con los coeficientes indicados en `x`.
"""
type trigPoly
	coefs::Vector{Complex128}
end

"""
    innerProduct(f,g)

Producto interior L² entre los trigPoly `f` y `g`.
"""
function innerProduct(f::trigPoly, g::trigPoly)
	ac = 0.0 + 0.0im
	for k = 1:min(length(f.coefs),length(g.coefs))
		ac += f.coefs[k] * conj(g.coefs[k])
	end
	2π*ac
end

# Sin cargar ningún paquete Julia ya tiene una función llamada norm; la nuestra
# será un método más
"""
    norm(f::trigPoly)

Norma L² de `f`.
"""
Base.norm(f::trigPoly) = sqrt(real(innerProduct(f,f)))

"""
    nodes(n)

Arreglo de `n` nodos de Fourier.
"""
function nodes(n::Integer)
	2π*(0:n-1)/n
end

# Función que transforma un vector del orden C al orden A; Julia reemplazará T
# con el tipo de número que contenga el vector con el cual se alimente a esta
# función
function reordering{T}(v::Vector{T})
	n = length(v)
	out = Vector{T}(n) # Reserva de memoria
	floorhalfn = n÷2
	for k = 1:floorhalfn
		out[2*k] = v[n+1-k]
	end
	for k = 1:n-floorhalfn
		out[2*k-1] = v[k]
	end
	out
end

# Función que transforma un vector del orden A al orden C
function reverse_reordering{T}(v::Vector{T})
	n = length(v)
	floorhalfn = n÷2
	out = Vector{T}(n)
	for k = 1:floorhalfn
		out[n+1-k] = v[2*k]
	end
	for k = 1:n-floorhalfn
		out[k] = v[2*k-1]
	end
	out
end

# Interpolador que asume que los valores recibidos corresponden a la evaluación
# de una función en los nodos de Fourier
function interp{T}(v::Vector{T})
	modes = fft(v)/length(v)
	trigPoly(reordering(modes))
end

"""
    trigPoly(f, n)

Polinomio trigonométrico que interpola a la función `f` en `n` nodos de Fourier.
"""
function trigPoly(f::Function, n::Integer)
	vals = map(f, nodes(n)) # Vector con evaluaciones de f
	interp(vals)
end

"""
    vals(f)

Evaluación del polinomio trigonométrico `f` en tantos nodos de Fourier como coeficientes almacena `f`.
"""
function vals(f::trigPoly)
	ifft(reverse_reordering(f.coefs*length(f.coefs)))
end

"""
    vals(f, x)

Evaluación del polinomio trigonométrico `f` en el arreglo o escalar real `x`.
"""
function vals(f::trigPoly, x)
	coefsplus = f.coefs[3:2:end]
	coefsminus = f.coefs[2:2:end]
	# Podría escribir exp(im*x), pero si x es escalar Julia emite una
	# advertencia (no un error) que invita a usar la versión vectorizada de
	# exp; usando map x puede ser escalar o vectorial indistintamente
	expplusix = map(exp, im*x)
	expminusix = conj(expplusix)
	# Arreglo con ceros o escalar cero según lo que sea x
	outplus = zero(expplusix)
	outminus = zero(expplusix)
	for val in reverse(coefsplus)
		outplus = expplusix.*(outplus + val)
	end
	for val in reverse(coefsminus)
		outminus = expminusix.*(outminus + val)
	end
	outplus + outminus + f.coefs[1]
end

# Evaluador por defecto de un trigPoly; permite evaluar un trigPoly f en x con
# la sintaxis clásica f(x)
(f::trigPoly)(x) = vals(f, x)

# Suma
function +(f::trigPoly, g::trigPoly)
	fl = length(f.coefs)
	gl = length(g.coefs)
	retl = max(fl, gl)
	retcoefs = zeros(Complex128, retl)
	retcoefs[1:fl] = f.coefs
	retcoefs[1:gl] += g.coefs
	trigPoly(retcoefs)
end

# Cambio de signo
-(f::trigPoly) = trigPoly(-f.coefs)

# Resta
-(f::trigPoly, g::trigPoly) = f + (-g)

# Suma, resta, multiplicación y división por escalar
function +(f::trigPoly, a::Number)
	retcoefs = f.coefs
	retcoefs[1] += a
	trigPoly(retcoefs)
end
-(f::trigPoly, a::Number) = f + (-a)
*(f::trigPoly, a::Number) = trigPoly(f.coefs*a)
/(f::trigPoly, a::Number) = trigPoly(f.coefs/a)
*(a::Number, f::trigPoly) = f*a
+(a::Number, f::trigPoly) = f+a
-(a::Number, f::trigPoly) = (-f) + a

# Multiplicación
function *(f::trigPoly, g::trigPoly)
	ndof = length(f.coefs) + length(g.coefs)
	paddedf = trigPoly([f.coefs; zeros(Complex128,ndof-length(f.coefs))])
	paddedg = trigPoly([g.coefs; zeros(Complex128,ndof-length(g.coefs))])
	productvals = vals(paddedf) .* vals(paddedg)
	interp(productvals)
end

"""
    diff(f::trigPoly)

Derivada de `f`.
"""
function Base.diff(f::trigPoly)
	retcoefs = Vector{Complex128}(length(f.coefs))
	for k = 1:length(f.coefs)
		#             ↓     modo     ↓
		retcoefs[k] = (-1)^(k+1)*(k÷2) * im * f.coefs[k]
	end
	trigPoly(retcoefs)
end

"""
    diff(f::trigPoly, k::Int64)

`k`-ésima derivada de `f`.
"""
function Base.diff(f::trigPoly, k::Int64)
	@assert k>=0
	if k == 0
		return f
	elseif k == 1
		return diff(f)
	else
		return diff(f, k-1)
	end
end

"""
    tpfplot(f)
    tpfplot(f1, f2, ...)
    tpfplot(f1, f2, ..., labels=[label1, label2, ...])

Grafica a un polinomio trigonométrico (trigPoly) o función (Function) `f` o a la colección de polinomios trigonométricos o funciones `f1`, `f2`, ...; el argumento opcional `labels`, de aparecer, debe especificar nombres para las curvas ingresadas como un arreglo de cadenas de tipo `String`.

La figura obtenida se puede guardar mediante la sintaxis

    fh = tpfplot(...)
    Gadfly.draw(Gadfly.SVG("myplot.svg", 15Gadfly.cm, 9Gadfly.cm), fh)
"""
function tpfplot(vf::Vararg{Union{Function,trigPoly}}; labels=[])
	defaultNamesFlag = isempty(labels)
	layersReal = []
	layersImag = []
	layerLabels = []
	layerColors = []
	for i = 1:length(vf)
		c = ["blue" "green" "red" "deepskyblue" "magenta" "gold" "lightcoral" "black"][(i-1)%8+1]
		if defaultNamesFlag
			push!(labels, "Curve $(i)")
		end
		if typeof(vf[i]) <: trigPoly
			nsamples = max(150, 6*length(vf[i].coefs))
			x = nodes(nsamples)
			vfix = vals(trigPoly([vf[i].coefs;zeros(Complex128,nsamples-length(vf[i].coefs))]))
			# Agrego el punto de abscisa 2π
			x = [x;2π]
			vfix = [vfix;vfix[1]]
			# Ahora los puntos con abscisas en la colección de
			# nodos de Fourier que corresponde a vf[i]
			xn = nodes(length(vf[i].coefs))
			vfin = vals(vf[i])
			push!(layersReal, Gadfly.layer(x=x, y=real(vfix), Gadfly.Geom.line, Gadfly.Theme(default_color=c)))
			push!(layersReal, Gadfly.layer(x=xn, y=real(vfin), Gadfly.Geom.point, Gadfly.Theme(default_color=c)))
			push!(layersImag, Gadfly.layer(x=x, y=imag(vfix), Gadfly.Geom.line, Gadfly.Theme(default_color=c)))
			push!(layersImag, Gadfly.layer(x=xn, y=imag(vfin), Gadfly.Geom.point, Gadfly.Theme(default_color=c)))
			push!(layerColors, c); push!(layerLabels, labels[i])
		else
			x = linspace(0,2π,2000)
			vfix = map(vf[i], x)
			push!(layersReal, Gadfly.layer(x=x, y=real(vfix), Gadfly.Geom.line, Gadfly.Theme(default_color=c)))
			push!(layersImag, Gadfly.layer(x=x, y=imag(vfix), Gadfly.Geom.line, Gadfly.Theme(default_color=c)))
			push!(layerColors, c); push!(layerLabels, labels[i])
		end
	end
	labelSpecification = Gadfly.Guide.manual_color_key("", layerLabels, layerColors)
	Gadfly.hstack(
	       Gadfly.plot(layersReal..., Gadfly.Guide.XLabel("θ"), Gadfly.Guide.YLabel(""), Gadfly.Guide.Title("Real part"), Gadfly.Coord.cartesian(xmin=0, xmax=2π), Gadfly.Theme(key_position=:bottom), labelSpecification),
	       Gadfly.plot(layersImag..., Gadfly.Guide.XLabel("θ"), Gadfly.Guide.YLabel(""), Gadfly.Guide.Title("Imaginary part"), Gadfly.Coord.cartesian(xmin=0, xmax=2π), Gadfly.Theme(key_position=:bottom), labelSpecification)
	       )
end

end # module

max_stars_repo_name: lfiguero/paqueteFourier

max_issues_repo_name: lfiguero/paqueteFourier

max_forks_repo_name: lfiguero/paqueteFourier


text:
# TimeType-Year arithmetic
function (+)(y::YearDate, dy::Year)
    return YearDate(UTY(Year(value(y) + Dates.value(dy))))
end

function (+)(dy::Year, y::YearDate)
    return YearDate(UTY(Year(value(y) + Dates.value(dy))))
end

function (-)(y::YearDate, dy::Year)
    return YearDate(UTY(Year(value(y) - Dates.value(dy))))
end

function (+)(y::SemesterDate, dy::Year)
    return SemesterDate(UTS(value(y) + 2 * value(dy)))
end

function (+)(dy::Year, y::SemesterDate)
    return SemesterDate(UTS(value(y) + 2 * value(dy)))
end

function (-)(y::SemesterDate, dy::Year)
    return SemesterDate(UTS(value(y) - 2 * value(dy)))
end

function (+)(y::SemesterDate, dy::Semester)
    return SemesterDate(UTS(value(y) + value(dy)))
end

function (+)(dy::Semester, y::SemesterDate)
    return SemesterDate(UTS(value(y) + value(dy)))
end

function (-)(y::SemesterDate, dy::Semester)
    return SemesterDate(UTS(value(y) - value(dy)))
end

function (+)(y::QuarterDate, dy::Year)
    return QuarterDate(UTQ(value(y) + 4 * value(dy)))
end

function (+)(dy::Year, y::QuarterDate)
    return QuarterDate(UTQ(value(y) + 4 * value(dy)))
end

function (-)(y::QuarterDate, dy::Year)
    return QuarterDate(UTQ(value(y) - 4 * value(dy)))
end

function (+)(y::QuarterDate, dy::Semester)
    return QuarterDate(UTQ(value(y) + 2 * value(dy)))
end

function (+)(dy::Semester, y::QuarterDate)
    return QuarterDate(UTQ(value(y) + 2 * value(dy)))
end

function (-)(y::QuarterDate, dy::Semester)
    return QuarterDate(UTQ(value(y) - 2 * value(dy)))
end

function (+)(y::QuarterDate, dy::Quarter)
    return QuarterDate(UTQ(value(y) + value(dy)))
end

function (+)(dy::Quarter, y::QuarterDate)
    return QuarterDate(UTQ(value(y) + value(dy)))
end

function (-)(y::QuarterDate, dy::Quarter)
    return QuarterDate(UTQ(value(y) - value(dy)))
end

function (+)(y::MonthDate, dy::Year)
    return MonthDate(UTM(value(y) + 12 * value(dy)))
end

function (+)(dy::Year, y::MonthDate)
    return MonthDate(UTM(value(y) + 12 * value(dy)))
end

function (-)(y::MonthDate, dy::Year)
    return MonthDate(UTM(value(y) - 12 * value(dy)))
end

function (+)(y::MonthDate, dy::Semester)
    return MonthDate(UTM(value(y) + 6 * value(dy)))
end

function (+)(dy::Semester, y::MonthDate)
    return MonthDate(UTM(value(y) + 6 * value(dy)))
end

function (-)(y::MonthDate, dy::Semester)
    return MonthDate(UTM(value(y) - 6 * value(dy)))
end

function (+)(y::MonthDate, dy::Quarter)
    return MonthDate(UTM(value(y) + 3 * value(dy)))
end

function (+)(dy::Quarter, y::MonthDate)
    return MonthDate(UTM(value(y) + 3 * value(dy)))
end

function (-)(y::MonthDate, dy::Quarter)
    return MonthDate(UTM(value(y) - 3 * value(dy)))
end

function (+)(y::MonthDate, dy::Month)
    return MonthDate(UTM(value(y) + value(dy)))
end

function (+)(dy::Month, y::MonthDate)
    return MonthDate(UTM(value(y) + value(dy)))
end

function (-)(y::MonthDate, dy::Month)
    return MonthDate(UTM(value(y) - value(dy)))
end

function (+)(w::WeekDate, dy::Year)
    y = year(w)
    return WeekDate(UTW(value(w) + weeksinyear(y, y + value(dy))))
end

function (+)(dy::Year, w::WeekDate)
    y = year(w)
    return WeekDate(UTW(value(w) + weeksinyear(y, y + value(dy))))
end

function (-)(w::WeekDate, dy::Year)
    y = year(w)
    return WeekDate(UTW(value(w) - weeksinyear(y - value(dy), y)))
end

function (+)(w::WeekDate, dy::Semester)
    d1 = DayDate(UTD(Dates.WEEKEPOCH + 7 * value(w)))
    d2 = d1 + dy
    if dayofweek(d2) > 4
        return WeekDate(year(d2), week(d2)) + Week(1)
    else
        return WeekDate(year(d2), week(d2))
    end
end

function (+)(dy::Semester, w::WeekDate)
    return w + dy
end

function (-)(w::WeekDate, dy::Semester)
    d1 = DayDate(UTD(Dates.WEEKEPOCH + 7 * value(w)))
    d2 = d1 - dy
    if dayofweek(d2) > 4
        return WeekDate(year(d2), week(d2)) + Week(1)
    else
        return WeekDate(year(d2), week(d2))
    end
end

function (+)(w::WeekDate, dy::Quarter)
    d1 = DayDate(UTD(Dates.WEEKEPOCH + 7 * value(w)))
    d2 = d1 + dy
    if dayofweek(d2) > 4
        return WeekDate(year(d2), week(d2)) + Week(1)
    else
        return WeekDate(year(d2), week(d2))
    end
end

function (+)(dy::Quarter, w::WeekDate)
    return w + dy
end

function (-)(w::WeekDate, dy::Quarter)
    return WeekDate(UTW(value(w) - value(dy)))
end

function (+)(w::WeekDate, dy::Month)
    d1 = DayDate(UTD(Dates.WEEKEPOCH + 7 * value(w)))
    d2 = d1 + dy
    if dayofweek(d2) > 4
        return WeekDate(year(d2), week(d2)) + Week(1)
    else
        return WeekDate(year(d2), week(d2))
    end
end

function (+)(dy::Month, w::WeekDate)
    return w + dy
end

function (-)(w::WeekDate, dy::Month)
    return WeekDate(UTW(value(w) - value(dy)))
end

function (+)(w::WeekDate, dy::Week)
    # Date of Thursday week w
    d1 = DayDate(UTD(7 * (value(w) - 1) + 4))
    d2 = d1 + dy
    y2 = year(d2)
    w2 = week(d2)
    if dayofweek(d2) > 4
        if w2 == weeksinyear(y2)
            return WeekDate(y2 + 1, 1)
        else
            return WeekDate(y2, w2 + 1)
        end
    else
        return WeekDate(y2, w2)
    end
end

function (+)(dy::Week, w::WeekDate)
    return w + dy
end

function (-)(w::WeekDate, dy::Week)
    return WeekDate(UTW(value(w) - value(dy)))
end

function (+)(y::DayDate, dy::Year)
    yy, mm, dd = yearmonthday(y)
    y2 = yy + value(dy)
    ld = daysinmonth(y2, mm)
    return DayDate(y2, mm, (dd > ld) ? ld : dd)
end

function (+)(dy::Year, y::DayDate)
    return y + dy
end

function (-)(y::DayDate, dy::Year)
    yy, mm, dd = yearmonthday(y)
    y2 = yy - value(dy)
    ld = daysinmonth(y2, mm)
    return DayDate(y2, mm, (dd > ld) ? ld : dd)
end

function (+)(d::DayDate, ds::Semester)
    y1, m1, d1 = yearmonthday(d)
    dm = 6 * value(ds)
    y2 = Dates.yearwrap(y1, m1, dm)
    m2 = Dates.monthwrap(m1, dm)
    ld = daysinmonth(y2, m2)
    return DayDate(y2, m2, (d1 > ld) ? ld : d1)
end

function (+)(dy::Semester, y::DayDate)
    return y + dy
end

function (-)(d::DayDate, ds::Semester)
    y1, m1, d1 = yearmonthday(d)
    dm = 6 * value(ds)
    y2 = Dates.yearwrap(y1, m1, -6 * dm)
    m2 = Dates.monthwrap(m1, dm)
    ld = daysinmonth(y2, m2)
    return DayDate(y2, m2, (d1 > ld) ? ld : d1)
end

function (+)(d::DayDate, dq::Quarter)
    y1, m1, d1 = yearmonthday(d)
    dm = 3 * value(dq)
    y2 = Dates.yearwrap(y1, m1, dm)
    m2 = Dates.monthwrap(m1, dm)
    ld = daysinmonth(y2, m2)
    return DayDate(y2, m2, (d1 > ld) ? ld : d1)
end

function (+)(dy::Quarter, y::DayDate)
    return y + dy
end

function (-)(d::DayDate, dq::Quarter)
    y1, m1, d1 = yearmonthday(d)
    dm = 6 * value(dq)
    y2 = Dates.yearwrap(y1, m1, -dm)
    m2 = Dates.monthwrap(m1, dm)
    ld = daysinmonth(y2, m2)
    return DayDate(y2, m2, (d1 > ld) ? ld : d1)
end

function (+)(d::DayDate, dm::Month)
    y1, m1, d1 = yearmonthday(d)
    y2 = Dates.yearwrap(y1, m1, value(dm))
    m2 = Dates.monthwrap(m1, value(dm))
    ld = daysinmonth(y2, m2)
    return DayDate(y2, m2, (d1 > ld) ? ld : d1)
end

function (+)(dy::Month, y::DayDate)
    return y + dy
end

function (-)(d::DayDate, dm::Month)
    y1, m1, d1 = yearmonthday(d)
    y2 = Dates.yearwrap(y1, m1, -value(dm))
    m2 = Dates.monthwrap(m1, -value(dm))
    ld = daysinmonth(y2, m2)
    return DayDate(y2, m2, (d1 > ld) ? ld : d1)
end

function (+)(y::DayDate, dy::Week)
    return DayDate(UTD(value(y) + 7 * value(dy)))
end

function (+)(dy::Week, y::DayDate)
    return y + dy
end

function (-)(y::DayDate, dy::Week)
    return DayDate(UTD(value(y) - 7 * value(dy)))
end

function (+)(y::DayDate, dy::Day)
    return DayDate(UTD(value(y) + value(dy)))
end

function (+)(dy::Day, y::DayDate)
    return DayDate(UTD(value(y) + value(dy)))
end

function (-)(y::DayDate, dy::Day)
    return DayDate(UTD(value(y) - value(dy)))
end

function (+)(y::UndatedDate, dy::Undated)
    return UndatedDate(UTU(value(y) + value(dy)))
end

function (+)(dy::Undated, y::UndatedDate)
    return UndatedDate(UTU(value(y) + value(dy)))
end

function (-)(y::UndatedDate, dy::Undated)
    return UndatedDate(UTU(value(y) - value(dy)))
end


text:
"""
    dissipator(n::Int, l)

Construct the `n`-site vectorised dissipator corresponding to a single Lindblad
operator `l` acting on a two-dimensional one-local Hilbert space represented by ℂ^2. The
Lindblad operator argument `l` defaults to the Pauli-lowering operator (σ^{-}).
"""
function dissipator(n::Int; l=σm)
    L = Dict([])

    Llist = [one(l) for _ in 1:n]
    Llist[1] = l

    Ds = []

    for i in 1:n
        if n == 1
            L[1] = Llist[1]
        else
            L[i] = ⊗(circshift(Llist, i - 1)...)
        end
        push!(
            Ds,
            2 * (cj(L[i]) ⊗ L[i]) - one(L[i]) ⊗ (L[i]' * L[i]) -
            (cj(L[i])' * cj(L[i])) ⊗ one(L[i]),
        )
    end
    return sum(Ds)
end

@doc raw"""
    coherent(H::AbstractTensorMap)

Construct the `two`-local and vectorised operator `H` corresponding to the coherent part of
the Liouvillian super-operator. That is, given ``H_{\alpha, \beta}``,  constructs
```math
\mathcal{H}_{\alpha,\beta} = - i (\mathbb{I}_{\alpha, \beta} \otimes H_{\alpha,\beta}
- H^{T}_{\alpha,\beta} \otimes \mathbb{I}_{\alpha, \beta})
```
where ``\alpha`` and ``\beta`` are lattice indices corresponding to any two adjacent
lattice sites.
"""
coherent(H::AbstractTensorMap) = -im * (one(H) ⊗ H - cj(H)' ⊗ one(H))

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
module ParameterGradients

include("../src/Jaynes.jl")
using .Jaynes
using Gen

model = @jaynes (a₀, μ₀, λ₀, N) -> begin
    b₀ = learnable(:b₀)
    τ ~ Gamma(a₀, b₀)
    μ ~ Normal(μ₀, 1 / (λ₀ * τ))
    x = [{:x => i} ~ Normal(μ, 1 / τ) for i in 1 : N]
    x
end
init_param!(model, (:b₀, ), 5.0)

tr = simulate(model, (0.5, 3.0, 0.5, 50))
accumulate_param_gradients!(tr)
display(get_params_grad(model))

end # module


text:
@testset "init_log_table vectors" begin
    table = init_log_table(
        [:open_nodes],
        ["#Open"],
        [20],
        [:center],
        [2]
    )
    @test table.setup.ids == [:open_nodes]
end

@testset "init_log_table list of columns everything specified" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open", width=20, alignment=:center),
        (id=:closed_nodes, name="#Closed", width=30, alignment=:left, precision=3),
    )
    @test table.setup.ids == [:open_nodes, :closed_nodes]
    @test table.setup.names == ["#Open", "#Closed"]
    @test table.setup.widths == [20,30]
    @test table.setup.alignments == [:center,:left]
    @test table.setup.precisions == [2,3]
end

@testset "init_log_table list use defaults" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open", width=20, alignment=:center),
        (id=:closed_nodes, name="#Closed");
        width = 30,
        alignment = :left
    )
    @test table.setup.ids == [:open_nodes, :closed_nodes]
    @test table.setup.names == ["#Open", "#Closed"]
    @test table.setup.widths == [20,30]
    @test table.setup.alignments == [:center,:left]
end


text:
using VoronoiCells
using Base.Test
using GeometricalPredicates

@test (cells = poissonvoronoicells(10^3) ; true)

cells = poissonvoronoicells(10^2)

@test typeof(sarea(cells,1)) == Float64
@test typeof(area(cells[1])) == Float64
@test typeof(random_cell(cells)) == VoronoiCells.VoronoiCell
@test isvalid(findindex(cells,1.5,1.5))

# Test that any point in the bulk is covered by a cell.
function testrand(cells,n)
    c = 0
    s = 1.2
    f = 0.6
    for i in 1:n
        (x,y) = (f*rand()+s,f*rand()+s)
        if VoronoiCells.isexternal(cells, x,y) c += 1 end
    end
    (c,c/n)
end

@test testrand(poissonvoronoicells(10^5),10^5) == (0,0.0)

# Test mean area of cells in the bulk
function meanarea(cells::Array{VoronoiCell,1}, cutoff)
    ma = 0.0
    cnt = 0
    lc = 1 + cutoff
    hc = 2 - cutoff
    for i in 1:length(cells)
        c = cells[i]
        g = getgenerator(c)
        getx(g) < lc && continue
        getx(g) > hc && continue
        gety(g) < lc && continue
        gety(g) > hc && continue
        cnt += 1
        ma += VoronoiCells.area(cells[i])
    end
    ma / cnt
end

@test abs(meanarea(poissonvoronoicells(10^5)._cells,.1) / 1e-5 -1) < 1e-2


text:
"""
A handle to a cursor structure for navigating through a database.
"""
type Cursor
    handle::Ptr{Void}
    Cursor(cur::Ptr{Void}) = new(cur)
end

"Check if cursor is open"
isopen(cur::Cursor) = cur.handle != C_NULL

"Create a cursor"
function open(txn::Transaction, dbi::DBI)
    cur_ptr_ref = Ref{Ptr{Void}}(C_NULL)
    ret = ccall((:mdb_cursor_open, liblmdb), Cint,
                 (Ptr{Void}, Cuint, Ptr{Ptr{Void}}),
                  txn.handle, dbi.handle, cur_ptr_ref)
    (ret != 0) && throw(LMDBError(ret))
    return Cursor(cur_ptr_ref[])
end

"Wrapper of Cursor `open` for `do` construct"
function open(f::Function, txn::Transaction, dbi::DBI)
    cur = open(txn, dbi)
    try
        f(cur)
    finally
        close(cur)
    end
end

"Close a cursor"
function close(cur::Cursor)
    if cur.handle == C_NULL
        warn("Cursor is already closed")
    end
    ccall((:mdb_cursor_close, liblmdb), Void, (Ptr{Void},), cur.handle)
    cur.handle = C_NULL
    return
end

"Renew a cursor"
function renew(txn::Transaction, cur::Cursor)
    ret = ccall((:mdb_cursor_renew, liblmdb), Cint,
                 (Ptr{Void}, Ptr{Void}), txn.handle, cur.handle)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Return the cursor's transaction"
function txn(cur::Cursor)
    txn_ptr = ccall((:mdb_cursor_txn, liblmdb), Ptr{Void}, (Ptr{Void},), cur.handle)
    (txn_ptr == C_NULL) && return nothing
    return Transaction(txn_ptr)
end

"Return the cursor's database"
function dbi(cur::Cursor)
    dbi = ccall((:mdb_cursor_dbi, liblmdb), Cuint, (Ptr{Void},), cur.handle)
    (dbi == 0) && return nothing
    return DBI(dbi, "")
end

"""Retrieve by cursor.

This function retrieves key/data pairs from the database.
"""
function get{T}(cur::Cursor, key, ::Type{T}, op::CursorOps=FIRST)
    # Setup parameters
    mdb_key_ref = Ref(MDBValue(key))
    mdb_val_ref = Ref(MDBValue())

    # Get value
    ret = ccall( (:mdb_cursor_get, liblmdb), Cint,
                  (Ptr{Void}, Ptr{MDBValue}, Ptr{MDBValue}, Cint),
                   cur.handle, mdb_key_ref, mdb_val_ref, Cint(op))
    (ret != 0) && throw(LMDBError(ret))

    # Convert to proper type
    mdb_val = mdb_val_ref[]
    if T <: AbstractString
        return bytestring(convert(Ptr{UInt8}, mdb_val.data), mdb_val.size)
    else
        nvals = floor(Int, mdb_val.size/sizeof(T))
        value = pointer_to_array(convert(Ptr{T}, mdb_val.data), nvals)
        return length(value) == 1 ? value[1] : value
    end
end

"""Store by cursor.

This function stores key/data pairs into the database. The cursor is positioned at the new item, or on failure usually near it.
"""
function put!(cur::Cursor, key, val; flags::Cuint = zero(Cuint))
    mdb_key_ref = Ref(MDBValue(key))
    mdb_val_ref = Ref(MDBValue(val))

    ret = ccall((:mdb_cursor_put, liblmdb), Cint,
                 (Ptr{Void}, Ptr{MDBValue}, Ptr{MDBValue}, Cuint),
                  cur.handle, mdb_key_ref, mdb_val_ref, flags)

    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Delete current key/data pair to which the cursor refers"
function delete!(cur::Cursor; flags::Cuint = zero(Cuint))
    ret = ccall((:mdb_cursor_del, liblmdb), Cint,
                 (Ptr{Void}, Cuint), cur.handle, flags)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Return count of duplicates for current key"
function count(cur::Cursor)
    countp = Csize_t[0]
    ret = ccall( (:mdb_cursor_count, liblmdb), Cint,
                  (Ptr{Void}, Csize_t), cur.handle, countp)
    (ret != 0) && throw(LMDBError(ret))
    return Int(countp[1])
end


text:
using Test

import SurfaceFluxes
const SF = SurfaceFluxes
const UF = SF.UniversalFunctions

import CLIMAParameters
const CP = CLIMAParameters

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# TODO: Right now, we test these functions for
# type stability and correctness in the asymptotic
# limit. We may want to extend correctness tests.

@testset "UniversalFunctions" begin
    @testset "Type stability" begin
        FT = Float32
        ζ = FT(-2):FT(0.01):FT(200)
        for L in (-FT(10), FT(10))
            args = (L, param_set)
            for uf in (UF.Gryanik(args...), UF.Grachev(args...), UF.Businger(args...))
                for transport in (UF.MomentumTransport(), UF.HeatTransport())
                    ϕ = UF.phi.(uf, ζ, transport)
                    @test eltype(ϕ) == FT
                    ψ = UF.psi.(uf, ζ, transport)
                    @test eltype(ψ) == FT
                end
            end
        end

        # More type stability (phi/psi):
        FT = Float32
        ζ = (-FT(1), FT(0.5) * eps(FT), 2 * eps(FT))
        for L in (-FT(10), FT(10))
            args = (L, param_set)
            for uf in (UF.Gryanik(args...), UF.Grachev(args...), UF.Businger(args...))
                for transport in (UF.MomentumTransport(), UF.HeatTransport())
                    ϕ = UF.phi.(uf, ζ, transport)
                    @test eltype(ϕ) == FT
                    ψ = UF.psi.(uf, ζ, transport)
                    @test eltype(ψ) == FT
                end
            end
        end

        # More type stability (Psi):
        FT = Float32
        ζ = (-FT(1), -FT(0.5) * eps(FT), FT(0.5) * eps(FT), 2 * eps(FT))
        for L in (-FT(10), FT(10))
            args = (L, param_set)
            uf = UF.Businger(args...)
            for transport in (UF.MomentumTransport(), UF.HeatTransport())
                Ψ = UF.Psi.(uf, ζ, transport)
                @test eltype(Ψ) == FT
            end
        end
    end
    @testset "Conversions" begin
        FT = Float32
        ζ = FT(10)
        L = FT(10)
        args = (L, param_set)

        uf = UF.Gryanik(args...)
        @test UF.Businger(uf) isa UF.Businger
        @test UF.Grachev(uf) isa UF.Grachev

        uf = UF.Grachev(args...)
        @test UF.Businger(uf) isa UF.Businger
        @test UF.Gryanik(uf) isa UF.Gryanik

        uf = UF.Businger(args...)
        @test UF.Grachev(uf) isa UF.Grachev
        @test UF.Gryanik(uf) isa UF.Gryanik
    end
    @testset "Asymptotic range" begin
        FT = Float32

        ϕ_h_ζ∞(uf::UF.Grachev) = 1 + FT(UF.b_h(uf))
        ϕ_m_ζ∞(uf::UF.Grachev, ζ) = FT(UF.a_m(uf)) / FT(UF.b_m(uf)) * ζ^FT(1 / 3)

        ϕ_h_ζ∞(uf::UF.Gryanik) = FT(UF.Pr_0(uf)) * (1 + FT(UF.a_h(uf) / UF.b_h(uf)))
        ϕ_m_ζ∞(uf::UF.Gryanik, ζ) = FT(UF.a_m(uf) / UF.b_m(uf)^FT(2 / 3)) * ζ^FT(1 / 3)

        for L in (-FT(10), FT(10))
            args = (L, param_set)
            for uf in (UF.Grachev(args...), UF.Gryanik(args...))
                for ζ in FT(10) .^ (4, 6, 8, 10)
                    ϕ_h = UF.phi(uf, ζ, UF.HeatTransport())
                    @test isapprox(ϕ_h, ϕ_h_ζ∞(uf))
                end
                for ζ in FT(10) .^ (8, 9, 10)
                    ϕ_m = UF.phi(uf, ζ, UF.MomentumTransport())
                    @test isapprox(ϕ_m, ϕ_m_ζ∞(uf, ζ))
                end
            end
        end

    end

end


text:
# Create a random image
img = rand(UInt8 , 3, 500, 500)

# Test input as AbstractArray and cvtColor
img_gray = OpenCV.cvtColor(img, OpenCV.COLOR_RGB2GRAY)

@test size(img_gray, 1) == 1 && size(img_gray, 2) == size(img, 2) && size(img_gray, 3) == size(img, 3)

# Exception test
try
    # This should throw an error
    OpenCV.cvtColor(img_gray, OpenCV.COLOR_RGB2GRAY)
    exit(1)
catch
    # Error caught so we can continue
end

ve = view(img, :,200:300, 200:300)

# Auto-conversion from-to OpenCV types
ve_gray = OpenCV.cvtColor(ve, OpenCV.COLOR_RGB2GRAY)

# Shape check
@test size(ve_gray)[1] == 1 && size(img_gray)[1] == 1

print("imgproc test passed\n")


text:
using Agents
using Serialization

function restore_snapshot(filename = "snapshot100kr1.dat")
    space = GridSpace((k1,k2))
    covidmodel = ABM(CovidAgent.GenericAgent,space; scheduler = fastest)
    (agents, lochh, locf, unemplist, shorttimelist, empcount, unempcount, shorttimecount, oldcount, unemp, firms, hh, tau, divperhh, weeklyconsumption) = open(filename,"r") do snapshot_file
        deserialize(snapshot_file)
    end
    for (id, agent) in agents
        add_agent_pos!(agent,covidmodel)
    end
    return covidmodel, lochh, locf, unemplist, shorttimelist, empcount, unempcount, shorttimecount, oldcount, unemp, firms, hh, tau, divperhh, weeklyconsumption
end


text:
# unary
@i function NEG(a!::GVar)
    NEG(value(a!))
    NEG(grad(a!))
end

@i function CONJ(a!::GVar)
    CONJ(value(a!))
    CONJ(grad(a!))
end

# +-
@i function ⊖(identity)(a!::GVar, b::GVar)
    value(a!) ⊖ value(b)
    grad(b) ⊕ grad(a!)
end
@nograd ⊖(identity)(a!, b::GVar)
@nograd ⊖(identity)(a!::GVar, b)

# NOTE: it will error on `SWAP(a!::GVar, b)` or `SWAP(a!, b:GVar)`
@i function SWAP(a!::GVar, b!::GVar)
    SWAP(value(a!), value(b!))
    SWAP(grad(a!), grad(b!))
end

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

@i function ⊖(*)(out!::GVar, x, y::GVar)
    value(out!) -= value(x) * value(y)
    grad(y) += value(x) * grad(out!)
end

@i function ⊖(*)(out!::GVar, x::GVar, y)
    value(out!) -= value(x) * value(y)
    grad(x) += grad(out!) * value(y)
end

@i function ⊖(/)(out!::GVar{T}, x::GVar, y::GVar) where T
    value(out!) -= value(x)/value(y)
    @anc a1 = zero(T)
    @anc a2 = zero(T)
    grad(x) += grad(out!)/value(y)
    a1 += value(x)*grad(out!)
    a2 += a1/value(y)
    grad(y) -= a2/value(y)
    a2 -= a1/value(y)
    a1 -= value(x)*grad(out!)
end

@i function ⊖(/)(out!::GVar{T}, x, y::GVar) where T
    value(out!) -= x/value(y)
    @anc a1 = zero(T)
    @anc a2 = zero(T)
    a1 += x*grad(out!)
    a2 += a1/value(y)
    grad(y) -= a2/value(y)
    a2 -= a1/value(y)
    a1 -= x*grad(out!)
end

@i function ⊖(/)(out!::GVar, x::GVar, y)
    value(out!) -= value(x)/y
    grad(x) += grad(out!)/y
end

@i function ⊖(^)(out!::GVar{T}, x::GVar, n::GVar) where T
    ⊖(^)(value(out!), value(x), value(n))
    @anc anc1 = zero(T)
    @anc anc2 = zero(T)
    @anc jac = zero(T)

    # grad x
    @routine getjac begin
        n ⊖ 1
        anc1 += value(x)^value(n)
        n ⊕ 1
        jac += anc1 * value(n)
    end
    grad(x) += grad(out!) * jac
    ~@routine getjac

    # get grad of n
    @routine getnjac begin
        anc1 += log(value(x))
        anc2 += value(x) ^ value(n)
        jac += anc1*anc2
    end
    grad(n) += grad(out!) * jac
    ~@routine getnjac
end

@i function ⊖(^)(out!::GVar{T}, x::GVar, n) where T
    ⊖(^)(value(out!), value(x), n)
    @anc anc1 = zero(T)
    @anc anc2 = zero(T)
    @anc jac = zero(T)

    @routine getjac begin
        n ⊖ 1
        anc1 += value(x)^n
        n ⊕ 1
        jac += anc1 * n
    end
    grad(x) += grad(out!) * jac
    ~@routine getjac
end

@i function ⊖(^)(out!::GVar{T}, x, n::GVar) where T
    ⊖(^)(value(out!), x, value(n))
    @anc anc1 = zero(T)
    @anc anc2 = zero(T)
    @anc jac = zero(T)

    # get jac of n
    @routine getnjac begin
        anc1 += log(x)
        anc2 += x ^ value(n)
        jac += anc1*anc2
    end
    grad(n) += grad(out!) * jac
    ~@routine getnjac
end

@i function ⊖(abs)(out!::GVar, x::GVar{T}) where T
    value(out!) -= abs(value(x))
    if (x > 0, ~)
        grad(x) ⊕ grad(out!)
    else
        grad(x) ⊖ grad(out!)
    end
end

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


text:
# MonoidAlternative
# =================

# we don't define a Identity-only TypeClasses.neutral
# this could interfere with the more generic neutral definition within `Either` (every Const is neutral towards Identity)

# for convenience we forward Monoid definitions from the wrapped type
TypeClasses.combine(a::Identity, b::Identity) = Identity(a.value ⊕ b.value)

# Identity denotes success, and first success wins
TypeClasses.orelse(a::Identity, b::Identity) = a

# FunctorApplicativeMonad
# =======================

TypeClasses.pure(::Type{<:Identity}, a) = Identity(a)
TypeClasses.ap(f::Identity, a::Identity) = Identity(f.value(a.value))
# for convenience, Identity does not use convert, whatever monad is returned is valid, providing maximum flexibility.
TypeClasses.flatmap(f, a::Identity) = f(a.value)

# FlipTypes
# =========

TypeClasses.flip_types(a::Identity) = TypeClasses.map(Identity, a.value)


text:
module Match

using Base.Meta

export @match, @ismatch

include("matchutils.jl")
include("matchmacro.jl")

## Uncomment for debugging
# export unapply, unapply_array, gen_match_expr, subslicedim, getvars, getvar, arg1isa, joinexprs, let_expr, array_type_of, isexpr

end # module


text:
using ImageTransformations: center

"""
    mask_circle!(::AbstractMatrix, npix; fill=NaN)

In-place version of [`mask_circle`](@ref)
"""
function mask_circle!(arr::AbstractMatrix{T}, npix; fill = NaN) where {T <: AbstractFloat}
    yy = axes(arr, 1)
    xx = axes(arr, 2)
    yc, xc = center(arr)
    d = @. sqrt((xx' - xc)^2 + (yy - yc)^2)
    @. arr[d < npix] = T(fill)
    return arr
end

"""
    mask_circle(::AbstractMatrix, npix; fill=NaN)

Mask the inner-circle of an image with radius `npix` with value `fill`. Note that the input type must be compatible with the fill value's type.

# See Also
* [`mask_circle!`](@ref)
"""
mask_circle(arr::AbstractMatrix{<:AbstractFloat}, npix; fill = NaN) = mask_circle!(deepcopy(arr), npix, fill = fill)


text:
using Compat, Compat.Test
### Taken from the RCall tests.
### https://github.com/JuliaInterop/RCall.jl/blob/master/test/runtests.jl
hd = homedir()
pd = Pkg.dir()

if Compat.Sys.iswindows()
    Rhome = if haskey(ENV,"R_HOME")
        ENV["R_HOME"]
    else
        using WinReg
        WinReg.querykey(WinReg.HKEY_LOCAL_MACHINE, "Software\\R-Core\\R","InstallPath")
    end
    Rscript = joinpath(Rhome,"bin",Sys.WORD_SIZE==64 ? "x64" : "i386", "Rscript")
else
    Rscript = "Rscript"
end

libpaths = readlines(`$Rscript -e "writeLines(.libPaths())"`)
if VERSION < v"0.6.0"
    libpaths = map(chomp, libpaths)
end
### end stolen from RCall.

using RStan
using Distributions

# write your own tests here

function coveragetest()
    X = Normal()
    n = 50
    xs = rand(X, n)
    β0 = 2
    β1 = 0.75
    λ = exp.(2 + 0.75 * xs)
    ys = rand.(Poisson.(λ))
    model = "
    data {
        int n;
	    real x[n];
	    int y[n];
    }
    parameters {
        vector[2] beta;
    }
    model {
        beta ~ normal(0, 5);
        for (i in 1:n){
            y[i] ~ poisson(exp(beta[1] + beta[2] * x[i]));
        }
    }
    "
    data = Dict("n" => n, "x" => xs, "y" => ys)

    res = rstan(model, data, 1000, nchains = 4, pars = "beta")
    mat = res[:beta]
    quants0 = quantile(mat[:, 1], [0.25, 0.75])
    quants1 = quantile(mat[:, 2], [0.25, 0.75])

    test1 = β0 > quants0[1] && β0 < quants0[2]
    test2 = β1 > quants1[1] && β1 < quants1[2]
    test1, test2
end

function gencoverages()
    N = 100
    resmat = Matrix{Bool}(2, N)
    for i in 1:N
        resmat[:, i] .= coveragetest()
    end
    resmat
end

resmat = gencoverages()
coverage = mean(resmat, 2)

@testset "Poisson Interval Coverage Test" begin
    @test coverage[1] > 0.4 && coverage[1] < 0.6
    @test coverage[2] > 0.4 && coverage[2] < 0.6
end


text:
using DataFrames
using SoilPlantAirContinuum
using Photosynthesis
using LazyArtifacts
using PkgUtility
using Test

include("test_SPAC.jl")


text:
module Flips

using Random
using LightGraphs
using MetaGraphs
using SparseArrays
using DataStructures
using JSON
using Memoize

abstract type AbstractStat end

export
    Flips,
    Flip,
    DummyFlip,  # TODO: remove?
    AbstractFlip,
    update!,  # TODO: remove?
    Plan,
    IndexedGraph,
    recom,
    recom_until_step,
    reversible_recom,
    reversible_recom_until_step,
    random_mst,
    balanced_cuts,
    cut_assignment,
    traverse_mst,
    pychain

include("./indexed_graph.jl")
include("./plan.jl")
include("./flip.jl")
include("./constraints.jl")
include("./chain.jl")
include("./recom.jl")
include("./pychain.jl")
end # module


text:
function compute()
  ans=""
	for n=2:9
		a=10^(floor(Int,9/n))

		for i =1:a
			s=""
			for j=1:n
				s=s*string(i*j)

			end
			if Set(s)==Set("123456789") && length(s)==9
				ans=max(s,ans)
			end
		end
	end

	return ans
end

println(compute())


text:
p = random_micrm_params(2,2,0.5)

#convert to micrm_params
p_struct = micrm_params(p)

u0 = ones(4)
tspan = (0.0,1.0)
prob = ODEProblem(MiCRM.dx!,u0,tspan,p_struct)
sol = solve(prob)


text:
# Train an normalizing flow for seismic denosing
# Authors: Ali Siahkoohi, alisk@gatech.edu
#          Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
#          Philipp Witte, pwitte3@gatech.edu
# Date: September 2020
"""
Example:
julia train_hint_denoising.jl --batchsize 256 --max_epoch 5
"""

using DrWatson
@quickactivate :FastApproximateInference

using InvertibleNetworks
using JLD
using JLD2
using Random
using Statistics
using ArgParse
using BSON
using Printf
using Flux: gpu, cpu
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, ExpDecay, update!, ADAM

# Random seed
Random.seed!(19)

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 50
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-3
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 1
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 4
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 64
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 12
    "--sigma"
        help = "noise std"
        arg_type = Float32
        default = 1f-1
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "seismic-denoising-256x256-0.3CS"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
sigma = parsed_args["sigma"]
sim_name = parsed_args["sim_name"]

# Define raw data directory
mkpath(datadir("training-data"))
data_path = datadir("training-data", "seismic_samples_256_by_256_num_10k.jld")
label_path = datadir("training-data", "seismic_samples_256_by_256_num_10k_labels.jld")

# Download the dataset into the data directory if it does not exist
if isfile(data_path) == false
    run(`wget https://www.dropbox.com/s/vapl62yhh8fxgwy/'
        'seismic_samples_256_by_256_num_10k.jld -q -O $data_path`)
end
if isfile(label_path) == false
    run(`wget https://www.dropbox.com/s/blxkh6tdszudlcq/'
        'seismic_samples_256_by_256_num_10k_labels.jld -q -O $label_path`)
end

# Load seismic images and create training and testing data
JLD.jldopen(data_path, "r") do file
    X_orig = read(file, "X")

    labels = JLD.jldopen(label_path, "r")["labels"][:]
    idx = findall(x -> x == 1, labels)

    X_orig = X_orig[:, :, :, idx]
    global nx, ny, nc, nsamples = size(X_orig)

    # Whiten data
    X_orig = whiten_input(X_orig)

    # Forward operator
    global A = RombergOP(nx, ny, nc, 2f0)

    # Create observed data
    Y_orig = A*X_orig + sigma*randn(Float32, size(X_orig))

    # Split in training/testing
    global ntrain = Int(floor((nsamples*.9)))
    global train_idx = randperm(nsamples)[1:ntrain]

    # Save data for testing phase
    save_dict = @strdict train_idx Y_orig A
    @tagsave(
        datadir("testing-data", "noisy_seismic_samples_256_by_256_num_10k.jld2"),
        save_dict
    )

    # Dimensions after wavelet squeeze to increase no. of channels
    global nx = Int(nx/2)
    global ny = Int(ny/2)
    global n_in = Int(nc*4)
end

# Create network
CH = NetworkConditionalHINT(
    n_in, n_hidden, depth, k1=3, k2=3, p1=1, p2=1
) |> gpu

# Training
# Batch extractor
num_batches = cld(ntrain, batchsize)

X_train = JLD.jldopen(data_path, "r")["X"]
labels = JLD.jldopen(label_path, "r")["labels"][:]
idx = findall(x -> x == 1, labels)
X_train = X_train[:, :, :, idx]

Y_train = JLD2.jldopen(
    datadir("testing-data", "noisy_seismic_samples_256_by_256_num_10k.jld2"),
    "r"
)["Y_orig"]

train_loader = DataLoader(train_idx, batchsize=batchsize, shuffle=true)

# Optimizer
opt = Optimiser(ExpDecay(lr, .9f0, num_batches*lr_step, 0f0), ADAM(lr))

# Training log keeper
fval = zeros(Float32, num_batches*max_epoch)

for epoch=1:max_epoch
    for (itr, idx) in enumerate(train_loader)

        # Apply wavelet squeeze (change dimensions to -> n1/2 x n2/2 x nc*4 x ntrain)
        X = wavelet_squeeze(whiten_input(X_train[:, :, :, idx]))
        Y = wavelet_squeeze(Y_train[:, :, :, idx])

        X = X |> gpu
        Y = Y |> gpu

        fval[(epoch-1)*num_batches + itr] = loss_supervised(CH, X, Y)[1]

        @printf(
            "Epoch: [%d/%d] | Itr: [%d/%d] | Average negative log-likelihood: %4.2f\r",
            epoch, max_epoch, itr, num_batches, fval[(epoch-1)*num_batches + itr]
        )

        # Update params
        for p in get_params(CH)
            update!(opt, p.data, p.grad)
        end
        clear_grad!(CH)
    end

    # Saving parameters and logs
    Params = get_params(CH) |> cpu
    save_dict = @strdict epoch max_epoch lr lr_step batchsize n_hidden depth sigma sim_name A Params fval
    @tagsave(
        datadir(sim_name, savename(save_dict, "jld2"; digits=6)),
        save_dict;
        safe=true
    )

end


text:
push!(LOAD_PATH, pwd())
using QDLDL

c_float = Float64
c_int   = Int64
c_bool  = Bool

function qdMatExample()

    i = [1, 7, 2, 3, 6, 10, 2, 3, 10, 4, 8, 5, 2, 6, 1, 7, 9, 4, 8, 7, 9, 2, 3, 10]
    j = [1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 5, 6, 6, 7, 7, 7, 8, 8, 9, 9, 10, 10, 10]
    s = [1.000000e+00, 3.503210e-01, 4.606408e-01, -1.211894e-01, -2.900576e-02, 1.824522e-01, -1.211894e-01, 4.179283e-01, -1.565056e+00, 1.778279e-01, -8.453948e-02, 1.000000e-01, -2.900576e-02, -1.000000e+00, 3.503210e-01, -4.410924e-01, 1.786632e-01, -8.453948e-02, -3.162278e-01, 1.786632e-01, -2.990768e-01, 1.824522e-01, -1.565056e+00, -1.000000e-01]

    A = sparse(i,j,s)
end

A = qdMatExample()
b = randn(c_float,A.n)
println("Example (backslash)                : ",  norm(qdldl(A)\b - A\b))

F = qdldl(A)
println("Example (with solve())             : ",  norm(solve(F,b) - A\b))

x = copy(b)
solve!(F,x)
println("Example (in place with solve!())   : ",  norm(x - A\b))

F = qdldl(A,perm=collect(A.n:-1:1))
println("Example (User permutation)         : ",  norm(F\b - A\b))

F = qdldl(A,perm=nothing)
println("Example (No permutation)           : ",  norm(F\b - A\b))

#compute a logical factorisation only
Flog = qdldl(A,logical = true)
Fnum = qdldl(A,logical = true)

logicalCheck = all((full(Fnum.L) .!= 0) .== (full(Flog.L) .!= 0))
println("Logical factorisation works?: ", logicalCheck)

max_issues_repo_name: WatsonZhouAnda/Cheetah-Software

max_forks_repo_name: WatsonZhouAnda/Cheetah-Software


text:
"""
    animate(vis::MechanismVisualizer,
            times::Vector{Float64},
            configurations::Vector{Vector{Float64}};
            fps::Float64=60, realtimerate::Float64=1.)

Animate the given mechanism passing through a time-coded series of
configurations by linearly interpolating the configuration vectors.
"""
function animate(vis::MechanismVisualizer,
                 times::Vector{Float64},
                 configurations::AbstractVector{<:AbstractVector{Float64}};
                 fps::Float64 = 60., realtimerate::Float64 = 1.)
    @assert fps > 0
    @assert 0 < realtimerate < Inf

    state = vis.state
    interpolated_configurations = interpolate((times,), configurations, Gridded(Linear()))
    t0, tf = first(times), last(times)
    framenum = 0
    walltime0 = time()
    @throttle framenum while true
        t = min(tf, t0 + (time() - walltime0) * realtimerate)
        q = interpolated_configurations(t)
        set_configuration!(state, q)
        rbd.normalize_configuration!(state)
        _render_state!(vis)
        framenum += 1
        t == tf && break
    end max_rate = fps
end

function MeshCat.setanimation!(mvis::MechanismVisualizer,
                      times::AbstractVector{<:Real},
                      configurations::AbstractVector{<:AbstractVector{<:Real}};
                      fps::Integer=30,
                      play::Bool=true,
                      repetitions::Integer=1)
    q0 = copy(configuration(state(mvis)))
    @assert axes(times) == axes(configurations)
    interpolated_configurations = interpolate((times,), configurations, Gridded(Linear()))
    animation = Animation()
    num_frames = floor(Int, (times[end] - first(times)) * fps)
    for frame in 0:num_frames
        time = first(times) + frame / fps
        set_configuration!(state(mvis), interpolated_configurations(time))
        atframe(animation, visualizer(mvis), frame) do frame_visualizer
            _render_state!(MechanismVisualizer(state(mvis), frame_visualizer))
        end
    end
    setanimation!(visualizer(mvis), animation, play=play, repetitions=repetitions)
    set_configuration!(state(mvis), q0)
end


text:
using Pkg; Pkg.activate(@__DIR__);
include("BoxTrees.jl")
using .BoxTrees
bt = BoxTrees

AbstractPlotting.inline!(true)

f( x :: Vector{ R } where R<:Real ) = [
    1 - 1.2 * x[1]^2 + x[2]/5;
    5 * 0.2 * x[1]
]

#%%
root = bt.init_box_tree(fill(-2.0, 2), fill(2.0,2) );
box_collection = [ root, ];

function progress!( box_collection, f :: F where F<:Function )

    # Subdivision Step
    new_collection = bt.Box[];
    for box ∈ box_collection
        bt.subdivide!(box; method=:cycle)
        push!( new_collection, bt.children(box)... )
    end

    # Selection Step

    ## 1) generate test sites and evaluate
    for box ∈ new_collection
        bt.generate_sites!( box; method = :edges_center, num_edge_points = 3 )
        bt.eval!( box, f )
    end

    ## 2) Discrete Hit Test
    good_box_indices = Int[];
    for (i,boxᵢ) ∈ enumerate(new_collection)
        for boxⱼ ∈ new_collection
            intersection_empty = true
            for fₓ ∈ bt.results( boxⱼ )
                if bt.contains( boxᵢ, fₓ )
                    intersection_empty = false;
                    break;
                end
            end
            # box was hit …
            if !intersection_empty
                push!(good_box_indices, i);
                break;
            end
        end
    end

    # TODO: Make more robust by reintroducing boxes

    empty!( box_collection );
    push!( box_collection, new_collection[good_box_indices]... );
end

#%%
for i = 1 : 10
    progress!( box_collection, f )
end

#%%
# Plotting …
# should work in 3D and 2D
using Makie;
using GeometryBasics;

function box2rect( box )
    rect = Rect( Vec(bt.lb( box )... ), Vec( bt.width(box)... ) );
    rect_v = decompose( Point{ bt.n_vars( box), Float64 }, rect )
    rect_f = decompose( TriangleFace{Int}, rect );
    return rect_v, rect_f
end

scene, layout = Makie.layoutscene();
ax = layout[1,1] = Axis(scene);
for box ∈ box_collection
    Makie.mesh!(ax, box2rect(box)...; color=:blue, strokecolor=:lightblue, strokewidth = 5.0, shading=false)
end
ax.aspect = AxisAspect(1)
scene

#%%


text:
mutable struct AttachmentDescription
    flags::vk.VkAttachmentDescriptionFlags
    format::vk.VkFormat
    samples::vk.VkSampleCountFlagBits
    loadOp::vk.VkAttachmentLoadOp
    storeOp::vk.VkAttachmentStoreOp
    stencilLoadOp::vk.VkAttachmentLoadOp
    stencilStoreOp::vk.VkAttachmentStoreOp
    initialLayout::vk.VkImageLayout
    finalLayout::vk.VkImageLayout

    mHandleRef::Ref{vk.VkAttachmentDescription}

    function AttachmentDescription()
        this = new()
        this.flags = 0
        return this
    end
end
@class AttachmentDescription

function createWithDepth24Stencil8(::Type{AttachmentDescription})::AttachmentDescription
    result = AttachmentDescription()
    result.format = vk.VK_FORMAT_D24_UNORM_S8_UINT #format::VkFormat
    result.samples = vk.VK_SAMPLE_COUNT_1_BIT #samples::VkSampleCountFlagBits
    result.loadOp = vk.VK_ATTACHMENT_LOAD_OP_LOAD #loadOp::VkAttachmentLoadOp
    result.storeOp = vk.VK_ATTACHMENT_STORE_OP_STORE #storeOp::VkAttachmentStoreOp
    result.stencilLoadOp = vk.VK_ATTACHMENT_LOAD_OP_LOAD #stencilLoadOp::VkAttachmentLoadOp
    result.stencilStoreOp = vk.VK_ATTACHMENT_STORE_OP_STORE #stencilStoreOp::VkAttachmentStoreOp
    result.initialLayout = vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL #initialLayout::VkImageLayout
    result.finalLayout = vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL #finalLayout::VkImageLayout
    return result;
end

function createWithDepth16(::Type{AttachmentDescription})::AttachmentDescription
    result = AttachmentDescription()
    result.format = vk.VK_FORMAT_D16_UNORM #format::VkFormat
    result.samples = vk.VK_SAMPLE_COUNT_1_BIT #samples::VkSampleCountFlagBits
    result.loadOp = vk.VK_ATTACHMENT_LOAD_OP_CLEAR #loadOp::VkAttachmentLoadOp
    result.storeOp = vk.VK_ATTACHMENT_STORE_OP_STORE #storeOp::VkAttachmentStoreOp
    result.stencilLoadOp = vk.VK_ATTACHMENT_LOAD_OP_DONT_CARE #stencilLoadOp::VkAttachmentLoadOp
    result.stencilStoreOp = vk.VK_ATTACHMENT_LOAD_OP_DONT_CARE #stencilStoreOp::VkAttachmentStoreOp
    result.initialLayout = vk.VK_IMAGE_LAYOUT_UNDEFINED #initialLayout::VkImageLayout
    result.finalLayout = vk.VK_IMAGE_LAYOUT_DEPTH_READ_ONLY_STENCIL_ATTACHMENT_OPTIMAL #finalLayout::VkImageLayout
    return result
end

function createWithDepth32float(::Type{AttachmentDescription})::AttachmentDescription
    result = AttachmentDescription()
    result.format = vk.VK_FORMAT_D32_SFLOAT #format::VkFormat
    result.samples = vk.VK_SAMPLE_COUNT_1_BIT #samples::VkSampleCountFlagBits
    result.loadOp = vk.VK_ATTACHMENT_LOAD_OP_LOAD #loadOp::VkAttachmentLoadOp
    result.storeOp = vk.VK_ATTACHMENT_STORE_OP_STORE #storeOp::VkAttachmentStoreOp
    result.stencilLoadOp = vk.VK_ATTACHMENT_LOAD_OP_LOAD #stencilLoadOp::VkAttachmentLoadOp
    result.stencilStoreOp = vk.VK_ATTACHMENT_STORE_OP_STORE #stencilStoreOp::VkAttachmentStoreOp
    result.initialLayout = vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL #initialLayout::VkImageLayout
    result.finalLayout = vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL #finalLayout::VkImageLayout
    return result
end

 function createWithColor(::Type{AttachmentDescription}, vkFormat::vk.VkFormat)::AttachmentDescription
    result = AttachmentDescription()
    result.format = vkFormat #format::VkFormat
    result.samples = vk.VK_SAMPLE_COUNT_1_BIT #samples::VkSampleCountFlagBits
    result.loadOp = vk.VK_ATTACHMENT_LOAD_OP_LOAD #loadOp::VkAttachmentLoadOp
    result.storeOp = vk.VK_ATTACHMENT_STORE_OP_STORE #storeOp::VkAttachmentStoreOp
    result.stencilLoadOp = vk.VK_ATTACHMENT_LOAD_OP_DONT_CARE #stencilLoadOp::VkAttachmentLoadOp
    result.stencilStoreOp = vk.VK_ATTACHMENT_STORE_OP_DONT_CARE #stencilStoreOp::VkAttachmentStoreOp
    result.initialLayout = vk.VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL #initialLayout::VkImageLayout
    result.finalLayout = vk.VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL #finalLayout::VkImageLayout
    return result
end

function clear(this::AttachmentDescription)
    this.loadOp = vk.VK_ATTACHMENT_LOAD_OP_CLEAR
    this.stencilLoadOp = vk.VK_ATTACHMENT_LOAD_OP_CLEAR
    this.initialLayout = vk.VK_IMAGE_LAYOUT_UNDEFINED
end

function discard(this::AttachmentDescription)
    this.storeOp = vk.VK_ATTACHMENT_STORE_OP_DONT_CARE
    this.stencilStoreOp = vk.VK_ATTACHMENT_STORE_OP_DONT_CARE
    if aspectsOf(this.format) & vk.VK_IMAGE_ASPECT_COLOR_BIT != 0
        this.finalLayout = vk.VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL
    elseif aspectsOf(this.format) & vk.VK_IMAGE_ASPECT_DEPTH_BIT != 0
        this.finalLayout = vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL
    end
end

function finalLayout(this::AttachmentDescription, layout::vk.VkImageLayout)
    this.finalLayout = layout
end

function getFormat(this::AttachmentDescription)::vk.VkFormat
    return this.mHandleRef[].format
end

function commit(this::AttachmentDescription)
    this.mHandleRef = Ref(vk.VkAttachmentDescription(
                            this.flags, #flags::VkAttachmentDescriptionFlags
                            this.format, #format::VkFormat
                            this.samples, #samples::VkSampleCountFlagBits
                            this.loadOp, #loadOp::VkAttachmentLoadOp
                            this.storeOp, #storeOp::VkAttachmentStoreOp
                            this.stencilLoadOp, #stencilLoadOp::VkAttachmentLoadOp
                            this.stencilStoreOp, #stencilStoreOp::VkAttachmentStoreOp
                            this.initialLayout, #initialLayout::VkImageLayout
                            this.finalLayout #finalLayout::VkImageLayout
    ))
end

function handleRef(this::AttachmentDescription)::Ref{vk.VkAttachmentDescription}
    return this.mHandleRef
end

# TODO
# finalLayout_Attachment
# finalLayout_ShaderRead
# finalLayout_TransferSrc
# finalLayout_DepthStencilRead
# finalLayout_PresentSrc