text:
############################################################
# joKron ###################################################
############################################################

############################################################
## outer constructors

"""
    julia> op = joKron(op1[,op2][,...];[name=...])

Kronecker product of serial JOLI operators

# Signature

    joKron(ops::joAbstractLinearOperator...;name::String="joKron")

# Arguments

- `op#`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `name`: custom name

# Notes

- the domain and range types of joKron are equal respectively to domain type of rightmost operator and range type of leftmost operator
- all operators in the chain must have consistent passing domain/range types, i.e. domain type of operator on the left have to be the same as range type of operator on the right

# Example

define operators

    a=rand(ComplexF64,6,4);
    A=joMatrix(a;name="A")

    b=rand(ComplexF64,6,8);
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64,name="B")

    c=rand(ComplexF64,6,4);
    C=joMatrix(c;DDT=ComplexF64,RDT=ComplexF32,name="C")

define Kronecker product

    K=joKron(A,B,C)

"""
function joKron(ops::joAbstractLinearOperator...;name::String="joKron")

    isempty(ops) && throw(joKronException("empty argument list"))
    l=length(ops)
    for i=2:l
        reltype(ops[i])==deltype(ops[i-1]) || throw(joKronException("domain/range type mismatch for $i operator"))
    end
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        ms[i]=ops[i].m
        ns[i]=ops[i].n
    end
    m=prod(ms)
    n=prod(ns)
    for i=1:l
        push!(fops,ops[i])
        push!(fops_T,transpose(ops[i]))
        push!(fops_A,adjoint(ops[i]))
        push!(fops_C,conj(ops[i]))
    end
    return joKron{deltype(fops[l]),reltype(fops[1])}(name*"($l)",m,n,l,ms,ns,false,
                 fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end

############################
## overloaded Base functions

# show(jo)
function show(A::joKron)
    println("Type: $(typeof(A).name)")
    println("Name: $(A.name)")
    println("Size: $(size(A))")
    println(" DDT: $(deltype(A))")
    println(" RDT: $(reltype(A))")
    println(" Flipped: $(A.flip)")
    println("# of ops: $(A.l)")
    for i=1:A.l
        a=A.fop[i]
        println(" $i Type: $(typeof(a).name)")
        println("    Name: $(a.name)")
        println("    Size: $(size(a))")
        println("     DDT: $(deltype(a))")
        println("     RDT: $(reltype(a))")
    end
end

# conj(jo)
conj(A::joKron{DDT,RDT}) where {DDT,RDT} =
    joKron{DDT,RDT}("(conj("*A.name*"))",
        A.m,A.n,A.l,A.ms,A.ns,A.flip,
        A.fop_C,A.fop_A,A.fop_T,A.fop,
        A.iop_C,A.iop_A,A.iop_T,A.iop)

# transpose(jo)
transpose(A::joKron{DDT,RDT}) where {DDT,RDT} =
    joKron{RDT,DDT}("(transpose("*A.name*"))",
        A.n,A.m,A.l,A.ns,A.ms,!A.flip,
        A.fop_T,A.fop,A.fop_C,A.fop_A,
        A.iop_T,A.iop,A.iop_C,A.iop_A)

# adjoint(jo)
adjoint(A::joKron{DDT,RDT}) where {DDT,RDT} =
    joKron{RDT,DDT}("(adjoint("*A.name*"))",
        A.n,A.m,A.l,A.ns,A.ms,!A.flip,
        A.fop_A,A.fop_C,A.fop,A.fop_T,
        A.iop_A,A.iop_C,A.iop,A.iop_T)

# *(jo,vec)
function *(A::joKron{ADDT,ARDT},v::LocalVector{ADDT}) where {ADDT,ARDT}
    size(A,2) == size(v,1) || throw(joKronException("shape mismatch"))
    ksz=reverse(A.ns)
    V=reshape(v,Tuple(ksz))
    p=[x for x in 1:A.l]
    if A.flip
        p=circshift(p,1)
        for i=1:1:A.l
            ksz=circshift(ksz,1)
            V=permutedims(V,p)
            V=reshape(V,Tuple([ksz[1],prod(ksz[2:length(ksz)])]))
            V=A.fop[i]*V
            ksz[1]=A.fop[i].m
            V=reshape(V,Tuple(ksz))
        end
    else
        p=circshift(p,-1)
        for i=A.l:-1:1
            V=reshape(V,Tuple([ksz[1],prod(ksz[2:length(ksz)])]))
            V=A.fop[i]*V
            ksz[1]=A.fop[i].m
            V=reshape(V,Tuple(ksz))
            V=permutedims(V,p)
            ksz=circshift(ksz,-1)
        end
    end
    return vec(V)
end

# *(jo,mvec)
function *(A::joKron{ADDT,ARDT},mv::LocalMatrix{ADDT}) where {ADDT,ARDT}
    size(A, 2) == size(mv, 1) || throw(joKronException("shape mismatch"))
    MV=Matrix{ARDT}(undef,A.m,size(mv,2))
    for i=1:size(mv,2)
        MV[:,i]=A*mv[:,i]
    end
    return MV
end

# -(jo)
function -(A::joKron{DDT,RDT}) where {DDT,RDT}
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:A.l
        push!(fops,-A.fop[i])
        push!(fops_T,-A.fop_T[i])
        push!(fops_A,-A.fop_A[i])
        push!(fops_C,-A.fop_C[i])
    end
    return joKron{DDT,RDT}("(-"*A.name*")",
        A.m,A.n,A.l,A.ms,A.ns,A.flip,
        fops,fops_T,fops_A,fops_C,
        A.iop,A.iop_T,A.iop_A,A.iop_C)
end


text:
using QuantumComputer
using Random
using StatsBase
using Test

@testset "Register" begin
  register = QuantumComputer.Register(1)
  @test register.qubits == [1 0]

  register = QuantumComputer.Register(1, 1)
  @test register.qubits == [0 1]

  register = QuantumComputer.Register(2, 2)
  @test register.qubits == [0 1; 1 0]

  register = QuantumComputer.Register(10)
  @test size(register.qubits) == (10, 2)
end

@testset "Superposition" begin
  register = QuantumComputer.Register(3, 7)
  superposition = QuantumComputer.Superposition(register.qubits)
  @test length(superposition.state) == 2^3
  @test superposition.state[8] == 1
  @test sum(superposition.state) == 1
end

@testset "Measurement" begin
  @test_throws DomainError QuantumComputer.Measurement([1, 1], [1])
  @test_throws DomainError QuantumComputer.Measurement([1], [1, 1])
  @test_throws DomainError QuantumComputer.Measurement([1, 1], [1, 2])
  @test_throws DomainError QuantumComputer.Measurement([1, 2], [1, 1])
end

@testset "qubits_operated_on_by_unitary_matrix" begin
  identity_matrix = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
  @test QuantumComputer.qubits_operated_on_by_unitary_matrix(identity_matrix) == 2
  @test_throws DomainError QuantumComputer.qubits_operated_on_by_unitary_matrix(identity_matrix[1:3, :])
  @test_throws DomainError QuantumComputer.qubits_operated_on_by_unitary_matrix(identity_matrix[:, 1:3])
end

@testset "Gate" begin
  identity_matrix = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
  gate::QuantumComputer.Gate = QuantumComputer.Gate(identity_matrix)
  @test gate.matrix == identity_matrix
  @test gate.superposed_qubits_required == 2
end

@testset "gate_x" begin
  state::Array{Complex{Float64}, 1} = [1, 0]
  @test QuantumComputer.gate_x.matrix * state == [0, 1]

  state = [0, 1]
  @test QuantumComputer.gate_x.matrix * state == [1, 0]
end

@testset "gate_cx" begin
  gate::QuantumComputer.Gate = QuantumComputer.gate_cx(1, 3, 3)

  qubits::Matrix{Complex{Float64}} = [1 0; 1 0; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [1 0; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [0 1; 1 0; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0, 0, 0, 0, 0, 1, 0, 0]

  qubits = [0 1; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0, 0, 0, 0, 1, 0, 0, 0]

  gate = QuantumComputer.gate_cx(1, 2, 3)

  qubits = [0 1; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0, 0, 0, 0, 0, 0, 0, 1]
end

@testset "gate_cnx" begin
  gate::QuantumComputer.Gate = QuantumComputer.gate_cnx([1, 2], 3, 3)

  qubits::Matrix{Complex{Float64}} = [0 1; 1 0; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [1 0; 0 1; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [1 0; 1 0; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [0 1; 0 1; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0, 0, 0, 0, 0, 0, 0, 1]
end

@testset "gate_swap" begin
  gate = QuantumComputer.gate_swap(1, 3, 3)

  qubits::Matrix{Complex{Float64}} = [1 0; 1 0; 1 0]
  state::Array{Complex{Float64}} = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [0 1; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [0 1; 1 0; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0, 1, 0, 0, 0, 0, 0, 0]

  qubits = [1 0; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0, 0, 0, 0, 1, 0, 0, 0]
end

@testset "gate_multi_control" begin
  inner_gate = QuantumComputer.gate_swap(1, 3, 3)
  gate = QuantumComputer.gate_multi_control(inner_gate, 2, 5)

  qubits::Matrix{Complex{Float64}} = [1 0; 0 1; 0 1; 1 0; 1 0]
  state::Array{Complex{Float64}} = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [0 1; 1 0; 1 0; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == state

  qubits = [0 1; 0 1; 0 1; 1 0; 1 0]
  state = QuantumComputer.qubit_tensor_product(qubits)
  expected_qubits::Matrix{Complex{Float64}} = [0 1; 0 1; 1 0; 1 0; 0 1]
  expected_state::Array{Complex{Float64}} = QuantumComputer.qubit_tensor_product(expected_qubits)
  @test gate.matrix * state == expected_state

  qubits = [0 1; 0 1; 1 0; 0 1; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  expected_qubits = [0 1; 0 1; 0 1; 0 1; 1 0]
  expected_state = QuantumComputer.qubit_tensor_product(expected_qubits)
  @test gate.matrix * state == expected_state

  extended_gate = QuantumComputer.gate_extension(gate, 2, 6)

  qubits = [1 0; 0 1; 0 1; 1 0; 0 1; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  expected_qubits = [1 0; 0 1; 0 1; 0 1; 0 1; 1 0]
  expected_state = QuantumComputer.qubit_tensor_product(expected_qubits)
  @test extended_gate.matrix * state == expected_state

  qubits = [1 0; 0 1; 1 0; 1 0; 1 0; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test extended_gate.matrix * state == state
end

@testset "gate_fourier_transform" begin
  gate = QuantumComputer.gate_fourier_transform(2)

  qubits::Matrix{Complex{Float64}} = [1 0; 1 0]
  state::Array{Complex{Float64}} = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state == [0.5, 0.5, 0.5, 0.5]

  qubits = [0 1; 0 1]
  state = QuantumComputer.qubit_tensor_product(qubits)
  @test gate.matrix * state ≈ [0.5, -0.5im, -0.5, 0.5im] atol=0.000000000000001
end

@testset "circuit_convert_to_gate" begin
  initial_value = 13
  constant = 12
  qubit_count = 5

  register = QuantumComputer.Register(qubit_count, initial_value)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count)
  superposition = QuantumComputer.Superposition(register.qubits)
  adder = QuantumComputer.Circuits.constant_adder(constant, qubit_count)
  gate = QuantumComputer.circuit_convert_to_gate(adder)
  measurement = QuantumComputer.Measurement(Array(1:qubit_count), Array(1:qubit_count))

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_gate_to_circuit!(circuit, gate)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  @test classical_register.value == (initial_value + constant) % 2^qubit_count
end

function test_constant_adder(initial_value, constant, qubit_count)
  register = QuantumComputer.Register(qubit_count, initial_value)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count)
  superposition = QuantumComputer.Superposition(register.qubits)
  adder = QuantumComputer.Circuits.constant_adder(constant, qubit_count)
  measurement = QuantumComputer.Measurement(Array(1:qubit_count), Array(1:qubit_count))

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, adder)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  classical_register.value == (initial_value + constant) % 2^qubit_count
end

@testset "constant_adder" begin
  @test test_constant_adder(0, 3, 2)
  @test test_constant_adder(1, 1, 2)
  @test test_constant_adder(0, 63, 6)
  @test test_constant_adder(31, 31, 6)
end

function test_constant_adder_core(initial_value, constant, qubit_count)
  register = QuantumComputer.Register(qubit_count, initial_value)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count)
  superposition = QuantumComputer.Superposition(register.qubits)
  qft = QuantumComputer.gate_fourier_transform(qubit_count)
  adder_core = QuantumComputer.Circuits.constant_adder_core(constant, qubit_count)
  inverse_qft = QuantumComputer.gate_invert(qft)
  measurement = QuantumComputer.Measurement(Array(1:qubit_count), Array(1:qubit_count))

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_gate_to_circuit!(circuit, qft)
  QuantumComputer.add_subcircuit_to_circuit!(circuit, adder_core)
  QuantumComputer.add_gate_to_circuit!(circuit, inverse_qft)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  classical_register.value == (initial_value + constant) % 2^qubit_count
end

@testset "constant_adder_core" begin
  @test test_constant_adder_core(0, 3, 2)
  @test test_constant_adder_core(1, 1, 2)
  @test test_constant_adder_core(0, 63, 6)
  @test test_constant_adder_core(31, 31, 6)
end

@testset "shor2n3_controlled_controlled_modular_adder" begin
  a = 5
  n = 9

  modular_adder = QuantumComputer.Circuits.shor2n3_controlled_controlled_modular_adder(n, a)
  gate_qft = QuantumComputer.gate_fourier_transform(5)
  gate_inverse_qft = QuantumComputer.gate_invert(gate_qft)
  qft = QuantumComputer.gate_extension(gate_qft, 7, 11)
  inverse_qft = QuantumComputer.gate_extension(gate_inverse_qft, 7, 11)
  measurement = QuantumComputer.Measurement(Array(8:11), Array(8:11))

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_gate_to_circuit!(circuit, qft)
  QuantumComputer.add_subcircuit_to_circuit!(circuit, modular_adder)
  QuantumComputer.add_gate_to_circuit!(circuit, inverse_qft)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  classical_register = QuantumComputer.ClassicalRegister(11)

  qubits::Matrix{Complex{Float64}} = [0 1; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 1

  qubits = [0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 5

  qubits = [1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 5
end

@testset "shor2n3_controlled_modular_multiplier" begin
  a = 5
  n = 9

  modular_multiplier = QuantumComputer.Circuits.shor2n3_controlled_modular_multiplier(n, a)
  measurement = QuantumComputer.Measurement(Array(8:11), Array(8:11), 1)
  circuit = QuantumComputer.Circuit()

  QuantumComputer.add_subcircuit_to_circuit!(circuit, modular_multiplier)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  classical_register = QuantumComputer.ClassicalRegister(11)

  qubits::Matrix{Complex{Float64}} = [1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 5

  qubits = [1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 5

  qubits = [0 1; 1 0; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 1

  qubits = [0 1; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 0 1]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 6
end

@testset "shor2n3_controlled_ua" begin
  a = 5
  n = 9

  classical_register = QuantumComputer.ClassicalRegister(11)
  modular_multiplier = QuantumComputer.Circuits.shor2n3_controlled_ua(n, a)

  # first measure that we take no action when the control is not active
  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, modular_multiplier)
  measurement = QuantumComputer.Measurement(Array(8:11), Array(2:5))
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  qubits::Matrix{Complex{Float64}} = [1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 1

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, modular_multiplier)
  measurement = QuantumComputer.Measurement(Array(8:11), Array(8:11))
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  qubits = [1 0; 1 0; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 0

  # next we test that b is reset to 0 after the operation
  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, modular_multiplier)
  measurement = QuantumComputer.Measurement(Array(8:11), Array(8:11))
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  qubits = [0 1; 1 0; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 0

  # now we test that that x is multiplied by a
  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, modular_multiplier)
  measurement = QuantumComputer.Measurement(Array(8:11), Array(2:5))
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  qubits = [0 1; 1 0; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 5

  # finally we test that the modulus is applied
  qubits = [0 1; 1 0; 1 0; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 1

  # one more for good measure
  qubits = [0 1; 1 0; 1 0; 0 1; 0 1; 1 0; 1 0; 1 0; 1 0; 1 0; 1 0]
  superposition = QuantumComputer.Superposition(qubits)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == 6
end

@testset "11^x mod 13 period finder" begin
  qubit_count = 5
  bit_count = 3
  sample_size = 256

  register = QuantumComputer.Register(qubit_count)
  classical_register = QuantumComputer.ClassicalRegister(bit_count)
  superposition = QuantumComputer.Superposition(register.qubits)
  period_finder = QuantumComputer.Circuits.period_finding_for_11x_mod_15()
  measurement = QuantumComputer.Measurement([1, 2, 3], [3, 4, 5], sample_size)

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, period_finder)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  Random.seed!(0)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  # this circuit doesn't find the period, but instead a related value (4).
  # here we test that over 256 samples, at least 128 resulted in binary '100'.
  # success of measurement depends on the random seed above
  @test countmap(measurement.samples)[4] > sample_size / 2
  @test classical_register.value == 4
  @test sum(map(x -> abs(x^2), superposition.state)) ≈ 1 atol=0.00000000001
end

@testset "grover" begin
  qubit_count = 10
  # should work for any valid input
  target_input = rand(0:(2^qubit_count - 1))

  circuit_oracle = QuantumComputer.Circuits.identification_oracle(qubit_count + 1, target_input)
  oracle = QuantumComputer.circuit_convert_to_gate(circuit_oracle)
  grover = QuantumComputer.Circuits.grover(oracle)

  register = QuantumComputer.Register(qubit_count + 1, 1)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count, 0)
  superposition = QuantumComputer.Superposition(register.qubits)
  measurement = QuantumComputer.Measurement(Array(1:qubit_count), Array(1:qubit_count))

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, grover)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)

  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)
  @test classical_register.value == target_input
end

@testset "deutsch-jozsa, constant oracle f(x) = 1" begin
  qubit_count = 10

  oracle_circuit = QuantumComputer.Circuits.constant_oracle(qubit_count, true)
  oracle = QuantumComputer.circuit_convert_to_gate(oracle_circuit)
  deutsch_jozsa = QuantumComputer.Circuits.deutsch_jozsa(oracle)
  measurement = QuantumComputer.Measurement(Array(1:(qubit_count-1)), Array(1:(qubit_count-1)))

  register = QuantumComputer.Register(qubit_count, 1)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count - 1, 0)
  superposition = QuantumComputer.Superposition(register.qubits)

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, deutsch_jozsa)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  @test classical_register.value == 0
end

@testset "deutsch-jozsa, constant oracle f(x) = 0" begin
  qubit_count = 10

  oracle_circuit = QuantumComputer.Circuits.constant_oracle(qubit_count, false)
  oracle = QuantumComputer.circuit_convert_to_gate(oracle_circuit)
  deutsch_jozsa = QuantumComputer.Circuits.deutsch_jozsa(oracle)
  measurement = QuantumComputer.Measurement(Array(1:(qubit_count-1)), Array(1:(qubit_count-1)))

  register = QuantumComputer.Register(qubit_count, 1)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count - 1, 0)
  superposition = QuantumComputer.Superposition(register.qubits)

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, deutsch_jozsa)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  @test classical_register.value == 0
end

@testset "deutsch-jozsa, balanced oracle" begin
  qubit_count = 10

  oracle_circuit = QuantumComputer.Circuits.cx_balanced_oracle(qubit_count)
  oracle = QuantumComputer.circuit_convert_to_gate(oracle_circuit)
  deutsch_jozsa = QuantumComputer.Circuits.deutsch_jozsa(oracle)
  measurement = QuantumComputer.Measurement(Array(1:(qubit_count-1)), Array(1:(qubit_count-1)))

  register = QuantumComputer.Register(qubit_count, 1)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count - 1, 0)
  superposition = QuantumComputer.Superposition(register.qubits)

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, deutsch_jozsa)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  @test classical_register.value != 0
end

@testset "bernstein-vazirani" begin
  qubit_count = 10
  secret::Int64 = rand(0:(2^(qubit_count-1)-1))

  oracle_circuit = QuantumComputer.Circuits.bitwise_product_oracle(secret, qubit_count)
  oracle = QuantumComputer.circuit_convert_to_gate(oracle_circuit)
  bernstein_vazirani = QuantumComputer.Circuits.bernstein_vazirani(oracle)
  measurement = QuantumComputer.Measurement(Array(1:(qubit_count-1)), Array(1:(qubit_count-1)))

  register = QuantumComputer.Register(qubit_count, 0)
  classical_register = QuantumComputer.ClassicalRegister(qubit_count - 1, 0)
  superposition = QuantumComputer.Superposition(register.qubits)

  circuit = QuantumComputer.Circuit()
  QuantumComputer.add_subcircuit_to_circuit!(circuit, bernstein_vazirani)
  QuantumComputer.add_measurement_to_circuit!(circuit, measurement)
  QuantumComputer.apply_circuit_to_superposition!(superposition, circuit, classical_register)

  @test classical_register.value == secret
end

ext: jl


text:
@testset "SimpleGrid" begin
    @polyvar x y z w
    f = x^2+y^2+1
    F = AmoebaFiber2D(x^2+y^2+1)
    grid = Grid2D((-2, 2), 50)
    generator = PolynomialAmoebas.startvalue_generator(F)
    B = PolynomialAmoebas.simple_grid(F, grid, generator)
    @test any(B.data)
    B = PolynomialAmoebas.amoeba(f, grid=grid, alg=Simple())
    @test any(B.data)

    I = PolynomialAmoebas.imaginary_projection(f, grid=grid, alg=Simple())
    @test any(I.data)

    f = x^2+y^2+z^2+1
    F = AmoebaFiber3D(f)
    grid = Grid3D(xlims=(-2, 2), ylims=(-2, 2), zlims=(-2, 2), res=(21, 21, 21))
    generator = PolynomialAmoebas.startvalue_generator(F)
    B = PolynomialAmoebas.simple_grid(F, grid, generator)
    @test any(B.data)
    B = PolynomialAmoebas.amoeba(f, grid=grid, alg=Simple())
    @test any(B.data)

    f = 1 + z * cis(0.8π) + z^3*w*cis(1.6π) + z*w^2*cis(0.4π) + w*cis(1.2π)
    C = CoamoebaFiber2D(f)
    grid = Grid2D((0.0, 2π), 50)
    generator = PolynomialAmoebas.DomainStartValueGenerator2D((-5, 5), (-5, 5))
    B = simple_grid(C, grid, generator)
    @test any(B.data)
    B = coamoeba(f, alg=Simple())
    @test any(B.data)
end


text:
abstract type Kernel end
struct Linear <: Kernel end
struct Gaussian <: Kernel end

struct KernelType{T<:Real, K}
    γ::T
    scale::Bool
    mlkernels::Bool
    precomputed::Bool

    function KernelType(
        K::Type{<:Kernel};
        γ = 1,
        scale = true,
        T = Float32,
        mlkernels = true,
        precomputed = true,
    )

        return new{T, K}(γ, scale, mlkernels, precomputed)
    end
end

function Base.show(io::IO, ker::KernelType{T, K}) where {T, K<:Linear}
    print(io, "Linear(pre = $(ker.precomputed), ml = $(ker.mlkernels))")
end

function Base.show(io::IO, ker::KernelType{T, K}) where {T, K<:Gaussian}
    print(io, "Gaussian($(ker.γ), scale = $(ker.scale), pre = $(ker.precomputed), ml = $(ker.mlkernels))")
end

function init(ker::KernelType{T, K}, d) where {T, K}
    γ = compute_gamma(ker, d)
    return ker.mlkernels ? init_ml(K, γ, T) : init_kf(K, γ, T)
end

function compute_gamma(ker::KernelType{T, K}, d) where {T, K}
    return ker.scale ? T(ker.γ/d) : T(ker.γ)
end

# KernelFunctions.jl
function init_kf(::Type{Gaussian}, γ, T)
    return KernelFunctions.TransformedKernel(
        KernelFunctions.SqExponentialKernel(),
        KernelFunctions.ScaleTransform(T(sqrt(2)*sqrt(γ))),
    )
end

init_kf(::Type{Linear}, γ, T) = KernelFunctions.LinearKernel()

function kernelmatrix(kernel::KernelFunctions.Kernel, X)
    return KernelFunctions.kernelmatrix(kernel, X; obsdim = 1)
end

function kernelmatrix(kernel::KernelFunctions.Kernel, X, Y)
    return KernelFunctions.kernelmatrix(kernel, X, Y; obsdim = 1)
end

function kernelmatrix_diag(kernel::KernelFunctions.Kernel, X)
    return KernelFunctions.kernelmatrix_diag(kernel, X; obsdim = 1)
end

# MLKernels.jl
init_ml(::Type{Gaussian}, γ, T) = MLKernels.GaussianKernel{T}(γ)
init_ml(::Type{Linear}, γ, T) = MLKernels.PolynomialKernel{T}(2, 0, 1)

function kernelmatrix(kernel::MLKernels.Kernel, X)
    return MLKernels.kernelmatrix(Val(:row), kernel, X, true)
end

function kernelmatrix(kernel::MLKernels.Kernel, X, Y)
    return MLKernels.kernelmatrix(Val(:row), kernel, X, Y)
end

function kernelmatrix_diag(kernel::MLKernels.Kernel, X)
    return map(x -> MLKernels.kernel(kernel, x, x), eachrow(X))
end

# KernelMatrix
struct KernelMatrix{M, T<:Real, D}
    X::Matrix{T}
    y::BitVector
    kernel::D

    # auxilliary vars
    n::Int
    nα::Int
    nβ::Int
    perm::Vector{Int}

    # kernel matrix
    matrix::M

    function KernelMatrix(model, X::Matrix{T}, y, kernel; precomputed = false) where {T}
        y = BitVector(y)
        nα, nβ, perm = permutation(model, y)

        if precomputed
            matrix = kernelmatrix(kernel, X[perm, :])
            matrix[1:nα, (nα+1):end] .*= -1
            matrix[(nα+1):end, 1:nα] .*= -1
        else
            diag = kernelmatrix_diag(kernel, X[perm, :])
            matrix = OnFlyMatrix{T}(0, diag, similar(diag))
        end
        D, M = typeof(kernel), typeof(matrix)
        return new{M, T, D}(X, y, kernel, nα + nβ, nα, nβ, perm, matrix)
    end
end

Base.show(io::IO, K::KernelMatrix) = print(io, "$(K.n)x$(K.n) kernel matrix")
Base.size(K::KernelMatrix) = (K.n, K.n)
Base.eltype(::KernelMatrix{M, T}) where {M, T} = T
Base.getindex(K::KernelMatrix, args...) = getindex(K.matrix, args...)
Base.:*(K::KernelMatrix, s::AbstractVector) = K.matrix * s

inds_α(K::KernelMatrix)  = 1:K.nα
inds_β(K::KernelMatrix)  = (K.nα + 1):K.n

# OnFlyMatrix
mutable struct OnFlyMatrix{T<:Real}
    row_id::Int
    diag::Vector{T}
    row::Vector{T}
end

function computerow!(K::KernelMatrix{<:OnFlyMatrix}, i::Int)
    row = kernelmatrix(K.kernel, K.X[K.perm, :], K.X[K.perm[i:i], :])[:]
    row[i <= K.nα ? inds_β(K) : inds_α(K)] .*= -1

    K.matrix.row .= row
    K.matrix.row_id = i
    return
end

function Base.getindex(K::KernelMatrix{<:OnFlyMatrix}, i::Int, j::Int)
    if i == j
        K.matrix.diag[i]
    else
        id = K.matrix.row_id
        if i == id
            K.matrix.row[j]
        elseif j == id
            K.matrix.row[i]
        else
            computerow!(K, i)
            K.matrix.row[j]
        end
    end
end

function Base.getindex(K::KernelMatrix{<:OnFlyMatrix}, i::Int, ::Colon)
    i == K.matrix.row_id || computerow!(K, i)
    return copy(K.matrix.row)
end

Base.getindex(K::KernelMatrix{<:OnFlyMatrix}, ::Colon, i::Int) = getindex(K, i, :)

function Base.:*(K::KernelMatrix{<:OnFlyMatrix}, s::AbstractVector)
    x = zero(s)
    for (i, si) in enumerate(s)
        iszero(si) && continue
        x .+= si .* K[i, :]
    end
    return x
end


text:
function triangle(n::Int)

end


text:
using DandelionWebSockets: AbstractClientProtocol, SendTextFrame, SendBinaryFrame
using DandelionWebSockets: ClientPingRequest, PongMissed, CloseRequest, SocketClosed
using DandelionWebSockets: FrameFromServer
import DandelionWebSockets: handle

export ClientProtocolProxy

struct ClientProtocolProxy <: AbstractClientProtocol
    channel::Channel{Any}
    clientlogic::AbstractClientProtocol

    function ClientProtocolProxy(clientlogic::AbstractClientProtocol)
        proxy = new(Channel{Any}(Inf), clientlogic)
        @async run_clientlogicproxy(proxy)
        proxy
    end
end

function run_clientlogicproxy(c::ClientProtocolProxy)
    for p in c.channel
        handle(c.clientlogic, p)
    end
end

handle(proxy::ClientProtocolProxy, s::SendTextFrame) = put!(proxy.channel, s)
handle(proxy::ClientProtocolProxy, s::SendBinaryFrame) = put!(proxy.channel, s)
handle(proxy::ClientProtocolProxy, s::ClientPingRequest) = put!(proxy.channel, s)
handle(proxy::ClientProtocolProxy, s::PongMissed) = put!(proxy.channel, s)
handle(proxy::ClientProtocolProxy, s::CloseRequest) = put!(proxy.channel, s)
handle(proxy::ClientProtocolProxy, s::SocketClosed) = put!(proxy.channel, s)
handle(proxy::ClientProtocolProxy, s::FrameFromServer) = put!(proxy.channel, s)

stopproxy(proxy::ClientProtocolProxy) = close(proxy.channel)


text:
@with_kw mutable struct InfoCollectorPOMDP{S, A, O} <: POMDP{S,A,O}
    pomdp::POMDP{S,A,O}
    extract_info
    dataset = []
end

POMDPs.initialstate(mdp::InfoCollectorPOMDP) = initialstate(mdp.pomdp)
POMDPs.initialobs(mdp::InfoCollectorPOMDP, s) = initialobs(mdp.pomdp, s)

POMDPs.actions(mdp::InfoCollectorPOMDP) = actions(mdp.pomdp)
POMDPs.actionindex(mdp::InfoCollectorPOMDP, a) = actionindex(mdp.pomdp)

POMDPs.isterminal(mdp::InfoCollectorPOMDP, s) = isterminal(mdp.pomdp, s)
POMDPs.discount(mdp::InfoCollectorPOMDP) = discount(mdp.pomdp)

POMDPs.observation(mdp::InfoCollectorPOMDP, s) = observation(mdp.pomdp, s)

render(mdp::InfoCollectorPOMDP, s, a = nothing; kwargs...) = render(mdp.pomdp, s, a; kwargs...)
render(mdp::InfoCollectorPOMDP; kwargs...) = render(mdp.pomdp; kwargs...)

function POMDPs.gen(mdp::InfoCollectorPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, o, r = gen(mdp.pomdp, s, a, rng, info=info)
    data_point = mdp.extract_info(info)
    if !isnothing(data_point) && !ismissing(data_point)
        push!(mdp.dataset, data_point)
    end
    return (sp=sp, o=o, r=r)
end

@with_kw mutable struct InfoCollectorMDP{S, A} <: MDP{S,A}
    mdp::MDP{S,A}
    extract_info
    dataset = []
end

POMDPs.initialstate(mdp::InfoCollectorMDP) = initialstate(mdp.mdp)

POMDPs.actions(mdp::InfoCollectorMDP) = actions(mdp.mdp)
POMDPs.actionindex(mdp::InfoCollectorMDP, a) = actionindex(mdp.mdp)

POMDPs.isterminal(mdp::InfoCollectorMDP, s) = isterminal(mdp.mdp, s)
POMDPs.discount(mdp::InfoCollectorMDP) = discount(mdp.mdp)

render(mdp::InfoCollectorMDP, s, a = nothing; kwargs...) = render(mdp.mdp, s, a; kwargs...)
render(mdp::InfoCollectorMDP; kwargs...) = render(mdp.mdp; kwargs...)

function POMDPs.gen(mdp::InfoCollectorMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, r = gen(mdp.mdp, s, a, rng, info=info)
    data_point = mdp.extract_info(info)
    if !isnothing(data_point) && !ismissing(data_point)
        push!(mdp.dataset, data_point)
    end
    return (sp=sp, r=r)
end

InfoCollector(mdp::T, extract_info) where T <: MDP = InfoCollectorMDP(mdp=mdp, extract_info=extract_info)
InfoCollector(pomdp::T, extract_info) where T <: POMDP = InfoCollectorPOMDP(pomdp=pomdp, extract_info=extract_info)

function clear_dataset(mdp::T) where T <: Union{InfoCollectorMDP, InfoCollectorPOMDP}
    mdp.dataset = []
    return mdp
end


text:
module MCMCChains

import Showoff: showoff
import StatsBase: autocor, autocov, countmap, counts, describe, predict,
       quantile, sample, sem, summarystats, sample, AbstractWeights
import LinearAlgebra: diag
import Serialization: serialize, deserialize
import Base: sort, range, names, get, hash, convert, show, display
import Statistics: cor
import Core.Array
import DataFrames: DataFrame, names, eachcol

using RecipesBase
import RecipesBase: plot

using Serialization
using Distributions
using KernelDensity
using SpecialFunctions
using AxisArrays
const axes = Base.axes

export Chains, getindex, setindex!, chains, setinfo, chainscat
export describe, set_section, get_params, sections
export sample, AbstractWeights
export Array, DataFrame, sort_sections, convert
export summarize, summarystats, ChainDataFrame
export hpd, ess

# export diagnostics functions
export discretediag, gelmandiag, gewekediag, heideldiag, rafterydiag
export autocor

abstract type AbstractChains end

"""
    Chains type

Parameters:

- `value`: An `AxisArray` object with axes `iter` × `var` × `chains`
- `logevidence` : A field containing the logevidence.
- `name_map` : A `NamedTuple` mapping each variable to a section.
- `info` : A `NamedTuple` containing miscellaneous information relevant to the chain.
The `info` field can be set using `setinfo(c::Chains, n::NamedTuple)`.
"""
struct Chains{A, T, K<:NamedTuple, L<:NamedTuple} <: AbstractChains
    value::AxisArray{A,3}
    logevidence::T
    name_map::K
    info::L
end

# imports
include("utils.jl")

include("chains.jl")
include("constructors.jl")
include("summarize.jl")
include("discretediag.jl")
include("fileio.jl")
include("gelmandiag.jl")
include("gewekediag.jl")
include("heideldiag.jl")
include("mcse.jl")
include("rafterydiag.jl")
include("sampling.jl")
include("stats.jl")
include("plot.jl")

end # module


text:
function datafile(filename)
	joinpath(dirname(@__FILE__),"..","data", filename)
end

###### datum parsers

"""
Parses as string of the form #RRGGBB
eg \"#32CD32\" (limegreen) into a color
"""
parse_hexcode(h::AbstractString) = reinterpret(RGB24, parse(UInt32,h[2:end],base=16))

"Parses three strings that are integers between 0-255, and converts to color"
function parse_decimal256(r::AbstractString,
                            g::AbstractString,
                            b::AbstractString)
    ret = 0x00_00_00_00
    ret+=parse(UInt8,r)
    ret<<=8
    ret+=parse(UInt8,g)
    ret<<=8
    ret+=parse(UInt8,b)
    reinterpret(RGB24, ret)
end

"Parses strings in CIE l*a*b, eg \"2.03/-0.03/9.89\""
function parse_CieLab(vals)
    l,a,b = (x->parse(Float64, x)).(split(vals,"/"))
    Lab(l,a,b)
end

######### Formats parsers
"""
Loads an x11 style color list.
Lines are of the form
`r g b name which may have spaces`
the `r`, `g` and `b` terms are numbers between 0 and 256 in decimal
Lines starting with `!` are comments and are skipped
"""
function load_x11_style(path)
	lines = eachline(path)
    data = Dict{String, RGB24}()
    for line in lines
        line[1]=='!'  && continue
        fields = split(line)
        length(fields)==-0 && continue
		r, g, b = fields[1:3]
        name = strip(join(fields[4:end]," "))
        data[name] = parse_decimal256(r,g,b)
    end
    data
end

"""
Load a XCMS color database
See: https://www.x.org/archive/X11R7.5/doc/man/man3/XcmsColor.3.html
"""
function load_xcms_colordb(filename)
    ret = Dict{String, RGB24}()
    lines = eachline(open(filename,"r"))
    for line in lines
        startswith(line, "XCMS_COLORDB_START") && break
    end

    #Now at the data
    for line in lines
        startswith(line, "XCMS_COLORDB_END") && break

        name, colorcode = split(line, '\t'; keepempty=false)
        spec, vals = split(colorcode, ":")
        @assert spec == "CIELab" "Only CIELab format supported right now"
        col = parse_CieLab(vals)

        ret[name] = convert(RGB24, col)
    end
    ret
end

########## Specific formats for specific instances

"""
Load the color list from the [XKCD color survey](https://blog.xkcd.com/2010/05/03/color-survey-results/).
These colors are suitable for use with (modern) monitors
"""
function load_xkcd()
	linefields = (split(line,"\t") for line in eachline(datafile("xkcd.txt")))
	Dict((name => parse_hexcode(code)) for (name,code) in linefields)
end

"""
Load the color list from [Paul Tol's Notes](https://personal.sron.nl/~pault/).
These colors are suitable for use with (modern) monitors
"""
function load_paul_tol()
    linefields = []
    for line in eachline(datafile("paul-tol.txt"))
        startswith(line,";") && continue
        push!(linefields, split(line,"\t"))
    end
	Dict((name => parse_hexcode(code)) for (name,code) in linefields)
end

"""
load the color list from Resene Paints Ltd
These colors are best used with paint.
"""
function load_resene()

    lines = eachline(datafile("resenecolours.txt"))
    for line in lines
        startswith(line, "\"Colour Name\"") && break
    end
    #get to the good part

    data = Dict{String, RGB24}()

    for line in lines
        namefield, rr, gg, bb = split(line,"\t")
        name = namefield[1+length("\"Resene ") : end-1] #skip the end \"
        data[name] = parse_decimal256(rr, gg, bb)
    end
    data
end

########  specific instances of common formats

"""
Loads the NBS-ISCC color dictionary from

    Kenneth L. Kelly and Deanne B. Judd.
    "Color: Universal Language and Dictionary of Names",
    National Bureau of Standards,
    Spec. Publ. 440, Dec. 1976, 189 pages

Featuring the corrections by John Foster
to the orignal computer entry of David Mundie

This is suitable for surface colors
"""
load_nbs() = load_x11_style(datafile("NBS-ISCC-rgb.txt"))

"""
Loads the ``Standard'' X11 color-lists
These colours are good for monitors.
To be explict, these colours are good for Paul Raveling's HP Monitor
"""
load_x11() = load_x11_style(datafile("x11rgb.txt"))

"""
Load the W3C defined "Basic Colors".
In theory these are as safe as it gets.
These should be good everywhere.
"""
load_css3() = load_x11_style(datafile("css3.txt"))

"""
Load the set of colors from Crayola Crayons.
This is not the complete list, but it is extensive.
These should be good for crayons
"""
load_crayola() = load_x11_style(datafile("crayola.txt"))

"""
Loads the Winsor-Newton Water-colours
These should be good for surface colors
"""
load_winsor_newton() = load_xcms_colordb(datafile("winsor-newton-lab.txt"))


text:
# we write a function that returns a root of func in the
# interval [low, high] with relative error tol, assuming such a root exists
function binarySearch(func, low, high, tol)

	# we repeat until the tolorance condition is met
	while abs(high - low)/2 > tol * abs(high + low)/2

		# we define mid as the midpoint of the interval [low, high]
		mid = (low + high)/2

		# we check if the value of func at high and mid have the same sign
		if func(high) * func(mid) > 0

			# and reassign [low, high] accordingly
		    high = mid
		else
		    low = mid
		end
	end

	# we return the avarage of low and high as the approximate root
	return (low + high)/2;
end

max_stars_repo_path: 1a/Root Finding/binarySearch.jl

max_issues_repo_path: 1a/Root Finding/binarySearch.jl

max_forks_repo_path: 1a/Root Finding/binarySearch.jl


text:
using SafeTestsets

@safetestset "Test Cell Maps" begin
    include("maps_and_mesh/test_cell_map.jl")
end

@safetestset "Test CG Mesh" begin
    include("maps_and_mesh/test_1d_cg_mesh.jl")
    include("maps_and_mesh/test_cg_mesh.jl")
end

@safetestset "Test DG Mesh" begin
    include("maps_and_mesh/test_1d_dg_mesh.jl")
    include("maps_and_mesh/test_dg_mesh.jl")
end

@safetestset "Test LevelSet" begin
    include("miscellaneous/test_levelset.jl")
end

@safetestset "Test Cut Mesh" begin
    include("maps_and_mesh/test_cut_mesh.jl")
end

@safetestset "Test Cell Quadratures" begin
     include("quadratures/test_cell_quadratures.jl")
end

@safetestset "Test Face Quadratures" begin
    include("quadratures/test_face_quadratures.jl")
end

@safetestset "Test Interface Quadratures" begin
    include("quadratures/test_interface_quadratures.jl")
end

@safetestset "Test Cell Merging" begin
    include("miscellaneous/test_cell_merging.jl")
end

@safetestset "Test Hooke Stiffness" begin
    include("miscellaneous/test_hooke_stiffness.jl")
end

@safetestset "Test Assembly" begin
    include("miscellaneous/test_assembly.jl")
end

@safetestset "Test Simple Tension" begin
    include("simple_tension/test_simple_tension.jl")
    include("simple_tension/test_cut_simple_tension.jl")
    include("simple_tension/test_merged_simple_tension.jl")
    include("simple_tension/test_transformation_strain_simple_tension.jl")
    include("simple_tension/test_incoherent_interface_simple_tension.jl")
    include("simple_tension/test_incoherent_interface_transformation_simple_tension.jl")
end

@safetestset "Test Vertical Plane Interface L2 Displacement Convergence" begin
    include("standard_convergence/test_vertical_plane_interface_convergence.jl")
end

@safetestset "Test Inclined Plane Interface Displacement Convergence" begin
    include("standard_convergence/test_inclined_plane_interface_convergence.jl")
end

@safetestset "Test Curved Interface Displacement Convergence" begin
    include("standard_convergence/test_curved_interface_convergence.jl")
end

@safetestset "Test Displacement + Traction BC Displacement Convergence" begin
    include(
        "standard_convergence/test_displacement_and_traction_bc_convergence.jl",
    )
end

@safetestset "Test Mixed BC Displacement Convergence" begin
    include("standard_convergence/test_mixed_bc_convergence.jl")
end

@safetestset "Test Stress Convergence" begin
    include("standard_convergence/test_stress_convergence.jl")
end

@safetestset "Test Periodic BC" begin
    include("periodic_bc/test_modify_cell_connectivity.jl")
    include("periodic_bc/test_vertical_periodic_bc_convergence.jl")
end

@safetestset "Test Coherent Interface Transformation Strain Displacement Convergence" begin
    include(
        "transformation_strain_convergence/test_transformation_strain_displacement_convergence.jl",
    )
end

@safetestset "Test Coherent Interface Transformation Strain Displacement + Traction BC Displacement Convergence" begin
    include(
        "transformation_strain_convergence/test_transformation_strain_displacement_traction_bc_convergence.jl",
    )
end

@safetestset "Test Coherent Interface Transformation Strain Mixed BC Displacement Convergence" begin
    include(
        "transformation_strain_convergence/test_transformation_strain_mixed_bc_convergence.jl",
    )
end

@safetestset "Test Coherent Interface Stress Convergence" begin
    include(
        "transformation_strain_convergence/test_transformation_strain_L2_stress_convergence.jl",
    )
end

@safetestset "Test Levelset Reinitialize" begin
    include("levelset_routines/test_seed_zero_levelset.jl")
    include("levelset_routines/test_plane_interface_closest_points.jl")
    include("levelset_routines/test_circular_interface_closest_points.jl")
    include("levelset_routines/test_levelset_reinitialize_convergence.jl")
    include("levelset_routines/test_closest_points_on_merged_mesh.jl")
end

@safetestset "Test Levelset Propagate" begin
    include("levelset_routines/test_levelset_propagate.jl")
    include("levelset_routines/test_levelset_propagate_convergence.jl")
end


text:
"""
Fallback for signals that are not visualizable
"""
function GLVisualize._default(s::Signal, style::GLAbstraction.Style, data)
    GLVisualize._default(map(string, s), style, data)
end

get_pairval(pairtuple, i) = pairtuple[i][2]

"""
Visualization of tuples with kw_args
"""
function GLVisualize.visualize{X, KW_ARGS<:Tuple}(
        v_kw_args::Signal{Tuple{X, KW_ARGS}},
        style::GLAbstraction.Style, data
    )
    to_visualize = map(first, v_kw_args)
    val = value(to_visualize)
    if isa(val, Tuple)
        to_visualize = map(first, to_visualize)
        style = Style{val[2]}()
    end
    kw_args = map(last, v_kw_args)
    for (i, (name, val)) in enumerate(value(kw_args))
        data[name] = const_lift(get_pairval, kw_args, i)
    end
    visualize(to_visualize, style, data)
end

"""
Visualization of a graph with Rects and Lines
"""
function GLVisualize.visualize{R<:SimpleRectangle, L<:LineSegment}(
        rects_lines::Tuple{Vector{R}, Vector{L}},
        style::GLAbstraction.Style, data
    )
    rects, lines = rects_lines
    positions = Point2f0[Point2f0(r.x, r.y) for r in rects]
    scales = Vec2f0[Vec2f0(r.w, r.h) for r in rects]

    rect_vis = visualize((SimpleRectangle, positions), scale=scales,
        color=RGBA{Float32}(0,0,0,0),
        stroke_width=2f0,
        stroke_color=GLVisualize.default(RGBA)
    )
    line_vis = visualize(lines; data...)

    Context(rect_vis, line_vis)
end

function view_signal(signal::Signal, window)
    adjecency_list, signal_nodes = to_adjency_list(signal)
    label_sizes = [SimpleRectangle(0f0, 0f0, 100f0, 100f0) for i=1:length(adjecency_list)]
    rects, lines = layout_tree(adjecency_list, label_sizes)

    graph = visualize((rects, lines))
    view(graph, window)
    for (area, s) in zip(rects, signal_nodes)
        s_area = map(window.cameras[:orthographic_pixel].projectionview, window.inputs[:framebuffer_size]) do pv, fs
            xy = Vec4f0(area.x+2, area.y+2, 0, 1)
            xywh = Vec4f0(xy[1]+area.w-4, xy[2]+area.h-4, 0, 1)
            _xy = round(((Vec2f0(pv * xy) + 1f0)/2f0) .* fs)
            _wh = round(((Vec2f0(pv * xywh) + 1f0)/2f0) .* fs)
            SimpleRectangle{Int}(_xy..., (_wh-_xy)...)
        end
        screen = Screen(window, area=s_area)
        if applicable(vizzedit, s, screen)
            _s, vis = vizzedit(s, screen)
        else
            vis = visualize(s)
        end

        view(vis, screen)
        center!(screen, vis.children[][:preferred_camera])
    end
end


text:
@testset "xlogx, xlogy, and xlog1py" begin
    @test iszero(xlogx(0))
    @test xlogx(2) ≈ 2.0 * log(2.0)
    @test_throws DomainError xlogx(-1)
    @test isnan(xlogx(NaN))

    @test iszero(xlogy(0, 1))
    @test xlogy(2, 3) ≈ 2.0 * log(3.0)
    @test_throws DomainError xlogy(1, -1)
    @test isnan(xlogy(NaN, 2))
    @test isnan(xlogy(2, NaN))
    @test isnan(xlogy(0, NaN))

    @test iszero(xlog1py(0, 0))
    @test xlog1py(2, 3) ≈ 2.0 * log1p(3.0)
    @test_throws DomainError xlog1py(1, -2)
    @test isnan(xlog1py(NaN, 2))
    @test isnan(xlog1py(2, NaN))
    @test isnan(xlog1py(0, NaN))

    # Since we allow complex/negative values, test for them. See comments in:
    # https://github.com/JuliaStats/StatsFuns.jl/pull/95

    @test xlogx(1 + im) == (1 + im) * log(1 + im)
    @test isnan(xlogx(NaN + im))
    @test isnan(xlogx(1 + NaN * im))

    @test xlogy(-2, 3) == -xlogy(2, 3)
    @test xlogy(1 + im, 3) == (1 + im) * log(3)
    @test xlogy(1 + im, 2 + im) == (1 + im) * log(2 + im)
    @test isnan(xlogy(1 + NaN * im, -1 + im))
    @test isnan(xlogy(0, -1 + NaN * im))
    @test isnan(xlogy(Inf + im * NaN, 1))
    @test isnan(xlogy(0 + im * 0, NaN))
    @test iszero(xlogy(0 + im * 0, 0 + im * Inf))

    @test xlog1py(-2, 3) == -xlog1py(2, 3)
    @test xlog1py(1 + im, 3) == (1 + im) * log1p(3)
    @test xlog1py(1 + im, 2 + im) == (1 + im) * log1p(2 + im)
    @test isnan(xlog1py(1 + NaN * im, -1 + im))
    @test isnan(xlog1py(0, -1 + NaN * im))
    @test isnan(xlog1py(Inf + im * NaN, 1))
    @test isnan(xlog1py(0 + im * 0, NaN))
    @test iszero(xlog1py(0 + im * 0, -1 + im * Inf))
end

@testset "xexpx" begin
    for x in (false, 0, 0.0, 0f0, -Inf, -Inf32)
        @test (@inferred xexpx(x)) === zero(exp(x))
    end
    for x in (NaN16, NaN32, NaN64, Inf16, Inf32, Inf64)
        @test (@inferred xexpx(x)) === x
    end
    for x in (1, true, 1.0, 1f0)
        @test (@inferred xexpx(x)) === exp(x)
    end
    for a in (2, 2f0, 2.0), x in -a:a
        @test (@inferred xexpx(x)) === x * exp(x)
    end
end

@testset "xexpy" begin
    for x in (0, 1, 1.0, 1f0, Inf, Inf32), y in (-Inf, -Inf32)
        @test (@inferred xexpy(x, y)) === zero(x * exp(y))
    end
    for x in (0, 1, 1.0, 1f0, Inf, Inf32, -Inf, -Inf32, NaN, NaN32), nan in (NaN, NaN32)
        @test (@inferred xexpy(x, nan)) === oftype(x * exp(nan), NaN)
        @test (@inferred xexpy(nan, x)) === oftype(nan * exp(x), NaN)
    end
    for x in (2, -2f0, 2.0), y in (1, -1f0, 1.0)
        @test (@inferred xexpy(x, y)) ≈ x * exp(y)
    end
    for x in (randn(), randn(Float32))
        @test xexpy(x, x) ≈ xexpx(x)
    end
end

@testset "logistic & logit" begin
    @test logistic(2) ≈ 1.0 / (1.0 + exp(-2.0))
    @test logistic(-750.0) === 0.0
    @test logistic(-740.0) > 0.0
    @test logistic(+36.0) < 1.0
    @test logistic(+750.0) === 1.0
    @test iszero(logit(0.5))
    @test logit(logistic(2)) ≈ 2.0
end

@testset "logcosh" begin
    for x in (randn(), randn(Float32))
        @test @inferred(logcosh(x)) isa typeof(x)
        @test logcosh(x) ≈ log(cosh(x))
        @test logcosh(-x) == logcosh(x)
    end

    # special values
    for x in (-Inf, Inf, -Inf32, Inf32)
        @test @inferred(logcosh(x)) === oftype(x, Inf)
    end
    for x in (NaN, NaN32)
        @test @inferred(logcosh(x)) === x
    end
end

@testset "log1psq" begin
    @test iszero(log1psq(0.0))
    @test log1psq(1.0) ≈ log1p(1.0)
    @test log1psq(2.0) ≈ log1p(4.0)
end

# log1pexp, log1mexp, log2mexp & logexpm1

@testset "log1pexp" begin
    for T in (Float16, Float32, Float64, BigFloat), x in 1:40
        @test (@inferred log1pexp(+log(T(x)))) ≈ T(log1p(big(x)))
        @test (@inferred log1pexp(-log(T(x)))) ≈ T(log1p(1/big(x)))
    end

    # special values
    @test (@inferred log1pexp(0)) ≈ log(2)
    @test (@inferred log1pexp(0f0)) ≈ log(2)
    @test (@inferred log1pexp(big(0))) ≈ log(2)
    @test (@inferred log1pexp(+1)) ≈ log1p(ℯ)
    @test (@inferred log1pexp(-1)) ≈ log1p(ℯ) - 1

    # large arguments
    @test (@inferred log1pexp(1e4)) ≈ 1e4
    @test (@inferred log1pexp(1f4)) ≈ 1f4
    @test iszero(@inferred log1pexp(-1e4))
    @test iszero(@inferred log1pexp(-1f4))

    # compare to accurate but slower implementation
    correct_log1pexp(x::Real) = x > 0 ? x + log1p(exp(-x)) : log1p(exp(x))
    # large range needed to cover all branches, for all floats (from Float16 to BigFloat)
    for T in (Int, Float16, Float32, Float64, BigInt, BigFloat), x in -300:300
        @test (@inferred log1pexp(T(x))) ≈ float(T)(correct_log1pexp(big(x)))
    end
    # test BigFloat with multiple precisions
    for prec in (10, 20, 50, 100), x in -300:300
        setprecision(prec) do
            y = big(float(x))
            @test @inferred(log1pexp(y)) ≈ correct_log1pexp(y)
        end
    end
end

@testset "log1mexp" begin
    @test log1mexp(-1.0)  ≈ log1p(- exp(-1.0))
    @test log1mexp(-10.0) ≈ log1p(- exp(-10.0))
end

@testset "log2mexp" begin
    @test log2mexp(0.0)  ≈ 0.0
    @test log2mexp(-1.0) ≈ log(2.0 - exp(-1.0))
end

@testset "logexpm1" begin
    @test logexpm1(2.0)            ≈  log(exp(2.0) - 1.0)
    @test logexpm1(log1pexp(2.0))  ≈  2.0
    @test logexpm1(log1pexp(-2.0)) ≈ -2.0

    @test logexpm1(2f0)            ≈  log(exp(2f0) - 1f0)
    @test logexpm1(log1pexp(2f0))  ≈  2f0
    @test logexpm1(log1pexp(-2f0)) ≈ -2f0
end

@testset "log1pmx" begin
    @test iszero(log1pmx(0.0))
    @test log1pmx(1.0) ≈ log(2.0) - 1.0
    @test log1pmx(2.0) ≈ log(3.0) - 2.0
end

@testset "logmxp1" begin
    @test iszero(logmxp1(1.0))
    @test logmxp1(2.0) ≈ log(2.0) - 1.0
    @test logmxp1(3.0) ≈ log(3.0) - 2.0
end

@testset "logsumexp" begin
    @test logaddexp(2.0, 3.0)     ≈ log(exp(2.0) + exp(3.0))
    @test logaddexp(10002, 10003) ≈ 10000 + logaddexp(2.0, 3.0)

    for x in ([1.0], Complex{Float64}[1.0])
        @test @inferred(logsumexp(x)) == 1.0
        @test @inferred(logsumexp((xi for xi in x))) == 1.0
    end

    for x in ([1.0, 2.0, 3.0], Complex{Float64}[1.0, 2.0, 3.0])
        @test @inferred(logsumexp(x)) ≈ 3.40760596444438
        @test logsumexp(x .+ 1000) ≈ 1003.40760596444438
    end

    for x in ((1.0, 2.0, 3.0), map(complex, (1.0, 2.0, 3.0)))
        @test @inferred(logsumexp(x)) ≈ 3.40760596444438
    end

    _x = [[1.0, 2.0, 3.0] [1.0, 2.0, 3.0] .+ 1000.]
    for x in (_x, complex(_x))
        @test @inferred(logsumexp(x; dims=1)) ≈ [3.40760596444438 1003.40760596444438]
        @test @inferred(logsumexp(x; dims=[1, 2])) ≈ [1003.4076059644444]
        y = copy(x')
        @test @inferred(logsumexp(y; dims=2)) ≈ [3.40760596444438, 1003.40760596444438]
    end

    # check underflow
    @test logsumexp([1e-20, log(1e-20)]) ≈ 2e-20
    @test logsumexp(Complex{Float64}[1e-20, log(1e-20)]) ≈ 2e-20

    let cases = [([-Inf, -Inf], -Inf),   # correct handling of all -Inf
                 ([-Inf, -Inf32], -Inf), # promotion
                 ([-Inf32, -Inf32], -Inf32), # Float32
                 ([-Inf, Inf], Inf),
                 ([-Inf, 9.0], 9.0),
                 ([Inf, 9.0], Inf),
                 ([0, 0], log(2.0))] # non-float arguments
        for (arguments, result) in cases
            @test logaddexp(arguments...) ≡ result
            @test logsumexp(arguments) ≡ result
            @test logsumexp(complex(arguments)) ≡ complex(result)
        end
    end

    @test isnan(logsubexp(Inf, Inf))
    @test logsubexp(-Inf, -Inf) ≡ -Inf
    @test logsubexp(Inf, -Inf) ≡ Inf
    @test logsubexp(-Inf, Inf) ≡ Inf
    @test logsubexp(Inf, 9.0) ≡ Inf
    @test logsubexp(9.0, Inf) ≡ Inf
    @test logsubexp(-Inf, 9.0) ≡ 9.0
    @test logsubexp(9.0, -Inf) ≡ 9.0
    @test logsubexp(1f2, 1f2) ≡ -Inf32
    @test logsubexp(0, 0) ≡ -Inf
    @test logsubexp(3, 2) ≈ 2.541324854612918108978

    # NaN propagation
    for f in (logaddexp, logsubexp)
        @test isnan(f(NaN, 9.0))
        @test isnan(f(NaN, Inf))
        @test isnan(f(NaN, -Inf))
        @test isnan(f(9.0, NaN))
        @test isnan(f(Inf, NaN))
        @test isnan(f(-Inf, NaN))
    end

    @test isnan(logsumexp([NaN, 9.0]))
    @test isnan(logsumexp([NaN, Inf]))
    @test isnan(logsumexp([NaN, -Inf]))
    @test isnan(logsumexp(Complex{Float64}[NaN, 9.0]))
    @test isnan(logsumexp(Complex{Float64}[NaN, Inf]))
    @test isnan(logsumexp(Complex{Float64}[NaN, -Inf]))
    @test isnan(logsumexp(Complex{Float64}[NaN * im, 9.0]))
    @test isnan(logsumexp(Complex{Float64}[NaN * im, Inf]))
    @test isnan(logsumexp(Complex{Float64}[NaN * im, -Inf]))

    # logsumexp with general iterables (issue #63)
    xs = range(-500, stop = 10, length = 1000)
    @test @inferred(logsumexp(x for x in xs)) == logsumexp(xs)
    xs = range(-500 + 0.5im, stop = 10 + 30im, length = 1000)
    @test @inferred(logsumexp(x for x in xs)) == logsumexp(xs)

    # complex numbers
    xs = randn(Complex{Float64}, 10, 5)
    @test @inferred(logsumexp(xs)) ≈ log(sum(exp.(xs)))
    @test @inferred(logsumexp(xs; dims=1)) ≈ log.(sum(exp.(xs); dims=1))
    @test @inferred(logsumexp(xs; dims=2)) ≈ log.(sum(exp.(xs); dims=2))
    @test @inferred(logsumexp(xs; dims=[1, 2])) ≈ log(sum(exp.(xs); dims=[1, 2]))
    @test @inferred(logsumexp(x for x in xs)) == logsumexp(xs)
end

@testset "softmax" begin
    x = [1.0, 2.0, 3.0]
    r = exp.(x) ./ sum(exp.(x))

    # in-place versions
    for T in (Float32, Float64)
        s = Vector{T}(undef, 3)
        softmax!(s, x)
        @test s ≈ r

        fill!(s, zero(T))
        softmax!(s, x; dims=1)
        @test s ≈ r

        s = Matrix{T}(undef, 1, 3)
        softmax!(s, x)
        @test s ≈ permutedims(r)

        @test_throws DimensionMismatch softmax!(s, x; dims=1)

        fill!(s, zero(T))
        softmax!(s, permutedims(x); dims=2)
        @test s ≈ permutedims(r)

        fill!(s, zero(T))
        softmax!(s, permutedims(x); dims=1:2)
        @test s ≈ permutedims(r)
    end
    softmax!(x)
    @test x ≈ r

    for (S, T) in ((Int, Float64), (Float64, Float64), (Float32, Float32))
        x = S[1, 2, 3]
        s = softmax(x)
        @test s ≈ r
        @test eltype(s) === T

        x = repeat(S[1, 2, 3], 1, 3)
        s = softmax(x; dims=1)
        @test s ≈ repeat(r, 1, 3)
        @test eltype(s) === T

        x = repeat(S[1 2 3], 3, 1)
        s = softmax(x; dims=2)
        @test s ≈ repeat(permutedims(r), 3, 1)
        @test eltype(s) === T

        x = S[1 2 3]
        s = softmax(x; dims=1:2)
        @test s ≈ permutedims(r)
        @test eltype(s) === T
    end

    x = [1//2, 2//3, 3//4]
    r = exp.(x) ./ sum(exp.(x))
    s = softmax(x)
    @test s ≈ r
    @test eltype(s) === Float64

    # non-standard indices: #12
    x = OffsetArray(1:3, -2:0)
    s = softmax(x)
    @test s isa OffsetArray{Float64}
    @test axes(s, 1) == OffsetArrays.IdOffsetRange(-2:0)
    @test collect(s) ≈ softmax(1:3)
end


text:
"""
Struct for LogNormalMixture

		lnmModel=LogNormalMixture(η::Array{num1},λ::Array{num2}) where {num1,num2<: Number}

Where:\n
		η  =	Array of volatilities.
		λ  = 	Array of weights.
"""
mutable struct LogNormalMixture{num <: Number, num2 <: Number, abstrUnderlying <: AbstractUnderlying, numtype <: Number} <: ItoProcess{numtype}
    η::Array{num, 1}
    λ::Array{num2, 1}
    underlying::abstrUnderlying
    function LogNormalMixture(η::Array{num, 1}, λ::Array{num2, 1}, underlying::abstrUnderlying) where {num <: Number, num2 <: Number, abstrUnderlying <: AbstractUnderlying}
        if minimum(η) <= 0.0
            error("Volatilities must be positive")
        elseif minimum(λ) <= 0.0
            error("weights must be positive")
        elseif sum(λ) > 1.0
            error("λs must be weights")
        elseif length(λ) != length(η) - 1
            error("Check vector lengths")
        else
            zero_typed = zero(num) + zero(num2)
            return new{num, num2, abstrUnderlying, typeof(zero_typed)}(η, λ, underlying)
        end
    end
end

export LogNormalMixture;

function simulate!(S, mcProcess::LogNormalMixture, rfCurve::AbstractZeroRateCurve, mcBaseData::AbstractMonteCarloConfiguration, T::Number)
    @assert T > 0.0
    r = rfCurve.r
    S0 = mcProcess.underlying.S0
    d = dividend(mcProcess)
    η_gbm = mcProcess.η
    λ_gmb = mcProcess.λ
    push!(λ_gmb, 1.0 - sum(mcProcess.λ))
    mu_gbm = r - d
    simulate!(S, GeometricBrownianMotion(η_gbm[1], mu_gbm, S0), mcBaseData, T)
    S .*= λ_gmb[1]
    tmp_ = similar(S)
    for (η_gbm_, λ_gmb_) in zip(η_gbm[2:end], λ_gmb[2:end])
        simulate!(tmp_, GeometricBrownianMotion(η_gbm_, mu_gbm, S0), mcBaseData, T)
        S .+= λ_gmb_ * tmp_
    end
    return nothing
end


text:
using JuMP
using Clp
using LinearAlgebra
using Plots

n = 100
grid = range(-1, 1, length=n)
f = ℯ.^grid
coeff = zeros(n, 5)
for degree in 0:4
    coeff[:, degree + 1] = grid .^ degree
end
errors = zeros(5)
polynomials = zeros(5, 5)
taylors = [1 0 0 0 0; 1 1 0 0 0; 1 1 1/2 0 0; 1 1 1/2 1/6 0; 1 1 1/2 1/6 1/24]

for d in 0:4

    model = Model(Clp.Optimizer)

    @variables(model, begin
        m >= 0
        a[1:5]
    end
    )

    @objective(model, Min, m)

    @constraints(model, begin
        [i=d+2:5], a[i] == 0
        m .>= f - coeff * a
        m .>= coeff * a - f
    end
    )

    optimize!(model)
    errors[d + 1] = objective_value(model)
    polynomials[d + 1, :] = value.(a)
end

plot(grid, f, label="e^x")
for d in 0:4
    plot!(grid, coeff * (polynomials[d + 1, :]), label="cheb $d")
end
savefig("polynomial.png")

plot(0:4, errors)
savefig("error.png")

for d in 0:4
    plot(grid, f, label="e^x")
    plot!(grid, coeff * polynomials[d + 1, :], label="cheb $d")
    plot!(grid, coeff * taylors[d + 1, :], label="tayl $d")
    savefig("compare-$d.png")
end

max_stars_repo_name: tansongchen/learn-optimization

max_issues_repo_name: tansongchen/learn-optimization

max_forks_repo_name: tansongchen/learn-optimization


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the position of the Moon.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/sun/moon_position.jl
# ================================

# Function moon_position_i
# ------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Example 5-3: Finding the Moon position vector [1, p. 289]
#
#   In April 28, 1994, 00:00 UTC, one gets:
#
#             | -134_240.626 |
#   r_J2000 = | -311_571.590 | km
#             | -126_693.785 |
#
#   Note: it assumes that UTC ≈ UT1 ≈ TDB.
#
################################################################################

@testset "Function moon_position_i" begin
    JD_UTC = DatetoJD(1994,4,28,0,0,0)
    JD_TBD = JD_UTC

    r_J2000 = moon_position_i(JD_TBD)

    @test r_J2000[1]/1000 ≈ -134_240.626 atol = 1e-3
    @test r_J2000[2]/1000 ≈ -311_571.590 atol = 1e-3
    @test r_J2000[3]/1000 ≈ -126_693.785 atol = 1e-3
end


text:
using Oceananigans
using JLD2
using Oceananigans.Architectures: arch_array
using Oceananigans: AbstractModel
using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize

function get_parameter(observation, parameter_path)
    file = jldopen(observation.path)
    parameter = file[parameter_path]
    close(file)
    return parameter
end

vectorize(observation) = [observation]
vectorize(observations::Vector) = observations

"""
    OneDimensionalEnsembleModel(observations::SyntheticObservationsBatch; architecture = CPU(), ensemble_size = 50, kwargs...)

Build an Oceananigans `HydrostaticFreeSurfaceModel` with many independent
columns. The model grid is given by the data in `observations`, and the
dynamics in each column is attached to its own `CATKEVerticalDiffusivity` closure stored in
an `(Nx, Ny)` Matrix of closures. `ensemble_size = 50` is the
desired count of ensemble members for calibration with Ensemble Kalman Inversion (EKI), and the
remaining keyword arguments `kwargs` define the default closure across all columns.

In the "many columns" configuration, we run the model on a 3D grid with `(Flat, Flat, Bounded)` boundary
conditions so that many independent columns can be evolved at once with much of the computational overhead
split among the columns. The `Nx` rows of vertical columns are each reserved for an "ensemble" member
whose attached parameter value (updated at each iteration of EKI) sets the diffusivity closure
used to predict the model solution for the `Ny` physical scenarios described by the simulation-specific
`SyntheticObservations` objects in `observations`.
"""
function OneDimensionalEnsembleModel(observations;
    architecture = CPU(),
    ensemble_size = 50,
    closure = CATKEVerticalDiffusivity(Float64; warning = false)
)
    observations = vectorize(observations)

    Lz = get_parameter(observations[1], "grid/Lz")
    g = get_parameter(observations[1], "buoyancy/gravitational_acceleration")
    α = get_parameter(observations[1], "buoyancy/equation_of_state/α")

    Nx = ensemble_size
    Ny = length(observations)
    Nz = observations[1].grid.Nz
    column_ensemble_size = ColumnEnsembleSize(Nz = Nz, ensemble = (Nx, Ny), Hz = 1)

    ensemble_grid = RectilinearGrid(size = column_ensemble_size, z = (-Lz, 0), topology = (Flat, Flat, Bounded))
    closure_ensemble = arch_array(architecture, [closure for i = 1:Nx, j = 1:Ny])

    get_Qᵇ(observation) = get_parameter(observation, "parameters/boundary_condition_θ_top") * α * g
    get_Qᵘ(observation) = get_parameter(observation, "parameters/boundary_condition_u_top")
    get_f₀(observation) = get_parameter(observation, "coriolis/f")
    get_dudz_bottom(observation) = get_parameter(observation, "parameters/boundary_condition_u_bottom")
    get_dbdz_bottom(observation) = get_parameter(observation, "parameters/boundary_condition_θ_bottom") * α * g

    to_arch_array(f) = arch_array(architecture, [f(observations[j]) for i = 1:Nx, j = 1:Ny])

    Qᵇ_ensemble = to_arch_array(get_Qᵇ)
    Qᵘ_ensemble = to_arch_array(get_Qᵘ)
    ensemble_u_bottom = to_arch_array(get_dudz_bottom)
    ensemble_b_bottom = to_arch_array(get_dbdz_bottom)

    ensemble_b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ_ensemble), bottom = GradientBoundaryCondition(ensemble_b_bottom))
    ensemble_u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ_ensemble), bottom = GradientBoundaryCondition(ensemble_u_bottom))
    coriolis_ensemble = arch_array(architecture, [FPlane(f = get_f₀(observations[j])) for i = 1:Nx, j = 1:Ny])

    ensemble_model = HydrostaticFreeSurfaceModel(grid = ensemble_grid,
        tracers = (:b, :e),
        buoyancy = BuoyancyTracer(),
        boundary_conditions = (; u = ensemble_u_bcs, b = ensemble_b_bcs),
        coriolis = coriolis_ensemble,
        closure = closure_ensemble)

    return ensemble_model
end


text:
# This file is a part of SimilaritySearch.jl

module SimilaritySearch
abstract type Index end
abstract type AbstractSearchContext end

using Parameters
import Distances: evaluate, PreMetric
export AbstractSearchContext, PreMetric, evaluate, search, searchbatch

include("utils/knnresult.jl")

const GlobalKnnResult = [KnnResult(30)]   # see __init__ function at the end of this file

@inline getknnresult(res=nothing) = res !== nothing ? res : @inbounds GlobalKnnResult[Threads.threadid()]

"""
    search(searchctx::AbstractSearchContext, q, k::Integer=maxlength(searchctx.res))
    search(searchctx::AbstractSearchContext, q)

This is the most generic search function. It calls almost all implementations whenever an integer k is given.

"""
function search(searchctx::AbstractSearchContext, q, k::Integer=10)
    res = getknnresult()
    empty!(res, k)
    search(searchctx, q, res)
end

function searchbatch(searchctx::AbstractSearchContext, Q, k::Integer=10; parallel=false)
    searchbatch(searchctx, Q, [KnnResult(k) for i in 1:length(Q)]; parallel)
end

function searchbatch(searchctx::AbstractSearchContext, Q, KNN; parallel=false)
    if parallel
        Threads.@threads for i in eachindex(Q, KNN)
            @inbounds search(searchctx, Q[i], KNN[i])
        end
    else
        @inbounds for i in eachindex(Q, KNN)
            search(searchctx, Q[i], KNN[i])
        end
    end

    KNN
end

@inline Base.length(searchctx::AbstractSearchContext) = length(searchctx.db)
@inline Base.getindex(searchctx::AbstractSearchContext, i) = searchctx.db[i]
@inline Base.eachindex(searchctx::AbstractSearchContext) = 1:length(searchctx)
@inline Base.eltype(searchctx::AbstractSearchContext) = eltype(searchctx.db)

include("distances/bits.jl")
include("distances/sets.jl")
include("distances/strings.jl")
include("distances/vectors.jl")
include("distances/cos.jl")

include("utils/perf.jl")

include("indexes/pivotselection.jl")
include("indexes/seq.jl")
include("indexes/pivottable.jl")
include("indexes/pivotselectiontables.jl")
include("indexes/kvp.jl")

include("graph/graph.jl")

function __init__()
    __init__visitedvertices()
    __init__beamsearch()
    __init__neighborhood()
    for i in 2:Threads.nthreads()
        push!(GlobalKnnResult, KnnResult(32))
    end
end

end


text:
"""
    find_global_maximum_complex(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}})

`find_global_maximum_complex` is similar to `find_global_maximum` instead is desigend
for handling graphs with weights.

# Arguments
- `graph::Dict{Int64,Array{Tuple{Int64,Int64},1}}`: Graph of the connected nodes with weights
"""
function find_global_maximum_complex(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}})
    global_maximum = 0
    # Going through the dictionary
    for (key, value) in graph

        # Check the keys as new global maximum
        if key > global_maximum
            global_maximum = key
        end

        # Check the values as new local maximum
        for local_maximum in value

            # Check if new local maximum becomes global maximum
            if local_maximum[1] > global_maximum
                global_maximum = local_maximum[1]
            end

        end

    end
    return global_maximum
end

"""
    find_parent_in_spanning_tree(parent::Array{Int64,1}, i::Int64, ref_key::Int64)

The find parent algorithm is a recursive function to find the subeset of an item i for
graph test of having a spanning tree.

# Arguments
- `parent::Array{Int64,1}`: Array of the subset of the items
- `i::Int64`: Index of the parents-item
"""
function find_parent_in_spanning_tree(parent::Array{Int64,1}, i::Int64)
    if parent[i] == i
        return i
    end
    return find_parent_in_spanning_tree(parent, parent[i])
end

"""
Provided by: [https://discourse.julialang.org/t/sort-matrix-based-on-the-elements-of-a-specific-column/23475/5](https://discourse.julialang.org/t/sort-matrix-based-on-the-elements-of-a-specific-column/23475/5)
"""
function sortrows(M, by = zeros(0))
    if by == zeros(0)
        order = copy(M)
    else
        order = copy(by)
    end
    if size(order, 2) > 1
        order = Int64.(order .- minimum(order, dims = 1))
        order = (order ./ maximum(order, dims = 1)) * (10) .^ (size(order, 2):-1:1)
    end
    order = sortperm(order[:, 1])
    return M[order, :]
end

"""
    minimum_spanning_tree(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}}))

The minimum spanning tree (MST) algorithm detects a subset of the edges of a connected,
edge-weighted undirected graph that connects all the vertices together. The MST algorithms
focus is a) to exclude any cycles and b) to find the minimum possible total edge weight,
which will create a spanning tree whose sum of edge weights is as small as possible.
The Kruskal's algorithm is used to find the minimum spanning forest of an undirected
edge-weighted graph.
For more information see: [https://en.wikipedia.org/wiki/Minimum_spanning_tree](https://en.wikipedia.org/wiki/Minimum_spanning_tree)
and [https://en.wikipedia.org/wiki/Kruskal%27s_algorithm](https://en.wikipedia.org/wiki/Kruskal%27s_algorithm)

# Arguments
- `graph::Dict{Int64,Array{Tuple{Int64,Int64},1}}`: Graph of the connected nodes with the weights

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> graph_with_spanning_tree = Dict(1 => [(2, 10), (3, 6), (4, 5)], 2 => [(4, 15)], 3 => [(4, 4)]
julia> ClassicAlgorithmsCollections.minimum_spanning_tree(graph_with_spanning_tree)
[3 4 4; 1 4 5; 1 2 10]
```
"""
function minimum_spanning_tree(graph::Dict{Int64,Array{Tuple{Int64,Int64},1}})
    # Get the total number of nodes
    vertex_size = find_global_maximum_complex(graph)

    # Intialize array list
    graph_array = zeros(Int64, false)

    # Append the dictionary to 1D array
    for (key, value) in graph
        for items in value
            append!(graph_array, [key, items[1], items[2]])
        end
    end

    # Retransform from 1D-array to 2D-array with three columns for node, connected node,
    # and weight.
    graph_array = adjoint(reshape(graph_array, (3, :)))
    # Because step on is sorting all the edges in non-decreasing order of their
    # weight, the dictionary has to be translated into an array to be able to be
    # sorted.
    graph_array = sortrows(graph_array, graph_array[:, [3]])

    # Create array-list for finding the parents and intialize all subsets to vertex_size
    # with an increasing range, because it has to be reordered by the ranking later.
    parent = collect(Int64, 1:vertex_size)

    # Create an array for ranking
    rank = zeros(Int64, vertex_size)

    # Create a array list for the results
    result = zeros(Int64, false)

    # Number of edges to be taken is equal to vertex_size
    i = 1 # An index variable, used for sorted edges
    e = 1 # An index variable, used for result[]
    while e < vertex_size
        # Step 2: Pick the smallest edge and increment the index for next iteration
        u, v, w = graph_array[i, :]
        i += 1
        x = find_parent_in_spanning_tree(parent, u)
        y = find_parent_in_spanning_tree(parent, v)

        # If including this edge does't cause cycle, include it in result and
        # increment the index of result for next edge
        if x != y
            e += 1
            append!(result, [u, v, w])
            #  union of two sets of x and y by ranking
            xroot = find_parent_in_spanning_tree(parent, x)
            yroot = find_parent_in_spanning_tree(parent, y)

            # Attach smaller rank tree under root of high rank tree (Union by Rank)
            if rank[xroot] < rank[yroot]
                parent[xroot] = yroot
            elseif rank[xroot] > rank[yroot]
                parent[yroot] = xroot

                # If ranks are same, then make one as root and increment its rank by one
            else
                parent[yroot] = xroot
                rank[xroot] += 1
            end
        end
        # Else discard the edge
    end
    # Reshape the 1D array to 3D array
    return adjoint(reshape(result, (3, :)))
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
@testset "55.jump-game.jl" begin
    @test can_jump([2, 3, 1, 1, 4]) == true
    @test can_jump([3, 2, 1, 0, 4]) == false
    @test can_jump([3, 2, 1, 0, 4]) == false
    @test can_jump([0]) == true
    @test can_jump([5, 9, 3, 2, 1, 0, 2, 3, 3, 1, 0, 0]) == true
end


text:
export align_all, distille_PSF, nn_mutual, select_rois, clear_border, border_dist, remove_border, precise_extract

"""
    align_all(rois; shifts=nothing, optimizer=LBFGS(), iterations=500, verbose=false)

aligns a series of regions of interest `rois` of equal sizes to another. This is achieved by fitting a Gaussian to each of the ROIs and applying a Fourier-based shift
to align all of them to the center. The aligned raw data is then summed.
#
`positions`: If provided, these will be added to the determined shifts for the valid fits.
"""
function align_all(rois; shifts=nothing, positions=nothing, iterations=500, verbose=false)
    res_shifts = []
    shifted = []
    if length(rois) == 0
        return [], [], []
    end
    mysum = zeros(Float32, size(rois[1]))
    for n=1:length(rois)
        myshift, valid = let
            if isnothing(shifts)
                to_center = rois[n] # gaussf(rois[n], 1.0)
                params, _, _ = gauss_fit(to_center, iterations=iterations, verbose=verbose)
                params[:μ], all(params[:FWHM] .< size(rois[n])) && all(params[:FWHM] .> 0.0)
            else
                shifts[n], true
            end
        end
        if !valid
            println("Warning fit failed (sigma .> roi_size). discarding bead $n")
        else
            if isnothing(positions)
                push!(res_shifts, myshift)
            else
                push!(res_shifts, myshift .+ positions[n])
            end
            ashifted = shift(rois[n], .-myshift)
            push!(shifted, ashifted)
            mysum .+= ashifted
        end
    end
    println("Averaged $(length(shifted)) beads.")
    return mysum./length(shifted), shifted, res_shifts
end

"""
    nn_mutual(points)
    returns the mutual nearest neighbor distances and the index of the neares neighbor, both in matrix form.
"""
function nn_mutual(points)
    kdtree = KDTree(points)
    nns, nndist = knn(kdtree, points, 2)
    nns = [mynn[1] for mynn in nns]
    nndist = [myd[1] for myd in nndist]
    return nndist, nns
end

"""
    select_rois(img::AbstractArray, roi_pos::Vector; valid=nothing, roi_size=(16,16))

extracts a series of ROIs based from the array `img` on the positions supplied in `Vector`.
The latter can be a vector of Integer positions or a vector of `CartesionCoord` as optained by
Image.local.

#arguments
+ `img`: image to extract ROIs from
+ `roi_pos`: a list of ND-indices where to extract. Alternatively also a matrix works, where the columns correspond to the positions.
+ `valid`: if provided a binary vector indicates which of the positions are valid and will be extracted
+ `roi_size`: the size of the region of interst to extract at each position via `NDTools.select_region()`

#returns
a tuple of a vector of extracted regions of interest and a list of cartesian indices at which they were extracted.
"""
function select_rois(img::AbstractArray, roi_pos::Vector; valid=nothing, roi_size=(16,16))
    rois = []
    cart_list = Vector{CartesianIndex{ndims(img)}}([])
    for n in 1:length(roi_pos)
        if isnothing(valid) || valid[n]
            pos = Tuple(round.(Int,roi_pos[n]))
            push!(rois, select_region(img, center=pos, new_size=roi_size))
            push!(cart_list, CartesianIndex(pos))
        end
    end
    return rois, cart_list
end

function select_rois(img::AbstractArray, roi_pos::Matrix; valid=nothing, roi_size=(16,16))
    select_rois(img, mat_to_vec(roi_pos); valid=valid, roi_size=roi_size)
end

"""
    clear_border(img, roi_size; pad_value=zero(eltype(img)))

inserts `pad_value` into the border region.

"""
function clear_border(img, roi_size; pad_value=zero(eltype(img)))
    sz = size(img)
    return select_region(select_region(img, new_size=sz .- roi_size, pad_value=pad_value), new_size=sz, pad_value=pad_value)
end

"""
    border_dist(vec, sz)

calculates the closest distance to the border along each dimension, given a vector of coordinate vectors.
"""
function border_dist(vec, sz, ndims=length(sz))
    [min.(Tuple(v[1:ndims]).-1, sz[1:ndims].-Tuple(v[1:ndims])) for v in vec]
end

"""
    remove_border(vec::Vector, roi_size)

removes all positions which are too close to the border in a vector of cordinates.
"""
function remove_border(vec::Vector, sz, roi_size; valid=nothing)
    res = []
    dist = border_dist(vec, sz, length(roi_size))
    for n=1:length(vec)
        if all(dist[n] .> roi_size) && (isnothing(valid) || valid[n])
            push!(res, vec[n])
        end
    end
    return res
end

function remove_border(mat::Matrix, sz, roi_size; valid=nothing)
    remove_border(mat_to_vec(mat), sz, roi_size; valid=valid)
end

"""
    precise_extract(img, positions, roi_size)

extracts a region of size `roi_size` from the ND-dataset `img` at each position centering each region precisely at the supixel cooredinates as specified by `positions`.

#returns
a tuple of `(extracted, cart_ids, mysum)` with `extracted` referring to the extracted and shifted regions and `mymean` being the average of all of these regions.
`cart_ids` refers to cartesian indices of the nearest pixel positions.
"""
function precise_extract(img, positions, roi_size)
    all_max_vec = [round.(Int, positions[n]) for n=1:length(positions)]
    shifts = [positions[n].- all_max_vec[n] for n=1:length(positions)]
    rois, cart_ids = select_rois(img, all_max_vec; roi_size=roi_size);
    mymean, extracted, shifts = align_all(rois, shifts=shifts);
    return extracted, cart_ids, mymean
end

"""
    destille_PSF(img, σ=1.3; positions=nothing, force_align=false, rel_thresh=0.1, min_dist=16.0, roi_size=(16,16), upper_thresh=nothing, pixelsize=1.0)

automatically extracts multiple PSFs from one dataset, alignes and averages them. The input image `img` should contain a sparse set of PSF measurements
obtained by imagin beads, QDots or NV centers.
If you want to apply this to multicolor or multimode datasets, run it first on one channel and then again on the other channels using the `positions` argument.

#arguments
+ `img`: ND-image to destill the PSF from
+ `σ`: the size of the filtering kernel usde in the preprocessing step before finding the local maxima. This may be noise-dependent.
+ `positions`: if a list of (sub-pixel precision) positions is provided, these will be used instead of aligning them.
+ `force_align`: If true, the subpixel-alignment will be done, even though positions are given.
+ `rel_trhesh`: the threshold spedifying which local maxima are valid. This is a relative value compared to the maximum of the Gauss-filtered data.
+ `min_dist`: The minimum distance in pixels to the nearest other maximum.
+ `roi_size`: The size of the region of interest to extract. The default is 2D but this should also work for higher dimensions assuming the size of `img` for the other dimensions.
+ `upper_thresh`: if provided, also an upper relative threshold will be applied eliminating the very bright particles. If you have clusters, try `upper_thresh = 0.45`.
+ `pixelsize`: size of a pixel. Scales the FWHMs and σ in the result parameters.

#returns
a tuple of  `(psf, shifted, shifts, cart_ids, selected)`
+ `psf`: the distilled PSF
+ `shifted`: the individual aligned beads
+ `shifts`: the subpixels shifts
+ `cart_ids`: a vector of cartesian indices where the valid beads were selected
+ `selected`: a Float32 image with selection rings around each bead
+ `params`: the result of the fit of the final PSF.
"""
function distille_PSF(img, σ=1.3; positions=nothing, force_align=false, rel_thresh=0.1, min_dist=nothing, roi_size=(16,16), verbose=true, upper_thresh=nothing, pixelsize=1.0, preferred_z=nothing)
    # may also upcast to Float32
    img, bg = remove_background(img,2 .*σ)
    println("Subtracted a background of $(bg)")
    if isnothing(min_dist)
        min_dist = maximum(roi_size)
    end
    all_max_vec=[]
    if isnothing(positions)
        gimg, roisz = let
            # decide whether to select only one slice for finding the beads
            if isnothing(preferred_z)
                gaussf(img, σ), roi_size
            else
                gaussf(img[:,:,preferred_z], σ), roi_size[1:2] #slice(img,3,preferred_z)
            end
        end
        gimg .*= (gimg .> maximum(gimg) .* rel_thresh)
        # gimg = clear_border(gimg, roi_size)
        all_max = findlocalmaxima(gimg)
        all_val = gimg[all_max]
        all_max_vec = cart_to_mat(all_max)
        println("Found $(length(all_max)) maxima to consider.")

        nndist, _ = nn_mutual(all_max_vec)
        valid_nn = nndist .> min_dist

        println("Found $(sum(valid_nn)) beads to consider with sufficient distance min_dist=$(min_dist).")
        if !isnothing(upper_thresh)
            max_b = maximum(all_val) .* upper_thresh
            valid_nn .*= (all_val .< max_b)
            println("Found $(sum(valid_nn)) beads to consider with correct brightness $(max_b).")
        end
        all_max_vec = remove_border(all_max_vec, size(gimg), roisz./2; valid=valid_nn)
        println("Removed border. $(length(all_max_vec)) beads remaining.")

        all_max_vec = let
            # decide whether to select only one slice for finding the beads
            if isnothing(preferred_z)
                all_max_vec
            else
                [vcat(all_max_vec[n][1:2],preferred_z) for n=1:length(all_max_vec)]
            end
        end
        rois, cart_ids = select_rois(img, all_max_vec; roi_size=roi_size);
        force_align=true
    else
        rois, cart_ids, psf = precise_extract(img, positions, roi_size)
        all_max_vec = positions
    end
    if force_align
        println("Averaging $(length(rois)) regions of interest.")
        psf, rois, positions = align_all(rois, positions=all_max_vec, verbose=verbose);
    end
    selected = make_rings(size(img), cart_ids, expand_size(roi_size,size(img)))

    psf, _ = remove_background(psf,σ)
    if !isempty(psf)
        params, fwd, allp = gauss_fit(psf, verbose=verbose, pixelsize=pixelsize);

        return (psf, rois, positions, selected, params, fwd)
    else
        return [],[],[],positions,[],[]
    end
end

function make_rings(sz, cart_ids, roi_size)
    res = zeros(sz)
    res[cart_ids] .= 1.0
    res = gaussf(res)
    tmp = Float32.((rr(sz, scale=2 ./roi_size) .> 0.9) .* (rr(sz, scale=2 ./roi_size) .< 1.1))
    return conv_psf(res,tmp)
end


text:
module Persa

export Preference, UserPreference, Rating, Dataset
export preference, users, items, machine, predict, train!

# package code goes here
include("preference.jl")
include("rating.jl")
include("dataset.jl")
include("learn.jl")
include("machine.jl")

end # module


text:
@doc raw"""
    GroupOperationAction(group::AbstractGroupManifold, AD::ActionDirection = LeftAction())

Action of a group upon itself via left or right translation.
"""
struct GroupOperationAction{G,AD} <: AbstractGroupAction{AD}
    group::G
end

function GroupOperationAction(
    G::AbstractGroupManifold,
    ::TAD=LeftAction(),
) where {TAD<:ActionDirection}
    return GroupOperationAction{typeof(G),TAD}(G)
end

function Base.show(io::IO, A::GroupOperationAction)
    return print(io, "GroupOperationAction($(A.group), $(direction(A)))")
end

base_group(A::GroupOperationAction) = A.group

g_manifold(A::GroupOperationAction) = A.group

function switch_direction(A::GroupOperationAction)
    return GroupOperationAction(A.group, switch_direction(direction(A)))
end

apply(A::GroupOperationAction, a, p) = translate(A.group, a, p, direction(A))

apply!(A::GroupOperationAction, q, a, p) = translate!(A.group, q, a, p, direction(A))

function inverse_apply(A::GroupOperationAction, a, p)
    return inverse_translate(A.group, a, p, direction(A))
end

function inverse_apply!(A::GroupOperationAction, q, a, p)
    return inverse_translate!(A.group, q, a, p, direction(A))
end

function apply_diff(A::GroupOperationAction, a, p, X)
    return translate_diff(A.group, a, p, X, direction(A))
end

function apply_diff!(A::GroupOperationAction, Y, a, p, X)
    return translate_diff!(A.group, Y, a, p, X, direction(A))
end

function inverse_apply_diff(A::GroupOperationAction, a, p, X)
    return inverse_translate_diff(A.group, a, p, X, direction(A))
end

function inverse_apply_diff!(A::GroupOperationAction, Y, a, p, X)
    return inverse_translate_diff!(A.group, Y, a, p, X, direction(A))
end

function optimal_alignment(A::GroupOperationAction, p, q)
    return inverse_apply(switch_direction(A), p, q)
end

function optimal_alignment!(A::GroupOperationAction, x, p, q)
    return inverse_apply!(switch_direction(A), x, p, q)
end

function center_of_orbit(
    A::GroupOperationAction,
    pts::AbstractVector,
    p,
    mean_method::AbstractEstimationMethod,
)
    μ = mean(A.group, pts, mean_method)
    return inverse_apply(switch_direction(A), p, μ)
end
function center_of_orbit(A::GroupOperationAction, pts::AbstractVector, p)
    μ = mean(A.group, pts)
    return inverse_apply(switch_direction(A), p, μ)
end


text:
Base.@propagate_inbounds Base.getindex(A::DSSASolution, i::Int) = [A.u[i], A.channel[i]]
Base.@propagate_inbounds Base.getindex(A::DSSASolution, i::Int, ::Colon) = [A.u[j][i] for j in 1:length(A.t)]

function (A::DSSASolution)(s::Symbol,i::Int)
    if s == :channel
        @assert i <= length(A.channel[1])
        return [A.channel[j][i] for j in 1:length(A.t)]
    elseif s ==:u
        @assert i <= length(A.u[1])
        return [A.u[j][i] for j in 1:length(A.t)]
    end
end

function (A::DSSASolution)(tval)
    @unpack odesol = A
    odesol(tval)
end

(integrator::DSSAIntegrator)(t) = copy(integrator.u)
(integrator::DSSAIntegrator)(out,t) = (out .= integrator.u)

function Base.show(io::IO, m::MIME"text/plain", A::DSSASolution)
    println(io,string("retcode: ",A.odesol.retcode))
    println(io,string("Interpolation: "),DiffEqBase.interp_summary(A.odesol.interp))
    print(io,"t: ")
    show(io,m,A.t)
    println(io)
    print(io,"u: ")
    show(io,m,A.u)
    println(io)
    print(io,"channel: ")
    show(io,m,A.channel)
    println(io,"\n===\nUse sol.u to check the state variable and sol.channel to check the delay channel solution.\n===")
end

# """
#     function get_reaction_idx(rn::ReactionSystem)

# Get the rearranged order for a reaction system, returns a matrix `m` with first column as the rearranged reactions, second column is the order in the orign `ReactionSystem`
# If you know the old idx for a certain reaction in `ReactionSystem`, you can use
# `m[:,2][idx]` to get the new idx.
# """

# function get_reaction_idx(rn::Catalyst.ReactionSystem)
#     massactionjump_order = Int64[]
#     constantjump_order = Int64[]
#     rxvars = []
#     for (i, rx) in enumerate(ModelingToolkit.get_eqs(rn))
#         empty!(rxvars)
#         (Catalyst.reactionrates(rn) isa ModelingToolkit.Symbolic) && ModelingToolkit.get_variables!(rxvars, Catalyst.reactionrates(rn))
#         @inbounds for j = 1:length(rxvars)
#             if isequal(rxvars[j], ModelingToolkit.get_iv(rn))
#                 error("Does not support VariableRateJump")
#             end
#         end
#         if Catalyst.ismassaction(rx, rn)
#             push!(massactionjump_order, i)
#         else
#             push!(constantjump_order, i)
#         end
#     end
#     order = vcat(massactionjump_order, constantjump_order)
#     # println("The rearranged order :")
#     vcat([[ModelingToolkit.get_eqs(rn)[order[i]] order[i]] for i in eachindex(order)]...)
# end
# export get_reaction_idx


text:
module S2
abstract type S2Region end

using Markdown, LinearAlgebra

import Base:(*), (+), (-), (/), (<), (≤)
import Base:min, max

S2_M_PI = Float64(π)

include("S1Angle.jl")
include("S1ChordAngle.jl")
include("S2Cap.jl")
include("S2LatLng.jl")
include("R1Interval.jl")
include("S1Interval.jl")
include("S2LatLngRect.jl")
include("S2Point.jl")
include("S2Cell.jl")

export S1Angle, radians, degrees, distance, distance2, normalizeRange
export S2Point, norm2, ortho, rotate, full, isfull, isinverted
export S1ChordAngle, toAngle
export R1Interval, S1Interval, center, complement, complementcenter
export interiorcontains, intersects, interiorintersects, union, intersection
export addPoint, clampPoint, fromPointPair, fromPoint, expanded, getDirectedHausdorffDistance
export S2Cap, getCapBound
export S2LatLng, S2LatLngRect, bound
export S2Cell

end