text:
@testset "Basic CHP operations" begin

			using Juqst

           t = Tableau(4)
           start = Tableau(4)
           @test size(t.state)==(8,9)
           @test t.state ==  [1  0  0  0  0  0  0  0  0
							 0  1  0  0  0  0  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  0
							 0  0  0  0  1  0  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  0  0  0  0  1  0  0
							 0  0  0  0  0  0  0  1  0]
           hadamard(t,2)
           @test t.state == [1  0  0  0  0  0  0  0  0
 							 0  0  0  0  0  1  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  0
							 0  0  0  0  1  0  0  0  0
							 0  1  0  0  0  0  0  0  0
							 0  0  0  0  0  0  1  0  0
							 0  0  0  0  0  0  0  1  0]
		   phase(t,1)
		   @test t.state == [1  0  0  0  1  0  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  0
							 0  0  0  0  1  0  0  0  0
							 0  1  0  0  0  0  0  0  0
							 0  0  0  0  0  0  1  0  0
							 0  0  0  0  0  0  0  1  0]
			cnot(t,2,1)
			@test t.state == [1  0  0  0  1  1  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  0
							 0  0  0  0  1  1  0  0  0
							 1  1  0  0  0  0  0  0  0
							 0  0  0  0  0  0  1  0  0
							 0  0  0  0  0  0  0  1  0]
			X(t,3)
			@test t.state == [ 1  0  0  0  1  1  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  0
							 0  0  0  0  1  1  0  0  0
							 1  1  0  0  0  0  0  0  0
							 0  0  0  0  0  0  1  0  1
							 0  0  0  0  0  0  0  1  0]
			Z(t,4)
			@test t.state == [1  0  0  0  1  1  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  1
							 0  0  0  0  1  1  0  0  0
							 1  1  0  0  0  0  0  0  0
							 0  0  0  0  0  0  1  0  1
							 0  0  0  0  0  0  0  1  0]
			Z(t,3)
			@test t.state == [1  0  0  0  1  1  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  1  0  0  0  0  0  1
							 0  0  0  1  0  0  0  0  1
							 0  0  0  0  1  1  0  0  0
							 1  1  0  0  0  0  0  0  0
							 0  0  0  0  0  0  1  0  1
							 0  0  0  0  0  0  0  1  0]
			@test kets(t) == "+|0010> \n+|1110> \n"
			@test t.commands == ["initialise(4)"
								 "hadamard(2)"
								 "phase(1)"
								 "cnot(2,1)"
								 "hadamard(3)"
								 "phase(3)"
								 "phase(3)"
								 "hadamard(3)"
								 "phase(4)"
								 "phase(4)"
								 "phase(3)"
								 "phase(3)" ]
			t.trackCommands=false
			hadamard(t,4)
			@test t.commands == ["initialise(4)"
								 "hadamard(2)"
								 "phase(1)"
								 "cnot(2,1)"
								 "hadamard(3)"
								 "phase(3)"
								 "phase(3)"
								 "hadamard(3)"
								 "phase(4)"
								 "phase(4)"
								 "phase(3)"
								 "phase(3)" ]

			initialise(t)
			@test t.state ==  [1  0  0  0  0  0  0  0  0
							 0  1  0  0  0  0  0  0  0
							 0  0  1  0  0  0  0  0  0
							 0  0  0  1  0  0  0  0  0
							 0  0  0  0  1  0  0  0  0
							 0  0  0  0  0  1  0  0  0
							 0  0  0  0  0  0  1  0  0
							 0  0  0  0  0  0  0  1  0]
           @test t.commands == ["initialise(4)"
								 "hadamard(2)"
								 "phase(1)"
								 "cnot(2,1)"
								 "hadamard(3)"
								 "phase(3)"
								 "phase(3)"
								 "hadamard(3)"
								 "phase(4)"
								 "phase(4)"
								 "phase(3)"
								 "phase(3)" ]
			reinitialise(t)
			@test t.commands == ["initialise(4)"]

    end;

    @testset "Basic Symplectic operations" begin

 			t = cliffordToTableau(3,3472,3)
 			@test t.state == [0  0  0  0  1  0  1
						 0  0  0  1  0  1  1
						 1  0  1  1  0  1  0
						 0  1  0  1  0  1  0
						 1  0  0  0  1  1  0
						 0  0  0  0  0  1  0]
			hadamard(t,3)
			cnot(t,2,1)
			cnot(t,1,2)
			cnot(t,1,3)
			@test tableauToClifford(t) == 1280169
			@test decomposeState(t)
			@test t.commands == ["initialise(3)"
									"cnot(3,2)"
									"phase(2)"
									"phase(2)"
									"phase(2)"
									"phase(3)"
									"phase(3)"
									"cnot(3,2)"
									"phase(3)"
									"phase(3)"
									"phase(3)"
									"phase(2)"
									"phase(2)"
									"phase(2)"
									"phase(1)"
									"phase(1)"
									"hadamard(3)"
									"hadamard(2)"
									"hadamard(1)"
									"cnot(2,1)"
									"phase(1)"
									"cnot(2,1)"
									"phase(2)"
									"phase(2)"
									"phase(2)"
									"phase(1)"
									"phase(1)"
									"phase(1)"
									"cnot(2,1)"]
			@test drawCircuit(t) == (Any[("CNOT", "q2", "q3"), ("P", "q2"), ("P", "q2"), ("P", "q2"), ("P", "q3"), ("P", "q3"), ("CNOT", "q2", "q3"), ("P", "q3"), ("P", "q3"), ("P", "q3"), ("P", "q2"), ("P", "q2"), ("P", "q2"), ("P", "q1"), ("P", "q1"), ("H", "q3"), ("H", "q2"), ("H", "q1"), ("CNOT", "q1", "q2"), ("P", "q1"), ("CNOT", "q1", "q2"), ("P", "q2"), ("P", "q2"), ("P", "q2"), ("P", "q1"), ("P", "q1"), ("P", "q1"), ("CNOT", "q1", "q2")], Any["q1", "q2", "q3"])
			@test qiskitCircuit(t) == "#Pass in the circuit and qubit register\ndef createCircuit(circs,qreg):\n    circs.cx(qreg[3],qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[3])\n    circs.s(qreg[3])\n    circs.cx(qreg[3],qreg[2])\n    circs.s(qreg[3])\n    circs.s(qreg[3])\n    circs.s(qreg[3])\n    circs.s(qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[1])\n    circs.s(qreg[1])\n    circs.h(qreg[3])\n    circs.h(qreg[2])\n    circs.h(qreg[1])\n    circs.cx(qreg[2],qreg[1])\n    circs.s(qreg[1])\n    circs.cx(qreg[2],qreg[1])\n    circs.s(qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[2])\n    circs.s(qreg[1])\n    circs.s(qreg[1])\n    circs.s(qreg[1])\n    circs.cx(qreg[2],qreg[1])"
			@test getNumberOfCliffords(3) == 92897280
			@test getNumberOfSymplecticCliffords(4) == 47377612800
			@test getNumberOfBitStringsCliffords(5) == 1024
			@test vec(generateRawCliffords(hadamardGate=1/sqrt(3).*([1 -1;-1 1]))[3] ) == vec([0.5773502692 + 0.0im 0.5773502692 + 0.0im 0.5773502692 + 0.0im  0.5773502692 + 0.0im])
			state = makeFromCommand(t.commands)

			end;

	@testset "Ket tests" begin
	# These tests are a bit 'output specific' but
		t = Tableau(3)
		phase(t,1)
		hadamard(t,1)
		phase(t,1)
		cnot(t,1,2)
		cnot(t,2,3)
		hadamard(t,3)
		phase(t,3)
		@test kets(t) == "+|000> \n+i|110> \n+i|001> \n+|111> \n"
		t = Tableau(3)
		phase(t,1)
		hadamard(t,1)
		phase(t,1)
		cnot(t,1,2)
		cnot(t,2,3)
		hadamard(t,3)
		phase(t,3)
		hadamard(t,2)
		phase(t,2)
		phase(t,2)
		hadamard(t,2)
		phase(t,1)
		phase(t,2)
		cnot(t,2,3)
		cnot(t,1,2)
		@test kets(t) == "+|010> \n-i|111> \n-i|011> \n+|110> \n"
		t = cliffordToTableau(3,34681,400)
		@test kets(t) == "+|010> \n-i|101> \n-i|001> \n+|110> \n"
	end

	@testset "Measure tests" begin
		for i in 1:5
			t = Tableau(3)
			phase(t,1)
			hadamard(t,1)
			phase(t,1)
			cnot(t,1,2)
			cnot(t,2,3)
			hadamard(t,3)
			phase(t,3)
			x = measure(t,2)
			if x == 0
				@test kets(t) == "+|000> \n+i|001> \n"
			else
				@test kets(t) == "+|110> \n-i|111> \n"
			end
		end
		for i in 1:5
			t = cliffordToTableau(3,34681,400)
			x = measure(t,1)
			if x == 1
				@test kets(t) == "+|110> \n-i|101> \n"
			else
				@test kets(t) == "+|010> \n-i|001> \n"
			end
		end
		for i in 1:10
			t = cliffordToTableau(4,434681,400)
			z_m = Int32[1 0 1 1 0 0 1 1 1; 0 0 1 0 0 1 0 1 1; 1 1 0 0 1 1 1 0 0; 0 0 1 0 1 0 0 0 0; 0 0 0 0 0 0 0 1 0; 0 1 1 0 0 0 0 1 1; 1 1 1 0 1 1 1 1 0; 1 0 0 0 1 0 0 1 1]
			o_m = Int32[1 0 1 1 0 0 1 1 1; 0 0 1 0 0 1 0 1 1; 1 1 0 0 1 1 1 0 0; 0 0 1 0 1 0 0 0 0; 0 0 0 0 0 0 0 1 1; 0 1 1 0 0 0 0 1 1; 1 1 1 0 1 1 1 1 0; 1 0 0 0 1 0 0 1 1]
			x = measure(t,4)
			if x == 0
				@test t.state == z_m
			else
				@test t.state == o_m
			y = measure(t,3)
			if x==1 && y == 0
				@test kets(t) == "+|0001> \n+i|1001> \n"
			elseif x==0 && y == 0
				@test kets(t) == "+|0100> \n-i|1100> \n"
			elseif x==0 && y == 1
				@test kets(t) == "+|0010> \n-i|1010> \n"
			else
				@test kets(t) == "+|0111> \n+i|1111> \n"
			end
		end
		for i in 1:5
			t = cliffordToTableau(3,401,4)
			x = measure(t,1)
			@test x == 0
			@test kets(t) == "+|000> \n+i|010> \n+i|001> \n-|011> \n"
		end
		for i in 1:5
			t = cliffordToTableau(3,401,4)
			hadamard(t,1)
			phase(t,1)
			phase(t,1)
			hadamard(t,1)
			x = measure(t,1)
			@test x == 1
			@test kets(t) == "+|100> \n+i|110> \n+i|101> \n-|111> \n"
		end
	end
end


text:
# Overflow checks

function all_minus_inf(x)
    for y in x
        if y > typemin(y)
            return false
        end
    end
    return true
end

function all_plus_inf(x)
    for y in x
        if y < typemax(y)
            return false
        end
    end
    return true
end

function all_zeros(x)
    for y in x
        if 1 / abs(x) < typemax(x)
            return false
        end
    end
    return true
end

function all_nan(x)
    for y in x
        if !isnan(y)
            return false
        end
    end
    return true
end

ext: jl


text:
using Flux
using NeuralVerification: ReLU, Id, Layer, network

@testset "Flux" begin

    DI = Dense(5,5)
    DR = Dense(5,5, relu)
    LI = Layer(rand(5,5), rand(5), Id())
    LR = Layer(rand(5,5), rand(5), ReLU())

    # just test that they don't error:
    @test network(Chain(DR,DR,DI)) isa Network
    @test Chain(Network([LR,LR,LI])) isa Chain

end


text:
@testset "refco" begin
    phpa = 800.0
    tc = 10.0
    rh = 0.9
    wl = 0.4
    refa, refb = refco(phpa, tc, rh, wl)
    @test refa ≈ 0.2264949956241415009e-3 atol=1e-15
    @test refb ≈ -0.2598658261729343970e-6 atol=1e-18
end

@testset "rm2v" begin
    w = rm2v([0.0 -0.8 -0.6;
                   0.8 -0.36 0.48;
                   0.6 0.48 -0.64])
    @test w[1] ≈ 0.0 atol=1e-12
    @test w[2] ≈ 1.413716694115406957 atol=1e-12
    @test w[3] ≈ -1.884955592153875943 atol=1e-12
    @test_throws ArgumentError rm2v([0.0 -0.8 -0.6])
end

@testset "rv2m" begin
    r = rv2m([0.0, 1.41371669, -1.88495559])
    @test r[1,1] ≈ -0.7071067782221119905 atol=1e-14
    @test r[1,2] ≈ -0.5656854276809129651 atol=1e-14
    @test r[1,3] ≈ -0.4242640700104211225 atol=1e-14
    @test r[2,1] ≈ 0.5656854276809129651 atol=1e-14
    @test r[2,2] ≈ -0.0925483394532274246 atol=1e-14
    @test r[2,3] ≈ -0.8194112531408833269 atol=1e-14
    @test r[3,1] ≈ 0.4242640700104211225 atol=1e-14
    @test r[3,2] ≈ -0.8194112531408833269 atol=1e-14
    @test r[3,3] ≈ 0.3854415612311154341 atol=1e-14
    @test_throws ArgumentError rv2m([1.41371669, -1.88495559])
end

