text:
 @printf("\nTesting linear operators syntax\n")

###### ' ######
n,m = 5,3
A = randn(n,m)
B = randn(n,m)
C = randn(n,m)
x1 = randn(m)
x2 = randn(n)
opA = MatrixOp(A)
opB = MatrixOp(B)
opC = MatrixOp(C)

opt = opA'
y1 = opt*x2
y2 = A'*x2
@test norm(y1-y2) < 1e-9

######+,-######
opp = +opA
y1 = opp*x1
y2 = A*x1
@test norm(y1-y2) < 1e-9

opm = -opA
y1 = opm*x1
y2 = -A*x1
@test norm(y1-y2) < 1e-9

ops = opA+opB
y1 = ops*x1
y2 = A*x1+B*x1
@test norm(y1-y2) < 1e-9

ops = opA-opB
y1 = ops*x1
y2 = A*x1-B*x1
@test norm(y1-y2) < 1e-9

ops = opA+opB
opss = ops+opC
y1 = opss*x1
y2 = A*x1+B*x1+C*x1
@test norm(y1-y2) < 1e-9

ops = opB+opC
opss = opA-ops
y1 = opss*x1
y2 = A*x1-B*x1-C*x1
@test norm(y1-y2) < 1e-9

ops = opB+opC
opss = ops-opA
y1 = opss*x1
y2 = -A*x1+B*x1+C*x1
@test norm(y1-y2) < 1e-9

###### * ######

n,m,l = 5,3,7
A = randn(n,m)
B = randn(l,n)
x1 = randn(m)
opA = MatrixOp(A)
opB = MatrixOp(B)
Id = Eye(m)
alpha = pi
beta  = 4

ops = alpha*opA
y1 = ops*x1
y2 = alpha*A*x1
@test norm(y1-y2) < 1e-9

ops = alpha*opA
opss = beta*ops
y1 = opss*x1
y2 = beta*alpha*A*x1
@test norm(y1-y2) < 1e-9

ops1 = alpha*opA
ops2 = beta*opB
opss = ops2*ops1
y1 = opss*x1
y2 = beta*B*alpha*A*x1
@test norm(y1-y2) < 1e-9

opc = opB*opA
y1 = opc*x1
y2 = B*A*x1
@test norm(y1-y2) < 1e-9

opc = Id*opA
y1 = opc*x1
y2 = A*x1
@test norm(y1-y2) < 1e-9

opc = opA*Id
y1 = opc*x1
y2 = A*x1
@test norm(y1-y2) < 1e-9

opc = Id*Id
y1 = opc*x1
y2 = x1
@test norm(y1-y2) < 1e-9

###### getindex ######

ops = opA[1:n-1]
y1 = ops*x1
y2 = (A*x1)[1:n-1]
@test norm(y1-y2) < 1e-9

opF = DCT(n,m,l)
x3 = randn(n,m,l)
ops = opF[1:n-1,:,2:l]
y1 = ops*x3
y2 = dct(x3)[1:n-1,:,2:l]
@test norm(y1-y2) < 1e-9

opF = DCT(n,m,l)
x3 = randn(n,m,l)
ops = opF[1:n-1,2:m,1]
y1 = ops*x3
y2 = dct(x3)[1:n-1,2:m,1]
@test norm(y1-y2) < 1e-9

opV = Variation(n,m,l)
ops = opV[1:4]
y1 = ops*x3
y2 = (opV*x3)[1:4]
@test norm(y1-y2) < 1e-9

ops = (opB*opA)[1:l-1]
y1 = ops*x1
y2 = (B*A*x1)[1:l-1]
@test norm(y1-y2) < 1e-9

ops = (10.0*opA)[1:n-1]
y1 = ops*x1
y2 = (10*A*x1)[1:n-1]
@test norm(y1-y2) < 1e-9

#slicing HCAT

n,m1,m2,m3 = 5,6,7,8
A = randn(n,m1)
B = randn(n,m2)
C = randn(n,m3)
x1 = randn(m1)
x2 = randn(m2)
x3 = randn(m3)
opA = MatrixOp(A)
opB = MatrixOp(B)
opC = MatrixOp(C)
opH = HCAT(opA,opB,opC)
opH2 = opH[1:2]
y1 = opH2*ArrayPartition(x1,x2)
y2 = A*x1+B*x2
@test all(norm.(y1 .- y2) .<= 1e-12)
opH3 = opH[3]
y1 = opH3*x3
y2 = C*x3
@test all(norm.(y1 .- y2) .<= 1e-12)

opHperm = opH[[3,1,2]]
@test norm(opH*ArrayPartition(x1,x2,x3) - opHperm*ArrayPartition(x3,x1,x2)) <1e-12

@test opHperm[1] == opC
@test opHperm[2] == opA
@test opHperm[3] == opB

opHperm = opH[[3,1]]
@test norm(opC*x3+opA*x1 - opHperm*ArrayPartition(x3,x1)) <1e-12

# slicing Affine add of HCAT
d = randn(n)
opHA = AffineAdd(opH,d)
@test norm(opHA[1]*x1-(A*x1+d)) < 1e-9
@test norm(opHA[2]*x2-(B*x2+d)) < 1e-9
@test norm(opHA[3]*x3-(C*x3+d)) < 1e-9

m4 = 9
x4 = randn(m4)
D = randn(n,n)
E = randn(n,m4)
opD = MatrixOp(D)
opE = MatrixOp(E)
opCH = opD*opH

opHCH = HCAT(opCH,opE)

opH4 = opHCH[4]
@test opH4 == opE
@test_throws ErrorException  opHCH[1]
@test_throws ErrorException  opHCH[1:2]

#slicing VCAT

n1,n2,n3,m = 5,6,7,8
A = randn(n1,m)
B = randn(n2,m)
C = randn(n3,m)
x1 = randn(m)
opA = MatrixOp(A)
opB = MatrixOp(B)
opC = MatrixOp(C)
opV = VCAT(opA,opB,opC)
opV2 = opV[1:2]
y1 = opV2*x1
y2 = ArrayPartition(A*x1,B*x1)
@test norm(y1 - y2) <= 1e-12
opV3 = opV[3]
y1 = opV3*x3
y2 = C*x3
@test norm(y1 - y2) <= 1e-12

###### hcat ######

n,m1,m2 = 5,6,7
A = randn(n,m1)
B = randn(n,m2)
opA = MatrixOp(A)
opB = MatrixOp(B)
opH = [opA opB]
x1 = randn(m1)
x2 = randn(m2)
y1 = opH*ArrayPartition(x1,x2)
y2 = [A B]*[x1;x2]
@test norm(y1 - y2) <= 1e-12

opHH = [opH opB]
y1 = opHH*ArrayPartition(x1, x2, x2)
y2 = [A B B]*[x1;x2;x2]
@test norm(y1 - y2) <= 1e-12

###### vcat ######

n,m1,m2 = 5,6,7
A = randn(m1,n)
B = randn(m2,n)
opA = MatrixOp(A)
opB = MatrixOp(B)
opH = [opA; opB]
x1 = randn(n)
y1 = opH*x1
y2 = ArrayPartition(A*x1,B*x1)
@test norm(y1 - y2) <= 1e-12

opVV = [opA; opH]
y1 = opVV*x1
y2 = ArrayPartition(A*x1, A*x1, B*x1)
@test norm(y1 - y2) <= 1e-12

###### reshape ######
n,m =  10,5
A = randn(n,m)
x1 = randn(m)
opA = MatrixOp(A)
opR = reshape(opA,2,5)
opR = reshape(opA,(2,5))
y1 = opR*x1
y2 = reshape(A*x1,2,5)
@test norm(y1 - y2) <= 1e-12

# testing ndims & ndoms
L = Variation((3,4,5))
@test ndims(L) == (2,3)
@test ndims(L,1) == 2
@test ndims(L,2) == 3
@test ndoms(L) == (1,1)
H = hcat(L,L)
@test ndims(H) == (2,(3,3))
@test ndims(H,1) == 2
@test ndims(H,2) == (3,3)
@test ndoms(H) == (1,2)
@test ndoms(H,1) == 1
@test ndoms(H,2) == 2
D = DCAT(L,L)
@test ndims(D) == ((2,2),(3,3))
@test ndoms(D) == (2,2)

###### jacobian ######
n,m =  10,5
A = MatrixOp(randn(n,m))
B = Sigmoid(Float64,(n,),100.)
op = B*A
J = jacobian(op,randn(m))

#### convert ####
L = Eye(10)
LL = convert(AbstractOperator, Float64, (10,), L)
@test LL == L

LL = convert(LinearOperator, Float64, (10,), L)
@test LL == L

@test_throws MethodError convert(NonLinearOperator, Float64, (10,), L)

### displacement ###
L = Eye(10)
@test displacement(L) == 0.0


text:
"""
Struct that holds the boundary arrays in grid-point space

    geopot_surf     ::Array{Complex{NF},2}  # spectral surface geopotential (orography * gravity) [m^2/s^2]
    landsea_mask    ::Array{NF,2}           # land-sea mask, grid-point
    albedo          ::Array{NF,2}           # annual mean surface albedo, grid-point
"""
struct Boundaries{NF<:AbstractFloat}        # number format NF
    geopot_surf     ::Array{Complex{NF},2}  # surface geopotential (i.e. orography) [m^2/s^2]
    landsea_mask    ::Array{NF,2}           # land-sea mask
    albedo          ::Array{NF,2}           # annual mean surface albedo
end

""" Generator function for a Boundaries struct. Loads the boundary conditions,
orography, land-sea mask and albedo from an netCDF file and stores the in a
`Boundaries` struct."""
function Boundaries(P::Parameters,
                    G::GeoSpectral)

    @unpack boundary_path, boundary_file, gravity, NF = P

    # LOAD NETCDF FILE
    if boundary_path == ""
        path = joinpath(@__DIR__,"../input_data",boundary_file)
    else
        path = joinpath(boundary_path,boundary_file)
    end
    ncfile = NetCDF.open(path)

    # READ, latitude is North to South in file
    orography = ncfile.vars["orog"][:,end:-1:1]         # height [m]
    landsea_mask = ncfile.vars["lsm"][:,end:-1:1]       # fraction of land [0-1]
    albedo = ncfile.vars["alb"][:,end:-1:1]             # annual-mean albedo [0-1]

    # GEOPOTENTIAL
    # transform to spectral space and truncate
    P_inputdata = Parameters(trunc=31,nlat=48,nlon=96)  # orography only avail on default grid
    G_inputdata = GeoSpectral(P_inputdata)              # spectral transform for default grid
    geopot_surf_lowres = spectral(gravity*orography,G_inputdata)
    spectral_truncation!(geopot_surf_lowres,G_inputdata)

    # interpolate orography by adding 0 coefficients to high resolution spectral geopotential
    geopot_surf = zeros(Complex{NF},G.spectral.mx,G.spectral.nx)    # highres with mx,nx from G
    @unpack mx,nx = G_inputdata.spectral                            # lowres mx,nx matching input data
    geopot_surf[1:mx,1:nx] .= geopot_surf_lowres

    Boundaries{NF}(geopot_surf,landsea_mask,albedo)
end

ext: jl


text:
module ExpansionItems

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export expansionitem

register_normal_element("q__expansion__item", context = @__MODULE__)

function expansionitem(args...; kwargs...)
  q__expansion__item(args...; kw(kwargs)...)
end

end


text:
atoms = Molecules.parse_string("""
C   -2.131551124300    2.286168823700    0.000000000000
H   -1.061551124300    2.286168823700    0.000000000000
H   -2.488213906200    1.408104616400    0.496683911300
H   -2.488218762100    2.295059432700   -1.008766153900
H   -2.488220057000    3.155340844300    0.512081313000""")

s = BasisFunction(0, [1/√2, 1/√2], [5.0, 1.2])
p = BasisFunction(1, [1/√2, 1/√2], [5.0, 1.2])
shells = [ [s,p], [s], [s], [s], [s]]
bset = BasisSet("Test", atoms, shells)

@testset "Misc" begin
    @test GaussianBasis.string_repr(s) == "S shell with 1 basis built from 2 primitive gaussians\n\nχ₀₀  =    0.7071067812⋅Y₀₀⋅exp(-5.0⋅r²)\n     +    0.7071067812⋅Y₀₀⋅exp(-1.2⋅r²)"
    @test GaussianBasis.string_repr(bset) == "Test Basis Set\nNumber of shells: 6\nNumber of basis:  8\n\nC: 1s 1p \nH: 1s \nH: 1s \nH: 1s \nH: 1s"
    @test bset[1] == [s,p]
    @test bset[1,2] == p
end


text:
module ErrorLineNumber

## insert_line_numbers
## Take an expression inprogram,
##
## Produce an outprogram, which is the rewritten inprogram
## with statements to update line number variable.
##
## There are three versions of routine; the correct one is
## selected by multiple dispatch on the first argument.

insert_line_numbers(inprogram::Any,
                             ::ASCIIString,
                             ::Symbol) = inprogram

insert_line_numbers(inprogram::LineNumberNode,
                             fnname::ASCIIString,
                             specialsym::Symbol) =
    Expr(:block,
         Expr(:(=), specialsym, fnname * "." * string(inprogram.line)),
         inprogram)

function insert_line_numbers(inprogram::Expr,
                             fnname::ASCIIString,
                             specialsym::Symbol)
    if inprogram.head == :line
        # If the expression is an expression of type :line, i.e.,
        # a new line number in the source code,
        # replace it with a block that consists of the
        # the line number and a tracking array incrementer
        outprogram =  Expr(:block,
                           Expr(:(=), specialsym, fnname * "." *
                                string(inprogram.args[1])),
                           inprogram)
    else
        outprogram = Expr(inprogram.head)
        for arg in inprogram.args
            push!(outprogram.args, insert_line_numbers(arg, fnname, specialsym))
        end
    end
    return outprogram
end

macro errorlinenumber(fn)
    fn.head != :function && error("errorlinenumber macro only available for functions")
    @assert fn.args[1].head == :call
    @assert fn.args[2].head == :block
    fnname = fn.args[1].args[1]
    g = gensym()
    g2 = gensym()
    newfn = Expr(:function,
                 fn.args[1],
                 Expr(:block,
                      :($g = ""),
                      Expr(:try,
                           insert_line_numbers(fn.args[2], string(fnname), g),
                           g2,
                           Expr(:block,
                                :(println(STDERR, "    !!!!!!!!!!!!!!!!!")),
                                :(println(STDERR, "ERROR LINE number = ", $g)),
                                :(println(STDERR, "    !!!!!!!!!!!!!!!!!")),
                                :(rethrow($g2))))))
    esc(newfn)
