text:
include("../optimizers/trustregion.jl")
include("../utils/functions.jl")
using Gadfly

xvals = trust_region([3.;4.], 6, 1, 0.1, fenton, fenton_g, fenton_h,
    2000, "dogleg");
cvals = trust_region([3.;4.], 6, 1, 0.1, fenton, fenton_g, fenton_h,
    2000, "cg_steihaug");

println(xvals[end])

nsamps = length(xvals)
nsamps2 = length(cvals)

fx = [fenton(xvals[i]) for i in 1:nsamps]
fx2 = [fenton(cvals[i]) for i in 1:nsamps2]

Gadfly.plot(layer(x=1:nsamps, y=fx, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=fx2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=fx3, Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("f(x)"), Guide.title("Value of function"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

nsamps = length(xvals)

grads = [norm(fenton_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(fenton_g(cvals[i]), 2) for i in 1:nsamps2]
#grads3 = [norm(fenton_g(cvals[i]), 2) for i in 1:nsamps3]

Gadfly.plot(
layer(x=1:nsamps-2, y=grads[2:nsamps-1,:]./grads[1:nsamps-2,:], Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2-2, y=grads2[2:nsamps2-1,:]./grads2[1:nsamps2-2,:], Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3-1, y=grads3[2:nsamps3,:]./grads3[1:nsamps3-1,:], Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("gradient norm ratios"), Guide.title("gradient norm ratios"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

Gadfly.plot(layer(x=1:nsamps, y=grads, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=grads2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=grads3, Geom.line, Theme(default_color=color("orange"))),
    Guide.xlabel("iteration"), Guide.ylabel("gradient norm"), Guide.title("gradient norms"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
    Scale.x_log10, Scale.y_log10)

xvals = trust_region(randn(100), 6, 3, 0.1, rosenbrock,
    rosenbrock_g, rosenbrock_h, 1000, "dogleg");
cvals = trust_region(randn(100), 6, 3, 0.1, rosenbrock,
    rosenbrock_g, rosenbrock_h, 1000, "cg_steihaug");

nsamps = length(xvals)
nsamps2 = length(cvals)

func = rosenbrock
func_g = rosenbrock_g

fx = [func(xvals[i]) for i in 1:nsamps]
fx2 = [func(cvals[i]) for i in 1:nsamps2]

Gadfly.plot(layer(x=1:nsamps, y=fx, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=fx2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=fx3, Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("f(x)"), Guide.title("Value of function"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"],
    ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

nsamps = length(xvals)
nsamps2 = length(cvals)

grads = [norm(func_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps2]
#grads3 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps3]

Gadfly.plot(
layer(x=1:nsamps-2, y=grads[2:nsamps-1,:]./grads[1:nsamps-2,:], Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2-2, y=grads2[2:nsamps2-1,:]./grads2[1:nsamps2-2,:], Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3-1, y=grads3[2:nsamps3,:]./grads3[1:nsamps3-1,:], Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("gradient norm ratios"), Guide.title("gradient norm ratios"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

Gadfly.plot(layer(x=1:nsamps, y=grads, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=grads2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=grads3, Geom.line, Theme(default_color=color("orange"))),
    Guide.xlabel("iteration"), Guide.ylabel("gradient norm"), Guide.title("gradient norms"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
    Scale.x_log10, Scale.y_log10)

xvals = trust_region(ones(50)*10, 6, 3, 0.1, cute, cute_g, cute_h, 2000, "dogleg");
cvals = trust_region(ones(50)*10, 6, 3, 0.1, cute, cute_g, cute_h, 2000, "cg_steihaug");

nsamps = length(xvals)
nsamps2 = length(cvals)

func = cute
func_g = cute_g

fx = [func(xvals[i]) for i in 1:nsamps]
fx2 = [func(cvals[i]) for i in 1:nsamps2]

Gadfly.plot(layer(x=1:nsamps, y=fx, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=fx2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=fx3, Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("f(x)"), Guide.title("Value of function"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

nsamps = length(xvals)

grads = [norm(func_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps2]
#grads3 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps3]

Gadfly.plot(
layer(x=1:nsamps-1, y=grads[2:nsamps,:]./grads[1:nsamps-1,:], Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2-1, y=grads2[2:nsamps2,:]./grads2[1:nsamps2-1,:], Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3-1, y=grads3[2:nsamps3,:]./grads3[1:nsamps3-1,:], Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("gradient norm ratios"), Guide.title("gradient norm ratios"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

Gadfly.plot(layer(x=1:nsamps, y=grads, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=grads2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=grads3, Geom.line, Theme(default_color=color("orange"))),
    Guide.xlabel("iteration"), Guide.ylabel("gradient norm"), Guide.title("gradient norms"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
    Scale.x_log10, Scale.y_log10)

xvals = trust_region(ones(50)*10, 6, 3, 0.1, cute2, cute2_g, cute2_h, 2000, "dogleg");
cvals = trust_region(ones(50)*10, 6, 3, 0.1, cute2, cute2_g, cute2_h, 2000, "cg_steihaug");

nsamps = length(xvals)
nsamps2 = length(cvals)

func = cute2
func_g = cute2_g

fx = [func(xvals[i]) for i in 1:nsamps]
fx2 = [func(cvals[i]) for i in 1:nsamps2]

Gadfly.plot(layer(x=1:nsamps, y=fx, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=fx2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=fx3, Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("f(x)"), Guide.title("Value of function"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

nsamps = length(xvals)

grads = [norm(func_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps2]
#grads3 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps3]

Gadfly.plot(
layer(x=1:nsamps-1, y=grads[2:nsamps,:]./grads[1:nsamps-1,:], Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2-1, y=grads2[2:nsamps2,:]./grads2[1:nsamps2-1,:], Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3-1, y=grads3[2:nsamps3,:]./grads3[1:nsamps3-1,:], Geom.line, Theme(default_color=color("orange"))),
Guide.xlabel("iteration"), Guide.ylabel("gradient norm ratios"), Guide.title("gradient norm ratios"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
Scale.x_log10, Scale.y_log10)

Gadfly.plot(layer(x=1:nsamps, y=grads, Geom.line, Theme(default_color=color("blue"))),
layer(x=1:nsamps2, y=grads2, Geom.line, Theme(default_color=color("red"))),
#layer(x=1:nsamps3, y=grads3, Geom.line, Theme(default_color=color("orange"))),
    Guide.xlabel("iteration"), Guide.ylabel("gradient norm"), Guide.title("gradient norms"),
Guide.manual_color_key("Legend", ["Newton dogleg", "Steihaug-CG", "quasi-Newton SR1"], ["blue", "red", "orange"]),
    Scale.x_log10, Scale.y_log10)


text:
# Calc max G(X \ {x}) for the still selected ones
function max_GXminus1(Xstrs::Vector{ASCIIString}, selected, selectedindices, compressor::Compressor)
  maxGxminus1 = -Inf
  largest = 0 # dummy value for now
  for i in selectedindices
    selected[i] = false
    newGxminus1 = clen(compressor, join(Xstrs[selected]))
    selected[i] = true
    if newGxminus1 > maxGxminus1
      maxGxminus1 = newGxminus1
      largest = i
    end
  end
  return maxGxminus1, largest
end

# Calculate the sequence of NCD1 values used in the heuristic approximation
# of NCDM. We stop the sequence when there is not enough string length left
# for the compressor to compress anything since we cannot trust results at that
# point.
function ncd1_sequence(Xstrs::Vector{ASCIIString}, compressor::Compressor = ZlibC)
  Gxs = map(x -> clen(compressor, x), Xstrs)
  n = length(Xstrs)

  # Output will be a list of tuples each saving the NCD1 of that step and the index
  # of the element eliminated in that step.
  result = (Float64,Int)[]
  maxncd1 = -Inf # We will also calculate the maximum ncd1 which can be used as approximate NCDM
  maxatindex = -1

  # We will deselect one in each step but start with all selected
  selected = collect(1:n) .> 0
  selectedindices = Set(collect(1:n))
  numselected = n

  # Concatenation of all remaining strings
  remainingXS = join(Xstrs[selected])

  while numselected >= 2 && length(remainingXS) > minlen(compressor)

    # Calc max G(X \ {x}) for the still selected ones
    maxGxminus1, largest = max_GXminus1(Xstrs, selected, selectedindices, compressor)

    # NCD1 = (G(X) - minimum(G(X))) / maxGxminus1
    ncd1 = (clen(compressor, remainingXS) - minimum(Gxs[selected])) / maxGxminus1

    push!(result, (ncd1, largest))

    if ncd1 > maxncd1
      maxncd1 = ncd1
      maxatindex = length(result)
    end

    # Now remove the string that gave max G(X \ {x})
    selected[largest] = false
    delete!(selectedindices, largest)
    numselected -= 1
    remainingXS = join(Xstrs[selected])

  end

  return maxncd1, result, maxatindex
end

# We can apply it for any array of objects by just mapping them to strings.
function ncd1_sequence{T <: Any}(X::Array{T, 1}, compressor::Compressor = ZlibC)
  ncd1_sequence(map(string, X), compressor)
end

function ncdm{T <: Any}(X::Array{T, 1}, compressor::Compressor = ZlibC)
  ncd1_sequence(X, compressor)[1]
end

# We can use an NCD1 sequence to sort the objects if we assume they should enter
# in the reverse order to how they were deleted from the sequence. The elements
# that could not be deleted are randomly added at the start of the sequence.
function sortperm_ncd1_order{T <: Any}(X::Array{T, 1}, ncd1seq = false)
  if ncd1seq == false
    maxncd1, ncd1seq, maxatindex = ncd1_sequence(X)
  end
  deletion_order = map(t -> t[2], ncd1seq)
  non_deleted = setdiff(Set(1:length(X)), deletion_order)
  indices = vcat(shuffle(non_deleted), reverse(deletion_order))
end

using Plotly
Plotly.signin("robertfeldt", "903bj0pymv")

# We plot in reverse order and insert zeroes for values were we did not have enough
# information.
function ncd1plot{T <: Any}(X::Array{T, 1}; compressor::Compressor = ZlibC, plotname = "ncd1plot")
  ncdmvalue, seq = ncd1_sequence(X, compressor)
  num_zeros = length(X) - length(seq)
  values = vcat(zeros(Float64, num_zeros), map(t -> t[1], reverse(seq)))
  x = 1:length(values)

  # Now plot it with plotly
  ncd1_sequence_data = [
    "x" => x,
    "y" => values,
    "type" => "scatter"
  ]
  data = [ncd1_sequence_data]
  response = Plotly.plot(data, ["filename" => plotname, "fileopt" => "overwrite"])
  response["url"]
end

ext: jl


text:
import LinearAlgebra.BLAS: gemm, gemm!

for (f,T) in (
    (:(:cublasDgemm),:Float64),
    (:(:cublasSgemm),:Float32))
    @eval begin
        function gemm!(tA::Char, tB::Char,
            alpha::$T, A::CuVecOrMat{$T}, B::CuVecOrMat{$T},
            beta::$T, C::CuVecOrMat{$T})

            m = size(A, tA == 'N' ? 1 : 2)
            k = size(A, tA == 'N' ? 2 : 1)
            n = size(B, tB == 'N' ? 2 : 1)
            if m != size(C,1) || n != size(C,2) || k != size(B, tB == 'N' ? 1 : 2)
                throw(DimensionMismatch())
            end
            @cublas($f, (
                Ptr{Cvoid},Cint,Cint,
                Cint,Cint,Cint,
                Ptr{$T},Ptr{$T},Cint,
                Ptr{$T},Cint,
                Ptr{$T},Ptr{$T},Cint),
                gethandle(), cublasop(tA), cublasop(tB), m, n, k,
                Ref(alpha), A, stride(A,2), B, stride(B,2), Ref(beta), C, stride(C,2))
            C
        end
        function gemm(tA::Char, tB::Char, alpha::$T, A::CuVecOrMat{$T}, B::CuVecOrMat{$T})
            C = similar(B, size(A, tA=='N' ? 1 : 2), size(B, tB=='N' ? 2 : 1))
            gemm!(tA, tB, alpha, A, B, zero($T), C)
        end
        gemm(tA::Char, tB::Char, A::CuVecOrMat{$T}, B::CuVecOrMat{$T}) = gemm(tA, tB, one($T), A, B)
        gemm(A::CuVecOrMat{$T}, B::CuVecOrMat{$T}; tA='N', tB='N', alpha=1) = gemm(tA, tB, $T(alpha), A, B)
    end
end

for (fname,elty) in ((:cublasDgemmBatched,:Float64), (:cublasSgemmBatched,:Float32))
    @eval begin
        function gemm_batched!(tA::Char, tB::Char,
            alpha::$elty, As::Vector{CuMatrix{$elty}}, Bs::Vector{CuMatrix{$elty}},
            beta::$elty, Cs::Vector{CuMatrix{$elty}})

            if (length(As) != length(Bs) || length(As) != length(Cs))
                throw(DimensionMismatch(""))
            end
            for i = 1:length(As)
                A, B, C = As[i], Bs[i], Cs[i]
                m = size(A, tA == 'N' ? 1 : 2)
                k = size(A, tA == 'N' ? 2 : 1)
                n = size(B, tB == 'N' ? 2 : 1)
                if m != size(C,1) || n != size(C,2) || k != size(B, tB == 'N' ? 1 : 2)
                    throw(DimensionMismatch(""))
                end
            end
            m = size(As[1], tA == 'N' ? 1 : 2)
            k = size(As[1], tA == 'N' ? 2 : 1)
            n = size(Bs[1], tB == 'N' ? 2 : 1)
            lda = max(1, stride(As[1],2))
            ldb = max(1, stride(Bs[1],2))
            ldc = max(1, stride(Cs[1],2))
            Aptrs = map(pointer, As)
            Bptrs = map(pointer, Bs)
            Cptrs = map(pointer, Cs)
            $fname(handle(C), cublasop(tA), cublasop(tB), m, n, k, Ref(alpha), pointer(Aptrs),
                lda, pointer(Bptrs), ldb, Ref(beta), pointer(Cptrs), ldc, length(As))
            Cs
        end
    end
end
function gemm_batched(tA::Char, tB::Char,
    alpha::T, A::Vector{CuVecOrMat{T}}, B::Vector{CuVecOrMat{T}}) where T
    C = CudaMatrix{T}[similar(B[1], (size(A[1], tA=='N' ? 1 : 2), size(B[1], tB=='N' ? 2 : 1))) for i in 1:length(A)]
    gemm_batched!(tA, tB, alpha, A, B, T(0), C)
end
function gemm_batched(tA::Char, tB::Char,
    A::Vector{CuVecOrMat{T}}, B::Vector{CuVecOrMat{T}}) where T
    gemm_batched(tA, tB, T(1), A, B)
end


text:
# Tests the GeoUnits
using Test
using GeoParams
#using Parameters

@testset "Units" begin

# test GeoUnits
a = GeoUnit(100km);
@test a.val==100
@test Unit(a)==km
@test isdimensional(a)==true

b = GeoUnit(100);
@test b.val==100
@test Unit(b)==NoUnits
@test isdimensional(b)==false

a = GeoUnit(100km:10km:200km);
@test a.val==100:10:200
@test Unit(a)==km
@test isdimensional(a)==true

b = GeoUnit(100:10:200);
@test b.val==100:10:200
@test Unit(b)==NoUnits
@test isdimensional(b)==false

a = GeoUnit([100km 200km; 300km 400km]);
@test a.val==[100 200; 300 400]
@test Unit(a)==km
@test isdimensional(a)==true

b = GeoUnit([100.0 200; 300 400]);
@test b.val==[100.0 200.0; 300.0 400.0]
@test Unit(b)==NoUnits
@test isdimensional(b)==false

# test creating structures
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=1000km);
@test CharUnits_GEO.length  ==  1000km
@test CharUnits_GEO.Pa      ==  10000000Pa
@test CharUnits_GEO.Mass    ==  1.0e37kg
@test CharUnits_GEO.Time    ==  1.0e12s
@test CharUnits_GEO.Length  ==  1000000m

CharUnits_SI   =   SI_units();
@test CharUnits_SI.length   ==  1000m
@test CharUnits_SI.Pa       ==  10Pa

CharUnits_NO   =   NO_units();
@test CharUnits_NO.length   ==  1
@test CharUnits_NO.Pa       ==  1

# test nondimensionization of various parameters
@test nondimensionalize(10cm/yr,CharUnits_GEO) ≈ 0.0031688087814028945   rtol=1e-10
@test nondimensionalize(10cm/yr,CharUnits_SI) ≈ 3.168808781402895e7   rtol=1e-10

@test nondimensionalize([10 1; 20 30.0],CharUnits_GEO) == [10 1; 20 30.0]

A = 10MPa*s^(-1)   # should give 1e12 in ND units
@test nondimensionalize(A,CharUnits_GEO)  ≈ 1e12

B = GeoUnit(10MPa*s^(-1))
@test nondimensionalize(B,CharUnits_GEO)  ≈ 1e12

A = 10Pa^1.2*s^(-1)
@test nondimensionalize(A,CharUnits_GEO)  ≈ 39810.71705534975

CharUnits   =   GEO_units(viscosity=1e23, length=100km, stress=0.00371MPa);
@test nondimensionalize(A,CharUnits)  ≈ 1.40403327e16

A = (1.58*10^(-25))*Pa^(-4.2)*s^(-1)        # calcite
@test nondimensionalize(A,CharUnits_GEO)  ≈ 3.968780561785161e16

R=8.314u"J/mol/K"
@test nondimensionalize(R,CharUnits_SI)  ≈ 8.314e-7

# test Dimensionalize in case we provide a number and units
v_ND      =   nondimensionalize(3cm/yr, CharUnits_GEO);
@test dimensionalize(v_ND, cm/yr, CharUnits_GEO) == 3.0cm/yr

# Test the GeoUnit struct
x=GeoUnit(8.1cm/yr)
@test  x.val==8.1
x_ND  = nondimensionalize(x,CharUnits_GEO)
@test  x_ND ≈ 0.002566735112936345 rtol=1e-8        # nondimensionalize a single value
x_dim = dimensionalize(x,CharUnits_GEO)
@test  x_dim.val == 8.1                              # Dimensionalize again

y   =   x+2cm/yr
@test  y.val==10.1                            # errors

xx  =GeoUnit([8.1cm/yr; 10cm/yr]);
@test  xx.val==[8.1; 10]
yy=xx/(1cm/yr);                                         # transfer to no-unt

z=GeoUnit([100.0km 1000km 11km; 10km 2km 1km]);       # array
@test z/1km == [100.0 1000.0 11.0; 10.0 2.0 1.0]    # The division by 1km transfer it to a GeoUnit structure with no units; the multiplying with a float creates a float array

x = Vector(1:10)*km
zz = GeoUnit(x)
@test zz/1km == 1:10

# Test non-dimensionalisation if z is an array
z_ND = nondimensionalize(z,CharUnits_GEO)
@test z_ND.val==[0.1   1.0    0.011; 0.01  0.002  0.001]
@test isdimensional(z_ND)==false

z_D = dimensionalize(z_ND,CharUnits_GEO)
@test z_D.val==[100 1000 11; 10 2 1]          # transform back
@test isdimensional(z_D)==true

# test extracting a value from a GeoUnit array
@test NumValue(z_D[2,2]) == 2.0

# test setting a new value
z_D[2,1]=3
@test z_D[2,1].val == 3.0

# Conversion to different units
@test convert(GeoUnit,Float64(10.1)).val==10.1
@test convert(GeoUnit,Vector(10.1:.1:20)).val == 10.1:.1:20
@test Unit(convert(GeoUnit, 10km/s))==km/s
@test convert(Float64,  GeoUnit(10.2)) == 10.2
@test convert(Float64,  GeoUnit([10.2 11.2])) == [10.2 11.2]

# test various calculations (using arrays with and without units)
T_vec   =   273K:10K:500K        # using units
T       =   400.0K               # Unitful quanity
T_nd    =   10:10:200            # no units
α       =   GeoUnit(3e-5/K)
T₀      =   GeoUnit(293K)
ρ₀      =   GeoUnit(3300kg/m^3)

T_geo   =   GeoUnit(400K)
T_f = 400

# Calculations with two GeoUnits
t = T_geo+T₀
@test isa(t, GeoUnit)    # addition
@test Unit(t) == K
@test t ≈ 693.0

t = T_geo-T₀
@test isa(t, GeoUnit)    # subtraction
@test t ≈ 107
@test Unit(t) == K

t = T_geo*T₀
@test isa(t, GeoUnit)    # multiplication
@test t ≈ 117200.0
@test Unit(t) == K^2

t = T_geo/T₀
@test isa(t, GeoUnit)    # division
@test t ≈ 1.3651877133105803
@test Unit(t) == NoUnits

# Case in which one of them is a Unitful quantity
t = T+T₀
@test isa(t, Quantity)    # addition
@test t==693K

t = T-T₀
@test isa(t, Quantity)    # subtraction
@test t ≈ 107K

t = T*T₀
@test isa(t, Quantity)    # multiplication
@test t ≈ 117200.0K^2

t = T/T₀
@test isa(t, Float64)    # division
@test t ≈ 1.3651877133105803

# case in which one of them is a Float
t = T_f+T₀
@test isa(t, Float64)    # addition
@test t==693

t = T_f-T₀
@test isa(t, Float64)    # subtraction
@test t ≈ 107

t = T_f*T₀
@test isa(t, Float64)    # multiplication
@test t ≈ 117200.0

t = T_f/T₀
@test isa(t, Float64)    # division
@test t ≈ 1.3651877133105803

# case in which one of them is an array with Floats
t = T_nd+T₀
@test isa(t, AbstractArray)    # addition
@test t==303.0:10.0:493.0

t = T_nd-T₀
@test isa(t, AbstractArray)    # subtraction
@test t == -283.0:10.0:-93.0

t = T_nd*T₀
@test isa(t, AbstractArray)    # multiplication
@test t == 2930.0:2930.0:58600.0

t = T_nd/T₀
@test isa(t, AbstractArray)    # division
@test t ≈ 0.034129692832764506:0.034129692832764506:0.6825938566552902

# case in which one of them is an array with Quantities
t = T_vec+T₀
@test isa(t, AbstractArray)    # addition
@test t==(566.0:10.0:786.0)K

t = T_vec-T₀
@test isa(t, AbstractArray)    # subtraction
@test t == (-20.0:10.0:200.0)K

t = T_vec*T₀
@test isa(t, AbstractArray)    # multiplication
@test t == (79989.0:2930.0:144449.0)*K^2

t = T_vec/T₀
@test isa(t, AbstractArray)    # division
@test t ≈ 0.931740614334471:0.034129692832764506:1.68259385665529

ρ = ρ₀*(1.0 .- α*(T-T₀))  # The second expression is turned into a GeoUnit, so ρ should be GeoUnits
@test ρ ≈ 3289.4069999999997

ρ = ρ₀*(1.0 .- α*(T_nd-T₀))  # The second expression is turned into a GeoUnit, so ρ should be GeoUnits
@test ρ ≈ 3328.017:-0.99:3309.207

# test conversion:
b = convert(GeoUnit{Float64}, 3300kg/m^3)
@test b.val == 3300.0

c = convert(GeoUnit{Float32}, 3300.0kg/m^3)
@test c.val == 3300.0f0

d = convert(GeoUnit{Float32},100)
@test d.val == 100.0f0

d = convert(GeoUnit{Float64},100)
@test d.val == 100.0

d = convert(GeoUnit{Float64},[100 200])
@test d.val == [100.0 200.0]

# test broadcasting with float and unit arrays
ar   = [1.0 2;  4 5]            # Float array
ar1  = [1.0km 2km;  4km 5km]    # Unit array
test = GeoUnit(10km)

@test ar .+ test == [11.0 12.0; 14.0 15.0]
@test test .- ar == [9.0 8.0; 6.0 5.0]
@test test .* ar1 == [10.0km^2 20.0km^2; 40.0km^2 50.0km^2]
@test ar ./ ar1 == [1/km 1/km ; 1/km 1/km ]

# The way we define different structures heavily affects the number of allocs
# that are done while computing with parameters defined within the struct
struct ConDensity   # 1 allocation
    ρ::Float64
end

# Define struct with types and dimensions
struct ConDensity1{T} <: AbstractMaterialParam  # 1 allocation
    test::String
    ρ::GeoUnit{T}
    v::GeoUnit{T}
end

struct ConDensity2 <: AbstractMaterialParam # 3001 allocation
    test::String
    ρ::GeoUnit
    v::GeoUnit
end

# No type info
mutable struct ConDensity3 <: AbstractMaterialParam  # 3001 allocation
    test::String
    ρ::GeoUnit
end

# add type info in name, but no default values
mutable struct ConDensity4{T <: AbstractFloat} <: AbstractMaterialParam  # 1 allocation
    test::String
    ρ::GeoUnit{T}
    v::GeoUnit{T}
end

# use keywords, but no info about type
Base.@kwdef struct ConDensity5 <: AbstractMaterialParam # 3001 allocation
    test::String =""
    ρ::GeoUnit = GeoUnit(3300.1kg/m^3)
    v::GeoUnit = GeoUnit(100/s)
end

# This is slow
Base.@kwdef struct ConDensity6 <: AbstractMaterialParam # 3001 allocation
    test::String = ""
    ρ::GeoUnit = GeoUnit(3300.1kg/m^3)
    v::GeoUnit = GeoUnit(100/s)
end

# Here, we indicate the type:
Base.@kwdef struct ConDensity7{T} <: AbstractMaterialParam # 1 allocation
    test::String =""
    ρ::GeoUnit{T} = GeoUnit(3300.1kg/m^3)
    v::GeoUnit{T}   = GeoUnit(100/s)
end

# This seems to work, but has 3001 allocations
Base.@kwdef struct ConDensity8 <: AbstractMaterialParam # 3001 allocation
    test::String =""
    ρ::GeoUnit = 3300.1kg/m^3
    v::GeoUnit   = 100/s
end

# This works with 1 allocation, and is the preferred way to use it within GeoParams
Base.@kwdef struct ConDensity9{T <: AbstractFloat} <: AbstractMaterialParam # 1 allocation
    test::String =""
    ρ::GeoUnit{T} = 3300.1kg/m^3
    v::GeoUnit{T} = 100.0m/s
end
ConDensity9(a...) = ConDensity9{Float64}(a...)  # in case we do not give a type

rho  = ConDensity(3300.1)
rho1 = ConDensity1("t",GeoUnit(3300.1kg/m^3), GeoUnit(100/s))
rho2 = ConDensity2("t",3300.1kg/m^3, GeoUnit(100km/s))
rho3 = ConDensity3("t",3300.1kg/m^3)
rho4 = ConDensity4("t",GeoUnit(3300.1kg/m^3), GeoUnit(100/s))
rho5 = ConDensity5(test="t")
rho6 = ConDensity6("t",GeoUnit(3300.1kg/m^3), GeoUnit(100/s))
rho7 = ConDensity7()
rho8 = ConDensity8()
rho9 = ConDensity9(ρ=2800kg/m^3)
rho9_ND = nondimensionalize(rho9, GEO_units())

# test automatic nondimensionalization of a MaterialsParam struct:
CD = GEO_units()
rho2_ND = nondimensionalize(rho2, CD)
@test rho2_ND.ρ ≈ 3.3000999999999995e-18
@test rho2_ND.v ≈ 9.999999999999999e11

b = 20.1
c = 10.1;
r = 0.0

# Simple function to test speed
function f!(r,x,y)
    for i=1:1000
        r += x.ρ * y          # compute
    end
    r
end

function g!(r,x,y)
    for i=1:1000
        r += x.Density[1].ρ * y          # compute
    end
    r
end

#Phase1 = SetMaterialParams(Name="test1", Phase=22, Density  = ConDensity9())

#=
# testing speed (# of allocs)
r = 0.0
@btime f!($r, $rho,  $c)    # 1 allocations
@btime f!($r, $rho1, $c)    # 1 allocation
@btime f!($r, $rho2, $c)    # 3001 allocations
@btime f!($r, $rho3, $c)    # 3001 allocations
@btime f!($r, $rho4, $c)    # 1 allocation
@btime f!($r, $rho5, $c)    # 3001 allocation
@btime f!($r, $rho6, $c)    # 3001 allocation (so also with keywords, it is crucial to indicate the type)
@btime f!($r, $rho7, $c)    # 1 allocation (shows that we need to encode the units)
@btime f!($r, $rho8, $c)    # 3001 allocation
@btime f!($r, $rho9, $c)    # 1 allocation (shows that we need to encode the units)
@btime f!($r, $rho9_ND, $c)
=#

end


text:
module JavaCoreNLP

export StanfordCoreNLP, Annotation, demo, timinginformation, get_properties, get_property, sentences

include("core.jl")

end


text:
if abspath(PROGRAM_FILE) == @__FILE__
    a = 10
    b = "test"
    c = 2 + 4
    str1 = "hello $(a + 1) world"
    @assert(str1 == "hello 11 world")
    str2 = "hello $(b) world $(a)"
    @assert(str2 == "hello test world 10")
    str3 = "hello $(b) world $(a) test $(c)"
    @assert(str3 == "hello test world 10 test 6")
    a = 10
    @assert("hello $(a + 1) world" == "hello 11 world")
    @assert("test $(a + 1) test $(a + 20) test $(a + 30)" == "test 11 test 30 test 40")
    a = 1
    str4 = "hello $(a + 1) world $(round(0.44444, digits=3))"
    @assert(str4 == "hello 2 world 0.444")
    println("OK")
end


text:
# Standard stanh layer: f(A,x,b)=stanh(A*x+b)
# Note that the bias (a column vector) gets expanded here to match the dimension of A*x
sf=2.5

FstanhAXplusBias(A,X,b)=begin; a=A*X+b*ones(1,size(X,2)); return (sf*tanh(a),[]); end

function FstanhAXplusBias_inplace(value,aux,A,X,b)
    a=A*X+b*ones(1,size(X,2))
    copy!(value,sf*tanh(a))
end

function DstanhAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X,b)

    if derivativeIDX==1
        axpy!(1.0,(sf*grad_c.*(1.0-(f_c/sf).^2))*X',grad_n)
    elseif derivativeIDX==2
        axpy!(1.0,sf*A'*(grad_c.*(1.0-(f_c/sf).^2)),grad_n)
    elseif derivativeIDX==3
        axpy!(1.0,sf*sum(grad_c.*(1.0-(f_c/sf).^2),2),grad_n)
    end
end

if PROC=="GPU"
    function FstanhAXplusBias(A::CudaArray,X::CudaArray,b::CudaArray)
        tmp=FAXplusBias(A,X,b)[1]
        stanh!(sf,tmp,tmp)
        return(tmp,[])
    end

    function FstanhAXplusBias_inplace(value,aux,A::CudaArray,X::CudaArray,b::CudaArray)
        stanh!(sf,FAXplusBias(A,X,b)[1],value)
    end

    function DstanhAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray,b::CudaArray)
        tmp=CudaArray(Float64,size(grad_c)); fill!(tmp,0.0)
        Dstanh!(sf,grad_c,f_c,tmp)
        if derivativeIDX==1
            gemm!('N','T',1.0,tmp,X,1.0,grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',1.0,A,tmp,1.0,grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float64,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',1.0,tmp,ons,1.0,grad_n)
            free(ons)
        end
        free(tmp)
    end
end

Derivative[FstanhAXplusBias]=DstanhAXplusBias
Inplace[FstanhAXplusBias]=FstanhAXplusBias_inplace

stanhAXplusBias(A,X,b)=ADnode(FstanhAXplusBias,[A X b])
export stanhAXplusBias

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using LinearAlgebra
using Distributions
using Pathfinder
using Test

@testset "single path pathfinder" begin
    @testset "IsoNormal" begin
        # here pathfinder finds the exact solution after 1 iteration
        logp(x) = -sum(abs2, x) / 2
        ∇logp(x) = -x
        ndraws = 100
        @testset for n in [1, 5, 10, 100]
            x0 = randn(n)
            q, ϕ, logqϕ = @inferred pathfinder(logp, ∇logp, x0, ndraws)
            @test q isa MvNormal
            @test q.μ ≈ zeros(n)
            @test q.Σ isa Pathfinder.WoodburyPDMat
            @test q.Σ ≈ I
            @test size(q.Σ.B) == (n, 2) # history contains only 1 iteration
            @test ϕ isa AbstractMatrix
            @test size(ϕ) == (n, ndraws)
            @test logqϕ ≈ logpdf(q, ϕ)

            q2, ϕ2, logqϕ2 = pathfinder(logp, ∇logp, x0, 2)
            @test size(ϕ2) == (n, 2)
        end
    end
    @testset "MvNormal" begin
        #! format: off
        Σ = [
            2.71   0.5    0.19   0.07   1.04
            0.5    1.11  -0.08  -0.17  -0.08
            0.19  -0.08   0.26   0.07  -0.7
            0.07  -0.17   0.07   0.11  -0.21
            1.04  -0.08  -0.7   -0.21   8.65
        ]
        #! format: on
        P = inv(Symmetric(Σ))
        logp(x) = -dot(x, P, x) / 2
        ∇logp(x) = -(P * x)
        x₀ = [2.08, 3.77, -1.26, -0.97, -3.91]
        rng = MersenneTwister(38)
        q, _, _ = pathfinder(logp, ∇logp, x₀, 10; rng=rng, ndraws_elbo=100)
        @test q.Σ ≈ Σ rtol = 1e-1
    end
end


text:
using StatsModels
using RegressionFormulae
using RegressionFormulae: fulldummy
using Test

include("dummymod.jl")

dat = (; y=zeros(10), a=["u","i","o"], b=["q","w","e"], c=["s","d","f"], x=1:10)

@testset "error checking" begin
    @test_throws ArgumentError fulldummy(term(:a))

    sch = schema(dat)
    @test_throws ArgumentError apply_schema(@formula(y ~ fulldummy(x)), sch, RegressionModel)
end


text:
#module C
#export cfunc

#include("a.jl")
#using ..A

function cfunc(::Astruct)
end

#end


text:
############################################################################################

"""
$(SIGNATURES)
Numerically stable version for log(mean(exp(vec))).

# Examples
```julia
```

"""
function logmeanexp(arr::Vector{T}) where {T<:Real}
    return log(1 / size(arr, 1)) + LogExpFunctions.logsumexp(arr)
end

############################################################################################
"""
$(SIGNATURES)
Check if `x` is larger than `threshold`, guaranteeing same type for arguments.

# Examples
```julia
```

"""
function issmaller(x::T, threshold::T) where {T<:Real}
    #!TODO: function name ambigous if threshold name not explicitly mentioned
    #!NOTE: "<=" so Resampling at each iteration if threshold == 1.0
    return x <= threshold
end

############################################################################################
"""
$(SIGNATURES)
Return allocation friendly index of array argument.

# Examples
```julia
```

"""
function grab(data, index, sorted)
    return nothing
end

@inline function grab(
    data::AbstractVector{R}, index::I, sorted::S
) where {R,I<:Integer,S<:Union{ByRows,ByCols}}
    @inbounds return data[index]
end
@inline function grab(
    data::AbstractVector{R}, range::Union{Vector{I},UnitRange{I}}, sorted::S
) where {R,I<:Integer,S<:Union{ByRows,ByCols}}
    @inbounds return view(data, range)
end

@inline function grabcols(
    data::AbstractArray{R,2}, index::Union{I,Vector{I},UnitRange{I}}
) where {R,I<:Integer}
    return view(data, index, :)
end
@inline function grabrows(
    data::AbstractArray{R,2}, index::Union{I,Vector{I},UnitRange{I}}
) where {R,I<:Integer}
    return view(data, :, index)
end
function grab(
    data::AbstractArray{R,2}, idx::Union{I,Vector{I},UnitRange{I}}, sorted::ByRows
) where {R,I<:Integer}
    return grabcols(data, idx)
end
function grab(
    data::AbstractArray{R,2}, idx::Union{I,Vector{I},UnitRange{I}}, sorted::ByCols
) where {R,I<:Integer}
    return grabrows(data, idx)
end

grab(data, idx, config::ArrayConfig) = grab(data, idx, config.sorted)

############################################################################################
# Export
export logsumexp, #import from LogExpFunctions
    logaddexp, #import from LogExpFunctions
    logmeanexp, #import from LogExpFunctions
    issmaller,
    grab


text:
export loss
export cost

# Internal function
function _sum(x::Variable{T}) where T
    y = Variable{T}([sum(ᵛ(x))], x.backprop)
    if x.backprop
        function _sumBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(y)
            end
            ifNotKeepδThenFreeδ!(y)
        end
        push!(graph.backward, _sumBackward)
    end
    return y
end

# Internal function
function _mean(x::Variable{T}) where T
    n = eltype(x)(1) / prod(size(x))
    μ = Variable{T}([sum(ᵛ(x)) * n], x.backprop)
    if x.backprop
        function _meanBackward()
            if need2computeδ!(x)
                δ(x) .+= δ(μ) .* n
            end
            ifNotKeepδThenFreeδ!(μ);
        end
        push!(graph.backward, _meanBackward)
    end
    return μ
end

"""
    loss(x::Variable{T}; reduction::String="sum") -> y::Variable{T}

Sums or takes mean over all elements in `value` of `x` as the loss `Variable`, i.e. ⤦\n
+ `y = Variable{T}([sum(ᵛ(x))], x.backprop)`, if reduction="sum"
+ `y = Variable{T}([sum(ᵛ(x))/length(x)], x.backprop)`, if reduction="mean"
This is very convenient for mutiple loss training. e.g. ⤦\n
    totalLoss = β*loss(mse(x₁, ̂x₁)) + (1 - β)*loss(crossEntropy(y₁, ̂y₁))
or in a lazy way:\n
    totalLoss = β*mseLoss(x₁, ̂x₁) + (1 - β)*crossEntropyLoss(y₁, ̂y₁)
where `β` is weight of mseLoss function.
"""
function loss(x::Variable{T}; reduction::String="sum") where T
    by = lowercase(reduction)
    by=="sum" && return _sum(x)
    by=="mean" && return _mean(x)
    @error "wrong reduction method" reduction=="sum" || reduction=="mean"
end

"""
    cost(x::Variable) -> y::eltype(x)

sums over all elements in `value` of `x` as the final loss value which is a scalar, i.e. ⤦\n
    y = sum(x.value), so δ(x) = ∂y/∂x = 𝟙
"""
function cost(x::Variable{T}) where T
    if x.backprop
        𝟙 = eltype(x)(1.0)
        function costBackward()
            if need2computeδ!(x)
                δ(x) .+= 𝟙
            end
        end
        push!(graph.backward, costBackward)
    end
    return sum(ᵛ(x))
end

include("./1-criterion-regression.jl")
include("./1-criterion-probabilistic.jl")

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using Parquet
using Test

function test_load(file, parcompat=joinpath(dirname(@__FILE__), "parquet-compatibility"))
    p = ParFile(joinpath(parcompat, file))
    println("loaded $file")
    @test isa(p.meta, Parquet.FileMetaData)

    rg = rowgroups(p)
    @test length(rg) > 0
    println("\tfound $(length(rg)) row groups")

    cc = columns(p, 1)
    println("\tfound $(length(cc)) column chunks in row group 1")

    cnames = colnames(rg[1])
    @test length(cnames) == length(cc)
    println("\tcolumns: $cnames")

    pg = pages(p, 1, 1)
    println("\tfound $(length(pg)) pages in column chunk 1 of row group 1")

    bpg = bytes(pg[1])
    println("\tread $(length(bpg)) bytes from page 1, column chunk 1 of row group 1")

    cc1 = columns(p, rg[1], cnames[1:2:3])
    @test length(cc) > length(cc1)
    @test length(rowgroups(p, [cnames[1], "xyz"])) == 0

    iob = IOBuffer()
    show(iob, p)
    sb = take!(iob)
    @test !isempty(sb)
    println("\t" * String(sb))

    println("\tsuccess")
end

function print_names(indent, t)
    for f in fieldnames(t)
        ftype = fieldtype(t, f)
        println("\t", indent, f, "::", ftype)
        if isa(ftype, DataType)
            print_names(indent * "    ", ftype)
        end
    end
end

function test_schema(file, schema_name::Symbol,  parcompat=joinpath(dirname(@__FILE__), "parquet-compatibility"))
    p = ParFile(joinpath(parcompat, file))
    println("loaded $file")

    #mod_name = string(schema_name) * "Module"
    #eval(Main, parse("module $mod_name end"))
    #mod = getfield(Main, Symbol(mod_name))

    mod = Main
    schema(JuliaConverter(mod), schema(p), schema_name)
    @test schema_name in names(mod, all=true)
    println("\tschema: \n\t", getfield(mod, schema_name))
    print_names("    ", getfield(mod, schema_name))

    io = IOBuffer()
    schema(ThriftConverter(io), schema(p), schema_name)
    thriftsch = strip(String(take!(io)))
    @test startswith(thriftsch, "struct " * string(schema_name))
    @test endswith(thriftsch, "}")
    println(thriftsch)

    io = IOBuffer()
    schema(ProtoConverter(io), schema(p), schema_name)
    protosch = strip(String(take!(io)))
    @test startswith(protosch, "message " * string(schema_name))
    @test endswith(protosch, "}")
    println(protosch)
end

function test_load_all_pages()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            test_load("parquet-testdata/impala/1.1.1-$encformat/$fname.impala.parquet")
            test_schema("parquet-testdata/impala/1.1.1-$encformat/$fname.impala.parquet", Symbol(fname * "_" * encformat))
        end
    end
end

test_load_all_pages()


text:
using BSON: @save, @load
using Nettle
using Sockets

int(s::String) = parse(Int, s)
bigint(s::String) = parse(BigInt, s)

"""
  `var_to_bytes(variable)::Bytes`

Makes Byte representation (Vector{UInt8}) of given variable
This allows sending a variable over a socket connection in byte form
"""
function var_to_bytes(variable)::Bytes
  io = IOBuffer(write=true, read=true, append=true)
  write(io, variable)
  return read(io)
end

"""
  `model_to_bytes(model)::Bytes`

Converts a Flux model to Bytes (Vector{UInt8})
This allows sending a Flux model over a socket connection in byte form
"""
function model_to_bytes(model)::Bytes
  io = IOBuffer(write=true, read=true, append=true)
  @save io model
  return read(io)
end

"""
  `var_from_bytes(bytes::Bytes, to_object)` -> variable

Writes byte data into already initialized object (to_object).
This allows sending a variable over a socket connection in byte form
"""
function var_from_bytes(bytes::Bytes, to_object)
  # write binary into buffer
  io = IOBuffer(read=true, write=true, append=true)
  write(io, bytes)

  # copy buffer data to object
  read!(io, to_object)
  return to_object
end

"""
  `model_from_bytes(bytes::Bytes)` -> Flux model

Loads a BSON model from Bytes
This allows receiving a Flux model over a socket connection in byte form
"""
function model_from_bytes(bytes::Bytes)
  #  write binary into buffer
  io = IOBuffer(read=true, write=true, append=true)
  write(io, bytes)

  # copy buffer data to object
  @load io model
  return model
end

"""
  `model_from_buffer(io::IOBuffer)` -> Flux model

Loads a BSON model from an IOBuffer (byte form)
This allows receiving a Flux model over a socket connection in byte form
"""
function model_from_buffer(io::IOBuffer)
  # copy buffer data to object
  seekstart(io)
  @load io model
  return model
end

"""
  `receive_multi_str(con::TCPSocket, start::String, stop::String)::String`

Receives multiple lines (as a String) from a TCP connection
Returns the string of everything between `start` and `stop`

`start`: prefix indicating start of desired message
`stop`: suffix indicating end of desired message
Both `start` and `stop` are removed from final resulting String
"""
function receive_multi_str(con::TCPSocket, start::String, stop::String)::String
  buf = IOBuffer(read=true, write=true, append=true)
  line = readline(con, keep=true)
  if !occursin(start, line)
    close(con)
    throw(ConnectionError("Server did not send expected start phrase: $start"))
  end
  line = split(line, start)[end]

  if occursin(stop, line)
    line = split(line, stop)[1]
    return line
  end

  while !occursin(stop, line)
    write(buf, line)
    line = readline(con, keep=true)
  end

  line = split(line, stop)[1]
  write(buf, line)
  content = read(buf)
  return to_str(content)
end

"""
  `send_int(con::TCPSocket, num::Integer)`

Sends integer to socket connection (First packs integer to bytes to reduce size)
"""
function send_int(con::TCPSocket, num::Integer)::Nothing
  num_bytes = pack_int(num)
  #  s = to_str(num_bytes) # <-- this becomes corrupted for large integers
  s = bytes2hex(num_bytes)
  prefix = "STARTINT"
  suffix = "STOPINT"
  res = write(con, prefix, s, suffix, "\n")
  if res != length(prefix) + length(s) + length(suffix) + 1
    close(con)
    throw(ConnectionError("Failed to send integer ($res, $(length(s) + 16))"))
  end
end

"""
  `receive_int(con::TCPSocket)::BigInt`

Receive BigInt from socket connection (Assumes number was packed into bytes)
"""
function receive_int(con::TCPSocket)::BigInt
  try
    prefix = "STARTINT"
    suffix = "STOPINT"
    #  msg = readline(con, keep=true)
    msg = receive_multi_str(con, prefix, suffix)
    #  msg_bytes = to_bytes(strip(msg)) <-- this becomes corrupted for large integers
    msg_bytes = hex2bytes(strip(msg))
    num = unpack_int(msg_bytes)
    return num
  catch err
    close(con)
    throw(ConnectionError("Failed to receive valid integer"))
  end
end

"""
  `send_bytes(con::TCPSocket, msg::Bytes)::Nothing`

Sends Bytes over TCP socket connection
"""
function send_bytes(con::TCPSocket, msg::Bytes)::Nothing
  s = bytes2hex(msg)
  res = write(con, s, "\n")
  if res != length(s) + 1
    close(con)
    throw(ConnectionError("Failed to send bytes"))
  end
end

"""
  `receive_bytes(con::TCPSocket)::Bytes`

Receive Bytes (Vector{UInt8}) from TCP socket connection
"""
function receive_bytes(con::TCPSocket)::Bytes
  msg = readline(con, keep=false)
  try
    msg_bytes = hex2bytes(msg)
    return msg_bytes
  catch err
    close(con)
    throw(ConnectionError("Failed to receive bytes as expected"))
  end
end

"""
  `receive_str(con::TCPSocket)::String`

Receive String from socket connection
"""
function receive_str(con::TCPSocket)::String
  msg_bytes = receive_bytes(con)
  return to_str(msg_bytes)
end

"""
  `handshake_client(con::TCPSocket, [clientID::String])::Bytes` -> symmetric key

Establishes shared symmetric key (in Bytes) with server. Performs key exchange securely with RSA.
"""
function handshake_client(con::TCPSocket, clientID = nothing)::Bytes
  # Request passphase that server uses to initiate conversation
  expected_phrase = "INIT"
  phrase = readline(con, keep=true)
  if !occursin(expected_phrase, phrase)
    close(con)
    throw(ConnectionError("Expected handshake to start with phrase $expected_phrase"))
  end

  # Provide client ID to let server know client is authorized (optional)
  if typeof(clientID) != Nothing
    try
      write(con, "ID ", string(clientID), "\n")
    catch err
      throw(ArgumentError("clientID should be a String"))
    end
  else
    write(con, "NO_ID\n")
  end

  # Request server public key for RSA
  #  try
  #  catch err
    #  throw(ConnectionError("Server did not provide valid public key"))
  #  end
  pk_e = receive_int(con)
  pk_n = receive_int(con)
  public_key = PublicRSAKey(pk_e, pk_n)

  # Generate 32-bit symmetric key (in Bytes)
  symmetric_key = gen_key32()
  enc = Encryptor("AES256", symmetric_key)
  dec = Decryptor("AES256", symmetric_key)

  # Encrypt symmetric key with server's public key using RSA
  symmetric_key_int = unpack_int(symmetric_key)
  hidden_key = powermod(symmetric_key_int, public_key.e, public_key.n)

  # Send symmetric key to server (for AES)
  send_int(con, hidden_key)

  # Receive random message from server, encrypted with symmetric key
  ciphertext_bytes = receive_bytes(con)

  # Decrypt random message (using AES)
  msg_bytes = decrypt(dec, ciphertext_bytes)

  # Encrypt message with server's public key (using RSA)
  msg_int = unpack_int(msg_bytes)
  enc_msg = powermod(msg_int, public_key.e, public_key.n)

  # Send back encrypted message
  send_int(con, enc_msg)

  return symmetric_key
end

"""
  `handshake_server(con::TCPSocket, [verifyID])::Bytes` -> symmetric key

Establishes shared symmetric key (in Bytes) with client. Performs key exchange securely with RSA.

`verifyID`: a callback function, used to check if the client ID is authorized (optional)
"""
function handshake_server(con::TCPSocket, verifyID = nothing)
  # Send expected phrase to start handshake
  expected_phrase = "INIT"
  write(con, expected_phrase, "\n")

  # Obtain Client ID if available
  line = readline(con, keep=false)
  if !occursin("ID", line)
    close(con)
    throw(ConnectionError("Client did not provide client ID, or indicate that no ID is available"))
  end
  if !occursin("NO_ID", line)
    clientID = split("ID ", line)[2]
  else
    clientID = nothing
  end

  # Verify that client ID is authorized using callback function (optional)
  if typeof(verifyID) != Nothing
    verifyID(clientID)
  end

  # Generate RSA public key and private key
  public_key, private_key = gen_RSAKeys()

  # Send client public key for RSA
  send_int(con, public_key.e)
  send_int(con, public_key.n)

  # Receive encrypted symmetric key for AES
  hidden_key = receive_int(con)

  # Decrypt key using RSA, and decode to Bytes
  symmetric_key_int = powermod(hidden_key, private_key.d, private_key.n)
  symmetric_key = pack_int(symmetric_key_int) # Bytes

  # Generate random message (64 bits)
  random_message = secure_randstring(64)
  random_message_bytes = Bytes(random_message)

  # Encrypt message with AES and send to client
  enc = Encryptor("AES256", symmetric_key)
  encrypted_message = encrypt(enc, random_message_bytes)
  send_bytes(con, encrypted_message)

  # Receive message encrypted with RSA public key
  hidden_msg = receive_int(con)

  # Decrypt message using RSA, and decode to String
  msg_int = powermod(hidden_msg, private_key.d, private_key.n)
  msg_bytes = pack_int(msg_int)
  msg = to_str(msg_bytes)

  # Check if handshake was successful
  if msg != random_message
    close(con)
    throw(ConnectionError("Failed to establish shared symmetric key with client"))
  end

  return symmetric_key
end

"""
  `secure_send(con::TCPSocket, msg::String, AES_key::Bytes)::Nothing`

Encrypts message (String) with AES256 and sends it over socket connection
"""
function secure_send(con::TCPSocket, msg::String, AES_key::Bytes)::Nothing
  enc = Encryptor("AES256", AES_key)
  pad_length = (16 - (length(msg) % 16)) % 16
  padded_msg = ("\n"^pad_length) * msg
  padded_msgmsg = to_bytes(msg)
  enc_msg = encrypt(enc, padded_msg)
  send_bytes(con, enc_msg)
end

"""
  `secure_receive(con::TCPSocket, AES_key::Bytes)::String`

Receives and decrypts message (String) from TCP Socket connection using AES256
(Note that white space padding will get stripped from the ends of the message)
"""
function secure_receive(con::TCPSocket, AES_key::Bytes)::String
  dec = Decryptor("AES256", AES_key)
  enc_msg = receive_bytes(con)
  msg_bytes =  decrypt(dec, enc_msg)
  return strip(to_str(msg_bytes))
end

"""
  `secure_send_bytes(con::TCPSocket, msg::Bytes, AES_key::Bytes)::Nothing`

Encrypts message (Bytes / Vector{UInt8}) with AES256 and sends it over socket connection
"""
function secure_send_bytes(con::TCPSocket, msg::Bytes, AES_key::Bytes)::Nothing
  msg_hex = bytes2hex(msg)
  secure_send(con, msg_hex, AES_key)
end

"""
  `secure_receive_bytes(con::TCPSocket, AES_key::Bytes)::Bytes`

Receives and decrypts message (Bytes / Vector{UInt8}) from TCP Socket connection using AES256
(Note that white space padding will get stripped from the ends of the message)
"""
function secure_receive_bytes(con::TCPSocket, AES_key::Bytes)::Bytes
  msg = secure_receive(con, AES_key)
  return hex2bytes(msg)
end

"""
  `secure_send_int(con::TCPSocket, msg::Integer, AES_key::Bytes)::Nothing`

Encrypts message (BigInt) with AES256 and sends it over socket connection
"""
function secure_send_int(con::TCPSocket, msg::Integer, AES_key::Bytes)::Nothing
  msg_bytes = pack_int(msg)
  secure_send_bytes(con, msg_bytes, AES_key)
end

"""
  `secure_receive_int(con::TCPSocket, AES_key::Bytes)::BigInt`

Receives and decrypts message (BigInt) from TCP Socket connection using AES256
(Note that white space padding will get stripped from the ends of the message)
"""
function secure_receive_int(con::TCPSocket, AES_key::Bytes)::BigInt
  msg_bytes = secure_receive_bytes(con, AES_key)
  return unpack_int(msg_bytes)
end


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

fxb = Symbol("f×ᵦg")

eqschema = @acset LabeledGraph begin
    V = 3
    E = 4
    vlabel = [:A,:B,fxb]
    elabel = [:f,:g, :π1, :π2]
    src = [1,1,3,3]
    tgt = [2,2,1,1]
end

eqconed = @acset LabeledGraph begin
    V = 3
    E = 2
    vlabel = [:A,:A,:B]
    elabel = [:f,:g]
    src = [1,2]
    tgt = [3,3]
end

eqfls = FLS(:equalizer, eqschema, [
    Cone(eqconed, fxb, [1=>:π1,2=>:π2])], []);


text:
## SMC initial proposal ##

@gen function init_proposal(obs)
    n_strokes ~ shifted_poisson(0.25, 1)
    init_x ~ trunc_normal(obs[1], 0.05, 0.0, 5.0)
    init_y ~ trunc_normal(obs[2], 0.05, 0.0, 5.0)
end

## SMC update kernels ##

@gen function fwd_stroke_proposal(trace, t::Int, obs)
    # Resample a larger number of strokes
    n_strokes ~ shifted_poisson(0.25, trace[:n_strokes]+1)
    # Disconnect stroke if new observation is far from previous one
    prev_obs = trace[:obs][t-1]
    p_disconnect = 1 - exp(-sum(((obs .- prev_obs)./1.5).^2))
    disconnect ~ bernoulli(p_disconnect)
    # Sample new stroke location based on observation
    if disconnect
        start_x ~ trunc_normal(obs[1], 0.05, 0.0, 5.0)
        start_y ~ trunc_normal(obs[2], 0.05, 0.0, 5.0)
    end
end

@gen function bwd_stroke_proposal(trace, t::Int, obs)
    n_strokes ~ poisson(t/10 + 1)
end

@transform stroke_transform (p_old, q_fwd) to (p_new, q_bwd) begin
    n_new = @read(q_fwd[:n_strokes], :discrete)
    n_old = @read(p_old[:n_strokes], :discrete)
    @write(p_new[:n_strokes], n_new, :discrete)
    @write(q_bwd[:n_strokes], n_old, :discrete)
    @copy(q_fwd[:disconnect], p_new[:strokes => n_old + 1 => :disconnect])
    if @read(q_fwd[:disconnect], :discrete)
        @copy(q_fwd[:start_x], p_new[:strokes => n_old + 1 => :start_x])
        @copy(q_fwd[:start_y], p_new[:strokes => n_old + 1 => :start_y])
    end
end

@gen function fwd_points_proposal(trace, t::Int, obs)
    n_points = sum(length.(trace[:points_per_stroke]))
    n_stroke_points_prev = trace[:strokes => trace[:n_strokes] => :n_points]
    n_stroke_points ~ shifted_poisson(2, (t - n_points) + n_stroke_points_prev)
end

@gen function bwd_points_proposal(trace, t::Int, obs)
    n_stroke_points ~ poisson(10)
end

@transform points_transform (p_old, q_fwd) to (p_new, q_bwd) begin
    n_strokes = @read(p_old[:n_strokes], :discrete)
    @copy(q_fwd[:n_stroke_points], p_new[:strokes => n_strokes => :n_points])
    @copy(p_old[:strokes => n_strokes => :n_points], q_bwd[:n_stroke_points])
end

@gen function split_proposal(trace, t::Int, obs, i_stroke::Int, offset::Int)
    # Sample new stroke parameters
    prev_stroke = trace[:strokes => i_stroke => :stroke]
    speed ~ gamma(2.0, 0.05)
    # Disconnect stroke if new observation is far from previous one
    prev_obs = trace[:obs][t-1]
    p_disconnect = 1 - exp(-sum(((obs .- prev_obs)./1.5).^2))
    disconnect ~ bernoulli(p_disconnect)
    # Sample new stroke location based on observation
    if disconnect
        stroke ~ categorical(type_prior)
        start_x ~ trunc_normal(obs[1], 0.05, 0.0, 5.0)
        start_y ~ trunc_normal(obs[2], 0.05, 0.0, 5.0)
    else
        est_dir = argmax(directions * (obs - prev_obs))
        dir_probs = ones(4) * 0.05
        dir_probs[est_dir] = 0.85
        stroke ~ categorical(dir_probs)
    end
    return (i_stroke, offset)
end

@transform split_transform (p_old, q_fwd) to (p_new, q_bwd) begin
    i_stroke, offset = @read(q_fwd[], :discrete)
    # Increment number of strokes
    n_strokes = @read(p_old[:n_strokes], :discrete)
    @write(p_new[:n_strokes], n_strokes + 1, :discrete)
    # Split points between old and new stroke
    n_points = @read(p_old[:strokes => i_stroke => :n_points], :discrete)
    points_1, points_2 = offset - 1, n_points - offset +1
    @write(p_new[:strokes => i_stroke => :n_points], points_1, :discrete)
    @write(p_new[:strokes => i_stroke + 1 => :n_points], points_2, :discrete)
    # Copy over sampled parameters for new stroke
    @copy(q_fwd[:stroke], p_new[:strokes => i_stroke + 1 => :stroke])
    @copy(q_fwd[:speed], p_new[:strokes => i_stroke + 1 => :speed])
    @copy(q_fwd[:disconnect], p_new[:strokes => i_stroke + 1 => :disconnect])
    if @read(q_fwd[:disconnect], :discrete)
        @copy(q_fwd[:start_x], p_new[:strokes => i_stroke + 1 => :start_x])
        @copy(q_fwd[:start_y], p_new[:strokes => i_stroke + 1 => :start_y])
    end
    # Shift indices forward for all subsequent strokes
    for j_stroke in i_stroke+1:n_strokes
        @copy(p_old[:strokes => j_stroke], p_new[:strokes => j_stroke + 1])
    end
end

@gen function merge_proposal(trace, t::Int, obs, i_stroke::Int)
    # Do nothing (no additional randomness needed)
    return i_stroke
end

@transform merge_transform (p_old, q_fwd) to (p_new, q_bwd) begin
    i_stroke = @read(q_fwd[], :discrete)
    # Decrement number of strokes
    n_strokes = @read(p_old[:n_strokes], :discrete)
    @write(p_new[:n_strokes], n_strokes - 1, :discrete)
    # Merge points from current stroke to previous stroke
    points_1 = @read(p_old[:strokes => i_stroke - 1 => :n_points], :discrete)
    points_2 = @read(p_old[:strokes => i_stroke => :n_points], :discrete)
    @write(p_new[:strokes => i_stroke - 1 => :n_points], points_1 + points_2, :discrete)
    # Copy deleted parameters to backward proposal
    @copy(p_old[:strokes => i_stroke => :stroke], q_bwd[:stroke])
    @copy(p_old[:strokes => i_stroke => :speed], q_bwd[:speed])
    @copy(p_old[:strokes => i_stroke => :disconnect], q_bwd[:disconnect])
    if @read(p_old[:strokes => i_stroke => :disconnect], :discrete)
        @copy(p_old[:strokes => i_stroke => :start_x], q_bwd[:start_x])
        @copy(p_old[:strokes => i_stroke => :start_y], q_bwd[:start_y])
    end
    # Shift indices backward for all subsequent strokes
    for j_stroke in i_stroke+1:n_strokes
        @copy(p_old[:strokes => j_stroke], p_new[:strokes => j_stroke - 1])
    end
end

@gen function fwd_proposal(trace, t::Int, obs)
    # Find total points, current stroke, and offset from stroke start
    stroke_lengths = length.(trace[:points_per_stroke])
    n_points = sum(stroke_lengths)
    i_stroke, offset = t > n_points ? (nothing, nothing) : segment_find(t, stroke_lengths)
    # Guess if new observation is part of the same stroke as previous one
    prev_obs = t == 1 ? obs : trace[:obs][t-1]
    p_connected = exp(-sum(((obs .- prev_obs)./1.5).^2))
    # Decide between sub-proposals
    if t > n_points # Extend with new stroke or points
        p_new_stroke, p_new_points = 1-p_connected, p_connected
        p_split, p_merge, p_default = 0.0, 0.0, 0.0
    elseif offset > 1 # Potentially split current stroke
        p_new_stroke, p_new_points, p_merge = 0.0, 0.0, 0.0
        p_split, p_default = 1-p_connected, p_connected
    elseif i_stroke > 1 # Potentially merge current and previous stroke
        p_new_stroke, p_new_points, p_split = 0.0, 0.0, 0.0
        p_merge, p_default = p_connected/2, 1-p_connected/2
    else # Use default proposal
        p_new_stroke, p_new_points = 0.0, 0.0
        p_split, p_merge, p_default = 0.0, 0.0, 1.0
    end
    # Sample a branch and call sub-proposal
    branch ~ categorical([p_new_stroke, p_new_points, p_split, p_merge, p_default])
    if branch == 1 # Extend with new stroke
        {*} ~ fwd_stroke_proposal(trace, t, obs)
    elseif branch == 2 # Extend with new points
        {*} ~ fwd_points_proposal(trace, t, obs)
    elseif branch == 3 # Split current stroke
        {*} ~ split_proposal(trace, t, obs, i_stroke, offset)
    elseif branch == 4 # Merge current and previous stroke
        {*} ~ merge_proposal(trace, t, obs, i_stroke)
    elseif branch == 5 # Use default proposal
        nothing
    end
end

@gen function bwd_proposal(trace, t::Int, obs)
    # Find current stroke, and offset from stroke start
    stroke_lengths = length.(trace[:points_per_stroke])
    n_points = sum(stroke_lengths)
    i_stroke, offset = t > n_points ? (nothing, nothing) : segment_find(t, stroke_lengths)
    # Sample a branch and call sub-proposal
    branch ~ categorical([0.1, 0.1, 0.05, 0.05, 0.7])
    if branch == 1 # Reverse stroke extension
        {*} ~ bwd_stroke_proposal(trace, t, obs)
    elseif branch == 2 # Reverse points extension
        {*} ~ bwd_points_proposal(trace, t, obs)
    elseif branch == 3 # Reverse split proposal
        {*} ~ merge_proposal(trace, t, obs, i_stroke)
    elseif branch == 4 # Reverse merge proposal
        {*} ~ split_proposal(trace, t, obs, i_stroke, offset)
    elseif branch == 5 # Do nothing for default proposal
        nothing
    end
end

@transform update_transform (p_old, q_fwd) to (p_new, q_bwd) begin
    branch = @read(q_fwd[:branch], :discrete)
    @write(q_bwd[:branch], branch, :discrete)
    if branch == 1
        @tcall(stroke_transform())
    elseif branch == 2
        @tcall(points_transform())
    elseif branch == 3
        @tcall(split_transform())
    elseif branch == 4
        @tcall(merge_transform())
    end
end

## SMC rejuvenation kernels ##

@gen function rejuv_speed(trace, t_prev::Int=0)
    i_stroke = trace[:n_strokes] - t_prev
    if i_stroke <= 0 return end
    prev_speed = trace[:strokes => i_stroke => :speed]
    {:strokes => i_stroke => :speed} ~ normal(prev_speed, 0.1)
end

@gen function rejuv_smart(trace, t_prev::Int=0)
    for t in t_prev:-1:0
        i_stroke = trace[:n_strokes] - t
        if i_stroke <= 0 continue end
        # Compute observation indices for stroke
        n_obs = get_args(trace)[1]
        n_points = trace[:strokes => i_stroke => :n_points]
        t_start = i_stroke == 1 ?
            1 : sum(length.(trace[:points_per_stroke][1:i_stroke-1])) + 1
        t_stop = min(t_start + n_points - 1, n_obs)
        if t_stop <= t_start return end
        # Estimate speed from first and last point
        start = trace[:obs][t_start]
        stop = trace[:obs][t_stop]
        est_speed = maximum(abs.(start - stop)) / (t_start - t_stop)
        {:strokes => i_stroke => :speed} ~ normal(est_speed, 0.1)
        # Estimate direction from first and last point
        est_dir = argmax(directions * (stop - start))
        dir_probs = ones(4) * 0.05
        dir_probs[est_dir] = 0.85
        {:strokes => i_stroke => :stroke} ~ categorical(dir_probs)
    end
end

@gen function rejuv_dir(trace)
    n_strokes = trace[:n_strokes]
    if trace[:strokes => n_strokes => :disconnect]
        {:strokes => n_strokes => :stroke} ~ categorical(type_prior)
    else
        prev_stroke = n_strokes == 1 ? 5 : trace[:strokes => n_strokes-1 => :stroke]
        {:strokes => n_strokes => :stroke} ~ categorical(type_transition[prev_stroke,:])
    end
end

@gen function rejuv_init(trace)
    init_x ~ normal(trace[:obs][1][1], 0.25)
    init_y ~ normal(trace[:obs][1][2], 0.25)
end

@gen function rejuv_connect(trace)
    T = get_args(trace)[1]
    i_stroke = trace[:n_strokes]
    points_before_stroke = i_stroke == 1 ?
        0 : sum(length.(trace[:points_per_stroke][1:i_stroke-1]))
    if i_stroke == 1 || T <= points_before_stroke return end
    new_start = trace[:obs][points_before_stroke + 1]
    prev_stop = trace[:obs][points_before_stroke]
    p_disconnect = 1 - exp(-sum(((new_start .- prev_stop)./1.5).^2))
    disconnect = {:strokes => i_stroke => :disconnect} ~ bernoulli(p_disconnect)
    # Sample new stroke location based on observation
    if disconnect
        {:strokes => i_stroke => :start_x} ~ trunc_normal(obs[1], 0.05, 0.0, 5.0)
        {:strokes => i_stroke => :start_y} ~ trunc_normal(obs[2], 0.05, 0.0, 5.0)
    end
end

@gen function rejuv_start(trace, t_prev::Int=0)
    i_stroke = trace[:n_strokes] - t_prev
    if i_stroke > 0 && trace[:strokes => i_stroke => :disconnect]
        {:strokes => i_stroke => :start_x} ~
            trunc_normal(trace[:strokes => i_stroke => :start_x], 0.1, 0.0, 5.0)
        {:strokes => i_stroke => :start_y} ~
            trunc_normal(trace[:strokes => i_stroke => :start_y], 0.1, 0.0, 5.0)
    end
end

## Top-level particle filter ##

function particle_filter(observations, n_particles, ess_thresh=0.5;
                         callback=render_cb!, callback_args...)
    # Concatenate observations to 2D array if necessary
    if observations isa AbstractArray && observations[1] isa AbstractVector
         observations = reduce(hcat, observations) end
    # Initialize particle filter with first observation
    n_obs = size(observations)[2]
    obs_choices = [choicemap((:obs => t => :pt, observations[:, t])) for t=1:n_obs]
    state = pf_initialize(curve_model, (0,), choicemap(),
                          init_proposal, (observations[:, 1],), n_particles)
    # Iterate across timesteps
    for t=1:n_obs
        # Resample and rejuvenate if the effective sample size is too low
        if effective_sample_size(state) < ess_thresh * n_particles
            # Perform residual resampling, pruning low-weight particles
            pf_resample!(state, :residual, priority_fn=x->x/2)
            if 2 <= t <= 10 pf_rejuvenate!(state, mh, (rejuv_init, ())) end
            # pf_rejuvenate!(state, mh, (rejuv_connect, ()))
            for k in 0:2
                pf_rejuvenate!(state, mh, (rejuv_start, (k,)))
                pf_rejuvenate!(state, mh, (rejuv_smart, (k,)))
            end
        end
        pf_rejuvenate!(state, mh, (rejuv_dir, ()), 2)
        pf_rejuvenate!(state, mh, (rejuv_speed, (0,)))
        pf_rejuvenate!(state, mh, (rejuv_speed, (1,)))
        pf_rejuvenate!(state, mh, (rejuv_speed, (2,)))
        # Smart update that proposes new strokes or new points
        obs = observations[:, t]
        pf_update!(state, (t,), (UnknownChange(),), obs_choices[t],
                   fwd_proposal, (t, obs), bwd_proposal, (t, obs), update_transform)
        # Run callback
        if !isnothing(callback)
            callback(t, state, observations; callback_args...)
        end
    end
    return state
end

export particle_filter


text:
using Ledgers
import Test; using Test: @testset, @test
using Assets: USD

function example(withnumbers=true)
    if withnumbers
        group = AccountGroup(Account{USD}, name="Account Group", number="0000000")
        assets = AccountGroup(Account{USD}, name="Assets", number="1000000", parent=group)
        liabilities = AccountGroup(Account{USD}, name="Liabilities", number="2000000", parent=group, isdebit=false)
        cash = Account(USD, name="Cash", number="1010000", parent=assets)
        payable = Account(USD, name="Accounts Payable", number="2010000", parent=liabilities, isdebit=false)
    else
        group = AccountGroup(Account{USD}, name="Account Group")
        assets = AccountGroup(Account{USD}, name="Assets", parent=group)
        liabilities = AccountGroup(Account{USD}, name="Liabilities", parent=group, isdebit=false)
        cash = Account(USD, name="Cash", parent=assets)
        payable = Account(USD, name="Accounts Payable", parent=liabilities, isdebit=false)
    end
    entry = Entry(cash, payable)
    group, assets, liabilities, cash, payable, entry
end

group, assets, liabilities, cash, payable, entry = example()

@testset "Account creation" begin
    @test id(group) isa AccountId
    @test id(assets) isa AccountId
    @test cash.number.value === "1010000"
    @test payable.name === "Accounts Payable"
    @test balance(group) === 0USD
    @test balance(assets) === 0USD
    @test balance(liabilities) === 0USD
end

@testset "Post entry" begin
    amt = 10USD
    post!(entry, amt)
    @test balance(group) === 0USD
    @test balance(assets) === amt
    @test balance(liabilities) === -amt
    @test balance(cash) === amt
    @test balance(payable) === -amt
end


text:
abstract type AbstractSimData end

abstract type SingleSimData{T,N,I} <: AbstractSimData end

@mix struct SimDataMixin{T,N,I<:AbstractArray{T,N},S,St,LSt,B,Ra,Ru,STi,CTi,CFr}
    init::I
    source::S
    dest::S
    sourcestatus::St
    deststatus::St
    localstatus::LSt
    buffers::B
    radius::Ra
    ruleset::Ru
    starttime::STi
    currenttime::CTi
    currentframe::CFr
end

SimDataOrReps = Union{AbstractSimData, Vector{<:AbstractSimData}}

(::Type{T})(data::AbstractSimData) where T <: AbstractSimData =
    T(init(data), source(data), dest(data), sourcestatus(data), deststatus(data),
      localstatus(data), buffers(data), radius(data), ruleset(data), starttime(data),
      currenttime(data), currentframe(data))

# Array interface
Base.length(d::SingleSimData) = length(source(d))
Base.firstindex(d::SingleSimData) = firstindex(source(d))
Base.lastindex(d::SingleSimData) = lastindex(source(d))
Base.size(d::SingleSimData) = size(source(d))
Base.broadcast(d::SingleSimData, args...) = broadcast(source(d), args...)
Base.broadcast!(d::SingleSimData, args...) = broadcast!(dest(d), args...)
Base.iterate(d::SingleSimData, args...) = iterate(source(d), args...)
Base.ndims(::SingleSimData{T,N}) where {T,N} = N
Base.eltype(::SingleSimData{T}) where T = T

# Getters
init(d::SingleSimData) = d.init
source(d::SingleSimData) = d.source
dest(d::SingleSimData) = d.dest
sourcestatus(d::SingleSimData) = d.sourcestatus
deststatus(d::SingleSimData) = d.deststatus
localstatus(d::SingleSimData) = d.localstatus
buffers(d::SingleSimData) = d.buffers
radius(d::SingleSimData) = d.radius
ruleset(d::SingleSimData) = d.ruleset
starttime(d::SingleSimData) = d.starttime
currenttime(d::SingleSimData) = d.currenttime
currentframe(d::SingleSimData) = d.currentframe

# Getters forwarded to ruleset
framesize(d::SingleSimData) = size(init(d))
rules(d::SingleSimData) = rules(ruleset(d))
mask(d::SingleSimData) = mask(ruleset(d))
overflow(d::SingleSimData) = overflow(ruleset(d))
timestep(d::SingleSimData) = timestep(ruleset(d))
cellsize(d::SingleSimData) = cellsize(ruleset(d))

"""
Get the actual current timestep, ie. not variable periods like Month
"""
currenttimestep(d::AbstractSimData) = currenttime(d) + timestep(d) - currenttime(d)

" Simulation data and storage is passed to rules for each timestep "
@SimDataMixin struct SimData{} <: SingleSimData{T,N,I} end

ConstructionBase.constructorof(::Type{SimData}) =
    (init, args...) -> SimData{eltype(init),ndims(init),typeof(init),typeof.(args)...}(init, args...)

"""
Generate simulation data to match a ruleset and init array.
"""
SimData(ruleset::Ruleset, init::AbstractArray, starttime, r=radius(ruleset)) = begin
    # We add one extra row/column so we dont have to worry about
    # special casing the last block
    if r > 0
        hoodsize = 2r + 1
        blocksize = 2r
        source = addpadding(init, r)
        nblocs = indtoblock.(size(source), blocksize) .+ 1
        sourcestatus = BitArray(zeros(Bool, nblocs))
        deststatus = deepcopy(sourcestatus)
        updatestatus!(parent(source), sourcestatus, deststatus, r)

        buffers = [zeros(eltype(init), hoodsize, hoodsize) for i in 1:blocksize]
        localstatus = zeros(Bool, 2, 2)
    else
        source = deepcopy(init)
        sourcestatus = deststatus = true
        buffers = nothing
        localstatus = nothing
    end
    dest = deepcopy(source)
    currentframe = 1
    currenttime = starttime

    SimData(init, source, dest, sourcestatus, deststatus, localstatus, buffers, r,
            ruleset, starttime, currenttime, currentframe)
end

"""
MultipleSimData is used for MultiRuleset models
"""
struct MultiSimData{I,D<:NamedTuple,Ru} <: AbstractSimData
    init::I
    data::D
    ruleset::Ru
end

data(d::MultiSimData) = d.data
ruleset(d::MultiSimData) = d.ruleset
framesize(d::MultiSimData) = framesize(first(data(d)))
interactions(d::MultiSimData) = interactions(ruleset(d))

Base.getindex(d::MultiSimData, key) = begin
    getindex(data(d), key)
end
Base.keys(d::MultiSimData) = keys(data(d))

# Getters
init(d::MultiSimData) = d.init
ruleset(d::MultiSimData) = d.ruleset
data(d::MultiSimData) = d.data

firstruleset(d::MultiSimData) = first(ruleset(ruleset(d)))
firstdata(d::MultiSimData) = first(data(d))

source(d::MultiSimData) = source(firstdata(d))
dest(d::MultiSimData) = dest(firstdata(d))
sourcestatus(d::MultiSimData) = sourcestatus(firstdata(d))
deststatus(d::MultiSimData) = deststatus(firstdata(d))
localstatus(d::MultiSimData) = localstatus(firstdata(d))
buffers(d::MultiSimData) = buffers(firstdata(d))
radius(d::MultiSimData) = radius(firstdata(d))
starttime(d::MultiSimData) = starttime(firstdata(d))
currenttime(d::MultiSimData) = currenttime(firstdata(d))
currentframe(d::MultiSimData) = currentframe(firstdata(d))

# Getters forwarded to ruleset
framesize(d::MultiSimData) = size(first(init(d)))
# rules(d::MultiSimData) = map(rules, ruleset(d))
mask(d::MultiSimData) = mask(firstruleset(d))
overflow(d::MultiSimData) = overflow(firstruleset(d))
timestep(d::MultiSimData) = timestep(firstruleset(d))
cellsize(d::MultiSimData) = cellsize(firstruleset(d))

"""
Swap source and dest arrays. Allways returns regular SimData.
"""
swapsource(data::SingleSimData) = begin
    src = data.source
    dst = data.dest
    @set! data.dest = src
    @set data.source = dst
end
swapsource(d::MultiSimData) = @set d.data = map(swapsource, d.data)

"""
Uptate timestamp
"""
updatetime(data::SimData, f::Integer) = begin
    @set! data.currentframe = f
    @set data.currenttime = timefromframe(data, f)
end
updatetime(data::AbstractVector{<:SimData}, f) = updatetime.(data, f)
updatetime(multidata::MultiSimData, f) =
    @set multidata.data = map(d -> updatetime(d, f), data(multidata))

timefromframe(data::AbstractSimData, f) = starttime(data) + (f - 1) * timestep(data)

"""
Find the maximum radius required by all rules
Add padding around the original init array, offset into the negative
So that the first real cell is still 1, 1
"""
addpadding(init::AbstractArray{T,N}, r) where {T,N} = begin
    sze = size(init)
    paddedsize = sze .+ 2r
    paddedindices = -r + 1:sze[1] + r, -r + 1:sze[2] + r
    sourceparent = similar(init, paddedsize...)
    source = OffsetArray(sourceparent, paddedindices...)
    source .= zero(eltype(source))
    # Copy the init array to he middle section of the source array
    for j in 1:size(init, 2), i in 1:size(init, 1)
        source[i, j] = init[i, j]
    end
    source
end

"""
Initialise the block status array.
This tracks whether anything has to be done in an area of the main array.
"""
updatestatus!(data) = updatestatus!(parent(source(data)), sourcestatus(data), deststatus(data), radius(data))
updatestatus!(source, sourcestatus::Bool, deststatus::Bool, r) = nothing
updatestatus!(source, sourcestatus, deststatus, r) = begin
    blocksize = 2r
    source = parent(source)
    for i in CartesianIndices(source)
        # Mark the status block if there is a non-zero value
        if source[i] != 0
            bi = indtoblock.(Tuple(i), blocksize)
            sourcestatus[bi...] = true
            deststatus[bi...] = true
        end
    end
end

# TODO document these behaviours
initdata!(data::AbstractSimData, ruleset, init::AbstractArray, starttime, nreplicates) =
    initdata!(data, ruleset, starttime)
initdata!(data::AbstractVector{<:AbstractSimData}, ruleset, init::AbstractArray, starttime, nreplicates::Integer) =
    map(d -> initdata!(d, ruleset, init, starttime, nreplicates), data)
initdata!(data::Nothing, ruleset::Ruleset, init, starttime, nreplicates::Nothing) =
    SimData(ruleset, init, starttime)
# initdata!(data::Nothing, ruleset::Ruleset, init, starttime, nreplicates::Integer) =
    # [SimData(ruleset, init, starttime) for r in 1:nreplicates]
initdata!(data::Nothing, multiruleset::MultiRuleset, init::NamedTuple, starttime, nreplicates::Nothing) = begin
    radii = NamedTuple{keys(init)}(radius(multiruleset))
    data = map((rs, ra, in) -> SimData(rs, in, starttime, ra), ruleset(multiruleset), radii, init)
    MultiSimData(init, data, multiruleset)
end
# initdata!(multidata::MultiSimData, multiruleset::MultiRuleset, init, starttime, nreplicates::Nothing) =
    # MultiSimData(map((d, rs) -> initdata!(d, rs, starttime), data(MultiSimData), ruleset(multiruleset))),
                 # interactions(multiruleset))
initdata!(data::AbstractSimData, ruleset::Ruleset, starttime) = begin
    for j in 1:framesize(data)[2], i in 1:framesize(data)[1]
        @inbounds source(data)[i, j] = dest(data)[i, j] = init(data)[i, j]
    end
    updatestatus!(data)
    @set! data.ruleset = ruleset
    @set! data.starttime = starttime
    data
end

indtoblock(x, blocksize) = (x - 1) ÷ blocksize + 1
blocktoind(x, blocksize) = (x - 1) * blocksize + 1

Base.@propagate_inbounds Base.getindex(d::SimData, I...) = getindex(source(d), I...)

" Simulation data and storage is passed to rules for each timestep "
@SimDataMixin struct WritableSimData{} <: SingleSimData{T,N,I} end

ConstructionBase.constructorof(::Type{WritableSimData}) =
    (init, args...) -> SimData{eltype(init),ndims(init),typeof(init),typeof.(args)...}(init, args...)

Base.@propagate_inbounds Base.setindex!(d::WritableSimData, x, I...) = begin
    r = radius(d)
    if r > 0
        bi = indtoblock.(I .+ r, 2r)
        deststatus(d)[bi...] = true
    end
    isnan(x) && error("NaN in setindex: ", (d, I))
    dest(d)[I...] = x
end
Base.@propagate_inbounds Base.getindex(d::WritableSimData, I...) = getindex(dest(d), I...)


text:
function cleansyntax(line)
    line = strip(line)
    if line == ""
        return "#=blank line=#"
    elseif line[1] == '#'
        return line
    elseif line[1] == ';'
        line = replace(line, r"^;[;]+" => ";")
    else # active option
        o, p = splitline(line)
        line = p == nothing ? uppercase(o) : uppercase(o) * " " * p
    end
    join(filter(c -> isascii(c[1]) && !iscntrl(c[1]), split(line, "")), "")
end

"""
    Run to clean up the configuration file lines prior to other function application.
"""
cleansyntax!(lines) = (for (i, li) in enumerate(lines) lines[i] = cleansyntax(li) end; lines)

isdisabled(line) = startswith(line, [';', '#'])
isenabled(line) = !isdisabled(line)

function splitline(line)
    arr = split(line, r"\s+", limit=2)
    if length(arr) < 2
        return (arr[1], nothing)
    end
    return (arr[1], arr[2])
end

function disable!(lines, opt)
    for (i, li) in enumerate(lines)
        if isenabled(li) && splitline(li)[1] == uppercase(opt)
            lines[i] = ";" * li
            break # note: only first one found is disabled
        end
    end
    lines
end

function enable!(lines, opt)
    for (i, li) in enumerate(lines)
        if isdisabled(li)
            s = li[2:end]
            if splitline(s)[1] == uppercase(opt)
                lines[i] = s
                break # note: only first one found is enabled
            end
        end
    end
    lines
end

function changeparam!(lines, opt, newparam)
    for (i, li) in enumerate(lines)
        if isenabled(li)
            o, p = splitline(li)
            if o == opt
                lines[i] = o * " " * string(newparam)
                break # note: only first one found is changed
            end
        end
    end
    lines
end

function activecfg(lines)
    cfgdict = Dict()
    for li in lines
        if isenabled(li)
            o, p = splitline(li)
            cfgdict[o] = p
        end
    end
    cfgdict
end

const filename = "fruit.cfg"
const cfh = open(filename)
const cfglines = cleansyntax!(readlines(cfh))
close(cfh)

const cfg = activecfg(cfglines)

disable!(cfglines, "NEEDSPEELING")
enable!(cfglines, "SEEDSREMOVED")
changeparam!(cfglines, "NUMBEROFBANANAS", 1024)

if !haskey(cfg, "NUMBEROFSTRAWBERRIES")
    push!(cfglines, "NUMBEROFSTRAWBERRIES 62000")
end
cfg["NUMBEROFSTRAWBERRIES"] = 62000
changeparam!(cfglines, "NUMBEROFSTRAWBERRIES", 62000)

const cfw = open(filename, "w")
for li in cfglines
    if li != ""
        if li == "#=blank line=#"
            li = ""
        end
        write(cfw, li * "\n")
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Metis,CSparse
using Base.Test

function specialize(A::SparseMatrixCSC)
    issym(A) && return Symmetric(triu!(A),:U)
    ishermitian(A) && return Hermitian(triu!(A),:U)
    istril(A) && return LowerTriangular(A)
    istriu(A) && return UpperTriangular(A)
    A
end

cd(joinpath(dirname(@__FILE__),"..","Matrix")) do
    global mats = [specialize(testmat(nm)) for nm in readdir()]
end

@test size(mats[1]) == (219,85)

@test size(mats[2]) == (48,48)
const o48 = ones(48)
@test_approx_eq mats[2]*(mats[2]\o48) o48
@test_approx_eq Ac_mul_B(mats[2],Ac_ldiv_B(mats[2],o48)) o48

const bc01t = mats[2]'
@test istriu(bc01t)
@test isa(bc01t,UpperTriangular)
@test size(bc01t) == (48,48)
@test_approx_eq bc01t*(bc01t\o48) o48
@test_approx_eq Ac_mul_B(bc01t,Ac_ldiv_B(bc01t,o48)) o48

@test size(mats[3]) == (4884,4884)
const o4884 = ones(4884)
@test_approx_eq mats[3]*(mats[3]\o4884) o4884
@test_approx_eq Ac_mul_B(mats[3],Ac_ldiv_B(mats[3],o4884)) o4884

const bc16t = mats[3]'
@test istriu(bc16t)
@test isa(bc16t,UpperTriangular)
@test size(bc16t) == (4884,4884)
@test_approx_eq bc16t*(bc16t\o4884) o4884
@test_approx_eq Ac_mul_B(bc16t,Ac_ldiv_B(bc16t,o4884)) o4884

const sc1 = Symmetric(triu(mats[2].data'mats[2].data),:U)
const perm1,iperm1 = nodeND(sc1)
const tr1 = etree(sc1)
const sc2 = symperm(sc1,iperm1)
const tr2,postperm = etree(sc2,true)
const sc3 = symperm(sc2,invperm(postperm))
const tr3 = etree(sc3)
@show(tr1)
@show(tr2)
@show(tr3)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

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

Sequential Gaussian simulation.

## Parameters

* `variogram` - Variogram model (default to `GaussianVariogram()`)
* `mean`      - Simple Kriging mean
* `degree`    - Universal Kriging degree
* `drifts`    - External Drift Kriging drift functions

Latter options override former options. For example, by specifying
`drifts`, the user is telling the algorithm to ignore `degree` and
`mean`. If no option is specified, Ordinary Kriging is used by
default with the `variogram` only.

* `neighborhood` - Neighborhood on which to search neighbors
* `maxneighbors` - Maximum number of neighbors (default to 10)
* `path`         - Simulation path (default to `LinearPath()`)
* `mapping`      - Data mapping method (default to `NearestMapping()`)

For each location in the simulation `path`, a maximum number of
neighbors `maxneighbors` is used to fit a Gaussian distribution.
The neighbors are searched according to a `neighborhood`.

### References

Gomez-Hernandez & Journel 1993. *Joint Sequential Simulation of
MultiGaussian Fields*
"""
@simsolver SGS begin
  @param variogram = GaussianVariogram()
  @param mean = nothing
  @param degree = nothing
  @param drifts = nothing
  @param neighborhood
  @param minneighbors = 1
  @param maxneighbors = 10
  @param path = nothing
  @param mapping = NearestMapping()
end

function preprocess(problem::SimulationProblem, solver::SGS)
  # retrieve problem info
  pdomain = domain(problem)

  params = []
  for covars in covariables(problem, solver)
    for var in covars.names
      # get user parameters
      varparams = covars.params[(var,)]

      # determine which Kriging variant to use
      if varparams.drifts ≠ nothing
        estimator = ExternalDriftKriging(varparams.variogram, varparams.drifts)
      elseif varparams.degree ≠ nothing
        estimator = UniversalKriging(varparams.variogram, varparams.degree, embeddim(pdomain))
      elseif varparams.mean ≠ nothing
        estimator = SimpleKriging(varparams.variogram, varparams.mean)
      else
        estimator = OrdinaryKriging(varparams.variogram)
      end

      # determine marginal distribution
      marginal = Normal()

      # determine simulation path
      path = varparams.path ≠ nothing ? varparams.path : RandomPath()

      # determine data mapping
      mapping = varparams.mapping

      # equivalent parameters for SeqSim solver
      param = (estimator=estimator,
               neighborhood=varparams.neighborhood,
               minneighbors=varparams.minneighbors,
               maxneighbors=varparams.maxneighbors,
               marginal=marginal, path=path,
               mapping=mapping)

      push!(params, var => param)
    end
  end

  preprocess(problem, SeqSim(params...))
end

solvesingle(problem::SimulationProblem, covars::NamedTuple, ::SGS, preproc) =
  solvesingle(problem, covars, SeqSim(), preproc)


text:
module VortexHelperKillLightningTrigger
using ..Ahorn, Maple

@mapdef Trigger "VortexHelper/KillLightningTrigger" KillLightningTrigger(
					x::Integer, y::Integer,
					width::Integer = Maple.defaultTriggerWidth,
					height::Integer = Maple.defaultTriggerHeight,
					permanent::Bool = false)

const placements = Ahorn.PlacementDict(
	"Kill Lightning Trigger (Vortex Helper)" => Ahorn.EntityPlacement(
		KillLightningTrigger,
		"rectangle"
	)
)

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
# =================================================================
# Bouncing-ball model
# See https://juliareach.github.io/SX.jl/latest/examples/bball.html
# =================================================================

using Reachability, HybridSystems, MathematicalSystems, LazySets, LinearAlgebra
using LazySets: HalfSpace  # resolve name-space conflicts with Polyhedra

function bouncing_ball()
    # automaton structure
    automaton = LightAutomaton(1)

    # mode 1
    A = [0.0 1.0; 0.0 0.0]
    B = reshape([0.0, -1.0], (2, 1))
    U = Singleton([1.0])
    inv = HalfSpace([-1.0, 0.0], 0.0) # x >= 0
    m1 = ConstrainedLinearControlContinuousSystem(A, B, inv, ConstantInput(U))

    # modes
    modes = [m1]

    # transition from m1 to m1 (self-loop)
    add_transition!(automaton, 1, 1, 1)
    A = [1.0 0.0; 0.0 -0.75]
    guard = HPolyhedron([HalfSpace([0.0, 1.0], 0.0),   # v ≤ 0
                         HalfSpace([-1.0, 0.0], 0.0),  # x ≥ 0
                         HalfSpace([1.0, 0.0], 0.0)])  # x ≤ 0
    t1 = ConstrainedLinearMap(A, guard)

    # transition annotations
    resetmaps = [t1]

    # switching
    switchings = [AutonomousSwitching()]

    ℋ = HybridSystem(automaton, modes, resetmaps, switchings)

    # initial condition in mode 1
    X0 = Hyperrectangle(low=[10.0, 0.0], high=[10.2, 0.0])
    initial_condition = [(1, X0)]

    system = InitialValueProblem(ℋ, initial_condition)

    options = Options(:mode=>"reach", :T=>5.0,
                      :plot_vars=>[1, 2], :max_jumps=>1, :verbosity=>1)

    return (system, options)
end


text:
println(Seq())
println(e"a" + E"b" | p"c" > tuple)

println("print ok")


text:
using LinearAlgebra, MixedModels, Random, Test

@testset "vectorterm" begin
    trm = MatrixTerm(ones(10))
    @test size(trm) == (10, 1)
    @test trm'trm == 10.0 * ones(1, 1)
    @test trm.cnames == [""]
    @test trm.x === trm.wtx
    wts = rand(MersenneTwister(1234321), 10)
    MixedModels.reweight!(trm, wts)
    @test !(trm.x === trm.wtx)
    @test trm'trm ≈ sum(abs2, wts) * ones(1, 1)
    @test mul!(Vector{eltype(trm)}(undef, 10), trm, ones(1)) == ones(10)
end

@testset "matrixterm" begin
    trm = MatrixTerm(hcat(ones(30), repeat(0:9, outer = 3)), ["(Intercept)", "U"])
    piv = trm.piv
    ipiv = invperm(piv)
    @test size(trm) == (30, 2)
    @test trm.x === trm.wtx
    prd = trm'trm
    @test typeof(prd) == Matrix{Float64}
    @test prd == [30.0 135.0; 135.0 855.0][piv, piv]
    wts = rand(MersenneTwister(123454321), 30)
    MixedModels.reweight!(trm, wts)
    @test mul!(prd, trm', trm)[ipiv[1], ipiv[1]] ≈ sum(abs2, wts)
end


text:
@testset "Contrast Limited Adaptive Histogram Equalisation" begin
    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        # Create a graylevel ramp
        img = Gray{Float32}.([i/255.0 for i = 64:127, j = 1:64])
        img = T.(img)
        for rblocks in 2:4
            for cblocks in 2:4
                # There should be no change in intensity along a row.
                algo = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1,
                                            rblocks = rblocks, cblocks = cblocks,
                                            clip = 0.0)
                imgeq = adjust_histogram(img, algo)
                target = repeat(imgeq[:,1], inner=(1, 1), outer=(1, 64))
                @test all(imgeq .≈ target)

                # There should be no change in intensity along a column
                imgeq = adjust_histogram(transpose(img), algo)
                target =  transpose(repeat(imgeq[1, :], inner=(1, 1), outer=(1, 64)))
                @test all(imgeq .≈ target)
            end
        end

    end

    #=
    When rblocks and cblocks equal one and clip is zero,  then the method boils
    down to canonical histogram equalization. Hence, we repeat the same tests
    here that we have in histgoram_equalization.jl.
    =#

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        #=
        Create an image that spans a narrow graylevel range. Then quantize
        the 256 bins down to 32 and determine how many bins have non-zero
        counts.
        =#

        img = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = T.(img)
        _, counts_before = build_histogram(img, 32, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Equalize the image histogram. Then quantize the 256 bins down to 32
        and verify that all 32 bins have non-zero counts. This will confirm
        that the dynamic range of the original image has been increased.
        =#
        algo = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1,
                                    rblocks = 1, cblocks = 1, clip = 0.0)
        imgeq = adjust_histogram(img, algo)
        edges, counts_after = build_histogram(imgeq, 32, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 32
    end

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float32}, RGB{Float64})
        #=
        Create a color image that spans a narrow graylevel range.  Then
        quantize the 256 bins down to 32 and determine how many bins have
        non-zero counts.
        =#

        imgg = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
        img = colorview(RGB,imgg,imgg,imgg)
        img = T.(img)
        _, counts_before = build_histogram(img, 32, minval = 0, maxval = 1)
        nonzero_before = sum(counts_before .!= 0)

        #=
        Equalize the histogram. Then quantize the 256 bins down to 32 and
        verify that all 32 bins have non-zero counts. This will confirm that
        the dynamic range of the original image has been increased.
        =#
        algo = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1,
                                    rblocks = 1, cblocks = 1, clip = 0.0)
        imgeq = adjust_histogram(img, algo)
        edges, counts_after = build_histogram(imgeq, 32, minval = 0, maxval = 1)
        nonzero_after = sum(counts_after .!= 0)
        @test nonzero_before < nonzero_after
        @test nonzero_after == 32
    end

    # Verify that the minimum and maximum values of the equalised image match the
    # specified minimum and maximum values, i.e. that the intensities of the equalised
    # image are in the interval [minvalue, maxvalue].
    algo = AdaptiveEqualization(nbins = 256, minval = 64, maxval = 128,
                                rblocks = 1, cblocks = 1, clip = 0.0)
    imgeq = adjust_histogram(hcat(collect(0:1:255), collect(0:1:255)), algo)
    @test all(imgeq[1:65, :] .== 64)
    @test all(imgeq[128+1:end, :] .== 128)

    algo = AdaptiveEqualization(nbins = 256, minval = 64/255, maxval = 128/255,
                                rblocks = 1, cblocks = 1, clip = 0.0)
    imgeq = adjust_histogram(hcat(collect(0:1/255:1), collect(0:1/255:1)), algo)
    @test all(imgeq[1:65, :] .== 64/255)
    @test all(imgeq[128+1:end, :] .== 128/255)

    # Verify that increasing the clip value reduces the strength of the contrast
    # adjustment
    img = Gray{Float32}.(testimage("mandril_gray"))
    # As we increase the clip factor the contrast adjustment is diminished and we
    # should approach the original image.
    algo₁ = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1, rblocks = 1, cblocks = 1, clip = 0.0)
    algo₂ = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1, rblocks = 1, cblocks = 1, clip = 0.25)
    algo₃ = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1, rblocks = 1, cblocks = 1, clip = 0.5)
    algo₄ = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1, rblocks = 1, cblocks = 1, clip = 0.75)
    algo₅ = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1, rblocks = 1, cblocks = 1, clip = 1.0)
    algos = [algo₁, algo₂, algo₃, algo₄, algo₅]
    differences = [norm(adjust_histogram(img, algos[i]) .- img) for i = 1:5]
    @test all([differences[i] > differences[i+1] for i = 1:4])
    @test last(differences) < 2.8

    # Verify that auotmatically converting a color image and storing the histogram
    # adjusted result "inplace" works correctly.
    algo = AdaptiveEqualization(nbins = 256, minval = 0, maxval = 1, rblocks = 1, cblocks = 1, clip = 0)
    imgg₁ = Gray{Float32}.([i/255.0 for i = 64:128, j = 1:10])
    imgg₂ = copy(imgg₁)
    img = colorview(RGB,imgg₁,imgg₁,imgg₁)
    imgeq₁ = adjust_histogram!(imgg₁,img, algo)
    imgeq₂ = adjust_histogram(imgg₂, algo)
    @test norm(imgeq₁ .- imgeq₂) ≈ 0.0
end


text:
using OpenQASM, YaoBlocksQASM, Yao
using Test

@testset "YaoBlocksQASM.jl" begin
    qasm = """
        OPENQASM 2.0;
        include "qelib1.inc";
        qreg q[3];
        creg c1[3];
        h q[0];
        CX q[1],q[2];
        cy q[1],q[0];
        cz q[0],q[2];
        x q[0];
        swap q[1],q[2];
        id q[0];
        t q[1];
        rz(0.7) q[2];
        z q[0];
        p(0.7) q[1];
        ry(0.7) q[2];
        y q[0];
        rx(0.7) q[1];
        measure q[0] -> c1[0];
        measure q[1] -> c1[1];
        measure q[2] -> c1[2];
    """

    ast = OpenQASM.parse(qasm)
    qc = chain(
        3,
        put(1 => H),
        control(2, 3 => X),
        control(2, 1 => Y),
        control(1, 3 => Z),
        put(1 => X),
        swap(2, 3),
        put(1 => I2),
        put(2 => T),
        put(3 => Rz(0.7)),
        put(1 => Z),
        put(2 => shift(0.7)),
        put(3 => Ry(0.7)),
        put(1 => Y),
        put(2 => Rx(0.7)),
        Yao.Measure(3),
    )
    ast1 = convert_to_qasm(qc, 1)

    @test ast.version == ast1.version

    for i = 1:length(ast.prog)
        @test "$(ast.prog[i])" == "$(ast1.prog[i])"
    end
end

max_stars_repo_name: Sov-trotter/YaoBlocksQASM

max_issues_repo_name: Sov-trotter/YaoBlocksQASM

max_forks_repo_name: Sov-trotter/YaoBlocksQASM


text:
"""
Julia の binomialtest() は，母比率 ≠ 0.5 で両側検定の場合に片側 p 値を 2 倍した結果を返す。

他の多くの統計プログラムと同じ結果を得るためのオプション設定がないので，以下のようなプログラムを使用すればよい。
"""

using Distributions
function binomial_pvalue(x, n, p=0.5; tail=:both)
    obj = Binomial(n, p)
    tail == :left  && return cdf(obj, x)
    tail == :right && return ccdf(obj, x - 1)
    p == 0 && return Float64(x == 0)
    p == 1 && return Float64(x == n)
    m = n * p
    x == m && return 1
    limit = 1.0000001pdf(obj, x)
    if x < m
        y = sum(pdf.(obj, ceil(Int, m):n) .<= limit)
        cdf(obj, x) + ccdf(obj, n - y)
    else
        y = sum(pdf.(obj, 0:floor(Int, m)) .<= limit)
        cdf(obj, y - 1) + ccdf(obj, x - 1)
    end
end;

"""
使用法

julia> using HypothesisTests

julia> BinomialTest(12, 20, 0.8)
Binomial test
-------------
Population details:
    parameter of interest:   Probability of success
    value under h_0:         0.8
    point estimate:          0.6
    95% confidence interval: (0.3605, 0.8088)

Test summary:
    outcome with 95% confidence: fail to reject h_0
    two-sided p-value:           0.0643

Details:
    number of observations: 20
    number of successes:    12

julia> binomial_pvalue(12, 20, 0.8)
0.04367187812694369

"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
@testset "query synthesis" begin
    @testset "initialization" begin
        epsilon = fill(0.1, 2)
        data = ones(2, 10)
        labels = fill(:U, 10)
        history = Array{Float64, 2}[]

        @testset "RandomQss" begin
            @test_throws ArgumentError RandomQss(epsilon=-0.1)
            @test_throws ArgumentError RandomQss(epsilon=-ones(2))
            @test RandomQss(epsilon=1.0) !== nothing
            @test RandomQss(epsilon=ones(2)) !== nothing

        end

        @testset "RandomOutlierQss" begin
            occ = SVDD.RandomOCClassifier(data)
            @test_throws ArgumentError RandomOutlierQss(occ, max_tries=0)
            @test_throws ArgumentError RandomOutlierQss(occ, epsilon=-0.1)
            @test_throws ArgumentError RandomOutlierQss(occ, epsilon=-ones(2))
            @test RandomOutlierQss(occ, epsilon=1.0) !== nothing
            @test RandomOutlierQss(occ, epsilon=ones(2)) !== nothing
        end
    end

    @testset "get query" begin
        data = rand(2, 10)
        labels = fill(:U, size(data, 2))
        history = Vector{Array{Float64, 2}}()
        optimizer = ParticleSwarmOptimization()

        occ = SVDD.VanillaSVDD(data)
        init_strategy = SVDD.SimpleCombinedStrategy(SVDD.FixedGammaStrategy(GaussianKernel(2.0)), SVDD.FixedCStrategy(1))
        SVDD.initialize!(occ, init_strategy)
        SVDD.fit!(occ, TEST_SOLVER)
        for qs_type in [RandomQss, RandomOutlierQss, TestQss, DecisionBoundaryQss, NaiveExplorativeMarginQss]
            @testset "$qs_type" begin
                qs = initialize_qs(qs_type, occ, data, Dict(:optimizer => optimizer))
                query = get_query_objects(qs, data, labels, history)
                @test size(query) == (size(data, 1), 1)
                if qs_type == NaiveExplorativeMarginQss
                    labels[end] = :Lout
                    query = get_query_objects(qs, data, labels, history)
                    @test size(query) == (size(data, 1), 1)
                end
            end
        end
        @testset "ExplorativeMarginQss" begin
            qs = initialize_qs(ExplorativeMarginQss, occ, data,
                               Dict(:solver => TEST_SOLVER,
                                    :optimizer => optimizer))
            query = get_query_objects(qs, data, labels, history)
            @test size(query) == (size(data, 1), 1)
            labels[end] = :Lout
            query = get_query_objects(qs, data, labels, history)
            @test size(query) == (size(data, 1), 1)
        end
    end
end


text:
@testset "RCF" begin

  @testset "abelian extensions && related examples" begin

    Qx, x = PolynomialRing(FlintQQ, "x")
    K, a = NumberField(x - 1, "a")
    l = Hecke.abelian_normal_extensions(K, Int[2,2], fmpz(10)^4)
    @test length(l)==47
    l1 = collect(Hecke.C22_extensions(10^4))
    @test length(l1)==47
    @test length(abelian_fields(FlintQQ, [3], fmpz(10)^3)) == 5

    K, a = number_field(x^2+1, "a")
    auts = small_generating_set(automorphisms(K, copy = false))
    l = Hecke.abelian_normal_extensions(K, Int[2], fmpz(10)^5)
    @test length(l) == 41
    l1 = typeof(l[1])[x for x in l if abs(discriminant(x, FlintQQ)) < 5*10^3]
    lnn = Hecke.abelian_extensions(K, Int[2], 5*fmpz(10)^3)
    ln = typeof(l[1])[x for x in lnn if isnormal(x)]
    @test length(ln) == length(l1)
    for x in l[1:5]
      L = number_field(x)
      autsL = Hecke.absolute_automorphism_group(x, auts)
      K, autos = Hecke._relative_to_absolute(L, autsL)
      @test length(autos)==2
      y = small_generating_set(closure(autos, *))
      @test length(y)==1 || length(y)==2
      if length(y) == 2
        @test y[1]*y[2] == y[2]*y[1]
      else
        @test y[1]^2 != id_hom(K)
      end
    end
  end

end


text:
# ---
# title: 421. Maximum XOR of Two Numbers in an Array
# id: problem421
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation, Trie
# link: <https://leetcode.com/problems/maximum-xor-of-two-numbers-in-an-array/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, return _the maximum result of`nums[i] XOR
# nums[j]`_, where `0 ≤ i ≤ j < n`.
#
# **Follow up:** Could you do this in `O(n)` runtime?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [3,10,5,25,2,8]
#     Output: 28
#     Explanation: The maximum result is 5 XOR 25 = 28.
#
# **Example 2:**
#
#
#
#     Input: nums = [0]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: nums = [2,4]
#     Output: 6
#
#
# **Example 4:**
#
#
#
#     Input: nums = [8,10,2]
#     Output: 10
#
#
# **Example 5:**
#
#
#
#     Input: nums = [14,70,53,83,49,91,36,80,92,51,66,70]
#     Output: 127
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2 * 104`
#   * `0 <= nums[i] <= 231 - 1`
#
#
## @lc code=start
using LeetCode

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


text:
using SerializeTables, Test, Tables, CodecZlib

"Randomly return `missing` instead if the argument 1% of the time, otherwise the argument."
randmissing(x) = rand() < 0.01 ? missing : x

@testset "write and read small dataset" begin
    tmp = tempname()
    table = collect((a = i, b = Float64(i), c = 'a' + (i - 1)) for i in 1:10)
    serialize_table_rows(tmp, table)
    rowsitr = deserialize_table_rows(tmp)
    @test Tables.schema(rowsitr) == Tables.schema(Tables.rows(table))
    @test rowtable(rowsitr) == table
end

@testset "write and read large dataset w/ missing data" begin
    N = 10^6
    a = randmissing.(rand(rand(Float64, 100), N))
    b = randmissing.(rand(1:10000, N))
    tmp = tempname()
    table = columntable((a = a, b = b))
    serialize_table_rows(GzipCompressor(; level = 1), tmp, table)
    @info("compression",
          rat1 = filesize(tmp) / (N*(2+sizeof(Int)+sizeof(Float64))), # crude, but fast
          rat2 = filesize(tmp) / (Base.summarysize(a) + Base.summarysize(b))) # slow
    rowsitr = deserialize_table_rows(GzipDecompressor(), tmp)
    @test Tables.schema(rowsitr) == Tables.schema(table)
    @test all(rowtable(rowsitr) .≡ rowtable(table))
end