@testset "rx" begin
    phi = 0.3456789
    r = [2.0 3.0 2.0;
         3.0 2.0 3.0;
         3.0 4.0 5.0]
    r = rx(phi, r)
    @test r[1,1] ≈ 2.0 atol=0.0
    @test r[1,2] ≈ 3.0 atol=0.0
    @test r[1,3] ≈ 2.0 atol=0.0
    @test r[2,1] ≈ 3.839043388235612460 atol=1e-12
    @test r[2,2] ≈ 3.237033249594111899 atol=1e-12
    @test r[2,3] ≈ 4.516714379005982719 atol=1e-12
    @test r[3,1] ≈ 1.806030415924501684 atol=1e-12
    @test r[3,2] ≈ 3.085711545336372503 atol=1e-12
    @test r[3,3] ≈ 3.687721683977873065 atol=1e-12
    @test_throws ArgumentError rx(phi, r[1:2,:])
end

@testset "rxp" begin
    r = [2.0 3.0 2.0;
         3.0 2.0 3.0;
         3.0 4.0 5.0]
    p = [0.2,1.5,0.1]
    rp = ERFA._rxp(r, p)
    @test rp[1] ≈ 5.1 atol=1e-12
    @test rp[2] ≈ 3.9 atol=1e-12
    @test rp[3] ≈ 7.1 atol=1e-12
    rp_exp = r * p
    @testset for i in eachindex(rp, rp_exp)
        @test rp[i] ≈ rp_exp[i]
    end
    for i in eachindex(rp, rp_exp)
        @test rp[i] ≈ rp_exp[i]
    end
    @test_throws ArgumentError ERFA._rxp(r[1:2,:], p)
    @test_throws ArgumentError ERFA._rxp(r, p[1:2])
end

@testset "rxpv" begin
    r = [2.0 3.0 2.0;
         3.0 2.0 3.0;
         3.0 4.0 5.0]
    pv = [[0.2,1.5,0.1],
          [1.5,0.2,0.1]]
    rpv = ERFA._rxpv(r, pv)
    @test rpv[1][1] ≈ 5.1 atol=1e-12
    @test rpv[1][2] ≈ 3.9 atol=1e-12
    @test rpv[1][3] ≈ 7.1 atol=1e-12
    @test rpv[2][1] ≈ 3.8 atol=1e-12
    @test rpv[2][2] ≈ 5.2 atol=1e-12
    @test rpv[2][3] ≈ 5.8 atol=1e-12
    @test_throws ArgumentError ERFA._rxpv(r[1:2,:], pv)
    pve = [[1.5,0.1],
          [1.5,0.2,0.1]]
    @test_throws ArgumentError ERFA._rxpv(r, pve)
    rpj = [r * pv[1], r * pv[2]]
    for (v, vj) in zip(rpv, rpj)
        for i in eachindex(v, vj)
            @test v[i] ≈ vj[i]
        end
    end
end

@testset "rxr" begin
    a = [2.0 3.0 2.0;
         3.0 2.0 3.0;
         3.0 4.0 5.0]
    b = [1.0 2.0 2.0;
         4.0 1.0 1.0;
         3.0 0.0 1.0]
    atb = ERFA._rxr(a, b)
    @test atb[1,1] ≈ 20.0 atol=1e-12
    @test atb[1,2] ≈ 7.0 atol=1e-12
    @test atb[1,3] ≈ 9.0 atol=1e-12
    @test atb[2,1] ≈ 20.0 atol=1e-12
    @test atb[2,2] ≈ 8.0 atol=1e-12
    @test atb[2,3] ≈ 11.0 atol=1e-12
    @test atb[3,1] ≈ 34.0 atol=1e-12
    @test atb[3,2] ≈ 10.0 atol=1e-12
    @test atb[3,3] ≈ 15.0 atol=1e-12
    @test_throws ArgumentError ERFA._rxr(a[1:2,:], b)
    @test_throws ArgumentError ERFA._rxr(a, b[1:2,:])
    @test atb == a * b
end

@testset "ry" begin
    theta = 0.3456789
    r = [2.0 3.0 2.0;
         3.0 2.0 3.0;
         3.0 4.0 5.0]
    r = ry(theta, r)
    @test r[1,1] ≈ 0.8651847818978159930 atol=1e-12
    @test r[1,2] ≈ 1.467194920539316554 atol=1e-12
    @test r[1,3] ≈ 0.1875137911274457342 atol=1e-12
    @test r[2,1] ≈ 3 atol=1e-12
    @test r[2,2] ≈ 2 atol=1e-12
    @test r[2,3] ≈ 3 atol=1e-12
    @test r[3,1] ≈ 3.500207892850427330 atol=1e-12
    @test r[3,2] ≈ 4.779889022262298150 atol=1e-12
    @test r[3,3] ≈ 5.381899160903798712 atol=1e-12
    @test_throws ArgumentError ry(theta, r[1:2,:])
end

@testset "rz" begin
    psi = 0.3456789
    r = [2.0 3.0 2.0;
         3.0 2.0 3.0;
         3.0 4.0 5.0]
    r = rz(psi, r)
    @test r[1,1] ≈ 2.898197754208926769 atol=1e-12
    @test r[1,2] ≈ 3.500207892850427330 atol=1e-12
    @test r[1,3] ≈ 2.898197754208926769 atol=1e-12
    @test r[2,1] ≈ 2.144865911309686813 atol=1e-12
    @test r[2,2] ≈ 0.865184781897815993 atol=1e-12
    @test r[2,3] ≈ 2.144865911309686813 atol=1e-12
    @test r[3,1] ≈ 3.0 atol=1e-12
    @test r[3,2] ≈ 4.0 atol=1e-12
    @test r[3,3] ≈ 5.0 atol=1e-12
    @test_throws ArgumentError rz(psi, r[1:2,:])
end


text:
export meta, show

"""
    meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet)::Dict{Any,Any}

Gets metadata about a spreadsheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
"""
function meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet)::Dict{Any,Any}
    return gsheet_api_spreadsheet_get(client; spreadsheetId=spreadsheet.id)
end

"""
    meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet, key::AbstractString, value)::Dict{Any,Any}

Gets metadata about a spreadsheet key-value pair.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
- `key:AbstractString`: key
- `value`: value
"""
function meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet, key::AbstractString, value)::Dict{Any,Any}
    metadata = meta(client, spreadsheet)
    sheets = metadata["sheets"]

    for sheet in sheets
        properties = sheet["properties"]

        if properties[key] == value
            return properties
        end
    end

    throw(KeyError(title))
end

"""
    meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet, title::AbstractString)::Dict{Any,Any}

Gets metadata about a spreadsheet sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
- `title::AbstractString`: sheet title
"""
meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet, title::AbstractString)::Dict{Any,Any} = meta(client, spreadsheet, "title", title)

"""
    meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet, sheet_id::Int64)::Dict{Any,Any}

Gets metadata about a spreadsheet sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
- `sheet_id::Int64`: sheet id
"""
meta(client::GoogleSheetsClient, spreadsheet::Spreadsheet, sheet_id::Int64)::Dict{Any,Any} = meta(client, spreadsheet, "sheetId", sheet_id)

"""
    meta(client::GoogleSheetsClient, sheet::Sheet)::Dict{Any,Any}

Gets metadata about a spreadsheet sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `sheet::Sheet`: sheet
"""
meta(client::GoogleSheetsClient, sheet::Sheet)::Dict{Any,Any} = meta(client, sheet.spreadsheet, sheet.id)

"""
    show(client::GoogleSheetsClient, spreadsheet::Spreadsheet)

Prints metadata about a spreadsheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
"""
function Base.show(client::GoogleSheetsClient, spreadsheet::Spreadsheet)
    m = meta(client, spreadsheet)
    println("Spreadsheet:")
    println(json(m, 4))
end

"""
    show(client::GoogleSheetsClient, spreadsheet::Spreadsheet, title::AbstractString)

Prints metadata about a spreadsheet sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
- `title::AbstractString`: sheet title
"""
function Base.show(client::GoogleSheetsClient, spreadsheet::Spreadsheet, title::AbstractString)
    m = meta(client, spreadsheet, title)
    println("Sheet:")
    println(json(m, 4))
end

"""
    show(client::GoogleSheetsClient, spreadsheet::Spreadsheet, sheet_id::Int64)

Prints metadata about a spreadsheet sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `spreadsheet::Spreadsheet`: spreadsheet
- `sheet_id::Int64`: sheet id
"""
function Base.show(client::GoogleSheetsClient, spreadsheet::Spreadsheet, sheet_id::Int64)
    m = meta(client, spreadsheet, sheet_id)
    println("Sheet:")
    println(json(m, 4))
end

"""
    show(client::GoogleSheetsClient, sheet::Sheet)

Prints metadata about a spreadsheet sheet.

# Arguments
- `client::GoogleSheetsClient`: client
- `sheet::Sheet`: sheet
"""
Base.show(client::GoogleSheetsClient, sheet::Sheet) = Base.show(client, sheet.spreadsheet, sheet.id)


text:
"""
    RollingBuffer{T}(n)
The RollingBuffer type implements a circular buffer of fixed capacity,
similar to CircularBuffer from DataStructures.jl, but with growing indexes.

New items are pushed to the back of the list, overwriting values
in a circular fashion, and increasing starting index by one.

Allocate a buffer of elements of type `T` with maximum capacity `n`.
"""
mutable struct RollingBuffer{T} <: AbstractVector{T}
    buffer::Vector{T}
    capacity::Int
    length::Int
    offset::Int # change to delta for faster _buffer_index? delta = buf.offset + buf.first - 1
    first::Int
    last::Int
    final::Int # !specific - number of elements that are already processed in some sense
end

function RollingBuffer{T}(maxlen::Int, offset::Int = 0) where T
    RollingBuffer{T}(Vector{T}(undef, maxlen), maxlen, 0, offset, 1, 0, -1)
end

RollingBuffer(capacity) = RollingBuffer{Any}(capacity)

function RollingBuffer(vec::Vector{T}, offset::Int = 0) where T
    capacity = length(vec)
    RollingBuffer{T}(vec, capacity, capacity, offset, 1, capacity, -1)
end

# no copy for native vector type
RollingBuffer(v::Vector{T}, inds::AbstractUnitRange) where {T} =
    RollingBuffer(v, first(inds)-1)

RollingBuffer(v::AbstractVector{T}, inds::AbstractUnitRange) where {T} =
    RollingBuffer(collect(v), first(inds)-1)

RollingBuffer(v::AbstractVector) = RollingBuffer(v, Base.axes1(v))

Base.@propagate_inbounds function _buffer_index_checked(buf::RollingBuffer, i::Int)
    @boundscheck if i < firstindex(buf) || i > lastindex(buf)
        throw(BoundsError(buf, i))
    end
    _buffer_index(buf, i)
end

@inline function _buffer_index(buf::RollingBuffer, i::Int)
    idx = i - buf.offset + buf.first - 1     # idx = mod1(buf.first + i - 1, buf.capacity) # idx = (buf.first + i - 2) % buf.capacity + 1
    if idx > buf.capacity
        idx -= buf.capacity
    end
    idx #idx = buf.first + i - 1 #idx > buf.capacity ? idx - buf.capacity : idx
end

@inline Base.@propagate_inbounds function Base.getindex(buf::RollingBuffer, i::Int)
    buf.buffer[_buffer_index_checked(buf, i)]
end

@inline Base.@propagate_inbounds function Base.setindex!(buf::RollingBuffer, data, i::Int)
    buf.buffer[_buffer_index_checked(buf, i)] = data
    buf
end

# pop_back
@inline function Base.pop!(buf::RollingBuffer)
    i = buf.last
    if buf.length > 0
        buf.length -= 1
        buf.last = (buf.last <= 1 ? buf.capacity : buf.last - 1)
    else
        throw(ArgumentError("array must be non-empty"))
    end
    buf.buffer[i]
end
@inline function Base.pop!(buf::RollingBuffer, Npop::Integer)
    if Npop > buf.length
        Npop = buf.length
    end
    if buf.length > 0
        buf.length -= Npop
        buf.last = (buf.last > Npop ? buf.last - Npop : buf.last - Npop + buf.capacity)
    end
    buf
end

# pop_front
function Base.popfirst!(buf::RollingBuffer)
    i = buf.first
    if buf.length > 0
        buf.length -= 1
        buf.first = (buf.first == buf.capacity ? 1 : buf.first + 1)
        buf.offset += 1
    else
        throw(ArgumentError("array must be non-empty"))
    end
    buf.buffer[i]
end
function Base.popfirst!(buf::RollingBuffer, Npop::Integer)
    if Npop > buf.length
        Npop = buf.length
    end
    if buf.length > 0
        buf.length -= Npop
        buf.first += Npop
        buf.first = (buf.first > buf.capacity ? buf.first - buf.capacity : buf.first)
        buf.offset += Npop
    end
    buf
end

# push_back
@inline function Base.push!(buf::RollingBuffer, data)
    if buf.length == buf.capacity # increase offset and first point, if full
        buf.first = (buf.first == buf.capacity ? 1 : buf.first + 1)
        buf.offset += 1
    else # increase capacity, if not full
        buf.length += 1
    end
    buf.last = (buf.last == buf.capacity ? 1 : buf.last + 1)
    buf.buffer[buf.last] = data
    buf
end

# push_front - additional function
function Base.pushfirst!(buf::RollingBuffer, data)
    if buf.length == buf.capacity
        buf.last = (buf.last <= 1 ? buf.capacity : buf.last - 1)
    else
        buf.length += 1
    end
    buf.offset -= 1
    buf.first = (buf.first <= 1 ? buf.capacity : buf.first - 1)
    buf.buffer[buf.first] = data
    buf
end

# push_back vector
@inline Base.push!(buf::RollingBuffer, datavec::AbstractVector) = Base.append!(buf::RollingBuffer, datavec::AbstractVector)