end

end # module


text:
using Test
using Qecsim.PauliTools

@testset "bsp" begin
    # III bsp III commute (i.e. return false)
    @test bsp(BitVector([0, 0, 0, 0, 0, 0])', BitVector([0, 0, 0, 0, 0, 0]))::Bool == 0
    # XII bsp XII commute
    @test bsp(BitVector([1, 0, 0, 0, 0, 0])', BitVector([1, 0, 0, 0, 0, 0]))::Bool == 0
    # XII bsp IZZ commute
    @test bsp(BitVector([1, 0, 0, 0, 0, 0])', BitVector([0, 0, 0, 0, 1, 1]))::Bool == 0
    # XII bsp ZII do not commute (i.e. return true)
    @test bsp(BitVector([1, 0, 0, 0, 0, 0])', BitVector([0, 0, 0, 1, 0, 0]))::Bool == 1
    # XXI bsp ZZI commute
    @test bsp(BitVector([1, 1, 0, 0, 0, 0])', BitVector([1, 1, 0, 0, 0, 0]))::Bool == 0
    # XXX bsp ZZZ do not commute
    @test bsp(BitVector([1, 1, 1, 0, 0, 0])', BitVector([0, 0, 0, 1, 1, 1]))::Bool == 1

    # XII bsp IYY commute
    @test bsp(BitVector([1, 0, 0, 0, 0, 0])', BitVector([0, 1, 1, 0, 1, 1]))::Bool == 0
    # XII bsp YII do not commute
    @test bsp(BitVector([1, 0, 0, 0, 0, 0])', BitVector([1, 0, 0, 1, 0, 0]))::Bool == 1
    # XXI bsp YYI commute
    @test bsp(BitVector([1, 1, 0, 0, 0, 0])', BitVector([1, 1, 0, 1, 1, 0]))::Bool == 0
    # XXX bsp YYY do not commute
    @test bsp(BitVector([1, 1, 1, 0, 0, 0])', BitVector([1, 1, 1, 1, 1, 1]))::Bool == 1

    stabilizers = BitMatrix(  # 5-qubit stabilizers
        [1 0 0 1 0 0 1 1 0 0   # XZZXI
         0 1 0 0 1 0 0 1 1 0   # IXZZX
         1 0 1 0 0 0 0 0 1 1   # XIXZZ
         0 1 0 1 0 1 0 0 0 1]  # ZXIXZ
    )
    error = BitVector([0, 0, 1, 1, 0, 0, 1, 0, 1, 0])  # IZXYI
    commutation = BitVector([0, 1, 1, 0])
    @test bsp(stabilizers, error)::BitVector == commutation
    @test bsp(error', stabilizers')::BitMatrix == commutation'
    errors = hcat(
        BitVector([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),  # XIIII
        BitVector([0, 0, 0, 0, 0, 0, 0, 1, 0, 0]),  # IIZII
        BitVector([0, 0, 0, 0, 1, 0, 0, 0, 0, 1]),  # IIIIY
        BitVector([0, 0, 1, 1, 0, 0, 1, 0, 1, 0]),  # IZXYI
    )
    commutation = BitMatrix([[0, 0, 0, 1] [0, 0, 1, 0] [0, 1, 1, 1] [0, 1, 1, 0]])
    @test bsp(stabilizers, errors)::BitMatrix == commutation
    @test bsp(errors', stabilizers')::BitMatrix == commutation'
end

@testset "pack_unpack" begin
    # specific test
    bsf = to_bsf("XIZIY")  # 1000100101
    packed_bsf = pack(bsf)
    @test packed_bsf == ("225", 10)
    @test unpack(packed_bsf) == bsf
    # long (>128) bit-vector (to test BigInt handles long binary vectors)
    bsf = BitVector(rand(0:1, 150))
    @test unpack(pack(bsf)) == bsf
end

@testset "to_bsf" begin
    # Single Paulis
    @test to_bsf("XIZIY") == BitVector([1, 0, 0, 0, 1, 0, 0, 1, 0, 1])
    @test to_bsf("IIIII") == BitVector([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
    @test to_bsf("XXXXX") == BitVector([1, 1, 1, 1, 1, 0, 0, 0, 0, 0])
    @test to_bsf("ZZZZZ") == BitVector([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
    @test to_bsf("YYYYY") == BitVector([1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
    # Multiple Paulis
    @test to_bsf(["XIZIY", "IXZYI"]) == BitMatrix(
        [1 0 0 0 1 0 0 1 0 1
         0 1 0 1 0 0 0 1 1 0]
    )
    @test to_bsf(["XXXXX", "ZZZZZ", "YYYYY"]) == BitMatrix(
        [1 1 1 1 1 0 0 0 0 0
         0 0 0 0 0 1 1 1 1 1
         1 1 1 1 1 1 1 1 1 1]
    )
end

@testset "to_pauli" begin
    # Single bsf
    @test to_pauli(BitVector([1, 0, 0, 0, 1, 0, 0, 1, 0, 1])) == "XIZIY"
    @test to_pauli(BitVector([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])) == "IIIII"
    @test to_pauli(BitVector([1, 1, 1, 1, 1, 0, 0, 0, 0, 0])) == "XXXXX"
    @test to_pauli(BitVector([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])) == "ZZZZZ"
    @test to_pauli(BitVector([1, 1, 1, 1, 1, 1, 1, 1, 1, 1])) == "YYYYY"
    # Multiple bsfs
    @test to_pauli(BitMatrix(
        [1 0 0 0 1 0 0 1 0 1
         0 1 0 1 0 0 0 1 1 0])
    ) == ["XIZIY", "IXZYI"]
    @test to_pauli(BitMatrix(
        [1 1 1 1 1 0 0 0 0 0
         0 0 0 0 0 1 1 1 1 1
         1 1 1 1 1 1 1 1 1 1])
    ) == ["XXXXX", "ZZZZZ", "YYYYY"]
    end

@testset "weight" begin
    for (pauli, wt) in (
        ("IIIII", 0),
        ("XIZIY", 3),
        ("XXXXX", 5),
        ("YYYYY", 5),
        ("ZZZZZ", 5),
        (["XIZIY", "XXXXX"], 8),
    )
        @test weight(pauli) == wt
        @test weight(to_bsf(pauli)) == wt
    end
end


text:
using BenchmarkTools, Test, CUDA

a = CUDA.zeros(1024)

function kernel(a)
    i = threadIdx().x
    a[i] += 1
    return
end

@cuda threads=length(a) kernel(a)

##

N = 2^20
x_d = CUDA.fill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = CUDA.fill(2.0f0, N)  # a vector stored on the GPU filled with 2.0
y_d .+= x_d

function add_broadcast!(y, x)
    CUDA.@sync y .+= x
    return
end

##
@btime add_broadcast!($y_d, $x_d)

##

function gpu_add1!(y, x)
    for i = 1:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y_d, 2)
@cuda gpu_add1!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

##

function bench_gpu1!(y, x)
    CUDA.@sync begin
        @cuda gpu_add1!(y, x)
    end
end
@btime bench_gpu1!($y_d, $x_d)

##

const nx = 1024  # do 1024 x 1024 2D FFT
xc = CuArray{ComplexF64}(CUDA.randn(Float64, nx, nx))
p = plan_fft!( xc )

##
@btime CUDA.@sync(p * x) setup=(
    x=CuArray{ComplexF64}(CUDA.randn(Float64, nx, nx)));

##
for device in CUDA.devices()
    @show capability(device)
end

##

using AbstractFFTs
using CUDA.CUFFT
##
b = CUDA.rand(ComplexF32,64,64,64)
# pa = plan_fft( a )
@btime fft(b);


text:
### Operations and Strategies
include("operations/models.jl")
include("operations/options.jl")
include("operations/portfolios.jl")


text:
cd(@__DIR__)
using Pkg
CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
CI && Pkg.activate(@__DIR__)
CI && Pkg.instantiate()
CI && (ENV["GKSwstype"] = "100")
using Documenter
using DocumenterTools
using DocumenterMarkdown

using UncertainData
using Distributions
using KernelDensity
using StatsBase
using Measurements
using Interpolations

PAGES = [
    "index.md",
    "Uncertain values" => [
        "uncertain_values/uncertainvalues_overview.md",
        "uncertain_values/uncertainvalues_theoreticaldistributions.md",
        "uncertain_values/uncertainvalues_kde.md",
        "uncertain_values/uncertainvalues_fitted.md",
        "uncertain_values/uncertainvalues_certainvalue.md",
        "uncertain_values/uncertainvalues_populations.md",
        "uncertain_values/uncertainvalues_Measurements.md",
        "uncertain_values/merging.md",
        "uncertain_values/uncertainvalues_examples.md",
    ],
	"Uncertain datasets" => [
        "uncertain_datasets/uncertain_datasets_overview.md",
        "uncertain_datasets/uncertain_index_dataset.md",
        "uncertain_datasets/uncertain_value_dataset.md",
        "uncertain_datasets/uncertain_indexvalue_dataset.md",
        "uncertain_datasets/uncertain_dataset.md",
	],
    "Uncertain statistics" => [
        "Core statistics" => [
            "uncertain_statistics/core_stats/core_statistics.md",
            "uncertain_statistics/core_stats/core_statistics_point_estimates.md",
            "uncertain_statistics/core_stats/core_statistics_pairwise_estimates.md",
            "uncertain_statistics/core_stats/core_statistics_datasets_single_dataset_estimates.md",
            "uncertain_statistics/core_stats/core_statistics_datasets_pairwise_estimates.md",
            "uncertain_statistics/core_stats/core_statistics_datasets.md"
        ],

        "Hypothesis tests" => [
			"uncertain_statistics/hypothesistests/hypothesis_tests_overview.md",
            "uncertain_statistics/hypothesistests/one_sample_t_test.md",
            "uncertain_statistics/hypothesistests/equal_variance_t_test.md",
            "uncertain_statistics/hypothesistests/unequal_variance_t_test.md",
            "uncertain_statistics/hypothesistests/exact_kolmogorov_smirnov_test.md",
            "uncertain_statistics/hypothesistests/approximate_twosample_kolmogorov_smirnov_test.md",
            "uncertain_statistics/hypothesistests/jarque_bera_test.md",
            "uncertain_statistics/hypothesistests/mann_whitney_u_test.md",
            "uncertain_statistics/hypothesistests/anderson_darling_test.md"
        ],
    ],
    "Sampling constraints" => [
        "sampling_constraints/available_constraints.md",
        "sampling_constraints/constrain_uncertain_values.md",
        "sampling_constraints/sequential_constraints.md"
    ],

    "Binning" => [
        "binning/bin.md"
    ],
    "Resampling" => [
        "resampling/resampling_overview.md",
        "resampling/resampling_uncertain_values.md",
        "resampling/resampling_uncertain_datasets.md",
        "resampling/resampling_uncertain_indexvalue_datasets.md",

        "resampling/sequential/resampling_uncertaindatasets_sequential.md",
        "resampling/sequential/resampling_indexvalue_sequential.md",
        "resampling/sequential/strictly_increasing.md",
        "resampling/sequential/strictly_decreasing.md",

        "resampling/interpolation/interpolation.md",
        "resampling/interpolation/gridded.md",
        "resampling/resampling_schemes/resampling_schemes_uncertain_value_collections.md",
        "resampling/resampling_schemes/resampling_schemes_uncertain_indexvalue_collections.md",
        "resampling/resampling_schemes/resampling_with_schemes_uncertain_value_collections.md",
        "resampling/resampling_schemes/resampling_with_schemes_uncertain_indexvalue_collections.md",
        "resampling/resampling_inplace.md"
        #"resampling/models/resampling_with_models.md"
    ],

    "Propagation of errors" => [
        "propagation_of_errors/propagation_of_errors.md"
    ],

    "Mathematics" => [
        "mathematics/elementary_operations.md",
        "mathematics/trig_functions.md"
    ],

    "Tutorials" => [
        "tutorials/tutorial_overview.md",
        "tutorials/tutorial_transforming_data_to_regular_grid.md"
    ],

    "implementing_algorithms_for_uncertaindata.md",

    "changelog.md",

    "publications.md",

    "citing.md"
]

# %% Build docs
#PyPlot.ioff()
cd(@__DIR__)
ENV["JULIA_DEBUG"] = "Documenter"

makedocs(
    modules = [UncertainData],
    sitename = "UncertainData.jl documentation",
    format = format = Documenter.HTML(
        prettyurls = CI,
        ),
    pages = PAGES
)

if CI
    deploydocs(
        repo = "github.com/kahaaga/UncertainData.jl.git",
        target = "build",
        push_preview = true
    )
end


text:
using Plots

Plots.test_examples(:gr, skip = Plots._backend_skips[:gr])
Plots.test_examples(:plotly, skip = Plots._backend_skips[:plotly], disp = false)


text:
using Mimi, MimiSSPs, DataFrames, CSVFiles, Query, Test

all_countries = load(joinpath(@__DIR__, "..", "data", "keys", "MimiSSPs_ISO.csv")) |> DataFrame

# BASIC API

m = Model()
set_dimension!(m, :time, 1750:2300)
set_dimension!(m, :countries, all_countries.ISO)
add_comp!(m, MimiSSPs.SSPs, first = 2010, last = 2300)
update_param!(m, :SSPs, :model, "IIASA GDP")
update_param!(m, :SSPs, :ssp, "SSP1")
update_param!(m, :SSPs, :rcp, "RCP1.9")
update_param!(m, :SSPs, :country_names, all_countries.ISO)

run(m)

m = Model()
set_dimension!(m, :time, 1750:2300)
set_dimension!(m, :countries, all_countries.ISO[1:10])
add_comp!(m, MimiSSPs.SSPs, first = 2010, last = 2300)
update_param!(m, :SSPs, :model, "IIASA GDP")
update_param!(m, :SSPs, :ssp, "SSP5")
update_param!(m, :SSPs, :rcp, "RCP8.5")
update_param!(m, :SSPs, :country_names, all_countries.ISO[1:10])

run(m)

# ERRORS

m = Model()
set_dimension!(m, :time, 1750:2300)
set_dimension!(m, :countries, all_countries.ISO)
add_comp!(m, MimiSSPs.SSPs)
update_param!(m, :SSPs, :model, "IIASA GDP")
update_param!(m, :SSPs, :ssp, "SSP1")
update_param!(m, :SSPs, :rcp, "RCP1.9")
update_param!(m, :SSPs, :country_names, all_countries.ISO)

error_msg = (try eval(run(m)) catch err err end).msg
@test occursin("Cannot run SSP component in year 1750", error_msg)

dummy_countries = ["Sun", "Rain", "Cloud"]

m = Model()
set_dimension!(m, :time, 1750:2300)
set_dimension!(m, :countries, dummy_countries)
add_comp!(m, MimiSSPs.SSPs)
update_param!(m, :SSPs, :model, "IIASA GDP")
update_param!(m, :SSPs, :ssp, "SSP1")
update_param!(m, :SSPs, :rcp, "RCP1.9")
update_param!(m, :SSPs, :country_names, dummy_countries) # error because countries aren't in SSP set

error_msg = (try eval(run(m)) catch err err end).msg
@test occursin("All countries in countries parameter must be found in SSPs component Socioeconomic Dataframe, the following were not found:", error_msg)

m = Model()
set_dimension!(m, :time, 1750:2300)
set_dimension!(m, :countries, all_countries.ISO)
add_comp!(m, MimiSSPs.SSPs, first = 2010, last = 2300)
update_param!(m, :SSPs, :model, "NOT A MODEL")
update_param!(m, :SSPs, :ssp, "SSP1")
update_param!(m, :SSPs, :rcp, "RCP1.9")
update_param!(m, :SSPs, :country_names, all_countries.ISO)

error_msg = (try eval(run(m)) catch err err end).msg
@test occursin("Model NOT A MODEL provided to SSPs component model parameter not found in available list:", error_msg)

update_param!(m, :SSPs, :model, "IIASA GDP")
update_param!(m, :SSPs, :ssp, "NOT A SSP")

error_msg = (try eval(run(m)) catch err err end).msg
@test occursin("SSP NOT A SSP provided to SSPs component ssp parameter not found in available list:", error_msg)

update_param!(m, :SSPs, :ssp, "SSP1")
update_param!(m, :SSPs, :rcp, "NOT A RCP")

error_msg = (try eval(run(m)) catch err err end).msg
@test occursin("RCP NOT A RCP provided to SSPs component ssp parameter not found in available list:", error_msg)

# DATA OUTPUT CHECK


text:
#! /usr/bin/env julia

import Pkg

(pwd() == @__DIR__) || cd(@__DIR__)

Pkg.activate(pwd())
Pkg.instantiate()

using Air

jobs = Symbol[]
if length(ARGS) == 0
    push!(jobs, :test)
    push!(jobs, :docs)
else
    for a in ARGS
        a = lowercase(a)
        push!(jobs, Symbol(a))
    end
end
for j in jobs
    if j == :test
        Pkg.test("Air", coverage=true)
    elseif j == :doctest
        Pkg.add("Documenter")
        import Documenter
        Documenter.doctest(Air)
    elseif j == :codecov
        Pkg.add("Coverage")
        using Coverage
        Codecov.submit(Codecov.process_folder())
    elseif j == :docs
        Pkg.add("Documenter")
        include("docs/make.jl")
    else
        error("Unknown make command: $j")
    end
end

exit(0)


text:
# To run benchmarks locally, BenchmarkCI should be added to root project.
# Then call:
# ```bash
# julia benchmark/run_benchmarks.jl
# ```
using BenchmarkCI
on_CI = haskey(ENV, "GITHUB_ACTIONS")

BenchmarkCI.judge()
on_CI ? BenchmarkCI.postjudge() : BenchmarkCI.displayjudgement()


text:
"""
    certifyindices(domain::SphericalTriangle{T},
                   λs::AbstractVector{arb}) where {T}

Attempt to determine the indices for the eigenvalues `λs` of the domain.
Returns a list consisting of tuples `(λ, minindex, maxindex)` where
`minindex` and `maxindex` are the minimum and maximum possible indices
for the eigenvalue λ. In particular if `minindex = maxindex` then that
is the exact index of λ.

The `minindex` is determined simply by determining how many
eigenvalues in `λs` are smaller then `λ`. For the upper bound it
computes eigenvalues of enclosing domains and uses those to lower
bound the eigenvalues of the domain.
"""
function certifyindices(
    domain::SphericalTriangle{T},
    λs::AbstractVector{arb};
    show_trace = false,
) where {T}
    indices = certifyindices(domain, λs, 1, show_trace = show_trace)

    if all(((λ, minindex, maxindex),) -> minindex == maxindex, indices)
        return indices
    end

    for vertex = 2:3
        indices2 = certifyindices(domain, λs, vertex, show_trace = show_trace)

        indices = [
            (λs[i], max(indices[i][2], indices2[i][2]), min(indices[i][3], indices2[i][3])) for i = 1:length(indices)
        ]

        if all(((λ, minindex, maxindex),) -> minindex == maxindex, indices)
            return indices
        end
    end

    return indices
end

function certifyindices(
    domain::SphericalTriangle{T},
    λs::AbstractVector{arb},
    vertex;
    show_trace = false,
) where {T}
    domain = SphericalTriangle(
        domain.angles[[mod1(i, 3) for i = vertex:vertex+2]],
        domain.parent,
    )

    # Define interval to work on and the parameters to use
    start = domain.parent(0)
    stop = -0.5 + sqrt(0.25 + getinterval(λs[end])[2] + 1)

    z = cos(theta_bound(domain, :upper))

    # Find all eigenvalues of the enclosing spherical cap sector in
    # the above interval
    eigenvalues = arb[]
    eigenvaluesflags = Int[]
    k = 1

    # Temporary for plotting
    νs = range(Float64(start), stop = Float64(stop), length = 2000)
    #λs = νs .* (νs .+ 1)
    res = []

    while true
        if T == fmpq
            μ = domain.parent(-k * inv(domain.angles[1]))
        else
            μ = -k * inv(domain.angles[1]) * domain.parent(π)
        end

        f = ν -> legendre_p(ν, μ, z)
        roots, flags = isolateroots(
            f,
            start,
            stop,
            evaltype = :taylor,
            maxevals = 5000,
            show_trace = false,
        )

        if !all(flags)
            @error "unable to isolate the roots of the legendre function"
        end

        push!(res, Float64.(f.(domain.parent.(νs))))

        if show_trace
            @show k
            @show length(roots)
        end

        if isempty(roots)
            break
        end

        append!(
            eigenvalues,
            [setinterval(l, u) * (setinterval(l, u) + 1) for (l, u) in roots],
        )
        append!(eigenvaluesflags, flags)
        k += 1
    end

    indices = Tuple{arb,Int,Int}[]
    for λ in λs
        minindex = length(filter(eig -> eig < λ, λs)) + 1

        # Find all eigenvalues of the enclosing domain that at least
        # as small as λ. Ensure that all of them have flag = 1
        idx = findall(eig -> !(eig > λ), eigenvalues)
        if !all(==(1), eigenvaluesflags[idx])
            @error "unable to upper bound the index for λ = $λ"
            maxindex = typemax(Int)
        else
            maxindex = length(idx)
        end
        push!(indices, (λ, minindex, maxindex))
    end

    return indices
    return [(λ, Int[length(filter(eig -> !(eig > λ), eigenvalues))]) for λ in λs]

    eigenvalues, λs, res
    # Compare with the given eigenvalues
end


text:
@testset "ColorJitter" begin
    @test ColorJitter <: Augmentor.ImageOperation

    @testset "constructor" begin
        @test_throws ArgumentError ColorJitter(1., 3.0:2.0)
        @test_throws ArgumentError ColorJitter(5.0:3.0, 1.)
    end
    @testset "randparam" begin
        img = testpattern()
        αs = [3.0, 1.0:0.2:5.0, [1.0, 3.0, 5.0]]
        βs = [1.0, 1.0:0.1:2.0, [1.0, 2.0]]
        for α in αs, β in βs
            op = ColorJitter(α, β)
            p = @inferred Augmentor.randparam(op, img)
            @test p[1] in α
            @test p[2] in β
        end
    end
    @testset "_get_M" begin
        imgs = [testpattern()]#, camera]
        for img in imgs
            T = eltype(img)
            op = ColorJitter(1., 0., usemax=true)
            M = @inferred Augmentor._get_M(op, img)
            @test M == T(gamutmax(T)...)
            op = ColorJitter(1., 0., usemax=false)
            M = @inferred Augmentor._get_M(op, img)
            @test M == convert(T, mean(img))
        end
    end
    @testset "eager" begin
        @test Augmentor.supports_eager(ColorJitter)

        α = 0.2
        β = 0.1

        imgs = [testpattern(), camera]

        # M = mean value
        for img in imgs
            img = convert.(color_type(eltype(img)), img)
            ref = clamp01.(α .* img .+ β * mean(img))
            res = Augmentor.applyeager(ColorJitter(α, β, usemax=false), img)
            @test assess_psnr(res, ref) > 25
        end
        # M = maximum value
        for img in imgs
            img = convert.(color_type(eltype(img)), img)
            T = eltype(img)
            ref = clamp01.(α .* img .+ β * T(gamutmax(T)...))
            res = Augmentor.applyeager(ColorJitter(α, β, usemax=true), img)
            @test assess_psnr(res, ref) > 25
        end
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(ColorJitter)

        α = 0.2
        β = 0.1

        imgs = [testpattern(), camera]

        # M = mean value
        for img in imgs
            img = convert.(color_type(eltype(img)), img)
            ref = clamp01.(α .* img .+ β * mean(img))
            res = Augmentor.applylazy(ColorJitter(α, β, usemax=false), img)
            @test assess_psnr(res, ref) > 25
        end
        # M = maximum value
        for img in imgs
            img = convert.(color_type(eltype(img)), img)
            T = eltype(img)
            ref = clamp01.(α .* img .+ β * T(gamutmax(T)...))
            res = Augmentor.applylazy(ColorJitter(α, β, usemax=true), img)
            @test assess_psnr(res, ref) > 25
        end
    end
end


text:
compute_moments(basis::Dictionary) = [moment(basis,i) for i in eachindex(basis)]

function apply_quad(w, x, f)
    z = w[1] * f(x[1])
    for i in 2:length(w)
        z += w[i] * f(x[i])
    end
    z
end

compute_weights(x, basis, B) = interpolation_matrix(basis, x)' \ B

"""
A container to collect information about a quadrature rule, mainly for use
in types describing exactness equations.
"""
struct QuadRuleData{B <: Dictionary,T}
    "The basis for which the rule is exact."
    basis   ::  B
    "The length of the quadrature rule"
    len     ::  Int
    "Moments of the basis with respect to some measure."
    moments ::  Vector{T}
end

eltype(::QuadRuleData{B,T}) where {B,T} = T

function QuadRuleData(basis::Dictionary{S,T}, len::Int, moments = compute_moments(basis)) where {S,T}
    QuadRuleData{typeof(basis),T}(basis, len, moments)
end

length(d::QuadRuleData) = d.len

basis(d::QuadRuleData) = d.basis

moments(d::QuadRuleData) = d.moments

"A type that describes the nonlinear equations of exactness of a quadrature rule, for use in Newton's method."
abstract type QuadRuleSystem{T} <: NonlinearSystem end

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

data(sys::QuadRuleSystem) = sys.data

for op in (:basis, :moments)
    @eval $op(sys::QuadRuleSystem) = $op(sys.data)
end

length(rule::QuadRuleSystem) = dofs(rule)

size(rule::QuadRuleSystem) = (length(rule), dofs(rule))

quadlength(rule::QuadRuleSystem) = length(rule.data)

# Evaluate the residual of the system of equations
function residual!(result, sys::QuadRuleSystem, newton_x)
    w,x = newton_to_quad(sys, newton_x)
    residual!(result, sys, w, x, basis(sys), moments(sys))
end

function residual!(result, sys::QuadRuleSystem, w, x, basis, moments)
    for (i,ϕ) in enumerate(basis)
        result[i] =  apply_quad(w, x, ϕ) - moments[i]
    end
    result
end

# Evaluate the Jacobian of the system of equations
function jacobian!(J, sys::QuadRuleSystem, newton_x)
    w,x = newton_to_quad(sys, newton_x)
    jacobian!(J, sys, w, x, basis(sys))
end

function newton_to_quad(sys::QuadRuleSystem, newton_x)
    T = eltype(newton_x)
    x = zeros(T, length(sys.data))
    w = zeros(T, length(sys.data))
    newton_to_quad!(sys, w, x, newton_x)
end

function quad_to_newton(sys::QuadRuleSystem, w, x)
    T = promote_type(eltype(w),eltype(x))
    newton_x = zeros(T, dofs(sys))
    quad_to_newton!(sys, w, x, newton_x)
end

"System of equations for a quadrature rule where all points are free."
struct QuadRuleFreePoints{B,T} <: QuadRuleSystem{T}
    data    ::  QuadRuleData{B,T}
end

"The total number of degrees of freedom in the system."
dofs(rule::QuadRuleFreePoints) = 2*length(rule.data)

"System of equations for a quadrature rule with some points fixed."
struct QuadRuleFixedPoints{B,T} <: QuadRuleSystem{T}
    "Data about the quadrature rule: basis, moments, etcetera."
    data        ::  QuadRuleData{B,T}
    "Indices of the fixed points."
    fixed_idxs  ::  Vector{Int}
    "Values of the fixed points."
    fixed_pts   ::  Vector{T}
    "Indices of the free points."
    free_idxs   ::  Vector{Int}

    function QuadRuleFixedPoints{B,T}(data, fixed_idxs, fixed_pts) where {B,T}
        P = length(fixed_idxs)
        @assert length(fixed_pts) == P
        @assert P > 0
        # Make sure indices and points are sorted.
        for i in 1:P-1
            @assert fixed_idxs[i] < fixed_idxs[i+1]
            @assert fixed_pts[i] < fixed_pts[i+1]
        end
        new(data, fixed_idxs, fixed_pts, setdiff(1:length(data), fixed_idxs))
    end
end

QuadRuleFixedPoints(data::QuadRuleData{B,T}, fixed_idxs, fixed_pts) where {B,T} =
    QuadRuleFixedPoints{B,T}(data, fixed_idxs, fixed_pts)

"The number of fixed points in the quadrature rule."
nbfixed(rule::QuadRuleFixedPoints) = length(rule.fixed_pts)

"The number of free points in the quadrature rule."
nbfree(rule::QuadRuleFixedPoints) = length(rule.free_idxs)

dofs(rule::QuadRuleFixedPoints) = 2*length(rule.data) - nbfixed(rule)

function setfixedpoint!(rule::QuadRuleFixedPoints, idx, xstar)
    rule.fixed_pts[rule.fixed_idxs[idx]] = xstar
end

# newton_x contains l weights, followed by l points
function newton_to_quad!(sys::QuadRuleFreePoints, w, x, newton_x)
    @assert length(newton_x) == dofs(sys)
    l = length(w)
    for i in 1:l
        w[i] = newton_x[i]
        x[i] = newton_x[l+i]
    end
    w,x
end

function quad_to_newton!(sys::QuadRuleFreePoints, w, x, newton_x)
    @assert length(newton_x) == dofs(sys)
    l = length(w)
    for i in 1:l
        newton_x[i] = w[i]
        newton_x[l+i] = x[i]
    end
    newton_x
end

# newton_x contains l weights, followed by l-nbfixed(rule) free points (the other points are fixed)
function newton_to_quad!(sys::QuadRuleFixedPoints, w, x, newton_x)
    @assert length(newton_x) == dofs(sys)

    l = length(w)
    P = nbfixed(sys)

    for i in 1:l
        w[i] = newton_x[i]
    end

    # Copy the free points into x
    for (j,i) in enumerate(sys.free_idxs)
        x[i] = newton_x[l+j]
    end
    # and copy the fixed points into x
    for (i,xi) in zip(sys.fixed_idxs, sys.fixed_pts)
        x[i] = xi
    end
    w,x
end

function quad_to_newton!(sys::QuadRuleFixedPoints, w, x, newton_x)
    @assert length(newton_x) == dofs(sys)

    l = length(w)
    P = nbfixed(sys)

    for i in 1:l
        newton_x[i] = w[i]
    end

    # Copy the free points
    for (j,i) in enumerate(sys.free_idxs)
        newton_x[l+j] = x[i]
    end
    newton_x
end

function jacobian!(J, sys::QuadRuleFreePoints, w, x, basis)
    l = length(w)

    # We are computing the derivative of \sum_{i=1}^l w_i ϕ_j(x_i) wrt w_i and x_i
    # - First the derivatives with respect to the weight w[i]
    for j in 1:length(basis)
        for i in 1:l
            J[j,i] = BasisFunctions.unsafe_eval_element(basis, j, x[i])
        end
    end
    # - Then the weight times the derivatives of the basis functions
    for j in 1:length(basis)
        for i in 1:l
            J[j,l+i] = w[i] * BasisFunctions.unsafe_eval_element_derivative(basis, j, x[i], 1)
        end
    end
    J
end

function jacobian!(J, sys::QuadRuleFixedPoints, w, x, basis)
    l = length(w)
    P = nbfixed(sys)

    # We are computing the derivative of \sum_{i=1}^l w_i ϕ_j(x_i) wrt w_i and x_i
    # as before, but now some of the x_i's are fixed and we have to skip them.
    # - First the derivatives with respect to the weight w[i]: same as before
    for j in 1:length(basis)
        for i in 1:l
            J[j,i] = BasisFunctions.unsafe_eval_element(basis, j, x[i])
        end
    end
    # - Then the weight times the derivatives of the basis functions
    for j in 1:length(basis)
        for (k,i) in enumerate(sys.free_idxs)
            J[j,l+k] = w[i] * BasisFunctions.unsafe_eval_element_derivative(basis, j, x[i], 1)
        end
    end
    J
end


text:
include("hydrogen_fem_julia.jl")
include("hydrogen_vh_julia.jl")
using Printf
using .Hydrogen_FEM
using .Hydrogen_Vh
#using Plots

function main()
    hfem_param, hfem_val = Hydrogen_FEM.construct()

    Hydrogen_FEM.make_wavefunction(hfem_param, hfem_val)
    vh_param, vh_val = Hydrogen_Vh.construct(hfem_param)
    Hydrogen_Vh.do_run(hfem_param, hfem_val, vh_val)

    @printf "計算が終わりました\n"
    Hydrogen_Vh.save_result(hfem_val, vh_param, vh_val)
    @printf "計算結果を%sに書き込みました\n" vh_param.RESULT_FILENAME
    #plot(val.node_r_glo, val.phi)
end
@time main()


text:
"""
    DataFrameRow{<:AbstractDataFrame}

A view of one row of an AbstractDataFrame.
"""
struct DataFrameRow{T<:AbstractDataFrame}
    df::T
    row::Int

    @inline function DataFrameRow{T}(df::T, row::Integer) where T<:AbstractDataFrame
        if row isa Bool
            throw(ArgumentError("invalid index: $row of type Bool"))
        end
        @boundscheck if !checkindex(Bool, axes(df, 1), row)
            throw(BoundsError("attempt to access a data frame with $(nrow(df)) " *
                              "rows at index $row"))
        end
        new{T}(df, row)
    end
end

@inline DataFrameRow(df::T, row::Integer) where {T<:AbstractDataFrame} =
    DataFrameRow{T}(df, row)

Base.parent(r::DataFrameRow) = getfield(r, :df)
Base.parentindices(r::DataFrameRow) = (row(r), axes(parent(r), 2))
row(r::DataFrameRow) = getfield(r, :row)

Base.view(adf::AbstractDataFrame, rowind::Integer, ::Colon) =
    DataFrameRow(adf, rowind)
Base.view(sdf::SubDataFrame, rowind::Integer, ::Colon) =
    DataFrameRow(parent(sdf), rows(sdf)[rowind])

# Here a corner case is when colinds=[] and we pass a valid rowind
# into adf. We will throw an error in this case.
# The consequence is that it is impossible to create a DataFrameRow without columns.
Base.view(adf::AbstractDataFrame, rowind::Integer, colinds::AbstractVector) =
    DataFrameRow(adf[colinds], rowind)
Base.view(sdf::SubDataFrame, rowind::Integer, colinds::AbstractVector) =
    DataFrameRow(parent(sdf)[colinds], rows(sdf)[rowind])

# Same here. It is impossible to create a DataFrameRow without columns.
Base.getindex(df::DataFrame, rowind::Integer, colinds::AbstractVector) =
    DataFrameRow(df[colinds], rowind)
Base.getindex(sdf::SubDataFrame, rowind::Integer, colinds::AbstractVector) =
    DataFrameRow(parent(sdf)[colinds], rows(sdf)[rowind])
Base.getindex(df::DataFrame, rowind::Integer, ::Colon) =
    DataFrameRow(df, rowind)
Base.getindex(sdf::SubDataFrame, rowind::Integer, ::Colon) =
    DataFrameRow(parent(sdf), rows(sdf)[rowind])

Base.getindex(r::DataFrameRow, idx::ColumnIndex) = parent(r)[row(r), idx]
Base.getindex(r::DataFrameRow, idxs::AbstractVector) =
    DataFrameRow(parent(r)[idxs], row(r))
Base.getindex(r::DataFrameRow, ::Colon) = r

Base.setindex!(r::DataFrameRow, value::Any, idx::Any) =
    setindex!(parent(r), value, row(r), idx)

Base.names(r::DataFrameRow) = names(parent(r))
_names(r::DataFrameRow) = _names(parent(r))

Base.haskey(r::DataFrameRow, key::Any) = haskey(parent(r), key)

Base.getproperty(r::DataFrameRow, idx::Symbol) = getindex(r, idx)
Base.setproperty!(r::DataFrameRow, idx::Symbol, x::Any) = setindex!(r, x, idx)
# Private fields are never exposed since they can conflict with column names
Base.propertynames(r::DataFrameRow, private::Bool=false) = names(r)

Base.view(r::DataFrameRow, col::ColumnIndex) = view(parent(r)[col], row(r))
Base.view(r::DataFrameRow, cols::AbstractVector) =
    DataFrameRow(parent(r)[cols], row(r))
Base.view(r::DataFrameRow, ::Colon) = r

index(r::DataFrameRow) = index(parent(r))
Base.size(r::DataFrameRow) = (size(parent(r), 2),)
Base.size(r::DataFrameRow, i) = size(r)[i]
Base.length(r::DataFrameRow) = size(parent(r), 2)
Base.ndims(r::DataFrameRow) = 1
Base.ndims(::Type{<:DataFrameRow}) = 1

Base.lastindex(r::DataFrameRow) = size(parent(r), 2)

Base.iterate(r::DataFrameRow) = iterate(r, 1)

function Base.iterate(r::DataFrameRow, st)
    st > length(r) && return nothing
    return (r[st], st + 1)
end

# Computing the element type requires going over all columns,
# so better let collect() do it only if necessary (widening)
Base.IteratorEltype(::DataFrameRow) = Base.EltypeUnknown()

function Base.convert(::Type{Vector}, dfr::DataFrameRow)
    T = reduce(promote_type, eltypes(parent(dfr)))
    convert(Vector{T}, dfr)
end
Base.convert(::Type{Vector{T}}, dfr::DataFrameRow) where T =
    T[dfr[i] for i in 1:length(dfr)]
Base.Vector(dfr::DataFrameRow) = convert(Vector, dfr)
Base.Vector{T}(dfr::DataFrameRow) where T = convert(Vector{T}, dfr)

Base.keys(r::DataFrameRow) = names(parent(r))
Base.values(r::DataFrameRow) = ntuple(col -> parent(r)[col][row(r)], length(r))

"""
    copy(dfr::DataFrameRow)

Convert a `DataFrameRow` to a `NamedTuple`.
"""
Base.copy(r::DataFrameRow) = NamedTuple{Tuple(keys(r))}(values(r))

# hash column element
Base.@propagate_inbounds hash_colel(v::AbstractArray, i, h::UInt = zero(UInt)) = hash(v[i], h)
Base.@propagate_inbounds function hash_colel(v::AbstractCategoricalArray, i, h::UInt = zero(UInt))
    ref = v.refs[i]
    if eltype(v) >: Missing && ref == 0
        hash(missing, h)
    else
        hash(CategoricalArrays.index(v.pool)[ref], h)
    end
end

# hash of DataFrame rows based on its values
# so that duplicate rows would have the same hash
# table columns are passed as a tuple of vectors to ensure type specialization
rowhash(cols::Tuple{AbstractVector}, r::Int, h::UInt = zero(UInt))::UInt =
    hash_colel(cols[1], r, h)
function rowhash(cols::Tuple{Vararg{AbstractVector}}, r::Int, h::UInt = zero(UInt))::UInt
    h = hash_colel(cols[1], r, h)
    rowhash(Base.tail(cols), r, h)
end

Base.hash(r::DataFrameRow, h::UInt = zero(UInt)) =
    rowhash(ntuple(i -> parent(r)[i], ncol(parent(r))), row(r), h)

# comparison of DataFrame rows
# only the rows of the same DataFrame could be compared
# rows are equal if they have the same values (while the row indices could differ)
# if all non-missing values are equal, but there are missings, returns missing
Base.:(==)(r1::DataFrameRow, r2::DataFrameRow) = isequal(r1, r2)

function Base.isequal(r1::DataFrameRow, r2::DataFrameRow)
    isequal_row(parent(r1), row(r1), parent(r2), row(r2))
end

# internal method for comparing the elements of the same data table column
isequal_colel(col::AbstractArray, r1::Int, r2::Int) =
    (r1 == r2) || isequal(Base.unsafe_getindex(col, r1), Base.unsafe_getindex(col, r2))

# table columns are passed as a tuple of vectors to ensure type specialization
isequal_row(cols::Tuple{AbstractVector}, r1::Int, r2::Int) =
    isequal(cols[1][r1], cols[1][r2])
isequal_row(cols::Tuple{Vararg{AbstractVector}}, r1::Int, r2::Int) =
    isequal(cols[1][r1], cols[1][r2]) && isequal_row(Base.tail(cols), r1, r2)

isequal_row(cols1::Tuple{AbstractVector}, r1::Int, cols2::Tuple{AbstractVector}, r2::Int) =
    isequal(cols1[1][r1], cols2[1][r2])
isequal_row(cols1::Tuple{Vararg{AbstractVector}}, r1::Int,
            cols2::Tuple{Vararg{AbstractVector}}, r2::Int) =
    isequal(cols1[1][r1], cols2[1][r2]) &&
        isequal_row(Base.tail(cols1), r1, Base.tail(cols2), r2)

function isequal_row(df1::AbstractDataFrame, r1::Int, df2::AbstractDataFrame, r2::Int)
    if df1 === df2
        if r1 == r2
            return true
        end
    elseif !(ncol(df1) == ncol(df2))
        throw(ArgumentError("Rows of the tables that have different number of columns cannot be compared. Got $(ncol(df1)) and $(ncol(df2)) columns"))
    end
    @inbounds for (col1, col2) in zip(columns(df1), columns(df2))
        isequal(col1[r1], col2[r2]) || return false
    end
    return true
end

# lexicographic ordering on DataFrame rows, missing > !missing
function Base.isless(r1::DataFrameRow, r2::DataFrameRow)
    (ncol(parent(r1)) == ncol(parent(r2))) ||
        throw(ArgumentError("Rows of the data tables that have different number of columns cannot be compared ($(ncol(df1)) and $(ncol(df2)))"))
    @inbounds for i in 1:ncol(parent(r1))
        if !isequal(parent(r1)[i][row(r1)], parent(r2)[i][row(r2)])
            return isless(parent(r1)[i][row(r1)], parent(r2)[i][row(r2)])
        end
    end
    return false
end


text:
module PermutedArrays

using Base.Order

export PermutedVector
struct PermutedVector{T, V} <: AbstractVector{T}
    data::V
    perm::Vector{Int}
    iperm::Vector{Int}
end

function PermutedVector(v::V) where {T, V <: AbstractVector{T}}
    n = length(v)
    perm = collect(1:n)
    iperm = copy(perm)
    PermutedVector{T,V}(v, perm, iperm)
end

Base.size(v::PermutedVector) = size(v.data)
Base.length(v::PermutedVector) = length(v.data)

function Base.permute!(v::PermutedVector, p::Vector{Int})
    permute!(v.perm, p)
    invpermute!(v.iperm, p)
    return v
end

Base.getindex(v::PermutedVector, i::Int) = getindex(v.data, v.perm[i])

function Base.setindex!(v::PermutedVector, x, i::Int)
    setindex!(v.data, x, v.perm[i])
end

function swap!(v::AbstractVector, i::Int, j::Int)
    (v[i], v[j]) = (v[j], v[i])
    return v
end

export swap!
function swap!(v::PermutedVector, i::Int, j::Int)
    p = v.perm
    q = v.iperm
    swap!(p, i, j)
    swap!(q, p[i], p[j])
    return v
end

swap!(::Nothing, i::Int, j::Int) = nothing

function Base.push!(pv::PermutedVector{T}, x::T) where {T}
    push!(pv.data, x)
    n = length(pv.data)
    push!(pv.perm, n)
    push!(pv.iperm, n)
    return pv
end

function Base.deleteat!(pv::PermutedVector, i::Int)
    n = pv.perm[i]
    deleteat!(pv.data, n)

    # Take out an element
    deleteat!(pv.perm, i)
    # Decrement anything above its value to fill in the gap
    for j in eachindex(pv.perm)
        if pv.perm[j] > n
            pv.perm[j] -= 1
        end
    end

    deleteat!(pv.iperm, n)
    for j in eachindex(pv.iperm)
        if pv.iperm[j] > i
            pv.iperm[j] -= 1
        end
    end

end

end # module

max_stars_repo_name: cscherrer/PermutedArrays

max_issues_repo_name: cscherrer/PermutedArrays

max_forks_repo_name: cscherrer/PermutedArrays


text:
# A draft for a "get some data" interface

const _Data = Dict(
  "lattices" => ("data/lattices", "https://github.com/thofma/Hecke.jl/releases/download/DBtest/lattice_test", 2),
  "quadratic_lattices" => ("data/quadratic_lattices", "https://github.com/thofma/Hecke.jl/releases/download/DBtest/qlat_small", 6),
  "hermitian_lattices" => ("data/hermitian_lattices", "https://github.com/thofma/Hecke.jl/releases/download/DBtest/hlat_small", 0.2),
  "small_groups_extended" => ("data/small_groups_extended", "https://github.com/thofma/Hecke.jl/releases/download/DBtest/small_groups_extended", 2.4),
 )

function download_data(;data)
  @assert haskey(_Data, data)
  data = _Data[data]
  println("This will download $(data[3]) MB of data. Proceed? [Y/n]")
  n = readline()
  if n == "" || n == "Y" || n == "yes" || n == "Yes" || n == "y"
    Base.download(data[2], joinpath(pkgdir, data[1]))
    return true
  end
  return false
end


text:
# This file contains code that was a part of Julia
# License is MIT: see LICENSE.md

const StringIndexError = UnicodeError

module Unicode
const normalize  = Base.UTF8proc.normalize_string
const graphemes  = Base.UTF8proc.graphemes
const isassigned = Base.UTF8proc.is_assigned_char
end

const pwc = print_with_color

@api public StringIndexError

Base.replace(str::String, pair::Pair{String,String}; count::Integer=0) =
    replace(str, pair.first, pair.second, count)

const IteratorSize = Base.iteratorsize

const is_letter = isalpha

@api public! IteratorSize

@static if !isdefined(module_parent(current_module()), :Compat)
## Start of code from operators.jl =================================================
##
## It is used to support the new string searching syntax on v0.6.2

@static if !isdefined(Base, :Fix2)
"""
    Fix2(f, x)

A type representing a partially-applied version of function `f`, with the second
argument fixed to the value "x".
In other words, `Fix2(f, x)` behaves similarly to `y->f(y, x)`.
"""
struct Fix2{F,T} <: Function
    f::F
    x::T

    Fix2(f::F, x::T) where {F,T} = new{F,T}(f, x)
    Fix2(f::Type{F}, x::T) where {F,T} = new{Type{F},T}(f, x)
end

(f::Fix2)(y) = f.f(y, f.x)
end

@static if !method_exists(isequal, (Any,))
"""
    isequal(x)

Create a function that compares its argument to `x` using [`isequal`](@ref), i.e.
a function equivalent to `y -> isequal(y, x)`.

The returned function is of type `Base.Fix2{typeof(isequal)}`, which can be
used to implement specialized methods.
"""
Base.isequal(x) = Fix2(isequal, x)

const EqualTo = Fix2{typeof(isequal)}
end

@static if !method_exists(==, (Any,))
"""
    ==(x)

Create a function that compares its argument to `x` using [`==`](@ref), i.e.
a function equivalent to `y -> y == x`.

The returned function is of type `Base.Fix2{typeof(==)}`, which can be
used to implement specialized methods.
"""
Base.:(==)(x) = Fix2(==, x)
end

@static if !method_exists(in, (Any,))
"""
    in(x)

Create a function that checks whether its argument is [`in`](@ref) `x`, i.e.
a function equivalent to `y -> y in x`.

The returned function is of type `Base.Fix2{typeof(in)}`, which can be
used to implement specialized methods.
"""
Base.in(x) = Fix2(in, x)
const OccursIn = Fix2{typeof(in)}
end
end

## end of code from operators.jl =================================================

## Start of codeunits support from basic.jl ======================================
##
##It is used for CodeUnit support in pre v0.7 versions of Julia

## code unit access ##

codeunit(str::AbstractString) = UInt8
codeunit(::Type{<:AbstractString}) = UInt8
ncodeunits(str::AbstractString) = sizeof(str)

"""
    CodeUnits(s::AbstractString)

Wrap a string (without copying) in an immutable vector-like object that accesses the code units
of the string's representation.
"""
struct CodeUnits{T,S<:AbstractString} <: DenseVector{T}
    s::S
    CodeUnits(s::S) where {S<:AbstractString} = new{codeunit(s),S}(s)
end

length(s::CodeUnits) = ncodeunits(s.s)
sizeof(s::CodeUnits{T}) where {T} = ncodeunits(s.s) * sizeof(T)
size(s::CodeUnits) = (length(s),)
strides(s::CodeUnits) = (1,)
getindex(s::CodeUnits, i::Int) = codeunit(s.s, i)
IndexStyle(::Type{<:CodeUnits}) = IndexLinear()
@static if NEW_ITERATE
    iterate(s::CodeUnits, i=1) = (s[i], i+1)
else
    start(s::CodeUnits) = 1
    next(s::CodeUnits, i) = (s[i], i+1)
    @inline done(s::CodeUnits, i) = i > length(s)
end

write(io::IO, s::CodeUnits) = write(io, s.s)

unsafe_convert(::Type{Ptr{T}},    s::CodeUnits{T}) where {T} = unsafe_convert(Ptr{T}, s.s)
unsafe_convert(::Type{Ptr{Int8}}, s::CodeUnits{UInt8}) = unsafe_convert(Ptr{Int8}, s.s)

"""
    codeunits(s::AbstractString)

Obtain a vector-like object containing the code units of a string.
Returns a `CodeUnits` wrapper by default, but `codeunits` may optionally be defined
for new string types if necessary.
"""
codeunits(s::AbstractString) = CodeUnits(s)

## end of codeunits support ============================================================

const UC = Base.UTF8proc

unsafe_crc32c(a, n, crc) = ccall(:jl_crc32c, UInt32, (UInt32, Ptr{UInt8}, Csize_t), crc, a, n)

const utf8crc = Base.crc32c

import Base: isalnum, isgraph, islower, isupper, lcfirst, ucfirst
const is_alphanumeric = isalnum
const is_graphic      = isgraph
const is_lowercase    = islower
const is_uppercase    = isupper
const lowercase_first = lcfirst
const uppercase_first = ucfirst

codepoint(v::Char) = v%UInt32

@noinline index_error(s::AbstractString, i::Integer) =
    strerror(StrErrors.INVALID_INDEX, Int(i), codeunit(s, i))

macro preserve(args...)
    syms = args[1:end-1]
    for x in syms
        isa(x, Symbol) || error("Preserved variable must be a symbol")
    end
    #=
    s, r = gensym(), gensym()
    esc(quote
        $s = $(Expr(:gc_preserve_begin, syms...))
        $r = $(args[end])
        $(Expr(:gc_preserve_end, s))
        $r
        $(args[end])
    end)
    =#
    esc(quote ; $(args[end]) ; end)
end

Base.SubString(str::AbstractString, rng::UnitRange) = SubString(str, first(rng), last(rng))

Base.checkbounds(::Type{Bool}, s::AbstractString, i::Integer) = 1 <= i <= ncodeunits(s)

# Handle some name changes between v0.6 and master
const copyto! = copy!
const unsafe_copyto! = unsafe_copy!
const Nothing = Void
const Cvoid = Void
abstract type AbstractChar end

# Handle changes in array allocation
create_vector(T, len) = Vector{T}(len)

# Add new short name for deprecated hex function
outhex(v, p=1) = hex(v,p)

function get_iobuffer(siz)
    out = IOBuffer(Base.StringVector(siz), true, true)
    out.size = 0
    out
end

@api public copyto!, unsafe_copyto!, Nothing, Cvoid, AbstractChar

import Base: find, ind2chr, chr2ind

const is_lowercase = islower
const is_uppercase = isupper


text:
using Test, NQCDynamics

ω = 3
β = 5
m = 10
mom = MomentumHarmonicWigner(ω, β)
pos = PositionHarmonicWigner(ω, β, m)
vel = VelocityHarmonicWigner(ω, β, m)

norm = tanh(β*ω/2) / π
@test norm ≈ 1 / (mom.σ * pos.σ * 2π) / sqrt(m)
@test norm ≈ 1 / (vel.σ * pos.σ * 2π) / sqrt(m) / m

vec = MomentumHarmonicWigner.([1, 2, 3, 4, 5], β)

d = DynamicalDistribution(vec, vec, (1, 1, 1))
@test_throws ErrorException rand(d)
d = DynamicalDistribution(vec, vec, (2, 5, 2))
rand(d)


text:
##
## Monkhorst-Pack algorithm
##

#gamma_shift( k::Vector ) = map(int,map(iseven,k))

mp_point( x::Integer, l::Integer ) = (2*x - 1 - l)/(2*l)

mp_point_shift( x::Integer, l::Integer ) = ( 2*x - 1 - l + Int(iseven(l)))/(2*l)

# Monkhorst-Pack algorithm in crystallographic coordinates
mp_crystal( n::Vector ) = collect([ [ mp_point(i, n[1]), mp_point(j, n[2]), mp_point(k, n[3]) ] for i=1:n[1], j=1:n[2], k=1:n[3] ])

mp_crystal_sort( k::Vector ) = sort( mp_crystal(k), by=norm )

gamma_crystal( n::Vector ) = collect([ [ mp_point_shift(i, n[1]), mp_point_shift(j, n[2]), mp_point_shift(k, n[3]) ] for i=1:n[1], j=1:n[2], k=1:n[3] ])

gamma_crystal_sort( k::Vector ) = sort( gamma_crystal(k), by=norm )

#function gamma_crystal
#end

# Counter of k-points in BZ
#  num_k( n::Vector ) = length( mp_crystal(n) )

# Converstion of k-points to cartesian coordinates
mp_cart( n::Vector ) = [ reciprocal_lattice'*mp_crystal(n)[i] for i=1:length(mp_crystal(n)) ]

gamma_cart( k::Vector ) = [ reciprocal_lattice'*gamma_crystal(k)[i] for i=1:length(gamma_crystal(k)) ]

# sorting k-points by norm

mp_cart_sort( k::Vector ) = sort( mp_cart(k), by=norm )

gamma_cart_sort( k::Vector ) = sort( gamma_cart(k), by=norm )


text:
using MonteCarloX
using Random
using HypothesisTests
using Distributions
import Distributions.pdf
import Distributions.cdf

""" Testing if sampling single event correctly"""
function test_poisson_single()
    LambdaMaxs = [10.0, 20.0] # testing for two different max-sample-rates
    Lambda(t) = 10.0
    pass = true

    # Regular sampler
    alg = InhomogeneousPoisson()
    for LambdaMax in LambdaMaxs
        rng = MersenneTwister(1000)
        nSamples = 1000
        samples = zeros(nSamples)

        for i in 1:nSamples
            t0 = next_time(alg, rng, Lambda, LambdaMax)
            samples[i] = t0
        end

        test = HypothesisTests.ExactOneSampleKSTest(samples, Exponential(1.0 / 10.0))
        pass &= pvalue(test) > 0.05
    end

    # Decreasing sampler
    alg = InhomogeneousPoissonPiecewiseDecreasing()
    rng = MersenneTwister(1000)
    nSamples = 1000
    samples = zeros(nSamples)

    for i in 1:nSamples
        t0 = next_time(alg, rng, Lambda)
        samples[i] = t0
    end

    test = HypothesisTests.ExactOneSampleKSTest(samples, Exponential(1.0 / 10.0))
    pass &= pvalue(test) > 0.05

    return pass
end

""" Testing if sampling the uniform case correctly"""
function test_poisson_constant()
    LambdaMaxs = [1.0, 2.0] # testing for two different max-sample-rates
    Lambda(t) = 1.0
    pass = true

    alg = InhomogeneousPoisson()
    for LambdaMax in LambdaMaxs
        rng = MersenneTwister(1000)
        nSamples = 1000
        samples = zeros(nSamples)

        t0 = 0
        for i in 1:nSamples
            t0 += next_time(alg, rng, t->Lambda(t + t0), LambdaMax)
            samples[i] = t0 % 1
        end

        test = HypothesisTests.ExactOneSampleKSTest(samples, Uniform(0, 1))
        pass &= pvalue(test) > 0.05
    end

    return pass
end

""" Testing if sampling correctly from a shifted sine-wave distribution."""
function test_poisson_sin_wave()
    LambdaMaxs = [2.0, 3.0] # testing for two different max-sample-rates
    Lambda(t) = sin(t) + 1.0
    pass = true

    alg = InhomogeneousPoisson()
    for LambdaMax in LambdaMaxs
        rng = MersenneTwister(1000)
        nSamples = 1000
        samples = zeros(nSamples)

        t0 = 0
        for i in 1:nSamples
            t0 += next_time(alg, rng, t->Lambda(t + t0), LambdaMax)
            samples[i] = t0 % (2 * pi)
        end

        test = HypothesisTests.ExactOneSampleKSTest(samples, SinDistribution())
        pass &= pvalue(test) > 0.05
    end

    return pass
end

""" Definition of the custom sine-wave distribution for hypothesis testing."""
struct SinDistribution <: ContinuousUnivariateDistribution end

function Distributions.minimum(d::SinDistribution)
    0.0
end

function Distributions.maximum(d::SinDistribution)
    2 * pi
end

function pdf(d::SinDistribution, x::Real)
    insupport(d, x) ? (sin(x) + 1) / (2 * pi) : 0.0
end

function cdf(d::SinDistribution, x::Real)
    insupport(d, x) ? (x - cos(x) + 1.0) / (2 * pi) : Float64(x > 0)
end


text:
"""
Polynomial smoothing with the Savitzky Golay filters. Adapted for julia >= 1.0.
Sources: https://github.com/BBN-Q/Qlab.jl/blob/master/src/SavitskyGolay.jl
         https://gist.github.com/lnacquaroli/c97fbc9a15488607e236b3472bcdf097
Requires LinearAlgebra and DSP modules loaded.
"""

# TODO: could probably improve performance with some @inbounds sprinkling
function savitzkyGolay(x::AbstractVector, windowSize::Int, polyOrder::Int; deriv::Int=0)

  isodd(windowSize) || throw("Window size must be an odd integer.")
  polyOrder < windowSize || throw("Polynomial order must me less than window size.")

  halfWindow = Int( ceil((windowSize-1)/2) )

  # Setup the S matrix of basis vectors
  S = zeros.(windowSize, polyOrder+1)
  for ct = 0:polyOrder
    S[:,ct+1] = (-halfWindow:halfWindow).^(ct)
  end

  ## Compute the filter coefficients for all orders

  # From the scipy code it seems pinv(S) and taking rows should be enough
  G = S * pinv(S' * S)

  # Slice out the derivative order we want
  filterCoeffs = G[:, deriv+1] * factorial(deriv)

  # Pad the signal with the endpoints and convolve with filter
  paddedX = vcat(x[1]*ones(halfWindow), x, x[end]*ones(halfWindow))
  y = conv(filterCoeffs[end:-1:1], paddedX)

  # Return the valid midsection
  return y[2*halfWindow+1:end-2*halfWindow]

end


text:
module NARMAX

include("narmax_node_description.jl")
include("narmax_rules_prototypes.jl")
include("narmax_update_rules.jl")

end

max_stars_repo_name: wmkouw/NARMAX

max_issues_repo_name: wmkouw/NARMAX

max_forks_repo_name: wmkouw/NARMAX


text:
_parallel_functions = [:time_modeling, :fwi_objective, :twri_objective, :extended_source_modeling, :lsrtm_objective]
# Shot and experiment parallelism
for func in _parallel_functions
    @eval $func(args...; kwargs...) = task_distributed($func, _worker_pool(), args...; kwargs...)
end

# Find number of experiments
"""
    _get_nexp(x)

Get number of experiments given a JUDI type. By default we have only one experiment unless we input
a Vector of judiType such as [model, model] to compute gradient for different cases at once.
"""
_get_nexp(x) = 1
for T in [judiVector, Model, judiWeights, judiWavefield, PhysicalParameter, Vector{Float32}]
    @eval _get_nexp(v::Vector{<:$T}) = length(v)
end

"""
    get_nexp(args...)

Get number of experiments given a list of arguments. By default we have only one experiment unless we input
a Vector of judiType such as [model, model] to compute gradient for different cases at once.
All arguments can only be either a Vector of size `n` or a basic type. Basic type are conisdered common
to all `n` experiments. For Example

    VALID: get_nexp(([model1, model2], dobs)...) = 2 (compute a function for common data and two different models)
    INVALID:  get_nexp(([model1, model2], [dobs, dobs, dobs])...)
"""
function get_nexp(args...)
    n_exp = setdiff([_get_nexp(a) for a ∈ args], 1)
    length(n_exp) < 2 || throw(ArgumentError("Incompatible number of experiments in arguments"))
    return max(n_exp..., 1)
end

# Filter arguments for given task
"""
    _get_exp(x, i)

Filter input `x`` for experiment number `i`. Returns `x` is a constant not depending on experiment.
"""
_get_exp(x, i) = x
for T in [judiVector, Model, judiWeights, judiWavefield, Vector{Float32}, PhysicalParameter]
    @eval _get_exp(v::Vector{<:$T}, i) = v[i]
end

"""
    get_exp(i, args...; kwargs...)

Filter input arguments and keyword arguments for experiment number `i`.
"""
get_exp(i, args...; kwargs...) = (_get_exp(a, i) for a in (args..., kwargs.data...))

# Find task iterator (number of sources and indices)
"""
    _get_iter(x)

Get source number iterator from `x`
"""
_get_iter(x::UnitRange) = x
_get_iter(x::judiVector) = 1:x.nsrc
_get_iter(x) = nothing

"""
    get_iter(args...)

Get source number iterator input argument list (number of sources in the data or provided UnitRange)
"""
function get_iter(args...)
    iter = try
        filter(x -> !isnothing(x), [_get_iter(a) for a ∈ args])[1]
    catch BoundsError
        throw(ArgumentError("No iterator found in argument list"))
    end
    return iter, filter(x -> x != iter, args)
end

###
subsample(x, i) = x

"""
    submit_task!(rfunc, pool, i, args...)

Submits a new task running the function `rfunc` on an available worker in the `pool`. The function `rfunc` taskes
for input each element of `args` at index `i`.
"""
submit_task!(rfunc, pool, i, args...) = remotecall(rfunc, pool, (subsample(a, i) for a ∈ args)...)

"""
    filter_exp(x, ::Val)

Filter output based on the number of experiments. Returns `x[1]` for a single experiment or
(sum(x[i][1] for i in 1:length(x)), [x[i][2] for i in 1:length(x)]) (sum of objective values and list of gradient)
for multiple experiments.
"""
filter_exp(x, ::Val{1}) = x[1]
filter_exp(x, ::Val) = gather_nexp(x)

gather_nexp(x) = (sum(x[i][1] for i in 1:length(x)), [x[i][2] for i in 1:length(x)])

"""
    as_vec(x, ::Val{Bool})

Vectorizes output when `return_array` is set to `true`.
"""
as_vec(x, ::Val) = x
as_vec(x::PhysicalParameter, ::Val{true}) = vec(x.data)
for T in [judiVector, judiWeights, judiWavefield]
    @eval as_vec(x::$T, ::Val{true}) = vcat([vec(x.data[i]) for i=1:length(x.data)]...)
end
####

Base.retry(f, pool; kwargs...) = retry(f; kwargs...)

"""
    task_distributed(func, args...; kwargs...)

Distribute the function `func` over the JUDI pool and reduce the result.
"""
function task_distributed(func, pool, args...; kwargs...)
    # Get number of experiments
    nexp = get_nexp(args...)
    judilog("Running $(func) for $(nexp) experiments")
    # Make it retry if fails
    rfunc = retry(func, pool, delays=ExponentialBackOff(n=3))
    # Allocate results
    out = Vector{Any}(undef, nexp)
    # Submit all tasks, asynchronously
    @sync begin
        for e ∈ 1:nexp
            # get local args
            args_e = get_exp(e, args...; kwargs...)
            iter, args_e = get_iter(args_e...)
            # Get options
            opt = findfirst(x -> typeof(x) <: Options, args_e)
            opt = args_e[opt]
            # Allocate Future results
            res_e = Vector{_TFuture}(undef, length(iter))
            judilog("Running experiment $e out of $(nexp) with $(length(iter)) sources")
            flush(stdout)
            for i ∈ iter
                res_e[i] = submit_task!(rfunc, pool, i, args_e...)
            end
            @async out[e] = as_vec(reduce!(res_e), Val(opt.return_array))
        end
    end

    return filter_exp(out, Val(nexp))
end

"""
    single_reduce!(x, y)

Inplace reduction of y into x.
"""
@inline single_reduce!(x, y) = begin x .+= y; nothing end

for T in [judiVector, judiWeights, judiWavefield]
    @eval @inline single_reduce!(x::$T, y) = begin push!(x, y); nothing end
end

@inline single_reduce!(x::T, y) where {T<:Number} = begin x += y; nothing end

@inline function single_reduce!(x::Tuple, y::Tuple)
    nx = length(x)
    ny = length(y)
    (nx == ny) || throw(ArgumentError("Incompatible lengths ($nx, $ny)"))
    for (xi, yi) in zip(x, y)
        single_reduce!(xi, yi)
    end
    nothing
end

"""
    safe_gc()

Generic GC, compatible with different julia versions of it.
"""
safe_gc() = try Base.GC.gc(); catch; gc() end

"""
    local_reduce!(future, other)

Reduce `other` future into local `future`. This is perform remotely on the julia worker
`future.where`.
"""
function local_reduce!(my_future::_TFuture, other_future::_TFuture)
    y = remotecall_fetch(fetch, other_future.where, other_future)
    x = fetch(my_future)
    single_reduce!(x, y)
    # Force GC since julia parallel still has issues with it
    y=1;safe_gc(); remotecall(safe_gc, other_future.where);
    nothing
end

"""
    reduce_level!(futures, nleaf)

Reduce one level of the reduction tree consisting of `nleaf` futures. Each `leaf`
reduces into itself `futures[i]`.
"""
function reduce_level!(futures::Vector{_TFuture}, nleaf::Int)
    nleaf == 0 && return
    @sync for i = 1:2:2*nleaf
        @async remotecall_fetch(local_reduce!, futures[i].where, futures[i], futures[i+1])
    end
    deleteat!(futures, 2:2:2*nleaf)
end

"""
    reduce!(x)
binary-tree reduction of a Vector of Futures `x`

Modified from `DistributedOperations.jl` (Striped from custom types in it) and extended to multiple outputs
with different reduction functions.

Copyright (c) 2020: Chevron U.S.A Inc.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
function reduce!(futures::Vector{_TFuture})
    # Get length a next power of two for binary reduction
    M = length(futures)
    L = round(Int,log2(prevpow(2,M)))
    m = 2^L
    # remainder
    R = M - m

    # Reduce the remainder
    reduce_level!(futures, R)

    # Binary tree reduction
    for l = L:-1:1
        reduce_level!(futures, 2^(l-1))
    end
    fetch(futures[myid()])
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
## probabilistic loss

export crossEntropy
export crossEntropyLoss
export crossEntropyCost

export binaryCrossEntropy
export binaryCrossEntropyLoss
export binaryCrossEntropyCost

"""
    crossEntropy(x::Variable{T}, label::Variable{T}) -> Variable{T}
cross entropy = - y * log(̂y) where y is target and ̂y is the output of the network.
"""
function crossEntropy(x::Variable{T}, label::Variable{T}) where T
    @assert (x.shape == label.shape)
    backprop = (x.backprop || label.backprop)
    ϵ = eltype(x)(1e-38)
    y = Variable{T}(- ᵛ(label) .* log.(ᵛ(x) .+ ϵ), backprop)
    if backprop
        function crossEntropyBackward()
            if need2computeδ!(x)
                δ(x) .-= δ(y) .* ᵛ(label) ./ (ᵛ(x) .+ ϵ)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, crossEntropyBackward)
    end
    return y
end

crossEntropyLoss(x::Variable{T}, label::Variable{T}; reduction::String="sum") where T = loss( crossEntropy(x, label), reduction=reduction )
crossEntropyCost(x::Variable{T}, label::Variable{T}) where T = cost( crossEntropy(x, label) )

"""
    binaryCrossEntropy(x::Variable{T}, l::Variable{T}) -> Variable{T}
binary cross entropy = - y * log(̂y) - (1 - y) * log(1-̂y)
"""
function binaryCrossEntropy(x::Variable{T}, label::Variable{T}) where T
    @assert (x.shape == label.shape)
    backprop = (x.backprop || label.backprop)
    TOO  = eltype(x)
    ϵ  = TOO(1e-38)
    𝟙  = TOO(1.0f0)
    tmp1 = - ᵛ(label) .* log.(ᵛ(x) .+ ϵ)
    tmp2 = - (𝟙 .- ᵛ(label)) .* log.(𝟙 .- ᵛ(x) .+ ϵ)
    y  = Variable{T}(tmp1 + tmp2, backprop)
    if backprop
        function binaryCrossEntropyBackward()
            if need2computeδ!(x)
                temp1 = (𝟙 .- ᵛ(label)) ./ (𝟙 .- ᵛ(x) .+ ϵ)
                temp2 = ᵛ(label) ./ (ᵛ(x) .+ ϵ)
                δ(x) .+= δ(y) .* (temp1 - temp2)
            end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, binaryCrossEntropyBackward)
    end
    return y
end

binaryCrossEntropyLoss(x::Variable{T}, label::Variable{T}; reduction::String="sum") where T = loss( binaryCrossEntropy(x, label), reduction=reduction)
binaryCrossEntropyCost(x::Variable{T}, label::Variable{T}) where T = cost( binaryCrossEntropy(x, label) )

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
"""

```
md_berendsen(sys::System{T},opt::Options) where T
```

Performs a MD simulation with a Berendsen bath.

"""
function md_berendsen(sys::System{T},opt::Options=Options()) where T

  println("""
  -----------------------------------------------------------------------------------------
  NVT - Berendsen-bath simulation
  -----------------------------------------------------------------------------------------
  """)

  # aliases to simplify code
  @unpack n, x0 = sys
  @unpack dt, nsteps, kT, iequil, tau = opt
  u(x) = potential(x,sys,opt)
  f!(f,x) = forces!(f,x,sys,opt)

  # Copy the initial point, to preserve it
  x = copy(x0)

  # Obtain initial velocities
  v = init_velocities(sys,opt)

  # Open trajectory file for writting
  trajectory_file = open(opt.trajectory_file,"w")
  velocity_file = open(opt.velocity_file,"w")

  # log matrix will contain potential, kinetic, total energies, and "temperature"
  out = zeros(nsteps,4)

  println(" Potential energy at initial point: " , u(x))
  println(" Kinetic energy at initial point: ", kinetic(v))
  println(" Total initial energy = ", u(x) + kinetic(v))

  # Write coordinates to trajectory file, and update log vectors
  printxyz(0.,x,v,sys,trajectory_file,velocity_file)
  out[1,:] .= (u(x), kinetic(v), u(x) + kinetic(v), kinetic(v)/n)

  # Initialize velocity vector and save first set of forces
  f = zeros(T,n)
  f!(f,x)
  flast = copy(f)

  # Running simulation
  time = 0.
  for istep = 1:nsteps

    # Write coordinates to trajectory file, and update log vectors

    # Updating positions
    @. x = x + v*dt + 0.5*f*dt^2

    # Updating the forces
    @. flast = f
    f!(f,x)

    # Updating velocities
    @. v = v + 0.5*(f + flast)*dt
    remove_drift!(v)

    # Update step and print data
    time += dt
    kstep = kinetic(v)
    ustep = u(x)
    energy = kstep + ustep
    kavg = kstep / n

    # Stop if simulation exploded
    if ustep > 1e10
      println(" Simulation exploded: Energy = ", energy)
      close(trajectory_file)
      return out[1:istep,:]
    end

    # Save point to log
    out[istep,:] .= (ustep,kstep,energy,kavg)
    if mod(istep,opt.iprint) == 0
      @printf(" TIME = %12.3f U = %12.5e K = %12.5e TOT = %12.5e \n",
        time, ustep, kstep, energy
      )
      out[istep,:] .= (ustep,kstep,energy,kavg)
    end
    if mod(istep,opt.iprintxyz) == 0
      printxyz(time,x,v,sys,trajectory_file,velocity_file)
    end

    #
    # Berendsen-bath
    #
    if istep <= iequil
      lambda = sqrt( 1 + (dt/tau)*((dim(T)*kT/2)/kavg-1) )
      @. v = v * lambda
    end

  end

  close(trajectory_file)
  close(velocity_file)
  println(" Wrote trajectory file: ", opt.trajectory_file)
  println(" Wrote velocities file: ", opt.velocity_file)
  return out
end


text:
export
    boostNet,   # BoostNet
    sigma,      # activition function
    Net2,       # two layer boost net model
    predict,    # output +-1
    normalize1  # normalize to be weight 1

"""
    BoostNet

arguments:

  features,     column vectors
  labels,
  m,            # of layers (2)
  delta,        # 1 - delta is the probability of getting a good result
  gamma,        # gamma separable
  T,            # of hidden nodes
  B,            [-B, B]^k cube (10)

note:
  delta, gamma, not used.

"""
function boostNet(features, labels; m=2, delta=.5, gamma=1, T=10, B=10)
  d, n = size(features) # feature dim
  size(labels) == (1, n) || error("boostNet: column number of features should equal to number of labels")

  net = Net2(1, zeros(T), zeros(T, d))  # init f0 = 0, b0 = 0

  for t = 1:T
    alpha = normalize1(exp(-labels .* sigma(evaluate(net, features))))

    assert(size(alpha) == (1, n))

    w = alg3(features, labels, m=1, k=10, T=1, B=10, weight=alpha[:])

    assert(size(w) == (1, d))

    G = sum(alpha .* sigma(-labels .* sigma(w * features)))

    mu = max(-.5, G)
    b = .5 * log((1 - mu) / (1 + mu))
    net.t = t
    net.b[t] = b
    net.W[t, :] = w
  end

  net
end

" normalize to be weight 1 "
normalize1(v) = v / sum(v)

" activition function "
sigma = tanh

" two layer boost net model "
type Net2
  t   # how many available
  b   # a column vector
  W   # w matrix, each row is a w
end

" output value "
evaluate(m::Net2, features) = sum(sigma(slice(m.W, 1:m.t, :) * features) .* m.b[1:m.t], 1)

" output +-1 "
predict(m::Net2, features) = (evaluate(m, features) .> 0) * 2 - 1

max_stars_repo_name: innerlee/Halfspace

max_issues_repo_name: innerlee/Halfspace

max_forks_repo_name: innerlee/Halfspace


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    Point2D, Point3D, Segment,
    toarray, x_components, y_components, z_components,
    distance, unitvector, midpoint,
    translate, trim_u, trim_v, trim_uv,
    cross2d, interiorangle, isclockwise,
    transformmatrix, rotation

struct Point2D <: Point
    x::Float64
    y::Float64
end

struct Point3D <: Point
    x::Float64
    y::Float64
    z::Float64
end

struct Segment{T<:Point}
    u::T
    v::T
end

Point2D(x::Real, y::Real) = Point2D(float(x), float(y))
Point2D(t::Tuple{T,T}) where {T<:Real} = Point2D(float(t[1]), float(t[2]))
Point2D(arr::AbstractArray{T}
    ) where {T<:Real} = Point2D(float(arr[1]), float(arr[2]))
Point2D(coords::AbstractMatrix{T}, i::Int
    ) where {T<:Real} = Point2D(coords[i, 1:2])

Point3D(x::Real, y::Real, z::Real) = Point3D(float(x), float(y), float(z))
Point3D(t::Tuple{T,T,T}
    ) where {T<:Real} = Point3D(float(t[1]), float(t[2]), float(t[3]))
Point3D(arr::AbstractArray{T}
    ) where {T<:Real} = Point3D(float(arr[1]), float(arr[2]), float(arr[3]))
Point3D(coords::AbstractMatrix{T}, i::Int
    ) where {T<:Real} = Point3D(coords[i, 1:3])

Segment{T}(coords::AbstractMatrix{Float64}, u::Int, v::Int
    ) where {T<:Point} = Segment(T(coords, u), T(coords, v))

Base.:+(u::Point2D, v::Point2D) = Point2D(u.x + v.x,  u.y + v.y)
Base.:+(p::Point2D, transl) = Point2D(p.x + transl[1],  p.y + transl[2])
Base.:+(u::Point3D, v::Point3D) = Point3D(u.x + v.x,  u.y + v.y,  u.z + v.z)
Base.:+(p::Point3D, transl
    ) = Point3D(p.x + transl[1],  p.y + transl[2],  p.z + transl[3])
Base.:+(s::Segment, transl) = Segment(s.u + transl, s.v + transl)

Base.:-(u::Point2D, v::Point2D) = Point2D(u.x - v.x,  u.y - v.y)
Base.:-(p::Point2D, transl) = Point2D(p.x - transl[1],  p.y - transl[2])
Base.:-(u::Point3D, v::Point3D) = Point3D(u.x - v.x,  u.y - v.y,  u.z - v.z)
Base.:-(p::Point3D, transl
    ) = Point3D(p.x - transl[1],  p.y - transl[2],  p.y - transl[3])
Base.:-(s::Segment, transl) = Segment(s.u - transl, s.v - transl)

Base.:*(p::Point2D, f::Float64) = Point2D(p.x * f, p.y * f)
Base.:*(f::Float64, p::Point2D) = Point2D(p.x * f, p.y * f)
Base.:*(p::Point3D, f::Float64) = Point3D(p.x * f, p.y * f, p.z * f)
Base.:*(f::Float64, p::Point3D) = Point3D(p.x * f, p.y * f, p.z * f)
Base.:*(s::Segment, f::Real) = Segment(s.u * f, s.v * f)
Base.:*(f::Real, s::Segment) = Segment(s.u * f, s.v * f)

toarray(p::Point2D) = [p.x, p.y]
toarray(p::Point3D) = [p.x, p.y, p.z]
toarray(u::Point, v::Point) = transpose(hcat(toarray(u), toarray(v)))
toarray(s::Segment) = toarray(s.u, s.v)
toarray(coords::AbstractMatrix{T}) where {T<:Real} = coords
toarray(coords::AbstractMatrix{T}, i::Int
    ) where {T<:Real} = @view coords[i, :]
toarray(coords::AbstractMatrix{T}, u::Int, v::Int
    ) where {T<:Real} = @view coords[[u, v], :]
toarray(coords::AbstractMatrix{T}, vertices::Vector{Int}
    ) where {T<:Real} = @view coords[vertices, :]

x_components(coords::AbstractMatrix{T}) where {T<:Real} = @view coords[:, 1]
y_components(coords::AbstractMatrix{T}) where {T<:Real} = @view coords[:, 2]
z_components(coords::AbstractMatrix{T}) where {T<:Real} = @view coords[:, 3]

LinearAlgebra.norm(p::Point) = norm(toarray(p))
LinearAlgebra.normalize(p::Point) = Point2D(normalize(toarray(p)))
LinearAlgebra.dot(u::Point, v::Point) = dot(toarray(u), toarray(v))

"""
    distance(u::Point, v::Point) -> Float64
    distance(s::Segment) -> Float64
    distance(coords::AbstractMatrix{T}) where {T<:Real} -> Float64

Return distance between two endpoints.
"""
distance(u::Point, v::Point) = norm(v - u)
distance(s::Segment) = distance(s.u, s.v)
distance(coords::AbstractMatrix{T}
    ) where {T<:Real} = norm(coords[2, :] - coords[1, :])

"""
    unitvector(u::Point, v::Point) -> Point
    unitvector(s::Segment) -> Point
    unitvector(coords::AbstractMatrix{T}) where {T<:Real} -> AbstractMatrix

Return u -> v vector of length 1.
"""
unitvector(u::Point, v::Point) = normalize(v - u)
unitvector(s::Segment) = unitvector(s.u, s.v)
unitvector(coords::AbstractMatrix{T}
    ) where {T<:Real} = normalize(coords[2, :] - coords[1, :])

"""
    midpoint(u::Point, v::Point) -> Point
    midpoint(s::Segment) -> Point
    midpoint(coords::AbstractMatrix{T}) where {T<:Real} -> AbstractMatrix

Return the midpoint of u and v.
"""
midpoint(u::Point, v::Point) = (u + v) * 0.5
midpoint(s::Segment) = midpoint(s.u, s.v)
midpoint(coords::AbstractMatrix{T}
    ) where {T<:Real} = (s[1, :] + s[2, :]) * 0.5

function translate(seg, ang, dist)
    rotation = [cos(ang) -sin(ang); sin(ang) cos(ang)]
    move = rotation * toarray(unitvector(seg)) * dist
    return seg + move
end

trim_u(s::Segment, k) = Segment(s.u + (s.v - s.u) * k, s.v)
trim_v(s::Segment, k) = Segment(s.u, s.v + (s.v - s.u) * -k)
trim_uv(s::Segment, k
    ) = Segment(s.u + (s.v - s.u) * 0.5k, s.v + (s.v - s.u) * -0.5k)

cross2d(u::Point2D, v::Point2D) = u.x * v.y - u.y * v.x
cross2d(u, v) = u[1] * v[2] - u[2] * v[1]

interiorangle(u, v) = acos(dot(u, v) / (norm(u) * norm(v)))

"""
    isclockwise(vertices::AbstractMatrix{Float64}) -> Union{Bool,Nothing}

Return true/false if given vertices of a polygon in 2D space are placed clockwise/anticlockwise. Return nothing if the polygon is self-intersecting or some vertices are overlapped.
"""
function isclockwise(vertices::AbstractMatrix{Float64})
    vlen = size(vertices, 1)
    clockwise = 0
    counter = 0
    cycpath = vcat(vertices, vertices)
    for i in 1:vlen
        p = cycpath[i, :]
        q = cycpath[i+1, :]
        r = cycpath[i+2, :]
        cp = cross2d(q - p, r - p)
        intangle = interiorangle(p - q, r - q)
        if cp < 0
            clockwise += intangle
            counter += 2pi - intangle
        else
            clockwise += 2pi - intangle
            counter += intangle
        end
    end
    if round(clockwise / pi) == vlen - 2
        return true
    elseif round(counter / pi) == vlen - 2
        return false
    end
    return # overlapped or self-intersecting
end

function transformmatrix(scale::Point2D, rotate::Point2D, transl::Point2D)
    s = Matrix{Float64}([scale.x 0 0; 0 scale.y 0; 0 0 1])
    r = Matrix{Float64}([rotate.x -rotate.y 0; rotate.y rotate.x 0; 0 0 1])
    t = Matrix{Float64}([1 0 transl.x; 0 1 transl.y; 0 0 1])
    return t * r * s
end

function rotation(axis::Point3D, angle)
    (x1, y1, z1) = (axis.x, axis.y, axis.z)
    c = cos(angle)
    s = sin(angle)
    a12 = x1 * y1 * (1 - c)
    a13 = x1 * z1 * (1 - c)
    a23 = y1 * z1 * (1 - c)
    return Matrix{Float64}([
        (c + x1^2 * (1 - c)) (a12 - z1 * s) (a13 + y1 * s);
        (a12 + z1 * s) (c + y1^2 * (1 - c)) (a23 - x1 * s);
        (a13 - y1 * s) (a23 + x1 * s) (c + z1^2 * (1 - c))
    ])
end


text:
ProjDir = @__DIR__
cd(ProjDir)

open(`$(nodejs_cmd())`, "w", stdout) do io
  for i in 1:4
    println(io, "console.log('Hi #$i');")
  end
  println(io, "console.log('Hi #5');")
  println(io, "a = 6")
  println(io, "console.log(a);")
  #close(io)
end

mutable struct Pipes
    gstdin :: Pipe
    gstdout :: Pipe
    gstderr :: Pipe
    Pipes() = new()
end

const P = Pipes()

# initialize gnuplot
function init_node()
  global P
  try
      success(`$(nodejs_cmd()) --version`)
  catch
      error("NodeJS cannot be initialized.")
  end
  gstdin = Pipe()
  gstdout = Pipe()
  gstderr = Pipe()
  gproc = run(pipeline(`$(nodejs_cmd())`,
    stdin = gstdin, stdout = gstdout, stderr = gstderr),
    wait = false)
  process_running(gproc) || error("There was a problem starting up nodejs.")
  close(gstdout.in)
  close(gstderr.in)
  close(gstdin.out)
  P.gstdin = gstdin
  P.gstdout = gstdout
  P.gstderr = gstderr

  return nothing
end

# write commands to gnuplot's pipe
function nodejs_send(s)
    config[:debug] && println(s)  # print gnuplot commands if debug enabled
    w = write(P.gstdin, string(s,"\n"))
    # check that data was accepted by the pipe
    if !(w > 0)
        println("Something went wrong writing to gnuplot STDIN.")
        return
    end
    flush(P.gstdin)
end

## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Asynchronously reads the IO and buffers the read content. When end marker
# (GastonDone, which is stored in `gmarker` to account for Unix/Windows
# variability in line endings) is found in the content, sends everything
# between start (GastonBegin) and end markers to the returned channel.
# In case of timeout sends :timeout; in case of end of file, sends :eof.

function async_reader(io::IO, timeout_sec)::Channel
    ch = Channel(1)
    task = @async begin
        reader_task = current_task()
        function timeout_cb(timer)
            put!(ch, :timeout)
            Base.throwto(reader_task, InterruptException())
        end

        buf = ""
        while (match_done = findfirst(gmarker_done, buf)) == nothing
            timeout = Timer(timeout_cb, timeout_sec)
            data = String(readavailable(io))
            if data == ""; put!(ch, :eof); return; end
            timeout_sec > 0 && close(timeout) # Cancel the timeout
            buf *= data
        end
        match_begin = findfirst(gmarker_start, buf)
        start = (match_begin != nothing) ? last(match_begin)+1 : 1
        put!(ch, buf[start:first(match_done)-1])
    end
    bind(ch, task)
    return ch
end


text:
using PopulationDynamics
using Test

SKIP = true
if !SKIP
    ### Load ###
    bsnm = "test"
    datadir = abspath("test")
    metadata = EPhys.getmetadata(bsnm, datadir)
    sessions = [2];

    spiketimes = loadspikes(SpikeTimes, metadata, sessions)
    spiketimes = filter(x->length(x)>100, spiketimes)
    Z = binZ(spiketimes)'
    t = loadtrigtimes("toneSucrose_pulse", metadata, sessions)
    Z3 = binZ3(spiketimes,t)

    ### Test ###
    function testmodel(model,name, Z)
        println("   $name")
        println("       fit")
        fit!(model, Z)
        @test typeof(weights(model)) == Array{Float64,2}
        println("       track")
        @test track(model, zeros(size(Z,1),100)) == zeros(model.k,100)
        println("       bootstrap")
        btstrp = bootstrap(model, Z, 3)
        @test typeof(btstrp) == Array{typeof(model),1}
        println("       confint")
        lo, hi = confint(weights, model, btstrp)
        @test sum(hi) > sum(lo)
    end

    println("Testing 2nd Order Models...")
    modelnames = ["PCA",
                  "PCAICA"]
    models = [PCAICA(),
              PCA()]
    testmodel.(models, modelnames, [Z])

    println("")
    println("Testing 3rd Order Models...")
    modelnames = ["Canonical Polyadic Decomposition (CPD)"]
    models = [CP(2)]
    testmodel.(models, modelnames, [Z3])
end


text:
module JuliaCon2020RANSAC

import Remark
import RANSAC
using RANSAC: IterLow1
import RANSACVisualizer
using ColorSchemes

using Logging
using LinearAlgebra
using ExtractMacro
using UnionFind
using NearestNeighbors: KDTree, inrange
using StaticArrays: SVector

export  makeslides,
        openslides

export  SLIDE_FOLDER

const SLIDE_FOLDER = joinpath(dirname(@__DIR__), "slides")

makeslides() = Remark.slideshow(SLIDE_FOLDER, options=Dict("ratio"=>"16:9"),
    title = "Efficient RANSAC in efficient Julia", documenter=false)

openslides() = Remark.open(SLIDE_FOLDER)

export  FittedTranslational,
        ExtractedTranslational

include("translational/measures.jl")
include("translational/profit.jl")
include("translational/translational.jl")

end # module


text:
"""
Created on Fri May 01 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Rank selection method in NMF based on average stability of latent factor across perturbed versions of the input matrix.
matrix V.
Alexandrov et al, Deciphering Signatures of Mutational Processes Operative in Human Cancer. 2013, Cell Reports.
"""

abstract type AbstractRSParams <: AbstractNMFParams end

"""
    RSParams([,K_min, K_max, n_iter, pert, seed])

Parameters passed to `rank_by_stability`.

# Arguments
- `K_min::Integer`     : Minimum rank evaluated.
- `K_max::Integer`     : Maximum rank evaluated.
- `n_iter::Integer`    : Total number of iterations, each iteration corresponding to one run of NMF on one perturbed
    version of the matrix `V`.
- `pert_meth::Symbol`  : Perturbation method of the matrix `V` at each iteration. Choose :multinomial for random draws
    from a Multinomial distribution or :noise for another distribution as specified by `pert_dist`.
- `pert_dist::Distribution` : Distribution from  which additive random noise will be drawn in case `pert_meth` is :noise.
- `seed::Integer`      : Seed for the random numbner generator.
"""
mutable struct RSParams <: AbstractRSParams
    K_min::Integer
    K_max::Integer
    n_iter::Integer
    pert_meth::Symbol
    pert_dist::Union{Nothing, Distribution}
    seed::Integer

    function RSParams(;
        K_min::Integer                         = 2,
        K_max::Integer                         = 10,
        n_iter::Integer                        = 250,
        pert_meth::Symbol                      = :multinomial,
        pert_dist::Union{Nothing,Distribution} = nothing,
        seed::Integer                          = 0)

        K_min > 0 || throw(ArgumentError("K_min must be superior or equal to 0"))
        K_max >= K_min || throw(ArgumentError("K_max must be superior or equal to K_min"))
        n_iter > 0 || throw(ArgumentError("n_iter must be superior or equal to 1"))

        new(K_min, K_max, n_iter, pert_meth, pert_dist, seed)
    end
end

abstract type AbstractRSNMFResults <: AbstractNMFResults end

"""
    RSNMFResults(grid_nmf_results, list_rank, list_clu_results, list_fid_results)

Struct for storing results returned by `rank_by_stability`.
"""
mutable struct RSNMFResults{T} <: AbstractRSNMFResults
    grid_nmf_results::Vector{NMFResults{T}}
    df_metrics::DataFrame
    list_clu_results::Vector{RSCluNMFResults{T}}
    list_fid_results::Vector{RSFidNMFResults{T}}

    function RSNMFResults{T}(
        grid_nmf_results::Vector{NMFResults{T}},
        df_metrics::DataFrame,
        list_clu_results::Vector{RSCluNMFResults{T}},
        list_fid_results::Vector{RSFidNMFResults{T}}) where T <: Real

        new{T}(grid_nmf_results, df_metrics, list_clu_results, list_fid_results)
    end
end

"""
    _rank_by_stability_one((V_pert, nmf))

Helper function for `rank_by_stability`.
"""
function _rank_by_stability_one((V, nmf)) where T <: Real
    nmf.solver(V, nmf.global_params, nmf.local_params)
end

"""
    rank_by_stability(V, rs_params, nmf)

Rank selection method in NMF based on average stability of latent factor across perturbed versions of the input matrix
V.

# Arguments
- `V::Matrix{T}` : matrix to be factorised
- `rs_params::RSParams` : see `RSParams`(@ref)
- `nmf::NMF` : see `NMF`(@ref)

# References
- Alexandrov et al, Deciphering Signatures of Mutational Processes Operative in Human Cancer. 2013, Cell Reports.
"""
function rank_by_stability(V::Matrix{T}, rs_params::RSParams, nmf::NMF) where T <: Real
    #### fix the rng
    rng = Random.seed!(rs_params.seed)

    K_min  = rs_params.K_min
    K_max  = rs_params.K_max
    n_iter = rs_params.n_iter
    pert_meth = rs_params.pert_meth
    pert_dist = rs_params.pert_dist

    #### fix seeds for W and H initializations
    seeds    = randperm(rng, n_iter * (K_max - K_min + 1))
    grid_nmf = Array{Tuple{Matrix{T}, NMF}, 1}()

    for (i, K) in enumerate(K_min:K_max)
        for j in 1:n_iter
            V_pert = _get_perturbed_matrix(V, pert_meth, pert_dist, rng)
            nmf_K = copy(nmf)
            nmf_K.global_params.rank = K
            nmf_K.global_params.seed = seeds[(i-1) * n_iter + j]
            push!(grid_nmf, (V_pert, nmf_K))
        end
    end

    #### run in parallel
    grid_nmf_results = @showprogress pmap(grid_nmf) do g
        _rank_by_stability_one(g)
    end

    #### get reconstruction and stability results by rank
    grid_V_pert  = (x -> x[1]).(grid_nmf)
    grid_rank    = (x -> x[2].global_params.rank).(grid_nmf)

    list_rank = Vector{Int64}()
    list_clu_results = Vector{RSCluNMFResults{T}}()
    list_fid_results = Vector{RSFidNMFResults{T}}()
    list_stab_avg    = Vector{Float64}()
    list_stab_std    = Vector{Float64}()
    list_fid_avg     = Vector{Float64}()
    list_fid_std     = Vector{Float64}()

    for rank = K_min:K_max
        clu_results = evaluate_stability(grid_nmf_results[grid_rank .== rank], dist_func=CosineDist(), rng=rng)
        fid_results = evaluate_fidelity(grid_nmf_results[grid_rank .== rank], grid_V_pert[grid_rank .== rank])

        push!(list_rank, rank)
        push!(list_clu_results, clu_results)
        push!(list_fid_results, fid_results)

        push!(list_stab_avg, mean(clu_results.stab_avg))
        push!(list_stab_std, std(clu_results.stab_avg))
        push!(list_fid_avg, mean(fid_results.list_fid))
        push!(list_fid_std, std(fid_results.list_fid))
    end

    df_metrics = DataFrame(
        rank     = list_rank,
        stab_avg = list_stab_avg,
        stab_std = list_stab_std,
        fid_avg  = list_fid_avg,
        fid_std  = list_fid_std,
       )

    RSNMFResults{T}(
        grid_nmf_results,
        df_metrics,
        list_clu_results,
        list_fid_results
    )
end


text:
module RLInterface

using POMDPs
using POMDPModelTools

# for the ZMQ part
using ZMQ
using JSON
using Random

export
    # Environment types
    AbstractEnvironment,
    POMDPEnvironment,
    MDPEnvironment,
    KMarkovEnvironment,
    # supporting methods
    reset!,
    step!,
    actions,
    sample_action,
    n_actions,
    obs_dimensions,
    render,
    # deprecated
    reset

abstract type AbstractEnvironment{OV} end

"""
    obsvector_type(::Union{MDP, POMDP})

Returns the type of the observation vector associated with a specific problem.
The `MDPEnvironment` and `POMDPEnvironment` wrappers will convert observations to an object of such type when `reset!` or `step!` is called.
"""
obsvector_type(::Union{MDP, POMDP}) = Vector{Float32}

mutable struct MDPEnvironment{OV, M<:MDP, S, R<:AbstractRNG} <: AbstractEnvironment{OV}
    problem::M
    state::S
    rng::R
end
function MDPEnvironment(problem::M,
                        ov::Type{A} = obsvector_type(problem);
                        rng::R=MersenneTwister(0)) where {A<:AbstractArray, M<:MDP, R<:AbstractRNG}
    S = statetype(problem)
    return MDPEnvironment{ov, M, S, R}(problem, rand(rng, initialstate(problem)), rng)
end

mutable struct POMDPEnvironment{OV, M<:POMDP, S, R<:AbstractRNG} <: AbstractEnvironment{OV}
    problem::M
    state::S
    rng::R
end
function POMDPEnvironment(problem::M,
                          ov::Type{A} = obsvector_type(problem);
                          rng::R=MersenneTwister(0)) where {A<:AbstractArray, M<:POMDP, R<:AbstractRNG}
    S = statetype(problem)
    return POMDPEnvironment{ov, M, S, R}(problem, rand(rng, initialstate(problem)), rng)
end

"""
    reset!(env::MDPEnvironment{OV})
Reset an MDP environment by sampling an initial state returning it.
"""
function reset!(env::MDPEnvironment{OV}) where OV
    s = rand(env.rng, initialstate(env.problem))
    env.state = s
    return convert_s(OV, s, env.problem)
end

"""
    reset!(env::POMDPEnvironment{OV})
Reset an POMDP environment by sampling an initial state,
generating an observation and returning it.
"""
function reset!(env::POMDPEnvironment{OV}) where OV
    s = rand(env.rng, initialstate(env.problem))
    env.state = s
    o = rand(env.rng, initialobs(env.problem, s))
    return convert_o(OV, o, env.problem)
end

"""
    step!(env::MDPEnvironment{OV}, a::A)
Take in an POMDP environment, and an action to execute, and
step the environment forward. Return the state, reward,
terminal flag and info
"""
function step!(env::MDPEnvironment{OV}, a::A) where {OV, A}
    s, r, info = @gen(:sp, :r, :info)(env.problem, env.state, a, env.rng)
    env.state = s
    t = isterminal(env.problem, s)
    obs = convert_s(OV, s, env.problem)
    return obs, r, t, info
end

"""
    step!(env::POMDPEnvironment{OV}, a::A)
Take in an MDP environment, and an action to execute, and
step the environment forward. Return the observation, reward,
terminal flag and info
"""
function step!(env::POMDPEnvironment{OV}, a::A) where {OV, A}
    s, o, r, info = @gen(:sp, :o, :r, :info)(env.problem, env.state, a, env.rng)
    env.state = s
    t = isterminal(env.problem, s)
    obs = convert_o(OV, o, env.problem)
    return obs, r, t, info
end

"""
    actions(env::Union{POMDPEnvironment, MDPEnvironment})
Return an action object that can be sampled with rand.
"""
function POMDPs.actions(env::Union{POMDPEnvironment, MDPEnvironment})
    return actions(env.problem)
end

"""
    sample_action(env::Union{POMDPEnvironment, MDPEnvironment})
Sample an action from the action space of the environment.
"""
function sample_action(env::Union{POMDPEnvironment, MDPEnvironment})
    return rand(env.rng, actions(env))
end

"""
    obs_dimensions(env::MDPEnvironment{OV}) where OV
returns the size of the observation vector.
It generates an initial state, converts it to an array and returns its size.
"""
function obs_dimensions(env::MDPEnvironment{OV}) where OV
    return size(convert_s(OV, rand(env.rng, initialstate(env.problem)), env.problem))
end

"""
    obs_dimensions(env::POMDPEnvironment{OV}) where OV
returns the size of the observation vector.
It generates an initial observation, converts it to an array and returns its size.
"""
function obs_dimensions(env::POMDPEnvironment{OV}) where OV
    s = rand(env.rng, initialstate(env.problem))
    o = rand(env.rng, initialobs(env.problem, s))
    return size(convert_o(OV, o, env.problem))
end

"""
    render(env::AbstractEnvironment)
Renders a graphic of the environment
"""
function render(env::AbstractEnvironment) end

include("ZMQServer.jl")
include("k_markov.jl")

# deprecations
import Base.reset
@deprecate reset(env::KMarkovEnvironment) reset!(env)
@deprecate reset(env::POMDPEnvironment) reset!(env)
@deprecate reset(env::MDPEnvironment) reset!(env)

end # module