@inbounds function Base.append!(buf::RollingBuffer, datavec::AbstractVector)
    Nadd = length(datavec)
    Nadd > buf.capacity && error("Added length larger than buffer length!")

    lenRight = buf.capacity - buf.last
    bshift = buf.last # dst = pointer(buf.buffer, buf.last + 1)
    shift = 0 #src = pointer(datavec, 1)
    if lenRight < Nadd # there is not enough space at the end of the buffer
        copyto!(buf.buffer, bshift+1, datavec, shift+1, lenRight) # unsafe_copyto!(dst, src, lenRight)
        buf.last = Nadd - lenRight
        bshift = 0 #dst = pointer(buf.buffer, 1)
        shift = lenRight #src = pointer(datavec, lenRight + 1)
        copyto!(buf.buffer, bshift+1, datavec, shift+1, buf.last) #unsafe_copyto!(dst, src, buf.last)
    else # buf.last < buf.capacity
        copyto!(buf.buffer, bshift+1, datavec, shift+1, Nadd) # unsafe_copyto!(dst, src, Nadd)
        buf.last += Nadd
    end

    if buf.length == buf.capacity # most frequent case = buffer is full
        buf.first += Nadd #_buffer_index(buf, buf.first + Nadd)
        buf.first = (buf.first > buf.capacity ? buf.first - buf.capacity : buf.first)
        buf.offset += Nadd
    else
        if buf.capacity < buf.length + Nadd
            Nadd -= buf.capacity - buf.length
            buf.first += Nadd
            buf.first = (buf.first > buf.capacity ? buf.first - buf.capacity : buf.first)
            buf.offset += Nadd
            buf.length = buf.capacity;
        else
            buf.length += Nadd;
        end
    end
    buf
end

# Preparing a place for the next append of a fixed size without actual append.
# One should assign `buf[nextindex:end]` after this call!
@inbounds function prepare_append!(buf::RollingBuffer, Nadd::Int)
    Nadd > buf.capacity && error("Added length larger than buffer length!")
    nextindex = buf.offset + buf.length + 1
    lenRight = buf.capacity - buf.last
    if lenRight < Nadd # there is not enough space at the end of the buffer
        buf.last = Nadd - lenRight
    else # buf.last < buf.capacity
        buf.last += Nadd
    end
    if buf.length == buf.capacity # most frequent case = buffer is full
        buf.first += Nadd #_buffer_index(buf, buf.first + Nadd)
        buf.first = (buf.first > buf.capacity ? buf.first - buf.capacity : buf.first)
        buf.offset += Nadd
    else
        if buf.capacity < buf.length + Nadd
            Nadd -= buf.capacity - buf.length
            buf.first += Nadd
            buf.first = (buf.first > buf.capacity ? buf.first - buf.capacity : buf.first)
            buf.offset += Nadd
            buf.length = buf.capacity;
        else
            buf.length += Nadd;
        end
    end
    nextindex
end

function Base.empty!(buf::RollingBuffer)
    buf.first = 1
    buf.last = buf.length = buf.offset = 0
    buf
end

function Base.fill!(cb::RollingBuffer, data)
    for i in 1:capacity(cb)-length(cb)
        push!(cb, data)
    end
    cb
end

Base.parent(A::RollingBuffer) = A.buffer

Base.@propagate_inbounds function Base.first(buf::RollingBuffer)
    @boundscheck if buf.length == 0
        throw(BoundsError(buf, 1))
    end
    buf.buffer[buf.first]
end
Base.@propagate_inbounds function Base.last(buf::RollingBuffer)
    @boundscheck if buf.length == 0
        throw(BoundsError(buf, 1))
    end
    buf.buffer[buf.last]
end

Base.length(buf::RollingBuffer) = buf.length
Base.size(buf::RollingBuffer) = (length(buf),)
Base.isempty(buf::RollingBuffer) = buf.length == 0
#Base.convert(::Type{Array}, buf::RollingBuffer{T}) where T = T[x for x in buf]
# convert buffer to a 1-based indexed array, dropping its offset
function Base.convert(::Type{Array}, buf::RollingBuffer{T}) where {T}
    len = length(buf)
    vec = Vector{T}(undef, len)
    for (i, j) in zip(1:len, eachindex(buf))
        @inbounds vec[i] = buf[j]
    end
    vec
end
# copy buffer to another array, dropping its offset
function Base.copyto!(dest::Vector, src::RollingBuffer)
    if length(dest) != length(src)
        if length(dest) < length(src)
            throw(BoundsError(dest, length(src)))
        else
            throw(BoundsError(src, firstindex(src)-1+length(dest)))
        end
    end

    for (i, j) in zip(eachindex(dest), eachindex(src))
        @inbounds dest[i] = src[j]
    end
    dest
end

Base.lastindex(buf::RollingBuffer) = buf.offset + buf.length #count
Base.firstindex(buf::RollingBuffer) = buf.offset + 1 #countfirst

# if final is not used, it is equal to the last element
finalindex(buf::RollingBuffer) = buf.final < 0 ? lastindex(buf) : buf.final # !specific
setfinalindex(buf::RollingBuffer, i::Int) = buf.final = i # !specific
addfinalindex(buf::RollingBuffer, i::Int) = buf.final = (buf.final < 0) ? i : buf.final + i # !specific
final(buf::RollingBuffer) = buf[finalindex(buf)]

# specific
offset(buf::RollingBuffer) = buf.offset
capacity(buf::RollingBuffer) = buf.capacity # maximum available length
isfull(buf::RollingBuffer) = buf.length == buf.capacity
#eachindex bounds(buf::RollingBuffer) = (buf.first, buf.last) # bound indexes

# clear and resize
function reset!(buf::RollingBuffer, capacity::Int = buf.capacity, offset::Int = 0) # clear and resize
    empty!(buf)
    buf.offset = offset
    if buf.capacity != capacity
        buf.capacity = capacity
        buf.buffer = Vector{T}(undef, capacity)
    end
end

# custom indexing
# @inline Base.axes(A::RollingBuffer) =  (_slice(A.offset + 1, A.length),)
#@inline _slice(start, capacity) = Base.IdentityUnitRange(Base._range(start, nothing, nothing, capacity))
Base.axes1(buf::RollingBuffer) = URange(1 + buf.offset, buf.length + buf.offset)
Base.axes(buf::RollingBuffer) = (Base.axes1(buf),)

function Base.similar(v::RollingBuffer, T::Type, inds::Tuple{URange})
    inds1 = inds[1]
    n = length(inds1)
    RollingBuffer(Array{T}(undef, n), first(inds1)-1)
end

function Base.similar(f::Union{Function,Type}, inds::Tuple{URange})
    inds1 = inds[1]
    n = length(inds1)
    RollingBuffer(f(Base.OneTo(n)), first(inds1)-1)
end

#=
function Base.similar(::Type{T}, inds::Tuple{URange}) where {T<:AbstractArray}
    inds1 = inds[1]
    n = length(inds1)
    RollingBuffer(Array{T}(undef, n), first(inds1)-1)
end
=#


text:
# coil_compress.jl

using MIRT: ir_mri_coil_compress
using MIRTjim: jim, prompt
using MIRT: ir_mri_sensemap_sim, ir_load_brainweb_t1_256
import MIRT: snr2sigma

using LinearAlgebra: svd, norm
using Plots: plot, plot!, scatter, scatter!
using Test: @inferred

xtrue = ir_load_brainweb_t1_256()
(nx, ny) = size(xtrue)
ncoil = 8
smap = ir_mri_sensemap_sim(dims = (nx, ny), dx=192/nx, ncoil=ncoil, rcoil=100)
idata = smap .* xtrue
σ = snr2sigma(50, idata)
idata = idata + σ * randn(ComplexF32, size(idata))
nkeep = 4
#@inferred ir_mri_coil_compress(idata, ncoil = nkeep) # todo: fails
(odata, σ, _) = ir_mri_coil_compress(idata, ncoil = nkeep)
for iz=1:nkeep # normalize for display
	odata[:,:,iz] = odata[:,:,iz] / maximum(abs.(odata[:,:,iz]))
end
p = scatter(1:nkeep, σ[1:nkeep], marker=:circle, label="")
scatter!(p, (nkeep+1):ncoil, σ[nkeep+1:ncoil], marker=:x, label="")
tmp = round((norm(σ[1:nkeep]) / norm(σ))^2*100, digits=2)
plot!(p, title="percent kept $tmp")
plot(jim(smap, "smap"), jim(idata, "idata"), jim(odata, "odata"), p)


text:
push!(LOAD_PATH, "../src/")

using Documenter, EarthEngineREST

pages = ["Home" => "index.md", "API" => "api.md"]

makedocs(;
    modules = [EarthEngineREST],
    authors = "Kel Markert",
    repo = "https://github.com/KMarkert/EarthEngineREST.jl/blob/{commit}{path}#L{line}",
    sitename = "EarthEngineREST.jl",
    pages = pages,
)

deploydocs(;
    repo = "github.com/KMarkert/EarthEngineREST.jl.git",
    devbranch = "main",
)


text:
@testset "test that synonyms work" begin
    @test precision_score === positive_predictive_value
    @test sensitivity === true_positive_rate
    @test recall === true_positive_rate
    @test specificity === true_negative_rate
end

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

y_true_p = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]
y_true_m = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1]
y_true_i = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2]
y_true_b = Vector{Bool}(y_true_p)
y_true_pf = Vector{Float64}(y_true_p)
y_true_mf = Vector{Float64}(y_true_m)
y_true_sym = convertlabel([:pos,:neg], y_true_p, LabelEnc.ZeroOne())
y_true_sym2 = [:a, :a, :a, :a, :a, :a, :a, :a, :a, :a, :b, :b, :b, :b, :c, :c, :c]
y_true_c = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3]

y_hat_p  = [1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0]
y_hat_m  = [1, 1,-1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, 1,-1]
y_hat_i  = [1, 1, 2, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2]
y_hat_b  = Vector{Bool}(y_hat_p)
y_hat_pf = Vector{Float64}(y_hat_p)
y_hat_mf = Vector{Float64}(y_hat_m)
y_hat_sym  = convertlabel([:pos,:neg], y_hat_p,  LabelEnc.ZeroOne())
y_hat_sym2  = [:a, :a, :b, :c, :a, :b, :a, :c, :c, :a, :c, :a, :c, :a, :b, :a, :c]
y_hat_c  = [1, 1, 2, 3, 1, 2, 1, 3, 3, 1, 3, 1, 3, 1, 2, 1, 3]

_accuracy_nonorm(args...) = accuracy(args...; normalize=false)
# (fun, binary, micro, macro)
const fun_refs = [
    (prevalence, 10/17, nothing, nothing),
    (f_score, 0.5555555555555556, 0.35294117647058826, 0.25925925925925924),
    (f1_score, 0.5555555555555556, 0.35294117647058826, 0.25925925925925924),
    (accuracy, 9/17, nothing, nothing),
    (_accuracy_nonorm, 9., nothing, nothing),
    (positive_predictive_value, 5/8, 0.35294117647058826, 0.2638888888888889),
    (false_discovery_rate, 3/8, nothing, nothing),
    (negative_predictive_value, 4/9, nothing, nothing),
    (false_omission_rate, 5/9, nothing, nothing),
    (true_positive_rate, 5/10, 0.35294117647058826, 0.27777777777777773),
    (false_positive_rate, 3/7, nothing, nothing),
    (false_negative_rate, 5/10, nothing, nothing),
    (true_negative_rate, 4/7, nothing, nothing),
    (positive_likelihood_ratio, 5/10 * 7/3, nothing, nothing),
    (negative_likelihood_ratio, 5/10 * 7/4, nothing, nothing),
    (diagnostic_odds_ratio, (5/10 * 7/3) / (5/10 * 7/4), nothing, nothing)
]

@testset "fuzzy with/without encoding" begin
    targets = [y_true_p, y_true_m, y_true_b, y_true_pf, y_true_mf]
    outputs = [y_hat_p, y_hat_m, y_hat_b, y_hat_pf, y_hat_mf, y_hat_b']
    for (fun, ref_binary, ref_micro, ref_macro) in fun_refs
        @testset "$fun: binary" begin
            for target in targets, output in outputs
                @testset "$(typeof(target)) against $(typeof(output))" begin
                    @test @inferred(fun(target, output, LabelEnc.FuzzyBinary())) === ref_binary
                    if eltype(target) <: Bool || eltype(output) <: Bool
                        @test @inferred(fun(target, output)) === ref_binary
                    end
                end
            end
        end
    end
end

@testset "arrays with/without encoding" begin
    binary_data = [
        (y_true_p,    y_hat_p,    LabelEnc.ZeroOne()),
        (y_true_p,    y_hat_pf,   LabelEnc.ZeroOne()),
        (y_true_pf,   y_hat_p,    LabelEnc.ZeroOne()),
        (y_true_pf,   y_hat_pf,   LabelEnc.ZeroOne()),
        (y_true_m,    y_hat_m,    LabelEnc.MarginBased()),
        (y_true_m,    y_hat_mf,   LabelEnc.MarginBased()),
        (y_true_mf,   y_hat_m,    LabelEnc.MarginBased()),
        (y_true_mf,   y_hat_mf,   LabelEnc.MarginBased()),
        (y_true_i,    y_hat_i,    LabelEnc.Indices(2)),
        (y_true_b,    y_hat_b,    LabelEnc.TrueFalse()),
        (y_true_sym,  y_hat_sym,  LabelEnc.NativeLabels([:pos,:neg])),
        (y_true_sym,  y_hat_sym,  LabelEnc.OneVsRest(:pos)),
        (y_true_sym2, y_hat_sym2, LabelEnc.OneVsRest(:a)),
    ]
    multiclass_data = [
        (y_true_sym2,    y_hat_sym2,  LabelEnc.NativeLabels([:a,:b,:c])),
        (y_true_c,          y_hat_c,  LabelEnc.Indices(3)),
        (Int32.(y_true_c),  y_hat_c,  LabelEnc.Indices(3)),
        (y_true_c, Float32.(y_hat_c), LabelEnc.Indices(3)),
    ]
    for (fun, ref_binary, ref_micro, ref_macro) in fun_refs
        @testset "$fun: binary" begin
            for (target, output, enc) in binary_data
                @testset "$(eltype(target)) against $(eltype(output)) (enc: $(typeof(enc)))" begin
                    @test @inferred(fun(target, output, enc)) === ref_binary
                    if nlabel(target) == 2
                        @test fun(target, output) === ref_binary
                    end
                end
            end
            # this is not type stable (number of labels not static)
            @test fun(y_true_sym, y_hat_sym, [:pos, :neg]) === ref_binary
        end
        @testset "$fun: multiclass" begin
            for (target, output, enc) in multiclass_data
                @testset "$(eltype(target)) against $(eltype(output)) (enc: $(typeof(enc)))" begin
                    if fun == _accuracy_nonorm
                        @test @inferred(fun(target, output, enc)) == 6
                        @test fun(target, output, label(enc)) == 6
                        @test fun(target, output) == 6
                    elseif fun == accuracy
                        @test @inferred(fun(target, output, enc)) == 6/17
                        @test fun(target, output, label(enc)) == 6/17
                        @test fun(target, output) == 6/17
                        @test fun(target, output, normalize=true) == 6/17
                    else
                        res = Dict(map(label(enc)) do l
                            l => fun(target, output, LabelEnc.OneVsRest(l))
                        end)
                        @test @inferred(fun(target, output, enc)) == res
                        @test @inferred(fun(target, output, enc, AvgMode.None())) == res
                        @test fun(target, output, enc, avgmode=:none) == res
                        @test fun(target, output) == res
                        @test fun(target, output, avgmode=:none) == res
                        if ref_micro != nothing
                            @test @inferred(fun(target, output, enc, AvgMode.Micro())) ≈ ref_micro
                            @test fun(target, output, label(enc), AvgMode.Micro()) ≈ ref_micro
                            @test fun(target, output, AvgMode.Micro()) ≈ ref_micro
                            @test fun(target, output, enc, avgmode=:micro) ≈ ref_micro
                            @test fun(target, output, label(enc), avgmode=:micro) ≈ ref_micro
                            @test fun(target, output, avgmode=:micro) ≈ ref_micro
                        end
                        if ref_macro != nothing
                            @test @inferred(fun(target, output, enc, AvgMode.Macro())) ≈ ref_macro
                            @test fun(target, output, label(enc), AvgMode.Macro()) ≈ ref_macro
                            @test fun(target, output, AvgMode.Macro()) ≈ ref_macro
                            @test fun(target, output, enc, avgmode=:macro) ≈ ref_macro
                            @test fun(target, output, label(enc), avgmode=:macro) ≈ ref_macro
                            @test fun(target, output, avgmode=:macro) ≈ ref_macro
                        end
                    end
                end
            end
        end
    end
end

@testset "f_score: extra tests for beta" begin
    y_true = [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2]
    y_hat  = [0, 1, 0, 1, 0, 1, 0, 1, 2, 2, 1]
    @test f_score(y_true, y_hat) == Dict{Int, Float64}(0 => 0.5714285714285715, 1 => 0.6, 2=>0.8)
    @test f_score(y_true, y_hat, 2.0) == Dict(0 => 0.35714285714285715, 2 => 0.5, 1 => 0.375)
    @test f_score(y_true, y_hat, unique(y_true), 2.0) == Dict(0 => 0.35714285714285715, 2 => 0.5, 1 => 0.375)
    @test f_score(y_true, y_hat, beta=2.0) == Dict(0 => 0.35714285714285715, 2 => 0.5, 1 => 0.375)
    @test f_score(y_true, y_hat, avgmode=:none, beta=2.0) == Dict(0 => 0.35714285714285715, 2 => 0.5, 1 => 0.375)
    @test f_score(y_true, y_hat, LabelEnc.NativeLabels(unique(y_true))) == Dict{Int, Float64}(0 => 0.5714285714285715, 1 => 0.6, 2=>0.8)
    @test f_score(y_true, y_hat, LabelEnc.NativeLabels(unique(y_true)), 1) == Dict{Int, Float64}(0 => 0.5714285714285715, 1 => 0.6, 2=>0.8)
    @test f_score(y_true, y_hat, AvgMode.Macro()) ≈ 0.6571428571428571
    @test f_score(y_true, y_hat, AvgMode.Macro(), 1) ≈ 0.6571428571428571
    @test f_score(y_true, y_hat, AvgMode.Micro()) ≈ 0.6363636363636364
    @test f_score(y_true, y_hat, avgmode=:micro) ≈ 0.6363636363636364
    @test f_score(y_true, y_hat, avgmode=:micro, beta=1) ≈ 0.6363636363636364
    @test f_score(y_true, y_hat, LabelEnc.NativeLabels(unique(y_true)), AvgMode.Micro()) ≈ 0.6363636363636364
    @test f_score(y_true, y_hat, LabelEnc.NativeLabels(unique(y_true)), AvgMode.Micro(), 1) ≈ 0.6363636363636364
    @test f_score(y_true, y_hat, LabelEnc.NativeLabels(unique(y_true)); avgmode=:micro, beta=1) ≈ 0.6363636363636364
    @test f_score(y_true, y_hat, unique(y_true); avgmode=:micro, beta=1) ≈ 0.6363636363636364
    @test sort(collect(keys(f_score(y_true, y_hat, [0, 1, 2])))) == [0, 1, 2]
end

#=
@testset "multiclass sanity check" begin
    # We count positive strickly positive matches as positives
    @test true_positives([1,2,3,4], [1,3,2,4]) === 2
    @test true_negatives([1,2,3,4], [4,3,2,1]) === 0
    @test accuracy([1,2,3,4], [1,3,2,4]) === .5
    @test accuracy([:a,:b,:b,:c], [:c,:b,:a,:a]) === .25
end
=#
#=
# matthews_corrcoef
@test matthews_corrcoef(y_true, y_hat) == 0.0
=#


text:
module Elfel

# Elfel (C) 2020, Petr Krysl

include("utilities.jl")
include("LocalAssemblers.jl")
include("Assemblers.jl")
include("RefShapes.jl")
include("FElements.jl")
include("FESpaces.jl")
include("QPIterators.jl")
include("FEIterators.jl")

# We can either use/import individual functions from Elfel like so:
# ```
# using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
# ```
# or we can bring into our context all exported symbols as
# ```
# using Elfel.Exports
# ```
include("Exports.jl")

end # module


text:
using Test

Test.@testset "estimations" begin

    include("test_emr_msm_estimate.jl")

end


text:
using OrdinaryDiffEq, Flux, DiffEqFlux, Plots
using KitBase
import KitML

# config
begin
    case = "homogeneous"
    maxTime = 3
    tlen = 16
    u0 = -5.0
    u1 = 5.0
    nu = 80
    nug = 0
    v0 = -5.0
    v1 = 5.0
    nv = 28
    nvg = 0
    w0 = -5.0
    w1 = 5.0
    nw = 28
    nwg = 0
    vMeshType = "rectangle"
    nm = 5
    knudsen = 1.0
    inK = 0
    alpha = 1.0
    omega = 0.5
    nh = 8
end

# dataset
begin
    tspan = (0.0, maxTime)
    tran = linspace(tspan[1], tspan[2], tlen)
    γ = heat_capacity_ratio(inK, 3)
    vSpace = VSpace3D(u0, u1, nu, v0, v1, nv, w0, w1, nw, vMeshType)

    f0 =
        Float32.(
            0.5 * (1 / π)^1.5 .*
            (exp.(-(vSpace.u .- 0.99) .^ 2) .+ exp.(-(vSpace.u .+ 0.99) .^ 2)) .*
            exp.(-vSpace.v .^ 2) .* exp.(-vSpace.w .^ 2),
        ) |> Array
    prim0 =
        conserve_prim(moments_conserve(f0, vSpace.u, vSpace.v, vSpace.w, vSpace.weights), γ)
    M0 = Float32.(maxwellian(vSpace.u, vSpace.v, vSpace.w, prim0)) |> Array

    mu_ref = ref_vhs_vis(knudsen, alpha, omega)
    kn_bzm = hs_boltz_kn(mu_ref, 1.0)
    τ0 = mu_ref * 2.0 * prim0[end]^(0.5) / prim0[1]

    phi, psi, phipsi = kernel_mode(
        nm,
        vSpace.u1,
        vSpace.v1,
        vSpace.w1,
        vSpace.du[1, 1, 1],
        vSpace.dv[1, 1, 1],
        vSpace.dw[1, 1, 1],
        vSpace.nu,
        vSpace.nv,
        vSpace.nw,
        alpha,
    )

    # Boltzmann
    prob = ODEProblem(boltzmann_ode!, f0, tspan, [kn_bzm, nm, phi, psi, phipsi])
    data_boltz = solve(prob, Tsit5(), saveat = tran) |> Array

    # BGK
    prob1 = ODEProblem(bgk_ode!, f0, tspan, [M0, τ0])
    data_bgk = solve(prob1, Tsit5(), saveat = tran) |> Array

    data_boltz_1D = zeros(Float32, axes(data_boltz, 1), axes(data_boltz, 4))
    data_bgk_1D = zeros(Float32, axes(data_bgk, 1), axes(data_bgk, 4))
    for j in axes(data_boltz_1D, 2)
        data_boltz_1D[:, j] .=
            reduce_distribution(data_boltz[:, :, :, j], vSpace.weights[1, :, :])
        data_bgk_1D[:, j] .=
            reduce_distribution(data_bgk[:, :, :, j], vSpace.weights[1, :, :])
    end
    f0_1D = reduce_distribution(f0, vSpace.weights[1, :, :])
    M0_1D = reduce_distribution(M0, vSpace.weights[1, :, :])

    X = Array{Float32}(undef, vSpace.nu, 1)
    for i in axes(X, 2)
        X[:, i] .= f0_1D
    end
    Y = Array{Float32}(undef, vSpace.nu, 1, tlen)
    for i in axes(Y, 2)
        Y[:, i, :] .= data_boltz_1D
    end
    M = Array{Float32}(undef, nu, size(X, 2))
    for i in axes(M, 2)
        M[:, i] .= M0_1D
    end
    τ = Array{Float32}(undef, 1, size(X, 2))
    for i in axes(τ, 2)
        τ[1, i] = τ0
    end
end

# neural model
begin
    model_univ = DiffEqFlux.FastChain(
        DiffEqFlux.FastDense(nu, nu * nh, tanh),
        DiffEqFlux.FastDense(nu * nh, nu),
    )
    p_model = DiffEqFlux.initial_params(model_univ)

    function dfdt(f, p, t)
        df = (M .- f) ./ τ .+ model_univ(M .- f, p)
    end
    prob_ube = ODEProblem(dfdt, X, tspan, p_model)

    function loss(p)
        sol_ube = solve(prob_ube, Midpoint(), u0 = X, p = p, saveat = tran)
        loss = sum(abs2, Array(sol_ube) .- Y)

        return loss
    end

    his = []
    cb = function (p, l)
        println("loss: $l")
        push!(his, l)
        return false
    end
end

# train
res = DiffEqFlux.sciml_train(loss, p_model, ADAM(), cb = cb, maxiters = 200)
res = DiffEqFlux.sciml_train(loss, res.minimizer, ADAM(), cb = cb, maxiters = 200)

# residual history
plot(log.(his))

# solution
ube = ODEProblem(KitML.ube_dfdt, f0_1D, tspan, [M0_1D, τ0, (model_univ, res.minimizer)]);
sol = solve(
    ube,
    Midpoint(),
    u0 = f0_1D,
    p = [M0_1D, τ0, (model_univ, res.minimizer)],
    saveat = tran,
);

# result
plot(
    vSpace.u[:, vSpace.nv÷2, vSpace.nw÷2],
    data_boltz_1D[:, 1],
    lw = 2,
    label = "Initial",
    color = :gray32,
    xlabel = "u",
    ylabel = "particle distribution",
)
plot!(
    vSpace.u[:, vSpace.nv÷2, vSpace.nw÷2],
    data_boltz_1D[:, 2],
    lw = 2,
    label = "Boltzmann",
    color = 1,
)
plot!(
    vSpace.u[:, vSpace.nv÷2, vSpace.nw÷2],
    data_bgk_1D[:, 2],
    lw = 2,
    line = :dash,
    label = "BGK",
    color = 2,
)
plot!(
    vSpace.u[:, vSpace.nv÷2, vSpace.nw÷2],
    M0_1D,
    lw = 2,
    label = "Maxwellian",
    color = 10,
)
scatter!(vSpace.u[:, vSpace.nv÷2, vSpace.nw÷2], sol.u[2], lw = 2, label = "UBE", color = 3)


text:
abstract type AbstractTimeStep end
abstract type AbstractScalarTimeStep{Adp,N} <: AbstractTimeStep end
abstract type AbstractScalarTimeStepWithIF{Adp,N} <: AbstractScalarTimeStep{Adp,N} end

@inline is_implicit(::Type{<:AbstractScalarTimeStep}) = false
@inline is_implicit(::Type{<:AbstractScalarTimeStepWithIF}) = true
@inline @par is_implicit(::Type{<:@par(AbstractSimulation)}) = is_implicit(VelocityTimeStepType)

has_variable_timestep(t::Type{<:AbstractScalarTimeStep{adp,N}}) where {adp,N} =
    adp

@inline has_variable_timestep(t::AbstractScalarTimeStep{adp,N}) where {adp,N} =
    has_variable_timestep(typeof(t))

set_dt!(t::Type{<:AbstractScalarTimeStep{false,N}},dt::Real) where {N} =
    nothing

@inline set_dt!(t::AbstractScalarTimeStep{false,N},dt::Real) where {N} = nothing

struct VectorTimeStep{cfl,Tx<:AbstractScalarTimeStep,Ty<:AbstractScalarTimeStep,Tz<:AbstractScalarTimeStep} <: AbstractTimeStep
    x::Tx
    y::Ty
    z::Tz
end

function VectorTimeStep{cfl}(tx,ty,tz) where {cfl}
    return VectorTimeStep{cfl,typeof(tx),typeof(ty),typeof(tz)}(tx,ty,tz)
end

@inline is_implicit(::Type{<:VectorTimeStep{cfl,Tx,Ty,Tz}}) where {cfl,Tx,Ty,Tz} = is_implicit(Tx)

get_cfl(t::Type{VectorTimeStep{cfl,Tx,Ty,Tz}}) where{cfl,Tx,Ty,Tz} =
    cfl

@inline get_cfl(t::VectorTimeStep) =
    get_cfl(typeof(t))

@par get_cfl(s::Union{@par(AbstractSimulation),Type{@par(AbstractSimulation)}}) =
    get_cfl(VelocityTimeStepType)

function initialize!(t::VectorTimeStep,rhs::VectorField,s::AbstractSimulation)
    vis = ν
    initialize!(t.x,rhs.rr.x,vis,s)
    initialize!(t.y,rhs.rr.y,vis,s)
    initialize!(t.z,rhs.rr.z,vis,s)
end

function set_dt!(t::VectorTimeStep,dt)
    set_dt!(t.x,dt)
    set_dt!(t.y,dt)
    set_dt!(t.z,dt)
end

@inline get_dt(t::A) where {A<:VectorTimeStep} =
    get_dt(t.x)

@inline get_dt(s::A) where {A<:AbstractSimulation} =
    get_dt(s.timestep)

has_variable_timestep(t::Type{VectorTimeStep{cfl,Tx,Ty,Tz}}) where {cfl,Tx,Ty,Tz} =
    has_variable_timestep(Tx)

@inline has_variable_timestep(t::VectorTimeStep) =
    has_variable_timestep(typeof(t))

@inline @par has_variable_timestep(s::Type{A}) where {A<:@par(AbstractSimulation)} =
    has_variable_timestep(VelocityTimeStepType)

@inline has_variable_timestep(s::A) where {A<:AbstractSimulation} =
    has_variable_timestep(A)

function (f::VectorTimeStep)(u::VectorField,rhs::VectorField,s::AbstractSimulation)
    f.x(u.c.x,rhs.c.x,s)
    f.y(u.c.y,rhs.c.y,s)
    f.z(u.c.z,rhs.c.z,s)

    add_forcing!(u,s.forcing)
    pressure_projection!(u,s)
    return nothing
end

include("Euller.jl")
include("Adams_Bashforth.jl")
include("ETD.jl")

@par function set_dt!(s::@par(AbstractSimulation))
    umax = maximum(s.reductionh)
    dx = 2π*LY/NY
    cfl = get_cfl(s)
    it = s.iteration[]

    if it == 0
        cfl *= 0.1
    elseif it == 1
        cfl *= 0.2
    elseif it == 2
        cfl *= 0.4
    elseif it == 3
        cfl *= 0.8
    end

    νt = ν
    newdt =cfl * dx/umax
    # is_implicit(s) || (newdt = min(newdt, (2νt/umax^2)/2))
    if hasdensity(s)
        ρmax = maximum(s.densitystratification.reduction)
        g = norm(gravity(s))
        #k = diffusivity(s)
        newdt = min(newdt,
        cfl * sqrt(dx/(ρmax*g)))#,
        #cfl * (dx^2)/(k)/2)#,
        #cfl * (((ρmax*g)^(-2/3))*(2k)^(1/3))/2,
        #cfl * (2k/umax^2)/2)
    end

    if it > 3 # do not allow dt to change abruptly
        dt = get_dt(s)
        if newdt < dt
            newdt = max(dt/1.25,newdt)
        else
            newdt = min(1.25*dt,newdt)
        end
    end

    set_dt!(s.timestep,newdt)
    hasdensity(s) && set_dt!(s.densitystratification.timestep,newdt)
    haspassivescalar(s) && set_dt!(s.passivescalar.timestep,newdt)
    return nothing
end


text:
module MaxHelpingHandShatterFlagSwitchGate

using ..Ahorn, Maple

@pardef ShatterFlagSwitchGate(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    sprite::String="block", persistent::Bool=false, flag::String="flag_touch_switch", icon::String="vanilla", inactiveColor::String="5FCDE4", activeColor::String="FFFFFF", finishColor::String="F141DF",
    shakeTime::Number=0.5, finishedSound::String="event:/game/general/touchswitch_gate_finish") =
    Entity("MaxHelpingHand/ShatterFlagSwitchGate", x=x, y=y, width=width, height=height,
    sprite=sprite, persistent=persistent, flag=flag, icon=icon, inactiveColor=inactiveColor, activeColor=activeColor, finishColor=finishColor,
    shakeTime=shakeTime, finishedSound=finishedSound)

const textures = String["block", "mirror", "temple", "stars"]
const bundledIcons = String["vanilla", "tall", "triangle", "circle", "diamond", "double", "heart", "square", "wide", "winged"]

const placements = Ahorn.PlacementDict(
    "Shatter Flag Switch Gate ($(uppercasefirst(texture)))\n(max480's Helping Hand + Vortex Helper)" => Ahorn.EntityPlacement(
        ShatterFlagSwitchGate,
        "rectangle",
        Dict{String, Any}(
            "sprite" => texture
        )
    ) for texture in textures
)

Ahorn.editingOrder(entity::ShatterFlagSwitchGate) = String["x", "y", "width", "height", "flag", "inactiveColor", "activeColor", "finishColor", "finishedSound", "shakeTime"]

Ahorn.editingOptions(entity::ShatterFlagSwitchGate) = Dict{String, Any}(
    "sprite" => textures,
    "icon" => bundledIcons
)

Ahorn.minimumSize(entity::ShatterFlagSwitchGate) = 16, 16
Ahorn.resizable(entity::ShatterFlagSwitchGate) = true, true

function Ahorn.selection(entity::ShatterFlagSwitchGate)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [Ahorn.Rectangle(x, y, width, height)]
end

function renderGateSwitch(ctx::Ahorn.Cairo.CairoContext, entity::ShatterFlagSwitchGate, x::Number, y::Number, width::Number, height::Number, sprite::String)
    icon = get(entity.data, "icon", "vanilla")

    iconResource = "objects/switchgate/icon00"
    if icon != "vanilla"
        iconResource = "objects/MaxHelpingHand/flagSwitchGate/$(icon)/icon00"
    end

    iconSprite = Ahorn.getSprite(iconResource, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    frame = "objects/switchgate/$sprite"

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x, y + (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, x + (i - 1) * 8, y + (j - 1) * 8, 8, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, x, y, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, x, y + height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, x + width - 8, y + height - 8, 16, 16, 8, 8)

    Ahorn.drawImage(ctx, iconSprite, x + div(width - iconSprite.width, 2), y + div(height - iconSprite.height, 2))
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ShatterFlagSwitchGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "block")

    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderGateSwitch(ctx, entity, x, y, width, height, sprite)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using Test

@testset "Benchmark" begin
    @info "Running benchmark tests..."

    include("../benchmark/benchmark_hooghoudt.jl")
end


text:
export @progress

using Logging: @logmsg, LogLevel

const PROGRESSLEVEL = LogLevel(-1)

"""
    progress(f::Function; name = "", msg = "")

Evaluates `f` with `id` as its argument and makes sure to destroy the progress
bar afterwards. To update the progress bar in `f` you can call a logging statement
like `@info` or even just `@logmsg` with `_id=id` and `progress` as arguments.

`progress` can take either of the following values:
  - `0 <= progress < 1`: create or update progress bar
  - `progress == nothing || progress = NaN`: set progress bar to indeterminate progress
  - `progress >= 1 || progress == "done"`: destroy progress bar

The logging message (e.g. `"foo"` in `@info "foo"`) will be used as the progress
bar's name.

```julia
Juno.progress() do id
    for i = 1:10
        sleep(0.5)
        @info "iterating" progress=i/10 _id=id
    end
end
```
"""
function progress(f; name = "")
  _id = gensym()
  @logmsg PROGRESSLEVEL name progress=NaN _id=_id
  try
    f(_id)
  finally
    @logmsg PROGRESSLEVEL name progress="done" _id=_id
  end
end

"""
    @progress [name="", threshold=0.005] for i = ..., j = ..., ...
    @progress [name="", threshold=0.005] x = [... for i = ..., j = ..., ...]

Show a progress meter named `name` for the given loop or array comprehension
if possible. Update frequency is limited by `threshold` (one update per 0.5% of
progress by default).
"""
macro progress(args...)
  _progress(args...)
end

_progress(ex) = _progress("", 0.005, ex)
_progress(name::Union{AbstractString, Expr}, ex) = _progress(name, 0.005, ex)
_progress(thresh::Real, ex) = _progress("", thresh, ex)

function _progress(name, thresh, ex)
  if ex.head == Symbol("=") &&
        ex.args[2].head == :comprehension &&
        ex.args[2].args[1].head == :generator
    # comprehension: <target> = [<body> for <iter_var> in <range>,...]
    loop = _comprehension
    target = esc(ex.args[1])
    result = target
    gen_ex = ex.args[2].args[1]
    body = esc(gen_ex.args[1])
    iter_exprs = gen_ex.args[2:end]
    iter_vars = [e.args[1] for e in iter_exprs]
    ranges = [e.args[2] for e in iter_exprs]
  elseif ex.head == :for &&
        ex.args[1].head == Symbol("=") &&
        ex.args[2].head == :block
    # single-variable for: for <iter_var> = <range>; <body> end
    loop = _for
    target = :_
    result = :nothing
    iter_vars = [ex.args[1].args[1]]
    ranges = [ex.args[1].args[2]]
    body = esc(ex.args[2])
  elseif ex.head == :for &&
        ex.args[1].head == :block &&
        ex.args[2].head == :block
    # multi-variable for: for <iter_var> = <range>,...; <body> end
    loop = _for
    target = :_
    result = :nothing
    # iter_vars and ranges are ordered from inner loop to outer loop, for
    # consistent computation of progress between for loops and comprehensions
    iter_vars = reverse([e.args[1] for e in ex.args[1].args])
    ranges = reverse([e.args[2] for e in ex.args[1].args])
    body = esc(ex.args[2])
  else
    error("@progress requires a for loop (for i in irange, j in jrange, ...; <body> end) " *
          "or array comprehension with assignment (x = [<body> for i in irange, j in jrange, ...])")
  end
  _progress(name, thresh, ex, target, result, loop, iter_vars, ranges, body)
end

function _progress(name, thresh, ex, target, result, loop, iter_vars, ranges, body)
  count_vars = [Symbol("i$k") for k=1:length(iter_vars)]
  iter_exprs = [:(($i,$(esc(v))) = enumerate($(esc(r))))
                  for (i,v,r) in zip(count_vars,iter_vars,ranges)]
  _id = "progress_$(gensym())"
  quote
    if isactive()
      @logmsg($PROGRESSLEVEL, $(esc(name)), progress=0.0, _id=Symbol($_id))
      $target = try
        ranges = $(Expr(:vect,esc.(ranges)...))
        nranges = length(ranges)
        lens = length.(ranges)
        n = prod(lens)
        strides = cumprod([1;lens[1:end-1]])
        _frac(i) = (sum((i-1)*s for (i,s) in zip(i,strides)) + 1) / n
        lastfrac = 0.0

        $(loop(iter_exprs,
            quote
                val = $body
                frac = _frac($(Expr(:vect, count_vars...)))
                if frac - lastfrac > $thresh
                    @logmsg($PROGRESSLEVEL, $(esc(name)), progress=frac, _id=Symbol($_id))
                    lastfrac = frac
                end
                val
            end
        ))

      finally
        @logmsg($PROGRESSLEVEL, $(esc(name)), progress="done", _id=Symbol($_id))
      end
      $result
    else
      $(esc(ex))
    end
  end
end

_comprehension(iter_exprs, body,) = Expr(:comprehension, Expr(:generator, body, iter_exprs...))
_for(iter_exprs, body) = Expr(:for, Expr(:block, reverse(iter_exprs)...), body)


text:
# included in client_serverWS_test.jl
# and in client_listen_test.jl

if !@isdefined SUBPROTOCOL
    const SUBPROTOCOL = "Server start the conversation"
    const SUBPROTOCOL_CLOSE = "Server start the conversation and close it from within websocket handler"
end
addsubproto(SUBPROTOCOL)
addsubproto(SUBPROTOCOL_CLOSE)
if !@isdefined(PORT)
    const PORT = 8000
    const SURL = "127.0.0.1"
    const EXTERNALWSURI = "ws://echo.websocket.org"
    const EXTERNALHTTP = "http://httpbin.org/ip"
    const MSGLENGTHS = [0 , 125, 126, 127, 2000]
end

"""
`test_handler` is called by WebSockets inner function `_servercoroutine` for all accepted http requests
that are not upgrades. We don't check what's actually requested.
"""
test_handler(req::HTTP.Request) = HTTP.Response(200, "OK")

"""
`test_wshandler` is called by WebSockets inner function
`_servercoroutine` for all http requests that
    1) qualify as an upgrade,
    2) request a subprotocol we claim to support
Based on the requested subprotocol, test_wshandler calls
    `initiatingws`
        or
    `echows`
"""
function test_wshandler(req::HTTP.Request, ws::WebSocket)
    WebSockets.origin(req) != "" && @error "test_wshandler, got origin header as from a browser."
    WebSockets.target(req) != "/" && @error "test_wshandler, got origin header as in a POST request."
    if WebSockets.subprotocol(req) == SUBPROTOCOL
        initiatingws(ws, msglengths = MSGLENGTHS)
    elseif WebSockets.subprotocol(req) == SUBPROTOCOL_CLOSE
        initiatingws(ws, msglengths = MSGLENGTHS,  closebeforeexit = true)
    else
        echows(ws)
    end
end

"""
`echows` is called by
    - `test_wshandler` (in which case ws will be a server side websocket)
    or
    - 'WebSockets.open' (in which case ws will be a client side websocket)

Takes an open websocket.
    1)  Reads a message
    2)  Echoes it
    3)  Repeats until the websocket closes, or a read fails.
The tests will be captured if the function is run on client side.
If started by the server side, this is called as part of a coroutine.
Therefore, test results will not propagate to the enclosing test scope.
"""
function echows(ws::WebSocket)
    while isopen(ws)
        data, ok = readguarded(ws)
        if ok
            if writeguarded(ws, data)
                @test true
            else
                break
            end
        else
            if !isopen(ws)
                break
            else
                break
            end
        end
    end
end

"""
`initiatingws` is called by
    - `test_wshandler` (in which case ws will be a server side websocket)
    or
    - 'WebSockets.open' (in which case ws will be a client side websocket)

Takes
    - an open websocket
    keyword arguments
    - msglengths = a vector of message lengths, defaults to MSGLENGTHS
    - closebeforeexit, defaults to false

1) Pings, but does not check for received pong. There will be console output from the pong side.
2) Send a message of specified length
3) Checks for an exact echo
4) Repeats 2-4 until no more message lenghts are specified.
"""
function initiatingws(ws::WebSocket; msglengths = MSGLENGTHS, closebeforeexit = false)
    send_ping(ws)
    # No ping check made, the above will just output some info message.

    # We need to yield since we are sharing the same process as the task on the
    # other side of the connection.
    # The other side must be reading in order to process the ping-pong.
    yield()
    for slen in msglengths
        test_str = Random.randstring(slen)
        forcecopy_str = test_str |> collect |> copy |> join
        if writeguarded(ws, test_str)
            yield()
            readback, ok = readguarded(ws)
            if ok
                # if run by the server side, this test won't be captured.
                if String(readback) == forcecopy_str
                    @test true
                else
                    if ws.server == true
                        @error "initatews, echoed string of length ", slen, " differs from sent "
                    else
                        @test false
                    end
                end
            else
                # if run by the server side, this test won't be captured.
                if ws.server == true
                    @error "initatews, couldn't read ", ws, " length sent is ", slen
                else
                    @test false
                end
            end
        else
            @test false
        end
    end
    closebeforeexit && close(ws, statusnumber = 1000)
end

test_serverws = WebSockets.ServerWS(
    HTTP.RequestHandlerFunction(test_handler),
    WebSockets.WSHandlerFunction(test_wshandler))

function startserver(serverws=test_serverws;url=SURL, port=PORT, verbose=false)
    servertask = @async WebSockets.serve(serverws,url,port,verbose)
    while !istaskstarted(servertask);yield();end
    if isready(serverws.out)
        # capture errors, if any were made during the definition.
        @error take!(serverws.out)
    end
    serverws
end


text:
ts = "2019-03-23T23:10:00"
te = "2019-03-23T23:17:00"
chans = [ "CC.JRO..BHZ",
          "CC.LON..BHZ",
          "CC.VALT..BHZ",
          "UW.HOOD..ENZ",
          "UW.HOOD..ENN",
          "UW.HOOD..ENE" ]

printstyled("  IRIS Web Services\n", color=:light_green)
printstyled("    IRISWS continuous data requests\n", color=:light_green)
printstyled("    SAC\n", color=:light_green)
for i = 1:length(chans)
  cha = chans[i]
  S = get_data("IRIS", cha, src="IRIS", s=ts, t=te, fmt="sacbl", v=0, w=true)
  sleep(5)
  T = get_data("IRIS", cha, src="IRIS", s=ts, t=te, fmt="mseed", v=0, w=true)

  if S.n == 0 || T.n == 0
    @warn(string(cha, " appears to be offline; trying next."))
    if i == lastindex(chans)
      error("No data for any station; failing test due to connection errors.")
    end
  else
    printstyled("     equivalence of SAC and MSEED requests\n", color=:light_green)
    sync!(S, s=ts, t=te)
    sync!(T, s=ts, t=te)
    @test S.x[1] ≈ T.x[1]

    # Only notes and src should be different
    for f in Symbol[:id, :name, :loc, :fs, :gain, :resp, :units, :misc , :t, :x]
      @test getfield(S,f) == getfield(T,f)
    end

    printstyled("    GeoCSV\n", color=:light_green)
    U = get_data("IRIS", cha, src="IRIS", s=ts, t=te, fmt="geocsv", v=0)
    break
  end
end

# Test a bum data format
printstyled("    bad data format\n", color=:light_green)
sta_matrix = vcat(["UW" "LON" "" "BHZ"],["UW" "LON" "" "BHE"])
test_sta = deepcopy(sta_matrix)
redirect_stdout(out) do
  T = get_data("IRIS", sta_matrix, s=-600, t=0, v=2, fmt="audio")
  T = get_data("IRIS", sta_matrix, s=-600, t=0, v=2, fmt="ascii")
end

# check that these aren't modified in-place by the request
@test sta_matrix == test_sta

# Test a bad request
printstyled("    bad request format\n", color=:light_green)
T = get_data("IRIS", "DE.NENA.99.LUFTBALLOONS", src="IRIS", s=ts, t=te, fmt="mseed", v=0)

printstyled("    complicated IRISWS request\n", color=:light_green)
chans = ["UW.TDH..EHZ", "UW.VLL..EHZ", "CC.JRO..BHZ"] # HOOD is either offline or not on IRISws right now
st = -86400.0
en = -86100.0
(d0,d1) = parsetimewin(st,en)
S = get_data("IRIS", chans, s=d0, t=d1, y=true)
if isempty(S)
  @warn(string("No data for channels ", join(chans, ", "), "; test skipped."))
else
  L = [length(S.x[i])/S.fs[i] for i = 1:S.n]
  t = [S.t[i][1,2] for i = 1:S.n]
  L_min = minimum(L)
  L_max = maximum(L)
  t_min = minimum(t)
  t_max = maximum(t)
  try
    @test(L_max - L_min <= maximum(2 ./ S.fs))
    @test(t_max - t_min <= round(Int64, 1.0e6 * 2.0/maximum(S.fs)))
  catch
    @warn(string("Unexpected request length; check for partial outage at IRIS for ", join(chans, ", "), "!"))
  end
end


text:
#=
Vector
=#

using .Libc: malloc, realloc

"""
Space complexity: O(n)
"""
mutable struct SimpleVector{T} <: AbstractVector{T}
    buffer::Ptr{T}
    size::Int
    capacity::Int

    # Only allow isbitstypes for simplicity
    function SimpleVector{T}(
        buffer::Ptr{T},
        size::Integer,
        capacity::Integer,
    ) where T
        isbitstype(T) || throw(ArgumentError("eltype of `SimpleVector` must be `isbitstype`"))
        new{T}(buffer, size, capacity)
    end
end

#=
Constructors
=#

function SimpleVector{T}(; size::Integer = 0, capacity::Integer = 2) where T
    size ≤ capacity || throw(ArgumentError("`size` must be ≤ `capacity`"))
    n_bits = sizeof(T) * 8 * capacity
    buffer = convert(Ptr{T}, malloc(n_bits))
    SimpleVector{T}(buffer, size, capacity)
end

function SimpleVector(xs::AbstractVector{T}) where T
    n = length(xs)
    sv = SimpleVector{T}(size = n, capacity = n)
    for i in 1:n
        sv[i] = xs[i]
    end
    return sv
end

#=
Core methods
=#

"""
Time complexity: O(1)
"""
function Base.getindex(sv::SimpleVector{T}, i::Integer)::T where T
    checkbounds(sv, i)
    return unsafe_load(sv.buffer, i)
end

"""
Time complexity: O(1)
"""
function Base.setindex!(sv::SimpleVector{T}, x, i::Integer) where T
    x_T = convert(T, x)
    checkbounds(sv, i)
    unsafe_store!(sv.buffer, x_T, i)
    return sv
end

"""
Time complexity: amortised O(1), O(n) worst case
"""
function Base.push!(sv::SimpleVector{T}, x) where T
    x_T = convert(T, x)
    sv.size += 1
    if sv.size > sv.capacity
        new_capacity = 2^ceil(Int, log2(sv.capacity))
        resize!(sv, new_capacity)
    end
    setindex!(sv, x_T, sv.size)
end

"""
Time complexity: O(n)
"""
function Base.resize!(sv::SimpleVector{T}, n::Integer) where T
    if sv.capacity == n
        return sv
    end
    n_bits = sizeof(T) * 8 * n
    sv.buffer = realloc(sv.buffer, n_bits)
    sv.capacity = n
    if n < sv.size
        sv.size = n
    end
    return sv
end

"""
Time complexity: O(n)
"""
function Base.insert!(sv::SimpleVector{T}, index::Integer, item::T) where T
    checkbounds(sv, index)
    if sv.size == sv.capacity
        resize!(sv, sv.capacity + 1)
    end
    sv.size += 1
    for i = sv.size-1:-1:index
        sv[i+1] = sv[i]
    end
    sv[index] = item
    return sv
end

"""
Time complexity: O(n)
"""
function Base.deleteat!(sv::SimpleVector, index::Integer)
    checkbounds(sv, index)
    for i = index:sv.size-1
        sv[i] = sv[i+1]
    end
    sv.size -= 1
    return sv
end

#=
Interface
=#

Base.size(sv::SimpleVector) = (sv.size,)

IndexStyle(::SimpleVector) = IndexLinear


text:
function backtracking_linesearch!{T}(d::Union{DifferentiableFunction,
                                              TwiceDifferentiableFunction},
                                     x::Vector{T},
                                     s::Vector,
                                     x_scratch::Vector,
                                     gr_scratch::Vector,
                                     lsr::LineSearchResults,
                                     alpha::Real = 1.0,
                                     mayterminate::Bool = false,
                                     c1::Real = 1e-4,
                                     c2::Real = 0.9,
                                     rho::Real = 0.9,
                                     iterations::Integer = 1_000)

    # Count the total number of iterations
    iteration = 0

    # Track calls to function and gradient
    f_calls = 0
    g_calls = 0

    # Count number of parameters
    n = length(x)

    # Store f(x) in f_x
    f_x = d.fg!(x, gr_scratch)
    f_calls += 1
    g_calls += 1

    # Store angle between search direction and gradient
    gxp = vecdot(gr_scratch, s)

    # Tentatively move a distance of alpha in the direction of s
    @simd for i in 1:n
        @inbounds x_scratch[i] = x[i] + alpha * s[i]
    end

    # Backtrack until we satisfy sufficient decrease condition
    f_x_scratch = d.f(x_scratch)
    f_calls += 1
    while f_x_scratch > f_x + c1 * alpha * gxp
        # Increment the number of steps we've had to perform
        iteration += 1

        # Ensure termination
        if iteration > iterations
            error("Too many iterations in backtracking_linesearch!")
        end

        # Shrink proposed step-size
        alpha *= rho

        # Update proposed position
        @simd for i in 1:n
            @inbounds x_scratch[i] = x[i] + alpha * s[i]
        end

        # Evaluate f(x) at proposed position
        f_x_scratch = d.f(x_scratch)
        f_calls += 1
    end

    return alpha, f_calls, g_calls
end


text:
include("help.jl")
function mark_rows!(r::Robot)
    num_ost = moves!(r, Ost)
    num_vert = moves!(r, Sud)
    num_hor = moves!(r, West)
    #УТВ: Робот - в Юго-Западном углу

    putmarkers!(r, num_hor)
    #УТВ: Все ряды поля промаркированы в соответствующем порядке

    moves!(r, Sud)
    moves!(r, West)
    moves!(r, Nord, num_vert)
    moves!(r, Ost, num_hor-num_ost)
    #УТВ: Робот - в исходном положении
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while move_if_possible!(r,side)
        num_steps+=1
    end
    return num_steps
end

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move_if_possible!(r,side)
    end
end

function put_num_markers!(r::Robot,side::HorizonSide,num_marks::Int)
    putmarker!(r)
    for _ in 1:num_marks
        move_if_possible!(r,side)
        putmarker!(r)
    end
end

function putmarkers!(r::Robot, num_marks::Int)
    while !isborder(r,Nord)
        put_num_markers!(r, Ost, num_marks)

        move!(r, Nord)
        moves!(r, West)
        num_marks -= 1
    end

    put_num_markers!(r, Ost, num_marks)
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a721d7c0-980a-11eb-1b4e-23bc8dbc5211
begin
	using Pluto, PlutoUI, HypertextLiteral
end

# ╔═╡ 71bbf1ec-c8c4-4e19-b544-0fd5752e2869
using ImageTransformations

# ╔═╡ 6a4b75d5-03f7-43cf-a9dc-e1a4781397da
using Images

# ╔═╡ cd8c0ee9-9c9e-4099-83a1-e208f569a455
md"# Classifying handwritten digits!"

# ╔═╡ 0e600c58-6156-41da-b0cb-c638c1ef6798
@htl("""

	$(@bind classify_canvas Button("Classify"))
	$(@bind clear_canvas Button("Clear"))

""")

# ╔═╡ 115b959e-4e44-42ca-a97d-73be8c3826de
@htl("""

Try writing a digit in the box!

<canvas id="canvas" width="140" height="140" style="border: 1px solid grey">HTML5 is required to run this notebook</canvas>

""")

# ╔═╡ 33652a07-96ef-4388-94f6-3b60e33df6ff
@htl("""

""")

# ╔═╡ 42743644-4984-4b84-a924-bb6e28f85fc8
@bind pixel_data_raw html"""<input type="hidden" id="hidden_imgdata"/>"""

# ╔═╡ 5e73ad1c-b2a5-40ee-ad6c-0a633f8974b1
pixel_data_flat = parse.(Int, split(pixel_data_raw, ","))

# ╔═╡ 68fd63ec-6181-44d0-a6bf-fa3684b12744
function get_pixel_data(x, y)
	return pixel_data_flat[((x - 1) * 140 + y) * 4]
end

# ╔═╡ 8f9c01ef-03c5-451b-a9c2-1a672b003f91
pixel_data = [get_pixel_data(x, y) for x ∈ 1:140, y ∈ 1:140]

# ╔═╡ fce4dd05-8b31-4d0e-9d4a-5dde003b90d0
image = Gray.(pixel_data ./ 255)

# ╔═╡ bb579671-5685-405a-a0b9-57cf367ef27d
small_image = imresize(image, 28, 28)

# ╔═╡ de39700a-7b6d-44bf-9296-c053fd997aea
input_images = permutedims(reshape(collect(channelview(small_image)), (28, 28, 1, 1)), (2, 1, 3, 4));

# ╔═╡ 636054d4-8c33-4995-8e59-e0c1cce1c4fe
classifier_notebook = PlutoNotebook("Connor Burns — MNIST classifier for WYSIWYR.jl")

# ╔═╡ 699f3821-41e3-4fdc-814a-f906be3fe2ce
output_labels = classifier_notebook(; input_images=input_images).output_labels

# ╔═╡ 27f58f6f-01d9-4d45-a876-04dc6d5e9d72
md"### Prediction: $(output_labels |> first)"

# ╔═╡ 2465eae0-f58e-4f6c-8571-ac5e816b71c6
md"### Guess: **$(output_labels |> first)**"

# ╔═╡ Cell order:
# ╟─cd8c0ee9-9c9e-4099-83a1-e208f569a455
# ╟─a721d7c0-980a-11eb-1b4e-23bc8dbc5211
# ╟─115b959e-4e44-42ca-a97d-73be8c3826de
# ╟─0e600c58-6156-41da-b0cb-c638c1ef6798
# ╟─27f58f6f-01d9-4d45-a876-04dc6d5e9d72
# ╟─33652a07-96ef-4388-94f6-3b60e33df6ff
# ╟─42743644-4984-4b84-a924-bb6e28f85fc8
# ╠═5e73ad1c-b2a5-40ee-ad6c-0a633f8974b1
# ╠═68fd63ec-6181-44d0-a6bf-fa3684b12744
# ╠═8f9c01ef-03c5-451b-a9c2-1a672b003f91
# ╠═71bbf1ec-c8c4-4e19-b544-0fd5752e2869
# ╠═6a4b75d5-03f7-43cf-a9dc-e1a4781397da
# ╠═fce4dd05-8b31-4d0e-9d4a-5dde003b90d0
# ╠═bb579671-5685-405a-a0b9-57cf367ef27d
# ╠═de39700a-7b6d-44bf-9296-c053fd997aea
# ╠═636054d4-8c33-4995-8e59-e0c1cce1c4fe
# ╠═699f3821-41e3-4fdc-814a-f906be3fe2ce
# ╟─2465eae0-f58e-4f6c-8571-ac5e816b71c6

max_stars_repo_path: MNIST Client.jl

max_stars_repo_name: ctrekker/pluto-static-test

max_issues_repo_path: MNIST Client.jl

max_issues_repo_name: ctrekker/pluto-static-test

max_forks_repo_path: MNIST Client.jl

max_forks_repo_name: ctrekker/pluto-static-test


text:
const MultiTaskLassoRegressor_ = sklm(:MultiTaskLasso)
@sk_reg mutable struct MultiTaskLassoRegressor <: MMI.Deterministic
    alpha::Float64      = 1.0::(_ ≥ 0)
    fit_intercept::Bool = true
    normalize::Bool     = false
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    copy_X::Bool        = true
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::MultiTaskLassoRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
    )
add_human_name_trait(MultiTaskLassoRegressor, "multi-target lasso regressor")

# ==============================================================================
const MultiTaskLassoCVRegressor_ = sklm(:MultiTaskLassoCV)
@sk_reg mutable struct MultiTaskLassoCVRegressor <: MMI.Deterministic
    eps::Float64        = 1e-3::(_ > 0)
    n_alphas::Int       = 100::(_ > 0)
    alphas::Any         = nothing::(_ === nothing || all(0 .≤ _ .≤ 1))
    fit_intercept::Bool = true
    normalize::Bool     = false
    max_iter::Int       = 300::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    copy_X::Bool        = true
    cv::Any             = 5
    verbose::Union{Bool, Int} = false
    n_jobs::Option{Int} = 1
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::MultiTaskLassoCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    mse_path  = fitresult.mse_path_,
    alphas    = fitresult.alphas_
    )
add_human_name_trait(MultiTaskLassoCVRegressor, "multi-target lasso regressor $CV")

# ==============================================================================
const MultiTaskElasticNetRegressor_ = sklm(:MultiTaskElasticNet)
@sk_reg mutable struct MultiTaskElasticNetRegressor <: MMI.Deterministic
    alpha::Float64      = 1.0::(_ ≥ 0)
    l1_ratio::Union{Float64, Vector{Float64}} = 0.5::(0 ≤ _ ≤ 1)
    fit_intercept::Bool = true
    normalize::Bool     = true
    copy_X::Bool        = true
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    warm_start::Bool    = false
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::MultiTaskElasticNetRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing)
    )
add_human_name_trait(MultiTaskElasticNetRegressor, "multi-target elastic net regressor")

# ==============================================================================
const MultiTaskElasticNetCVRegressor_ = sklm(:MultiTaskElasticNetCV)
@sk_reg mutable struct MultiTaskElasticNetCVRegressor <: MMI.Deterministic
    l1_ratio::Union{Float64, Vector{Float64}} = 0.5::(0 ≤ _ ≤ 1)
    eps::Float64        = 1e-3::(_ > 0)
    n_alphas::Int       = 100::(_ > 0)
    alphas::Any         = nothing::(_ === nothing || all(0 .≤ _ .≤ 1))
    fit_intercept::Bool = true
    normalize::Bool     = false
    max_iter::Int       = 1_000::(_ > 0)
    tol::Float64        = 1e-4::(_ > 0)
    cv::Any             = 5
    copy_X::Bool        = true
    verbose::Union{Bool,Int} = 0
    n_jobs::Option{Int} = nothing
    random_state::Any   = nothing
    selection::String   = "cyclic"::(_ in ("cyclic","random"))
end
MMI.fitted_params(model::MultiTaskElasticNetCVRegressor, (fitresult, _, _)) = (
    coef      = fitresult.coef_,
    intercept = ifelse(model.fit_intercept, fitresult.intercept_, nothing),
    alpha     = fitresult.alpha_,
    mse_path  = fitresult.mse_path_,
    l1_ratio  = fitresult.l1_ratio_
    )
add_human_name_trait(MultiTaskElasticNetCVRegressor, "multi-target elastic "*
                     "net regressor $CV")

const SKL_REGS_MULTI = Union{
       Type{<:MultiTaskLassoRegressor},
       Type{<:MultiTaskLassoCVRegressor},
       Type{<:MultiTaskElasticNetRegressor},
       Type{<:MultiTaskElasticNetCVRegressor}
    }

MMI.input_scitype(::SKL_REGS_MULTI)  = Table(Continuous)
MMI.target_scitype(::SKL_REGS_MULTI) = Table(Continuous)


text:
using DataDeps

register(DataDep(
  "libras-movement",
  "http://archive.ics.uci.edu/ml/datasets/Libras+Movement",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/libras/movement_libras.data",
  "97ebdaa6a9b28ab4a2cdd84b14f19a95a7456a46137c362b65a0669eca3c3c4d",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=91,
        feature_indices=1:90,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module Zm_analysis

using Printf, Plots, JSON, Glob, Statistics
gr(size=(800,500), html_output_format=:png)

export read_Zm, Zm_E, Zm_T, Zm_PT

function read_Zm(filename)
    return JSON.parsefile(filename)
end

function Zm_E(D; return_plot=true)
    p = plot(legend =:false)

    for (index, value) in enumerate(keys(D))
        E_tot, E_kin, E_pot = D[value]["E_tot"], D[value]["E_kin"], D[value]["E_pot"]

        plot!(p,  E_tot)
        plot!(p, E_kin)
        plot!(p, E_pot)
    end

    title!(D["1"]["info"])
    xlabel!("MD step")
    ylabel!("Energy (eV)")

    if return_plot == false
        savefig(@sprintf("E_%s.png", D["1"]["info"]))
    else
        display(p)
    end
end

function Zm_T(D; return_plot=true)
    p = plot(legend =:false)

    for (index, value) in enumerate(keys(D))
        T = D[value]["T"]

        plot!(p, T)
    end

    title!(D["1"]["info"])
    xlabel!("MD step")
    ylabel!("Temperature (K)")

    if return_plot == false
        savefig(@sprintf("T_%s.png", D["1"]["info"]))
    else
        display(p)
    end
end

function Zm_PT(D, delta_step=2, large_data_points=2, eq_length=10; GPa=true, return_plot=true)
    p = scatter(legend =:false)

    colors = ["blue", "red", "green"]

    for (index, value) in enumerate(keys(D))
        T, P = D[value]["T"], D[value]["P"]

        # println(T, P)

        #scatter!(P[eq_length:delta_step:end], T[eq_length:delta_step:end], markersize=1, markerstrokewidth=0, color=colors[index], alpha=0.5)

        large_data_point_step = (length(T)-eq_length)/large_data_points

        mean_T_list, mean_P_list = [], []
        error_T_list, error_P_list = [], []

        for j in 1:large_data_points
            st = round(Int,((j-1)*large_data_point_step + eq_length))
            en = round(Int,(j*large_data_point_step))

            T_slice = Float64.(T[ st : en])
            P_slice = Float64.(P[ st : en])

            mean_T = mean(T_slice)
            mean_P = mean(P_slice)

            error_T = std(T_slice)
            error_P = std(P_slice)

            push!(mean_T_list, mean_T)
            push!(error_T_list, error_T)
            push!(mean_P_list, mean_P)
            push!(error_P_list, error_P)
        end

        if GPa
            mean_P_list = mean_P_list .* 160.21766208
            error_P_list = error_P_list .* 160.21766208
        end

        scatter!(mean_P_list, mean_T_list, xerr=error_P_list, yerr=error_T_list, markersize=5, color=colors[index])
    end

    title!(D["1"]["info"])
    xlabel!("Pressure [GPa]")
    ylabel!("Temperature [K]")

    if return_plot == false
        savefig(@sprintf("PT_%s.png", D["1"]["info"]))
    else
        display(p)
    end
end

end


text:
#Code written by Skipper af Rosenborg and Esben Bager
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found

#Load packages and support functions
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")

#Set path
#Esben's path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexData"

#Set path for dataloading
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

#Set parameters
trainingSize = 240
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

#Define the input and transformations
modelMatrix = zeros(industriesTotal, possibilities)
noDurModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

#Create and fill X and Y arrays
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)
XArrays, YArrays = generateXandYs(industries, modelMatrix)

#Set
nGammas = 5
standY = YArrays[1]
nRows = size(standY)[1]
amountOfModels = nGammas^4

#Initialization of parameters
gammaArr = linspace(0,4,16) #2.4 #risk aversion
w1N = repeat([1/10], outer = 10) #1/N weights
validationPeriod = 0 # Number of periods to do parameter training
PMatrix = zeros(nRows-trainingSize, amountOfModels)

bestModelAmount = 1
bestModelConfigs = zeros(bestModelAmount, 4)
bestModelIndexes = zeros(bestModelAmount)

return1NMatrix = zeros(nRows-trainingSize)
returnCEOMatrix = zeros(nRows-trainingSize, bestModelAmount)
returnPerfectMatrix = zeros(nRows-trainingSize)

weightsPerfect = zeros(nRows-trainingSize, 10)
weightsCEO     = zeros(nRows-trainingSize, 10, bestModelAmount)
expectedReturnMatrix = zeros(nRows-trainingSize, 10)
forecastErrors = zeros(nRows-trainingSize, 10)

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/VIXTimeTA2.4/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"
weightsPerfect = Array{Float64}(CSV.read(path*"weightsPerfect.csv",header=false, datarow=1, nullable=false))
returnPerfectMatrix = Array{Float64}(CSV.read(path*"returnPerfectMatrix.csv", header=false, datarow=1, nullable=false))

#Set model configurations
bestModelConfigs[1,:] = [1 0.5  0.25  100] #to change!
#bestModelConfigs[2,:] = [1 0.5  0.25  100] #to change!

for i = 1:bestModelAmount
    bestModelIndexes[i] = i
end

for l = 1:16
    #Set GAMMA
    gamma gammaArr[l]

    inputArg1 = 0 #83 breaks the trainingSize
    inputArg1 = parse(Int64,ARGS[1])

    if 10+validationPeriod+(10*inputArg1) <= nRows-trainingSize-2
        println("Starting CEO Validation loop")
        for t=validationPeriod+1+(10*inputArg1):10+validationPeriod+(10*inputArg1)
            #Split data
            trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

            for m = 1:bestModelAmount
                #Set expectedReturns
                expectedReturns = zeros(10)

                #Build and solve CEO
                betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma))

                expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
                if m == 1
                    expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100
                end

                valY = zeros(10)
                for i = 1:10
                    valY[i] = validationY[i][1]
                end

                #Create sigma for MV
                trainX = trainingXArrays[1][:,1:10]
                Sigma =  cov(trainX)
                F = lufact(Sigma)
                U = F[:U]  #Cholesky factorization of Sigma

                return1N, returnCEO, wStar, forecastRow = performMVOptimization(expectedReturns, U, gamma, valY, valY)
                weightsCEO[t, 1:10, Int64(m)] = wStar
                return1NMatrix[t]      = return1N
                returnCEOMatrix[t, Int64(bestModelIndexes[m])] = returnCEO
                returnPerfect = returnPerfectMatrix[t]

                if m == 1
                    forecastErrors[t, 1:10]  = forecastRow-expectedReturnMatrix[t,:]
                end
                println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
                PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
            end
        end
    elseif validationPeriod+(10*inputArg1) <= nRows-trainingSize-2
        for t=validationPeriod+1+(10*inputArg1):(nRows-trainingSize-2)
            println("Time $t/50")
            trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

            for m = 1:bestModelAmount
                expectedReturns = zeros(10)

                betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma))

                expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
                if m == 1
                    expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100
                end

                valY = zeros(10)
                for i = 1:10
                    valY[i] = validationY[i][1]
                end

                trainX = trainingXArrays[1][:,1:10]
                Sigma =  cov(trainX)
                F = lufact(Sigma)
                U = F[:U]  #Cholesky factorization of Sigma

                return1N, returnCEO, wStar, forecastRow = performMVOptimization(expectedReturns, U, gamma, valY, valY)
                weightsCEO[t, 1:10, Int64(m)] = wStar
                return1NMatrix[t]      = return1N
                returnCEOMatrix[t, Int64(bestModelIndexes[m])] = returnCEO
                returnPerfect = returnPerfectMatrix[t]
                if m == 1
                    forecastErrors[t, 1:10]  = forecastRow-expectedReturnMatrix[t,:]
                end
                println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
                PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
            end
        end
    end

    #Output data
    combinedPortfolios = hcat(returnPerfectMatrix[1:nRows-trainingSize-2, 1], return1NMatrix[1:nRows-trainingSize-2, 1],
        returnCEOMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)], PMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)])
    #path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
    path = "/zhome/9f/d/88706/SpecialeCode/Results/CEO/"*string(l)*"/"
    writedlm(path*string(inputArg1)*"_returnPvalueOutcome1to200.csv", combinedPortfolios, ",")
    for i = 1:bestModelAmount
        writedlm(path*"wightsCEO_Model"*string(i)*"_"*string(inputArg1)*".csv",weightsCEO[:,:,i], ",")
    end

    combinedPortfolios2 = hcat(expectedReturnMatrix, forecastErrors)
    writedlm(path*string(inputArg1)*"_forecastsAndErrors1to200.csv", combinedPortfolios2, ",")
end
println("Finished Everything")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/CEO.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/CEO.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
"""
Translate mouse or modifier state into the corresponding XCB value.
"""
function state_xcb end

state_xcb(s::MouseState) = .&(UInt[s.left, s.middle, s.right, s.scroll_up, s.scroll_down, s.any] .* UInt[XCB_BUTTON_MASK_1, XCB_BUTTON_MASK_2, XCB_BUTTON_MASK_3, XCB_BUTTON_MASK_4, XCB_BUTTON_MASK_5, XCB_BUTTON_MASK_ANY]...)
state_xcb(s::KeyModifierState) = sum(2 .^ (0:3) .* Int[s.shift, s.ctrl, s.alt, s.super])
state_xcb(::Any) = 0

state_xcb(e::MouseEvent) = state_xcb(e.state)
state_xcb(e::KeyEvent) = state_xcb(e.modifiers)

"""
Translate an action into its corresponding XCB value.
"""
function response_type_xcb end

response_type_xcb(::ButtonPressed) = XCB_BUTTON_PRESS
response_type_xcb(::ButtonReleased) = XCB_BUTTON_RELEASE
response_type_xcb(::KeyPressed) = XCB_KEY_PRESS
response_type_xcb(::KeyReleased) = XCB_KEY_RELEASE
response_type_xcb(::PointerMoves) = XCB_MOTION_NOTIFY
response_type_xcb(::PointerEntersWindow) = XCB_ENTER_NOTIFY
response_type_xcb(::PointerLeavesWindow) = XCB_LEAVE_NOTIFY
response_type_xcb(::Expose) = XCB_EXPOSE
response_type_xcb(::Resize) = XCB_CONFIGURE_NOTIFY

event_type_xcb(::Union{ButtonPressed, ButtonReleased}) = xcb_button_press_event_t
event_type_xcb(::Union{KeyPressed, KeyReleased}) = xcb_key_press_event_t
event_type_xcb(::PointerMoves) = xcb_motion_notify_event_t
event_type_xcb(::Union{PointerEntersWindow, PointerLeavesWindow}) = xcb_enter_notify_event_t
event_type_xcb(::ResizeEvent) = xcb_configure_notify_event_t
event_type_xcb(::ExposeEvent) = xcb_expose_event_t

button_xcb(::ButtonLeft) = XCB_BUTTON_INDEX_1
button_xcb(::ButtonMiddle) = XCB_BUTTON_INDEX_2
button_xcb(::ButtonRight) = XCB_BUTTON_INDEX_3
button_xcb(::ButtonScrollUp) = XCB_BUTTON_INDEX_4
button_xcb(::ButtonScrollDown) = XCB_BUTTON_INDEX_5

detail_xcb(wh::XWindowHandler, e::MouseEvent) = button_xcb(e.button)
detail_xcb(wh::XWindowHandler, e::PointerMovesEvent) = XCB_MOTION_NORMAL
detail_xcb(wh::XWindowHandler, e::PointerEntersWindowEvent) = XCB_ENTER_NOTIFY
detail_xcb(wh::XWindowHandler, e::PointerLeavesWindowEvent) = XCB_LEAVE_NOTIFY
detail_xcb(wh::XWindowHandler, e::EventData) = 0
detail_xcb(wh::XWindowHandler, e::EventDetails) = detail_xcb(wh, e.data)
detail_xcb(wh::XWindowHandler, e::EventDetails{<:KeyEvent}) = keycode(wh, e)

event_xcb(wh::XWindowHandler, e::EventDetails) = event_type_xcb(action(e))(
    response_type_xcb(action(e)),
    detail_xcb(wh, e),
    0,
    e.time,
    e.win.parent_id,
    e.win.id,
    0,
    0,
    0,
    e.location...,
    state_xcb(e.data),
    true,
    false
)

event_xcb(wh::XWindowHandler, e::EventDetails{<:ExposeEvent}) = event_type_xcb(action(e))(
    response_type_xcb(action(e)),
    0,
    0,
    e.win.id,
    e.location...,
    extent(e.win)...,
    0,
    0
)

event_xcb(wh::XWindowHandler, e::EventDetails{<:ResizeEvent}) = event_type_xcb(action(e))(
    response_type_xcb(action(e)),
    0,
    0,
    e.win.id,
    e.win.id,
    0,
    e.location...,
    extent(e.win)...,
    0,
    0,
    0
)

send_event(wh::XWindowHandler, e::EventDetails) = send_event(e.win, event_xcb(wh, e))

function send_event(win::XCBWindow, event)
    ref = Ref(event)
    GC.@preserve ref begin
        event_ptr = Base.reinterpret(Cstring, Base.unsafe_convert(Ptr{typeof(event)}, ref))
        @flush @check :error xcb_send_event(win.conn, false, win.id, 0, event_ptr)
    end
end

hex(x) = "0x$(string(x, base=16))"

function send(wh::XWindowHandler, win::XCBWindow; location=(0, 0))
    (event; location=location) -> send_event(wh, EventDetails(event, location, floor(time()), win))
end


text:
abstract type AbstractMetric end

apply(::Type{M}, args...; kwargs...) where {M <: AbstractMetric} =
    apply(M, ConfusionMatrix(args...); kwargs...)

apply(::Type{M}, x::AbstractArray{<:ConfusionMatrix}; kwargs...) where {M <: AbstractMetric} =
    apply.(M, x; kwargs...)

macro metric(name)
    name_lw = Symbol(lowercase(string(name)))

    quote
        abstract type $(esc(name)) <: AbstractMetric end

        Base.@__doc__  function $(esc(name_lw))(args...; kwargs...)
            apply($(esc(name)), args...; kwargs...)
        end
    end
end

"""
    $(SIGNATURES)

Returns # true positive samples.
"""
@metric True_positive
apply(::Type{True_positive}, x::ConfusionMatrix) = x.tp

"""
    $(SIGNATURES)

Returns # true negative samples.
"""
@metric True_negative
apply(::Type{True_negative}, x::ConfusionMatrix) = x.tn

"""
    $(SIGNATURES)

Returns # false positive samples.
"""
@metric False_positive
apply(::Type{False_positive}, x::ConfusionMatrix) = x.fp

"""
    $(SIGNATURES)

Returns # false negative samples.
"""
@metric False_negative
apply(::Type{False_negative}, x::ConfusionMatrix) = x.fn

"""
    $(SIGNATURES)

Returns true positive rate `tp/p`.
Aliases: `sensitivity`,  `recall`, `hit_rate`.
"""
@metric True_positive_rate
apply(::Type{True_positive_rate}, x::ConfusionMatrix) = x.tp/x.p

const sensitivity = true_positive_rate
const recall      = true_positive_rate
const hit_rate    = true_positive_rate

"""
    $(SIGNATURES)

Returns true negative rate `tn/n`.
Aliases: `specificity`,  `selectivity`.
"""
@metric True_negative_rate
apply(::Type{True_negative_rate}, x::ConfusionMatrix) = x.tn/x.n

const specificity = true_negative_rate
const selectivity = true_negative_rate

"""
    $(SIGNATURES)

Returns false positive rate `fp/n`.
Aliases: `fall_out`, `type_I_error`.
"""
@metric False_positive_rate
apply(::Type{False_positive_rate}, x::ConfusionMatrix) = x.fp/x.n

const fall_out     = false_positive_rate
const type_I_error = false_positive_rate

"""
    $(SIGNATURES)

Returns false negative rate `fn/p`.
Aliases: `miss_rate`, `type_II_error`.
"""
@metric False_negative_rate
apply(::Type{False_negative_rate}, x::ConfusionMatrix) = x.fn/x.p

const miss_rate     = false_negative_rate
const type_II_error = false_negative_rate

"""
    $(SIGNATURES)

Returns precision `tp/(tp + fp)`.
Aliases: `positive_predictive_value`.
"""
@metric Precision
apply(::Type{Precision}, x::ConfusionMatrix) = (val = x.tp/(x.tp + x.fp); isnan(val) ? one(val) : val)

const positive_predictive_value = precision

"""
    $(SIGNATURES)

Returns negative predictive value `tn/(tn + fn)`.
"""
@metric Negative_predictive_value
apply(::Type{Negative_predictive_value}, x::ConfusionMatrix) = x.tn/(x.tn + x.fn)

"""
    $(SIGNATURES)

Returns false discovery rate `fp/(fp + tp)`.
"""
@metric False_discovery_rate
apply(::Type{False_discovery_rate}, x::ConfusionMatrix) = x.fp/(x.fp + x.tp)

"""
    $(SIGNATURES)

Returns false omission rate `fn/(fn + tn)`.
"""
@metric False_omission_rate
apply(::Type{False_omission_rate}, x::ConfusionMatrix) = x.fn/(x.fn + x.tn)

"""
    $(SIGNATURES)

Returns threat score `tp/(tp + fn + fp)`.
Aliases: `critical_success_index`.
"""
@metric Threat_score
apply(::Type{Threat_score}, x::ConfusionMatrix) = x.tp/(x.tp + x.fn + x.fp)

const critical_success_index = threat_score

"""
    $(SIGNATURES)

Returns accuracy `(tp + tn)/(p + n).
"""
@metric Accuracy
apply(::Type{Accuracy}, x::ConfusionMatrix) = (x.tp + x.tn)/(x.p + x.n)

"""
    $(SIGNATURES)

Returns balanced accuracy `(tpr + tnr)/2`.
"""
@metric Balanced_accuracy
apply(::Type{Balanced_accuracy}, x::ConfusionMatrix) =
    (true_positive_rate(x) + true_negative_rate(x))/2

"""
    $(SIGNATURES)

Returns error rate `1 - accuracy`.
"""
@metric Error_rate
apply(::Type{Error_rate}, x::ConfusionMatrix) = 1 - accuracy(x)

"""
    $(SIGNATURES)

Returns balanced error rate `1 - balanced_accuracy`.
"""
@metric Balanced_error_rate
apply(::Type{Balanced_error_rate}, x::ConfusionMatrix) = 1 - balanced_accuracy(x)

"""
    $(SIGNATURES)

Returns f1 score `2*precision*recall/(precision + recall)`.
"""
@metric F1_score
apply(::Type{F1_score}, x::ConfusionMatrix) =
    2*precision(x)*recall(x)/(precision(x) + recall(x))

"""
    $(SIGNATURES)

Returns fβ score `(1 + β^2)*precision*recall/(β^2*precision + recall)`.
"""
@metric Fβ_score
apply(::Type{Fβ_score}, x::ConfusionMatrix; β::Real = 1) =
    (1 + β^2)*precision(x)*recall(x)/(β^2*precision(x) + recall(x))

"""
    $(SIGNATURES)

Returns Matthews correlation coefficient `(tp*tn - fp*fn)/sqrt((tp+fp)(tp+fn)(tn+fp)(tn+fn))`.
Aliases: ` mcc`.
"""
@metric Matthews_correlation_coefficient
apply(::Type{Matthews_correlation_coefficient}, x::ConfusionMatrix) =
    (x.tp*x.tn - x.fp*x.fn)/sqrt((x.tp + x.fp)*(x.tp + x.fn)*(x.tn + x.fp)*(x.tn + x.fn))

const mcc = matthews_correlation_coefficient

"""
    $(SIGNATURES)

Returns quant `(fn + tn)/(p + n)`.
"""
@metric Quant
apply(::Type{Quant}, x::ConfusionMatrix) = (x.fn + x.tn)/(x.p + x.n)

"""
    $(SIGNATURES)

Returns topquant `1 - quant`.
"""
@metric Topquant
apply(::Type{Topquant}, x::ConfusionMatrix) = 1 - quant(x)

"""
    $(SIGNATURES)

Returns positive likelihood ratio `tpr/fpr`.
"""
@metric Positive_likelihood_ratio
apply(::Type{Positive_likelihood_ratio}, x::ConfusionMatrix) =
    true_positive_rate(x)/false_positive_rate(x)

"""
    $(SIGNATURES)

Returns negative likelihood ratio `fnr/tnr`.
"""
@metric Negative_likelihood_ratio
apply(::Type{Negative_likelihood_ratio}, x::ConfusionMatrix) =
    false_negative_rate(x)/true_negative_rate(x)

"""
    $(SIGNATURES)

Returns diagnostic odds ratio `tpr*tnr/(fpr*fnr)`.
"""
@metric Diagnostic_odds_ratio
apply(::Type{Diagnostic_odds_ratio}, x::ConfusionMatrix) =
    true_positive_rate(x)*true_negative_rate(x)/(false_positive_rate(x)*false_negative_rate(x))

"""
    $(SIGNATURES)

Returns prevalence `p/(p + n)`.
"""
@metric Prevalence
apply(::Type{Prevalence}, x::ConfusionMatrix) = x.p/(x.p + x.n)