text:
"""
    using QuantumComputer

a quantum computer simulator. built out of boredom and intrigue. if you are new to julia, be warned - array indices start at 1. throughout this module, qubit index 1 is considered the most significant qubit. in typical circuit diagrams, it represents the top line.

i am fully aware that i 'over-typed' this module. i am new to julia and was curious about performance differences but i assume the compiler will sort everything out in many of the cases i am being explicit.

what this module can do:
- create registers, superpositions, gates and circuits
- apply circuits to superpositions
- add gates to circuits
- add measurements to circuits
- add circuits to other circuits as subcircuits
- add `HybridComponents` to a circuit, which are components where there is a classical interaction with the quantum portion of the simulator
- measure arbitrary qubits from a superposition
- generate a variety of single- and multi-qubit static and dynamic gates
- generate some circuits
- construct a gate from a circuit
- cache gates on disk

what this module cannot do:
- decompose a gate into a circuit of smaller gates (planned, long term)
- print a graphical representation of the circuit (planned, relatively simple but also tedious so not planned for near future)
- use multiple registers in one circuit (right now you get one quantum and one classical)
- output bloch sphere animations of algorithms running (planned)

# Example
```julia
using QuantumComputer

qubit_count = 6
initial_value = 12
constant = 13

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)

println(classical_register.value)
# 25
```
"""
module QuantumComputer

using LinearAlgebra
using Serialization
using StatsBase

# base components (registers and superposition)
###############################################

"""
    ClassicalRegister(width, value)

a classical register is used to store the results of a measurement

# Arguments
- `width`: the width, in bits, of the register
- `value`: the intitial value of the register (optional, default = 0)
"""
mutable struct ClassicalRegister
    width::Int64
    value::Int64

    function ClassicalRegister(width::Int64, value::Int64 = 0)
        value < 2^width || throw(DomainError(value, "value does not fit in register"))

        new(width, value)
    end
end

"""
    Register(qubit_count, value)

a quantum register represents a vector of qubits

# Arguments
- `qubit_count`: the width of the register
- `value`: the initial value of the qubits (optional, default = 0)
"""
mutable struct Register
    qubits::Matrix{Complex{Float64}}

    function Register(qubit_count::Int64, value::Int64 = 0)
        value < 2^qubit_count || throw(DomainError(value, "value does not fit in register"))

        qubits::Matrix{Int64} = hcat(ones(Int64, qubit_count), zeros(Int64, qubit_count))

        for i::Int64 = 1:qubit_count
            if (value & 2^(i - 1)) != 0
                k::Int64 = (qubit_count - i + 1)
                qubits[k:k, :] = [0 1]
            end
        end

        new(qubits)
    end
end

"""
    Superposition(qubits)
    Superposition(state)

a superposition is a representation of the probabilities of all possible states of a vector of qubits, typically a register (but there is no hard requirement for a register as input in this simulator)
"""
mutable struct Superposition
    state::Array{Complex{Float64},1}

    function Superposition(qubits::Matrix{Complex{Float64}})
        new(qubit_tensor_product(qubits))
    end

    function Superposition(state::Array{Complex{Float64},1})
        new(state)
    end
end

"""
    qubit_tensor_product(qubits)

this custom function computes a tensor product on the n rows of a matrix of 2-column qubits. this function is threaded, and the core of the algorithm follows in the next function. this works by iteratively applying the tensor product qubit by qubit.

# Arguments
- `qubits`: a matrix of qubits, in the same arrangement as a `Register`
"""
function qubit_tensor_product(qubits::Matrix{Complex{Float64}})::Array{Complex{Float64},1}
    qubit_count::Int64 = size(qubits, 1)
    thread_count::Int64 = Threads.nthreads()

    product::Array{Complex{Float64},1} =
        qubit_count > 1 ? qubits[1, :] :
        throw(DomainError(qubits, "at least 2 qubits required"))
    for i::Int64 = 2:qubit_count
        basis::Array{Complex{Float64},1} = qubits[i, :]
        product_size::Int64 = length(product)
        next_product::Array{Complex{Float64},1} =
            Array{Complex{Float64}}(undef, 2 * product_size)

        # coming back to this code, i was confused why we're iterating to
        # min(thread_count, product_size) but passing the child function thread_count
        # turns out that since we need to do some computations with product size in
        # the child function, we end up handling this there.
        Threads.@threads for thread_number::Int64 = 1:min(thread_count, product_size)
            qubit_tensor_product_thread(
                next_product,
                basis,
                product,
                product_size,
                thread_number,
                thread_count,
            )
        end

        product = next_product
    end

    product
end

"""
    qubit_tensor_product_thread(next_product, basis, product, product_size, thread_number, thread_count)

the threaded portion from our qubit tensor product algorithm. this will compute a slice of the current tensor product expansion from the parent algorithm.
this fills in a portion of `next_product` by computing a partial tensor product of `basis` and `product`.

# Arguments
- `next_product`: the product (state) that is currently being computed
- `basis`: the qubit being applied to the previous product
- `product`: the previous product
- `product_size`: the number of states represented by `product`
- `thread_number`: the algorithm-assigned thread number of this thread
- `thread_count`: the total number of threads
"""
function qubit_tensor_product_thread(
    next_product::Array{Complex{Float64},1},
    basis::Array{Complex{Float64},1},
    product::Array{Complex{Float64},1},
    product_size::Int64,
    thread_number::Int64,
    thread_count::Int64,
)
    width::Int64, padded_count::Int64 = divrem(product_size, thread_count)

    if padded_count < thread_number
        initial_value =
            (width + 1) * padded_count + width * (thread_number - padded_count - 1) + 1
    else
        width += 1
        initial_value = width * (thread_number - 1) + 1
    end

    for i::Int64 = initial_value:min(initial_value + width - 1, product_size)
        next_product[2i-1] = basis[1] * product[i]
        next_product[2i] = basis[2] * product[i]
    end
end

# gates
#######

"""
    Gate(matrix)

a gate operating on a state of size 2^n is really just a 2^n x 2^n unitary matrix. this can be added to a `Circuit`.

# Arguments
- `matrix`: a complex unitary matrix of size 2^n, n a positive integer
"""
struct Gate
    matrix::Matrix{Complex{Float64}}
    superposed_qubits_required::Int64

    function Gate(matrix::Matrix)
        new(matrix, qubits_operated_on_by_unitary_matrix(matrix))
    end
end

"""
    qubits_operated_on_by_unitary_matrix(matrix)

a helper function to compute the number of qubits a gate operates on

# Arguments
- `matrix`: a unitary matrix
"""
function qubits_operated_on_by_unitary_matrix(matrix::Matrix)::Int64
    dimensions = size(matrix)
    dimensions[1] == dimensions[2] || throw(DomainError(matrix, "matrix must be square"))
    log2(dimensions[1])
end

# single qubit gates
####################

# https://qiskit.org/textbook/ch-states/single-qubit-gates.html

# pauli
gate_x = Gate([0 1; 1 0])
gate_y = Gate([0.0 -1.0im; 1.0im 0.0])
gate_z = Gate([1 0; 0 -1])
# hadamard
gate_h = Gate([1 1; 1 -1] / √2)
# phase
gate_p(ϕ::Float64) = Gate([1 0; 0 ℯ^(ϕ*im)])
# identity
gate_i = Gate([1 0; 0 1])
# pi / 2 and pi / 4 phase rotations
gate_s = gate_p(pi / 2)
gate_sdg = gate_p(-pi / 2)
gate_t = gate_p(pi / 4)
gate_tdg = gate_p(-pi / 4)
# u, the most general gate
gate_u(θ::Float64, ϕ::Float64, λ::Float64) =
    Gate([cos(θ / 2) -sin(θ / 2)*ℯ^(λ*im); sin(θ / 2)*ℯ^(ϕ*im) cos(θ / 2)*ℯ^((ϕ+λ)*im)])

# multi-qubit gates
###################

"""
[wikipedia](https://en.wikipedia.org/wiki/Quantum_logic_gate#Square_root_of_swap_gate)
"""
gate_root_swap =
    Gate([1 0 0 0; 0 (1.0+im)/2 (1.0-im)/2 0; 0 (1.0-im)/2 (1.0+im)/2 0; 0 0 0 1])

"""
    gate_swap(qubit_a_index, qubit_b_index, qubit_count)

a gate that swaps two arbitrary qubits in a superposition [explanation](https://quantumcomputing.stackexchange.com/questions/9181/swap-gate-on-2-qubits-in-3-entangled-qubit-system)

# Arguments
- `qubit_a_index`: the index of the first qubit to swap
- `qubit_b_index`: the index of the second qubit to swap
- `qubit_count`: the total number of qubits in the superposition
this gate will operate on
"""
function gate_swap(qubit_a_index::Int64, qubit_b_index::Int64, qubit_count::Int64)
    identity = ((1.0 + 0.0im) * I)(2)

    matrix_size = 2^qubit_count
    matrix::Matrix{Complex{Float64}} = zeros(Complex{Float64}, (matrix_size, matrix_size))
    matrices::Array{Matrix{Complex{Float64}},1} = Array{Matrix{Complex{Float64}}}(undef, 0)

    for i = 1:qubit_count
        push!(matrices, identity)
    end

    for i = 0:1
        for j = 0:1
            matrices[qubit_a_index] = ket_bra(i, j)
            matrices[qubit_b_index] = ket_bra(j, i)
            partial_sum = matrices[1]
            for k = 2:qubit_count
                partial_sum = kron(partial_sum, matrices[k])
            end
            matrix += partial_sum
        end
    end

    Gate(matrix)
end

"""
    ket_bra(i, j)

check the stackexchange post on the previous function to understand what is happening here.

# Arguments
- `i`: first value
- `j`: second value
"""
function ket_bra(i::Int64, j::Int64)::Matrix{Complex{Float64}}
    (i > 1 || i < 0 || j > 1 || j < 0) &&
        throw(DomainError((i, j), "arguments must be in range [0, 1]"))

    matrix::Matrix{Complex{Float64}} = zeros(Complex{Float64}, (2, 2))
    matrix[[(j << 1) + i + 1]] = [1.0 + 0.0im]
    matrix
end

function gate_expansion(gate::Gate, qubit_count::Int64)
    matrix_single = gate.matrix
    matrix = matrix_single
    for i in 2:qubit_count
         matrix = kron(matrix, matrix_single)
    end
    QuantumComputer.Gate(matrix)
end

"""
    gate_extension(gate, qubit_index, qubit_count)

a gate that acts on a subset of qubits by identifying a range of
qubits in a superposition and applying a smaller gate to that range,
while applying identity transformations to the qubits outside the
range. effectively, extends a gate smaller than required by the
superposition size to operate on a subset of qubits in that superposition.

# Arguments
- `gate`: the underlying gate
- `qubit_index`: an integer selecting the qubit where the underlying gate begins
- `qubit_count`: the number of qubits in the superposition this gate will act on
"""
function gate_extension(gate::Gate, qubit_index::Int64, qubit_count::Int64)
    qubit_index > qubit_count &&
        throw(DomainError(qubit_index, "index cannot be greater than count"))
    matrix::Matrix{Complex{Float64}} = ((1.0 + 0.0im) * I)(2^(qubit_index - 1))
    matrix = kron(matrix, gate.matrix)
    identity = ((1.0 + 0.0im) * I)(
        2^(qubit_count - gate.superposed_qubits_required - qubit_index + 1),
    )
    Gate(kron(matrix, identity))
end

"""
    gate_control(gate, control_qubits, qubit_index, qubit_count)

creates an n-controlled single qubit gate. acts on an arbitrary qubit in a superposition and is controlled by n arbitrary qubits.

# Arguments
- `gate`: the underlying single qubit gate
- `control_qubits`: an array of integers selecting control qubits
- `qubit_index`: an integer selecting the qubit to which the single-qubit gate is applied
- `qubit_count`: the number of qubits in the superposition this gate will act on
"""
function gate_control(
    gate::Gate,
    control_qubits::Array{Int64},
    qubit_index::Int64,
    qubit_count::Int64,
)
    qubit_index in control_qubits && throw(
        DomainError((control_qubits, qubit_index), "control and target qubits must differ"),
    )

    outer_size = 2^qubit_count
    matrix::Matrix{Complex{Float64}} =
        Matrix{Complex{Float64}}(undef, (outer_size, outer_size))

    mask = 0
    for control_qubit in control_qubits
        mask += (1 << (qubit_count - control_qubit))
    end
    output_mask = (1 << (qubit_count - qubit_index))

    qubits::Matrix{Complex{Float64}} = Matrix{Complex{Float64}}(undef, (qubit_count, 2))
    for n = 1:outer_size
        for i = 1:qubit_count
            if i == qubit_index && ((n - 1) & mask) == mask
                qubits[i:i, :] =
                    (((n - 1) & output_mask) == 0) ? gate.matrix[:, 1:1]' :
                    gate.matrix[:, 2:2]'
            else
                qubits[i:i, :] =
                    ((n - 1) & 2^(qubit_count - i)) == 0 ? gate_i.matrix[:, 1:1]' :
                    gate_i.matrix[:, 2:2]'
            end
        end
        matrix[:, n:n] = qubit_tensor_product(qubits)
    end

    Gate(matrix)
end

"""
    gate_cx(control_qubit, qubit_index, qubit_count)

controlled pauli x gate.

# Arguments
- `control_qubit`: the index of the qubit controlling the pauli x operation
- `qubit_index`: the index of the qubit upon which to perform the pauli x operation
- `qubit_count`: the total number of qubits the resultant gate acts on
"""
function gate_cx(control_qubit::Int64, qubit_index::Int64, qubit_count::Int64)
    gate_control(gate_x, [control_qubit], qubit_index, qubit_count)
end

"""
    gate_cnx(control_qubits, qubit_index, qubit_count)

n-controlled pauli x gate.

# Arguments
- `control_qubits`: the indices of the qubits controlling the pauli x operation
- `qubit_index`: the index of the qubit upon which to perform the pauli x operation
- `qubit_count`: the total number of qubits the resultant gate acts on
"""
function gate_cnx(control_qubits::Array{Int64}, qubit_index::Int64, qubit_count::Int64)
    gate_control(gate_x, control_qubits, qubit_index, qubit_count)
end

"""
    gate_multi_control(gate, control_count, qubit_count)

this special function creates a controlled-n gate efficiently, by utilizing the top n qubits

# Arguments
- `gate`
"""
function gate_multi_control(gate::Gate, control_count::Int64, qubit_count::Int64)
    gate_qubit_count::Int64 = gate.superposed_qubits_required
    gate_qubit_count + control_count == qubit_count ||
        throw(DomainError(qubit_count, "incorrect qubit count for gate + controls"))

    matrix::Matrix{Complex{Float64}} = ((1.0 + 0.0im) * I)(2^qubit_count)

    range = (2^qubit_count-2^gate_qubit_count+1):(2^qubit_count)
    matrix[range, range] = gate.matrix

    Gate(matrix)
end

"""
    gate_fourier_transform(qubit_count)

the quantum fourier transform. [wikipedia](https://en.wikipedia.org/wiki/Quantum_Fourier_transform)

# Arguments
- `qubit_count`: the number of qubits the gate operates on
"""
function gate_fourier_transform(qubit_count::Int64)
    gate_size = 2^qubit_count
    ω = ℯ^(2im * pi / gate_size)

    matrix::Matrix{Complex{Float64}} =
        Matrix{Complex{Float64}}(undef, (gate_size, gate_size))
    for i = 1:gate_size
        for j = i:gate_size
            value = (ω^((i - 1) * (j - 1)))
            matrix[i, j] = value
            matrix[j, i] = value
        end
    end

    Gate(matrix / √gate_size)
end

"""
    gate_invert(gate)

inverts a unitary matrix quickly, taking advantage of the fact that the conjugate transpose is the inverse of a unitary matrix

# Arguments
- `gate`: the gate to invert
"""
function gate_invert(gate::Gate)
    Gate(Array(gate.matrix'))
end

"""
    Measurement(bits_to_output, qubits_to_measure, sample_size)

creates a measurement component for incorporation in a `Circuit`. `bits_to_output` and `qubits_to_measure` must be equal in length and contain no duplicates. it is fine for there to be an intersection of the to arrays, as long as each array has unique values.

# Arguments
- `bits_to_output`: the bits to write in the classical output register. order must correspond to `qubits_to_measure` order.
- `qubits_to_measure`: the qubits to measure, in the order of measurement.
- `sample_size`: the number of samples to perform when measuring
"""
struct Measurement
    qubits_to_measure::Array{Int64,1}
    bits_to_output::Array{Int64,1}
    sample_size::Int64
    samples::Array{Int64,1}

    function Measurement(
        bits_to_output::Array{Int64,1},
        qubits_to_measure::Array{Int64,1},
        sample_size::Int64 = 1024,
    )
        length(qubits_to_measure) == length(bits_to_output) || throw(
            DomainError(
                length(bits_to_output),
                "tuples must have same number of unique elements",
            ),
        )
        qubits_to_measure = unique(qubits_to_measure)
        length(qubits_to_measure) == length(bits_to_output) || throw(
            DomainError(
                length(bits_to_output),
                "tuples must have same number of unique elements",
            ),
        )
        bits_to_output = unique(bits_to_output)
        length(qubits_to_measure) == length(bits_to_output) || throw(
            DomainError(
                length(bits_to_output),
                "tuples must have same number of unique elements",
            ),
        )

        new(
            qubits_to_measure,
            bits_to_output,
            sample_size,
            Array{Int64,1}(undef, sample_size),
        )
    end
end

# using Printf

"""
    measure_superposition(superposition, classical_register, measurement)

measure some qubits and store the result in a classical register

# Arguments
- `superposition`: the superposition being measured
- `classical_register`: the output register
- `measurement`: details of the measurement
"""
function measure_superposition(
    superposition::Superposition,
    classical_register::ClassicalRegister,
    measurement::Measurement,
)
    measurement_qubit_count::Int64 = length(measurement.qubits_to_measure)
    register_qubit_count::Int64 = log2(length(superposition.state))
    probability_of_ones::Array{Float64} = zeros(Float64, measurement_qubit_count)

    for value = 0:(length(superposition.state)-1)
        for i = 1:measurement_qubit_count
            qubit_to_measure = measurement.qubits_to_measure[i]
            exponent = register_qubit_count - qubit_to_measure
            if 2^exponent & value != 0
                probability_of_ones[i] += abs(superposition.state[value+1])^2
            end
        end
    end

    initial_value = classical_register.value
    for j = 1:measurement.sample_size
        measurement.samples[j] = initial_value
        for i = 1:measurement_qubit_count
            bit_to_output = measurement.bits_to_output[i]
            mask = 2^(classical_register.width - bit_to_output)
            # should this be >= or >? it surely doesn't matter but i'd like to be correct
            if probability_of_ones[i] >= rand()
                measurement.samples[j] |= mask
            else
                measurement.samples[j] &= ~mask
            end
        end
    end

    counts = countmap(measurement.samples)
    most_common_value = first(counts)[1]
    max_count = first(counts)[2]
    for (value, count) in counts
        if count > max_count
            most_common_value = value
            max_count = count
        end
    end

    classical_register.value = most_common_value
    # @printf(
    #     "classical register value is now %s\n",
    #     string(classical_register.value, base = 2, pad = classical_register.width)
    # )

    for i = 1:measurement_qubit_count
        bit_output = measurement.bits_to_output[i]
        output_mask::Int64 = 2^(classical_register.width - bit_output)
        bit_set = (classical_register.value & output_mask != 0)
        for value::Int64 = 0:(length(superposition.state)-1)
            qubit_to_measure = measurement.qubits_to_measure[i]
            input_mask::Int64 = 2^(register_qubit_count - qubit_to_measure)
            if bit_set
                if value & input_mask == 0
                    index = (value | input_mask) + 1
                    if superposition.state[index] == 0
                        # do we need to worry about phase here?
                        superposition.state[index] = superposition.state[value+1]
                    else
                        # the resulting phase ignores the component being zeroed out, is this correct?
                        scale =
                            √(
                                abs(superposition.state[index])^2 +
                                abs(superposition.state[value+1])^2,
                            ) / abs(superposition.state[index])
                        superposition.state[index] *= scale
                    end
                    superposition.state[value+1] = 0
                end
            else
                if value & input_mask != 0
                    index = (value & ~input_mask) + 1
                    if superposition.state[index] == 0
                        # do we need to worry about phase here?
                        superposition.state[index] = superposition.state[value+1]
                    else
                        # the resulting phase ignores the component being zeroed out, is this correct?
                        scale =
                            √(
                                abs(superposition.state[index])^2 +
                                abs(superposition.state[value+1])^2,
                            ) / abs(superposition.state[index])
                        superposition.state[index] *= scale
                    end
                    superposition.state[value+1] = 0
                end
            end
        end
    end
end

"""
    HybridComponent(executor, arguments)

a component that can be added to a `Circuit` that is capable of controlling quantum gates and subcircuits with classical logic.

# Arguments
- `executor`: a function that accepts circuit application arguments (see below for signature) and runs code that applies quantum `Gates` and `Circuits` intelligently
- `arguments`: fixed arguments known at the time of `Circuit` composition
"""
struct HybridComponent
    execute::Any
    arguments::Array

    # executor should be a function with the signature
    # (superposition::QuantumComputer.Superposition, classical_register::QuantumComputer.ClassicalRegister, arguments...)
    function HybridComponent(executor, arguments)
        new(executor, arguments)
    end
end

"""
    Circuit()

a quantum circuit. technically this is just an ordered collection of components.
"""
struct Circuit
    components::Array{Union{Circuit,Gate,HybridComponent,Measurement},1}

    function Circuit()
        new(Array{Union{Circuit,Gate,HybridComponent,Measurement}}(undef, 0))
    end
end

"""
    add_gate_to_circuit!(circuit, gate)

adds a `Gate` to a `Circuit`.

# Arguments
- `circuit`: the `Circuit`
- `gate`: the `Gate`
"""
function add_gate_to_circuit!(circuit::Circuit, gate::Gate)
    push!(circuit.components, gate)
end

"""
    add_measurement_to_circuit!(circuit, measurement)

adds a `Measurement` to a `Circuit`.

# Arguments
- `circuit`: the `Circuit`
- `measurement`: the `Measurement`
"""
function add_measurement_to_circuit!(circuit::Circuit, measurement::Measurement)
    push!(circuit.components, measurement)
end

"""
    add_hybrid_component_to_circuit!(circuit, component)

adds a `HybridComponent` to a `Circuit`.

# Arguments
- `circuit`: the `Circuit`
- `component`: the `HybridComponent`
"""
function add_hybrid_component_to_circuit!(circuit::Circuit, component::HybridComponent)
    push!(circuit.components, component)
end

"""
    add_subcircuit_to_circuit!(circuit, subcircuit)

adds a `Circuit` to a another `Circuit` as a component.

# Arguments
- `circuit`: the parent `Circuit`
- `subcircuit`: the child `Circuit`
"""
function add_subcircuit_to_circuit!(circuit::Circuit, subcircuit::Circuit)
    push!(circuit.components, subcircuit)
end

cache_base = nothing
cached_gates = Dict{String,Gate}()

function gate_set_cache_base(base)
    global cache_base = [base]
end

function gate_load_from_cache(cache_path)
    key = join(cache_path, "/")
    key in keys(cached_gates) && return cached_gates[key]
    cache_base == nothing && return nothing

    # after this pop, we need to make it to the push at the end of the function
    # so don't return early for some reason
    filename = pop!(cache_path)

    path = join(vcat(cache_base, cache_path), "/")
    run(`mkdir -p $path`)
    cd(path)

    gate = nothing
    try
        gate = deserialize(string(filename, ".qg"))
        cached_gates[key] = gate
    catch
    end

    ancestor_path = [".." for _ in vcat(cache_base, cache_path)]
    path = join(ancestor_path, "/")
    cd(path)

    push!(cache_path, filename)

    return gate
end

function gate_save_to_cache(cache_path, gate)
    if cache_base != nothing
        filename = pop!(cache_path)

        path = join(vcat(cache_base, cache_path), "/")
        run(`mkdir -p $path`)
        cd(path)

        serialize(string(filename, ".qg"), gate)

        ancestor_path = [".." for _ in vcat(cache_base, cache_path)]
        path = join(ancestor_path, "/")
        cd(path)

        push!(cache_path, filename)
    end

    key = join(cache_path, "/")
    cached_gates[key] = gate
end

function gate_remove_from_cache(cache_path)
    if cache_base != nothing
        filename = pop!(cache_path)

        path = join(vcat(cache_base, cache_path), "/")
        run(`mkdir -p $path`)
        try
            cd(path)
            rm(string(filename, ".qg"))
        catch
        finally
            # this is pretty fragile but the mkdir above kind of saves us from edge cases
            ancestor_path = [".." for _ in vcat(cache_base, cache_path)]
            path = join(ancestor_path, "/")
            cd(path)
        end

        push!(cache_path, filename)
    end

    key = join(cache_path, "/")
    delete!(cached_gates, key)
end

function circuit_convert_to_gate(circuit::Circuit, cache_path::Array{String,1} = Array{String,1}(undef, 0))
    if length(cache_path) > 0
        gate = gate_load_from_cache(cache_path)
        if typeof(gate) != Nothing
            return gate
        end
    end

    gates = filter(component -> typeof(component) == Gate, circuit.components)
    rank = size(gates[1].matrix, 1)

    # initialize with identity
    matrix = ((1.0 + 0.0im) * I)(rank)

    for component in circuit.components
        if typeof(component) == Circuit
            # if we are recursing too deep our circuit is likely pretty poor
            # so this is reasonable
            matrix = circuit_convert_to_gate(component).matrix * matrix
        elseif typeof(component) == Measurement
            throw(
                DomainError(circuit, "cannot contain measurements when converting to gate"),
            )
        elseif typeof(component) == HybridComponent
            throw(
                DomainError(
                    circuit,
                    "cannot contain hybrid components when converting to gate",
                ),
            )
        else
            matrix = component.matrix * matrix
        end
    end

    gate = Gate(matrix)

    if length(cache_path) > 0
      gate_save_to_cache(cache_path, gate)
    end

    return gate
end

"""
    apply_circuit_to_superposition!(superposition, circuit, classical_register)

applies the gates and measurements in a circuit to the superposition provided. a classical register is required for measurement outputs.

# Arguments
- `superposition`: the `Superposition` to operate on
- `circuit`: the `Circuit` to apply
- `classical_register`: a `ClassicalRegister` capable of holding all circuit measurement output bits
"""
function apply_circuit_to_superposition!(
    superposition::Superposition,
    circuit::Circuit,
    classical_register::ClassicalRegister = undef,
)
    for component in circuit.components
        if typeof(component) == Circuit
            # if we are recursing too deep our circuit is likely pretty poor
            # so this is reasonable
            apply_circuit_to_superposition!(superposition, component, classical_register)
        elseif typeof(component) == Measurement
            measure_superposition(superposition, classical_register, component)
        elseif typeof(component) == HybridComponent
            component.execute(superposition, classical_register, component.arguments...)
        else
            superposition.state = component.matrix * superposition.state
        end
    end
end

"""
    QuantumComputer.Circuits

a collection of simple quantum circuits.

# Example
```julia
using QuantumComputer

qubit_count = 6
initial_value = 12
constant = 13

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)

println(classical_register.value)
# 25
```
"""
module Circuits

using LinearAlgebra
using ..QuantumComputer

"""
    constant_adder(n, qubit_count)

a quantum circuit that adds `n` to the superposition's value (`mod 2^qubit_count`). based on Draper's circuit.

# Arguments:
- `n`: the constant to add
- `qubit_count`: the number of qubits in the superposition
"""
function constant_adder(n::Int64, qubit_count::Int64)
    0 <= n < 2^qubit_count || throw(DomainError(n, "must be less than 2^qubit_count"))

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()
    qft::QuantumComputer.Gate = QuantumComputer.gate_fourier_transform(qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, qft)
    subcircuit::QuantumComputer.Circuit = constant_adder_core(n, qubit_count)
    QuantumComputer.add_subcircuit_to_circuit!(circuit, subcircuit)
    inverse_qft::QuantumComputer.Gate = QuantumComputer.gate_invert(qft)
    QuantumComputer.add_gate_to_circuit!(circuit, inverse_qft)

    circuit
end

"""
    constant_adder_core(n, qubit_count)

a quantum circuit that adds `n` to the superposition's value (`mod 2^qubit_count`) without applying pre and post fourier transforms

# Arguments:
- `n`: the constant to add
- `qubit_count`: the number of qubits in the superposition
"""
function constant_adder_core(n::Int64, qubit_count::Int64)
    0 <= n < 2^qubit_count || throw(DomainError(n, "must be less than 2^qubit_count"))

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()
    for i::Int64 = qubit_count:-1:1
        for k::Int64 = 1:i
            if (1 << (i - k)) & n != 0
                gate::QuantumComputer.Gate = QuantumComputer.gate_extension(
                    QuantumComputer.gate_p(pi / 2^(k - 1)),
                    i,
                    qubit_count,
                )
                QuantumComputer.add_gate_to_circuit!(circuit, gate)
            end
        end
    end

    circuit
end

"""
    shor2n3_controlled_controlled_modular_adder(n, a, qubit_count)

a quantum circuit that adds `a` to the superposition's value (`mod n`). see the modular adder circuit in [this paper](https://arxiv.org/pdf/quant-ph/0205095.pdf).

# Arguments:
- `n`: the modulus
- `a`: the constant to add
- `no_cache`: if true, caching will not be employed when constructing gates
- `rebuild`: if true, this will destroy the relevant cached gates before regenerating them
"""
function shor2n3_controlled_controlled_modular_adder(n::Int64, a::Int64, no_cache = true, rebuild = false)
    n_qubit_count::Int64 = ceil(log2(n))
    qubit_count::Int64 = 2 * n_qubit_count + 3

    cache_path = no_cache ? Array{String,1}(undef, 0) : ["constant_adder_core", string(n_qubit_count + 1), string(a)]
    rebuild && !no_cache && gate_remove_from_cache(cache_path)
    gate_add_a = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
    if typeof(gate_add_a) == Nothing
        circuit_add_a::QuantumComputer.Circuit = constant_adder_core(a, n_qubit_count + 1)
        gate_add_a = QuantumComputer.circuit_convert_to_gate(circuit_add_a, cache_path)
    end
    gate_subtract_a::QuantumComputer.Gate = QuantumComputer.gate_invert(gate_add_a)

    offset::Int64 = qubit_count - gate_add_a.superposed_qubits_required - 1
    add_a::QuantumComputer.Gate =
        QuantumComputer.gate_extension(gate_add_a, offset, qubit_count - 2)
    subtract_a::QuantumComputer.Gate =
        QuantumComputer.gate_extension(gate_subtract_a, offset, qubit_count - 2)

    cc_add_a::QuantumComputer.Gate =
        QuantumComputer.gate_multi_control(add_a, 2, qubit_count)
    cc_subtract_a::QuantumComputer.Gate =
        QuantumComputer.gate_multi_control(subtract_a, 2, qubit_count)

    cache_path = no_cache ? Array{String,1}(undef, 0) : ["constant_adder_core", string(n_qubit_count + 1), string(n)]
    rebuild && !no_cache && gate_remove_from_cache(cache_path)
    gate_add_n = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
    if typeof(gate_add_n) == Nothing
      circuit_add_n::QuantumComputer.Circuit = constant_adder_core(n, n_qubit_count + 1)
      gate_add_n = QuantumComputer.circuit_convert_to_gate(circuit_add_n, cache_path)
    end
    gate_subtract_n::QuantumComputer.Gate = QuantumComputer.gate_invert(gate_add_n)
    add_n::QuantumComputer.Gate =
        QuantumComputer.gate_multi_control(gate_add_n, 1, n_qubit_count + 2)

    c_add_n::QuantumComputer.Gate = QuantumComputer.gate_extension(
        add_n,
        qubit_count - add_n.superposed_qubits_required + 1,
        qubit_count,
    )
    subtract_n::QuantumComputer.Gate = QuantumComputer.gate_extension(
        gate_subtract_n,
        qubit_count - gate_subtract_n.superposed_qubits_required + 1,
        qubit_count,
    )

    gate_qft::QuantumComputer.Gate =
        QuantumComputer.gate_fourier_transform(n_qubit_count + 1)
    gate_inverse_qft::QuantumComputer.Gate = QuantumComputer.gate_invert(gate_qft)

    offset = qubit_count - gate_qft.superposed_qubits_required + 1
    qft::QuantumComputer.Gate =
        QuantumComputer.gate_extension(gate_qft, offset, qubit_count)
    inverse_qft::QuantumComputer.Gate =
        QuantumComputer.gate_extension(gate_inverse_qft, offset, qubit_count)

    gate_cx::QuantumComputer.Gate = QuantumComputer.gate_cx(offset, offset - 1, qubit_count)
    gate_x::QuantumComputer.Gate =
        QuantumComputer.gate_extension(QuantumComputer.gate_x, offset, qubit_count)

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    QuantumComputer.add_gate_to_circuit!(circuit, cc_add_a)
    QuantumComputer.add_gate_to_circuit!(circuit, subtract_n)
    QuantumComputer.add_gate_to_circuit!(circuit, inverse_qft)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    QuantumComputer.add_gate_to_circuit!(circuit, qft)
    QuantumComputer.add_gate_to_circuit!(circuit, c_add_n)
    QuantumComputer.add_gate_to_circuit!(circuit, cc_subtract_a)
    QuantumComputer.add_gate_to_circuit!(circuit, inverse_qft)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_x)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_x)
    QuantumComputer.add_gate_to_circuit!(circuit, qft)
    QuantumComputer.add_gate_to_circuit!(circuit, cc_add_a)

    circuit
end

"""
    shor2n3_controlled_modular_multiplier(n, a)

see the modular multiplier circuit in [this paper](https://arxiv.org/pdf/quant-ph/0205095.pdf).

# Arguments:
- `n`: the modulus
- `a`: the constant to multiply by
- `no_cache`: if true, caching will not be employed when constructing gates
- `rebuild`: if true, this will destroy the relevant cached gates before regenerating them
"""
function shor2n3_controlled_modular_multiplier(n::Int64, a::Int64, no_cache = true, rebuild = false)
    n_qubit_count::Int64 = ceil(log2(n))
    qubit_count::Int64 = 2 * n_qubit_count + 3

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    gate_qft::QuantumComputer.Gate =
        QuantumComputer.gate_fourier_transform(n_qubit_count + 1)
    gate_inverse_qft::QuantumComputer.Gate = QuantumComputer.gate_invert(gate_qft)

    qft::QuantumComputer.Gate =
        QuantumComputer.gate_extension(gate_qft, n_qubit_count + 3, qubit_count)
    inverse_qft::QuantumComputer.Gate =
        QuantumComputer.gate_extension(gate_inverse_qft, n_qubit_count + 3, qubit_count)

    QuantumComputer.add_gate_to_circuit!(circuit, qft)
    for i = 1:n_qubit_count
        control_qubit = n_qubit_count - i + 2
        if control_qubit != 2
            gate_swap::QuantumComputer.Gate =
                QuantumComputer.gate_swap(control_qubit, 2, qubit_count)
            QuantumComputer.add_gate_to_circuit!(circuit, gate_swap)
        end

        cache_path = no_cache ? Array{String,1}(undef, 0) : ["shor2n3", "controlled_controlled_modular_adder", string(qubit_count), string(n), string((a * 2^(i - 1)) % n)]
        rebuild && !no_cache && gate_remove_from_cache(cache_path)
        modular_adder = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
        if typeof(modular_adder) == Nothing
            circuit_modular_adder = shor2n3_controlled_controlled_modular_adder(n, (a * 2^(i - 1)) % n, no_cache, rebuild)
            modular_adder = QuantumComputer.circuit_convert_to_gate(circuit_modular_adder, cache_path)
        end

        QuantumComputer.add_gate_to_circuit!(circuit, modular_adder)
        if control_qubit != 2
            QuantumComputer.add_gate_to_circuit!(circuit, gate_swap)
        end
    end
    QuantumComputer.add_gate_to_circuit!(circuit, inverse_qft)

    circuit
end

"""
    shor2n3_controlled_ua(n, a)

the controlled-Ua gate from [this paper](https://arxiv.org/pdf/quant-ph/0205095.pdf).

# Arguments
- `n`: the modulus
- `a`: the constant to multiply `x` by
- `no_cache`: if true, caching will not be employed when constructing gates
- `rebuild`: if true, this will destroy the relevant cached gates before regenerating them
"""
function shor2n3_controlled_ua(n::Int64, a::Int64, no_cache = true, rebuild = false)
    n_qubit_count::Int64 = ceil(log2(n))
    qubit_count::Int64 = 2 * n_qubit_count + 3

    a_inverse::Int64 = invmod(a, n)

    cache_path = no_cache ? Array{String,1}(undef, 0) : ["shor2n3", "controlled_modular_multiplier", string(qubit_count), string(n), string(a)]
    rebuild && !no_cache && gate_remove_from_cache(cache_path)
    modular_multiplier = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
    if typeof(modular_multiplier) == Nothing
        circuit_modular_multiplier::QuantumComputer.Circuit =
            QuantumComputer.Circuits.shor2n3_controlled_modular_multiplier(n, a, no_cache, rebuild)
        modular_multiplier = QuantumComputer.circuit_convert_to_gate(circuit_modular_multiplier, cache_path)
    end

    cache_path = no_cache ? Array{String,1}(undef, 0) : ["shor2n3", "controlled_modular_multiplier", string(qubit_count), string(n), string(a_inverse)]
    rebuild && !no_cache && gate_remove_from_cache(cache_path)
    gate_inverse_modular_divider = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
    if typeof(gate_inverse_modular_divider) == Nothing
        circuit_inverse_modular_divider::QuantumComputer.Circuit =
            QuantumComputer.Circuits.shor2n3_controlled_modular_multiplier(n, a_inverse, no_cache, rebuild)
        gate_inverse_modular_divider =
            QuantumComputer.circuit_convert_to_gate(circuit_inverse_modular_divider, cache_path)
    end
    modular_divider::QuantumComputer.Gate =
        QuantumComputer.gate_invert(gate_inverse_modular_divider)

    circuit = QuantumComputer.Circuit()

    QuantumComputer.add_gate_to_circuit!(circuit, modular_multiplier)

    # controlled register swap
    swap_matrix::Matrix{Complex{Float64}} = (1.0 + 0.0im) * I(2^(qubit_count - 1))
    for i = 1:n_qubit_count
        gate_swap::QuantumComputer.Gate =
            QuantumComputer.gate_swap(i, n_qubit_count + i + 2, qubit_count - 1)
        # these should all commute
        swap_matrix *= gate_swap.matrix
    end
    gate_swap_registers::QuantumComputer.Gate = QuantumComputer.Gate(swap_matrix)
    c_swap::QuantumComputer.Gate =
        QuantumComputer.gate_multi_control(gate_swap_registers, 1, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, c_swap)

    QuantumComputer.add_gate_to_circuit!(circuit, modular_divider)

    circuit
end

function shor2n3_execute_hybrid_logic(
    superposition::QuantumComputer.Superposition,
    classical_register::QuantumComputer.ClassicalRegister,
    gates_inverse_qft::Array{QuantumComputer.Gate},
    gate_ua::QuantumComputer.Gate,
    i::Int64,
)
    qubit_count::Int64 = gate_ua.superposed_qubits_required

    # this happens to be the hadamard matrix
    superposition.state = gates_inverse_qft[1].matrix * superposition.state

    # the core gate
    superposition.state = gate_ua.matrix * superposition.state

    # inverse qft
    if i != 1
        for j = 0:(i-2)
            if (classical_register.value & 2^(classical_register.width - j)) != 0
            # if (classical_register.value & 2^j) != 0
                superposition.state = gates_inverse_qft[i-j].matrix * superposition.state
            end
        end
    end
    superposition.state = gates_inverse_qft[1].matrix * superposition.state

    # measure
    measurement = QuantumComputer.Measurement([i], [1], 1)
    # measurement = QuantumComputer.Measurement([classical_register.width - i + 1], [1], 1)
    QuantumComputer.measure_superposition(superposition, classical_register, measurement)

    # conditional pauli x
    if (classical_register.value & 2^(classical_register.width - i + 1)) != 0
    # if (classical_register.value & 2^(i-1)) != 0
        x::QuantumComputer.Gate =
            QuantumComputer.gate_extension(QuantumComputer.gate_x, 1, qubit_count)
        superposition.state = x.matrix * superposition.state
    end
end

"""
    shor2n3_period_finding(n::Int64, a::Int64)

Beauregard's circuit for finding the period of a^x mod n.

# Arguments
- `n`: the modulus
- `a`: the base
- `no_cache`: if true, caching will not be employed when constructing gates
- `rebuild`: if true, this will destroy the relevant cached gates before regenerating them
"""
function shor2n3_period_finding(n::Int64, a::Int64, no_cache = true, rebuild = false)
    n_qubit_count::Int64 = ceil(log2(n))
    qubit_count::Int64 = 2 * n_qubit_count + 3

    gates_ua::Array{Union{QuantumComputer.Gate,Nothing}} =
        Array{Union{QuantumComputer.Gate,Nothing}}(undef, 2 * n_qubit_count)
    gates_inverse_qft::Array{QuantumComputer.Gate} =
        Array{QuantumComputer.Gate}(undef, 2 * n_qubit_count)
    for i = 1:(2*n_qubit_count)
        if i == 1
            gates_inverse_qft[i] =
                QuantumComputer.gate_extension(QuantumComputer.gate_h, 1, qubit_count)
        else
            gates_inverse_qft[i] = QuantumComputer.gate_extension(
                QuantumComputer.gate_p(-pi / 2^(i - 1)),
                1,
                qubit_count,
            )
        end

        local_a = powermod(a, 2^(i - 1), n)
        cache_path = no_cache ? Array{String,1}(undef, 0) : ["shor2n3", "controlled_ua", string(qubit_count), string(n), string(local_a)]
        rebuild && !no_cache && gate_remove_from_cache(cache_path)
        gates_ua[i] = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
        if typeof(gates_ua[i]) == Nothing
            circuit_ua::QuantumComputer.Circuit =
                shor2n3_controlled_ua(n, local_a, no_cache, rebuild)
            gates_ua[i] = QuantumComputer.circuit_convert_to_gate(circuit_ua, cache_path)
        end
    end

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    for i = 1:(2*n_qubit_count)
        component = QuantumComputer.HybridComponent(
            shor2n3_execute_hybrid_logic,
            [gates_inverse_qft, gates_ua[i], i],
        )
        QuantumComputer.add_hybrid_component_to_circuit!(circuit, component)
    end

    circuit
end

"""
    period_finding_for_11x_mod_15()

a 5 qubit circuit that returns a value related to the period of f(x) = 11^x mod 15. this effectively allows implementation of Shor's algorithm for this specific case.
"""
function period_finding_for_11x_mod_15()
    qubit_count = 5

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()
    gate_h::QuantumComputer.Gate = QuantumComputer.Gate(
        kron(((1.0 + 0.0im) * I)(4), QuantumComputer.gate_expansion(QuantumComputer.gate_h, 3).matrix),
    )
    QuantumComputer.add_gate_to_circuit!(circuit, gate_h)
    gate_cx::QuantumComputer.Gate = QuantumComputer.gate_cx(3, 4, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    gate_cx = QuantumComputer.gate_cx(3, 5, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    gate_h = QuantumComputer.gate_extension(QuantumComputer.gate_h, 4, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_h)
    gate_cp =
        QuantumComputer.gate_control(QuantumComputer.gate_p(pi / 2), [4], 5, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cp)
    gate_h = QuantumComputer.gate_extension(QuantumComputer.gate_h, 5, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_h)

    # uncompute ancilla
    gate_cx = QuantumComputer.gate_cx(3, 2, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    gate_cx = QuantumComputer.gate_cx(3, 1, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    # done uncomputing

    gate_cp =
        QuantumComputer.gate_control(QuantumComputer.gate_p(pi / 4), [4], 3, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cp)
    gate_cp =
        QuantumComputer.gate_control(QuantumComputer.gate_p(pi / 2), [5], 3, qubit_count)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_cp)

    circuit
end

function identification_oracle(qubit_count::Int64, target_input::Int64)
    circuit = QuantumComputer.Circuit()

    matrix_x = QuantumComputer.gate_x.matrix
    matrix_i::Matrix{Complex{Float64}} = ((1.0 + 0.0im)I)(2)

    matrix = matrix_i

    for i in 0:(qubit_count - 2)
        operator = (target_input & (2^i) != 0) ? matrix_i : matrix_x
        matrix = kron(operator, matrix)
    end

    gate_outer::QuantumComputer.Gate = QuantumComputer.Gate(matrix)
    gate_generalized_toffoli::QuantumComputer.Gate = QuantumComputer.gate_multi_control(QuantumComputer.gate_x, qubit_count - 1, qubit_count)

    QuantumComputer.add_gate_to_circuit!(circuit, gate_outer)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_generalized_toffoli)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_outer)

    circuit
end

function grover_diffusion_operator(qubit_count::Int64, gate_h::QuantumComputer.Gate = nothing)
    circuit = QuantumComputer.Circuit()

    if gate_h == nothing
        matrix_h = QuantumComputer.gate_h.matrix
        matrix = matrix_h
        for i in 2:qubit_count
             matrix = kron(matrix, matrix_h)
        end
        gate_h::QuantumComputer.Gate = QuantumComputer.Gate(matrix)
    end

    matrix_x = QuantumComputer.gate_x.matrix
    matrix = matrix_x
    for i in 2:qubit_count
            matrix = kron(matrix, matrix_x)
    end
    gate_x::QuantumComputer.Gate = QuantumComputer.Gate(matrix)

    matrix::Matrix{Complex{Float64}} = ((1.0im) * I)(2)
    identity::Matrix{Complex{Float64}} = ((1.0 + 0.0im) * I)(2^(qubit_count - 2))
    matrix = kron(matrix, identity)
    matrix = kron(matrix, QuantumComputer.gate_h.matrix)
    gate_inner::QuantumComputer.Gate = QuantumComputer.Gate(matrix)

    gate_generalized_toffoli = QuantumComputer.gate_multi_control(QuantumComputer.gate_x, qubit_count - 1, qubit_count)

    QuantumComputer.add_gate_to_circuit!(circuit, gate_h)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_x)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_inner)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_generalized_toffoli)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_inner)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_x)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_h)

    circuit
end

# https://arxiv.org/pdf/quant-ph/0301079.pdf
function grover(oracle::QuantumComputer.Gate, no_cache = true, rebuild = false)
    qubit_count::Int64 = oracle.superposed_qubits_required - 1

    circuit = QuantumComputer.Circuit()

    small_h::QuantumComputer.Gate = QuantumComputer.gate_expansion(QuantumComputer.gate_h, qubit_count)
    gate_h::QuantumComputer.Gate = QuantumComputer.Gate(kron(small_h.matrix, QuantumComputer.gate_h.matrix))
    QuantumComputer.add_gate_to_circuit!(circuit, gate_h)

    cache_path = no_cache ? Array{String,1}(undef, 0) : ["grover", "diffusor", string(qubit_count)]
    rebuild && !no_cache && gate_remove_from_cache(cache_path)
    gate_diffusor = no_cache ? nothing : QuantumComputer.gate_load_from_cache(cache_path)
    if typeof(gate_diffusor) == Nothing
        circuit_diffusor = grover_diffusion_operator(qubit_count, small_h)
        gate_diffusor = QuantumComputer.circuit_convert_to_gate(circuit_diffusor, cache_path)
    end
    gate_extended_diffusor = QuantumComputer.gate_extension(gate_diffusor, 1, qubit_count + 1)

    circuit_grover = QuantumComputer.Circuit()
    QuantumComputer.add_gate_to_circuit!(circuit_grover, oracle)
    QuantumComputer.add_gate_to_circuit!(circuit_grover, gate_extended_diffusor)
    gate_grover = QuantumComputer.circuit_convert_to_gate(circuit_grover)

    iterations = qubit_count < 3 ? 1 : convert(Int64, ceil(pi * √(2^qubit_count) / 4))
    for _ in 1:iterations
        QuantumComputer.add_gate_to_circuit!(circuit, gate_grover)
    end

    circuit
end

function cx_balanced_oracle(qubit_count::Int64)
    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    for i in 1:(qubit_count - 1)
        gate_cx = QuantumComputer.gate_cx(i, qubit_count, qubit_count)
        QuantumComputer.add_gate_to_circuit!(circuit, gate_cx)
    end

    circuit
end

function constant_oracle(qubit_count::Int64, result::Bool)
    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    gate = result ? QuantumComputer.gate_extension(QuantumComputer.gate_x, qubit_count, qubit_count) : QuantumComputer.Gate(convert(Matrix{Complex{Float64}}, ((1.0 + 0.0im) * I)(2^qubit_count)))
    QuantumComputer.add_gate_to_circuit!(circuit, gate)

    circuit
end

function deutsch_jozsa(oracle::QuantumComputer.Gate)
    qubit_count::Int64 = oracle.superposed_qubits_required - 1

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    postfix_h = QuantumComputer.gate_expansion(QuantumComputer.gate_h, qubit_count)
    prefix_h = QuantumComputer.Gate(kron(postfix_h.matrix, QuantumComputer.gate_h.matrix))
    postfix_h_extended = QuantumComputer.gate_extension(postfix_h, 1, qubit_count + 1)

    QuantumComputer.add_gate_to_circuit!(circuit, prefix_h)
    QuantumComputer.add_gate_to_circuit!(circuit, oracle)
    QuantumComputer.add_gate_to_circuit!(circuit, postfix_h_extended)

    circuit
end

function bitwise_product_oracle(secret::Int64, qubit_count::Int64)
    circuit = QuantumComputer.Circuit()

    for i in 1:(qubit_count - 1)
        gate = ((secret & 2^(qubit_count - i - 1)) == 0) ? nothing : QuantumComputer.gate_cx(i, qubit_count, qubit_count)
        gate != nothing && QuantumComputer.add_gate_to_circuit!(circuit, gate)
    end

    if length(circuit.components) == 0
        identity = QuantumComputer.Gate(convert(Matrix{Complex{Float64}}, ((1.0 + 0.0im)I)(2^qubit_count)))
        QuantumComputer.add_gate_to_circuit!(circuit, identity)
    end

    circuit
end

function bernstein_vazirani(oracle)
    qubit_count::Int64 = oracle.superposed_qubits_required - 1

    circuit::QuantumComputer.Circuit = QuantumComputer.Circuit()

    gate_h::QuantumComputer.Gate = QuantumComputer.gate_expansion(QuantumComputer.gate_h, qubit_count)
    gate_prefix::QuantumComputer.Gate = QuantumComputer.Gate(kron(gate_h.matrix, QuantumComputer.gate_z.matrix * QuantumComputer.gate_h.matrix))
    gate_postfix::QuantumComputer.Gate = QuantumComputer.gate_extension(gate_h, 1, qubit_count + 1)

    QuantumComputer.add_gate_to_circuit!(circuit, gate_prefix)
    QuantumComputer.add_gate_to_circuit!(circuit, oracle)
    QuantumComputer.add_gate_to_circuit!(circuit, gate_postfix)

    circuit
end

end # module Circuits

end # module QuantumComputer


text:
# Tests for Cost DEA Models
@testset "CostDEAModel" begin

    ## Test Cost DEA Model with Cooper et al. (2007)
    # Test agains book results
    X = [3 2; 1 3; 4 6]
    Y = [3; 5; 6]
    W = [4 2; 4 2; 4 2]

    deacostcooper = deacost(X, Y, W, rts = :CRS)
    @test efficiency(deacostcooper, :Economic)   ≈ [0.375; 1; 0.429] atol = 1e-3
    @test efficiency(deacostcooper, :Technical)  ≈ [0.9  ; 1; 0.6  ] atol = 1e-3
    @test efficiency(deacostcooper, :Allocative) ≈ [0.417; 1; 0.714] atol = 1e-3

    ## Test Cost DEA Model with Zofío and Prieto (2006) data.
    # Test agains results with R
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]
    W = [2 1; 2 1; 2 1; 2 1; 2 1.0]

    # Cost CRS
    deacostcrs = deacost(X, Y, W, rts = :CRS)
    @test efficiency(deacostcrs) ≈ [0.4307692308;
                                    1.000;
                                    1.000;
                                    0.250;
                                    0.2608695652]
    @test efficiency(deacostcrs, :Technical) ≈ [0.6364;
                                1.000;
                                1.000;
                                0.250;
                                0.2609]  atol = 1e-3
    @test efficiency(deacostcrs, :Allocative) ≈ [0.6769230769;
                                1.000;
                                1.000;
                                1.000;
                                1.000]

    # Cost VRS
    deacostvrs = deacost(X, Y, W, rts = :VRS)
    @test efficiency(deacostvrs) ≈ [0.6153846154;
                                    1.000;
                                    1.000;
                                    0.500;
                                    0.3478260870]
    @test efficiency(deacostvrs, :Technical) ≈ [0.750;
                                1.000;
                                1.000;
                                0.500 ;
                                0.375]  atol = 1e-3
    @test efficiency(deacostvrs, :Allocative) ≈ [0.8205128205;
                                1.000;
                                1.000;
                                1.000;
                                0.9275362319]

    # Check defaults
    @test efficiency(deacost(X, Y, W)) == efficiency(deacostvrs)
    @test efficiency(deacostvrs, :Economic) == efficiency(deacostvrs)

    # Print
    show(IOBuffer(), deacostcooper)

    # Test errors
    @test_throws ErrorException deacost([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws ErrorException deacost([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale
    @test_throws ErrorException deacost([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws ErrorException deacost([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3]) # Different number of input prices and inputs

end

ext: jl


text:
append!(empty!(LOAD_PATH), Base.DEFAULT_LOAD_PATH)
using Pkg

################################################################################
# Generate notebooks
################################################################################
exampledir = joinpath(@__DIR__, "..", "examples")
Pkg.activate(exampledir)
Pkg.instantiate()
include(joinpath(exampledir, "generate_notebooks.jl"))


text:
function euclidian(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T,nsz,nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = (x[i]-x[j])^2
      yd = (y[i]-y[j])^2
      dist[i,j] = dist[j,i] = round(sqrt(xd+yd),RoundNearestTiesUp)
    end
  end
  return dist
end

function att_euclidian(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T,nsz,nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = (x[i]-x[j])^2
      yd = (y[i]-y[j])^2
      dist[i,j] = dist[j,i] = ceil(sqrt((xd+yd)/10.0))
    end
  end
  return dist
end

function ceil_euclidian(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T,nsz,nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = (x[i]-x[j])^2
      yd = (y[i]-y[j])^2
      dist[i,j] = dist[j,i] = ceil(sqrt(xd+yd))
    end
  end
  return dist
end

function geo(x::Vector{T},y::Vector{T}) where T<:Real
  PI = 3.141592
  RRR = 6378.388
  nsz = length(x)
  dist = zeros(T,nsz,nsz)
  degs = trunc.(hcat(x,y))
  mins = hcat(x,y).-degs
  coords = PI.*(degs.+(5.0.*(mins./3.0)))./180.0
  lat = coords[:,1]
  lon = coords[:,2]

  for i in 1:nsz, j in 1:nsz
    if i<=j
      q1 = cos(lon[i]-lon[j])
      q2 = cos(lat[i]-lat[j])
      q3 = cos(lat[i]+lat[j])
      dij = RRR.*acos(0.5.*((1.0.+q1).*q2.-(1.0.-q1).*q3)).+1.0
      dist[i,j] = dist[j,i] = round(dij,RoundNearestTiesUp)
    end
  end
  return dist
end


text:
using Test
using SparseArrays
using LinearAlgebra
using Distributed

println("test MUMPS")
@testset "MUMPS" begin
    @testset "DivGrad" begin include("testDivGrad.jl") end
    @testset "Two systems" begin include("testTwoSystem.jl") end
    @testset "Two systems parallel" begin include("testTwoSystemParallel.jl") end
    @testset "Sparse RHS" begin include("testSparseRHS.jl") end
    println("Done!")
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, ImageBinarization, ColorTypes, ColorVectorSpace
makedocs(sitename="Documentation",
            Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"))
deploydocs(repo = "github.com/zygmuntszpak/ImageBinarization.jl.git")


text:
module FedApp

include("ml/model.jl")
include("ml/mnist.jl")

include("common.jl")

include("server.jl")
include("client.jl")

include("ml/train.jl")
include("ml/evaluate.jl")
include("ml/permuter.jl")

end # module


text:
# ---
# title: 735. Asteroid Collision
# id: problem735
# author: Indigo
# date: 2021-06-29
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/asteroid-collision/description/>
# hidden: true
# ---
#
# We are given an array `asteroids` of integers representing asteroids in a row.
#
# For each asteroid, the absolute value represents its size, and the sign
# represents its direction (positive meaning right, negative meaning left). Each
# asteroid moves at the same speed.
#
# Find out the state of the asteroids after all collisions. If two asteroids
# meet, the smaller one will explode. If both are the same size, both will
# explode. Two asteroids moving in the same direction will never meet.
#
#
#
# **Example 1:**
#
#
#
#     Input: asteroids = [5,10,-5]
#     Output: [5,10]
#     Explanation: The 10 and -5 collide resulting in 10.  The 5 and 10 never collide.
#
#
# **Example 2:**
#
#
#
#     Input: asteroids = [8,-8]
#     Output: []
#     Explanation: The 8 and -8 collide exploding each other.
#
#
# **Example 3:**
#
#
#
#     Input: asteroids = [10,2,-5]
#     Output: [10]
#     Explanation: The 2 and -5 collide resulting in -5. The 10 and -5 collide resulting in 10.
#
#
# **Example 4:**
#
#
#
#     Input: asteroids = [-2,-1,1,2]
#     Output: [-2,-1,1,2]
#     Explanation: The -2 and -1 are moving left, while the 1 and 2 are moving right. Asteroids moving the same direction never meet, so no asteroids will meet each other.
#
#
#
#
# **Constraints:**
#
#   * `1 <= asteroids <= 104`
#   * `-1000 <= asteroids[i] <= 1000`
#   * `asteroids[i] != 0`
#
#
## @lc code=start
using LeetCode

function asteroid_collision(asteroids::Vector{Int})
    res = Int[]
    for ast in asteroids
        ast > 0 && (push!(res, ast); continue)
        flg = true
        while !isempty(res) && res[end] > 0
            cmp = res[end] + ast
            if cmp < 0
                pop!(res)
            elseif cmp == 0
                pop!(res)
                flg = false
            else
                flg = false
                break
            end
        end
        flg && push!(res, ast)
    end
    return res
end
## @lc code=end


text:
using Plots, Test

@testset "Series Attributes" begin
    pl = plot([[1, 2, 3], [2, 3, 4]], lw = 5)
    @test hline!(deepcopy(pl), [1.75])[1].series_list[3][:label] ==
          hline!(deepcopy(pl), [1.75], z_order = :front)[1].series_list[3][:label] ==
          "y3"
    @test hline!(deepcopy(pl), [1.75], z_order = :back)[1].series_list[1][:label] == "y3"
    @test hline!(deepcopy(pl), [1.75], z_order = 2)[1].series_list[2][:label] == "y3"
end

@testset "Axis Attributes" begin
    pl = @test_nowarn plot(; tickfont = font(10, "Times"))
    for axis in (:xaxis, :yaxis, :zaxis)
        @test pl[1][axis][:tickfontsize] == 10
        @test pl[1][axis][:tickfontfamily] == "Times"
    end
end

@testset "Permute recipes" begin
    pl = bar(["a", "b", "c"], [1, 2, 3])
    ppl = bar(["a", "b", "c"], [1, 2, 3], permute = (:x, :y))
    @test xticks(ppl) == yticks(pl)
    @test yticks(ppl) == xticks(pl)
    @test filter(isfinite, pl[1][1][:x]) == filter(isfinite, ppl[1][1][:y])
    @test filter(isfinite, pl[1][1][:y]) == filter(isfinite, ppl[1][1][:x])
end


text:
using Clp
using JuMP
using CoordinateTransformations

import Base.copy

# redefine CoordinateTransformations map evaluation funtion (is this allowed? is this bad?)
function (trans::CoordinateTransformations.AffineMap{M, V})(x::Array{Float64,2}) where {M, V}
    broadcast( + , trans.m * x , trans.v )
end

mutable struct Polyhedron
    G::Array{Float64,2}
    h::Array{Float64,2}
    A::Array{Float64,2}
    b::Array{Float64,2}
end

function Polyhedron(G::Array{Float64,2},h::Array{Float64,2})
    return Polyhedron(G,h,Array{Float64,2}(0,size(G,2)),Array{Float64,2}(0,1))
end

# intersection of two polyhedra
function Polyhedron(pa::Polyhedron,pb::Polyhedron)
    return Polyhedron( [pa.G;pb.G] , [pa.h;pb.h] , [pa.A;pb.A] , [pa.b;pb.b] )
end

function Base.copy(p::Polyhedron)
    return Polyhedron(copy(p.G),copy(p.h),copy(p.A),copy(p.b))
end

function isinside(p::Polyhedron , x::Array{Float64,2} ; eq_tol=1e-4)
    dx =  broadcast( - , p.h , p.G * x )

    if( ~isempty(p.A) )
        eq_viol = broadcast( - , p.A * x , b )
        return all( dx[:] .> 0 ) && all(eq_viol .< eq_tol)
    else
        return all( dx[:] .> 0 )
    end
end

function checksamples(p::Polyhedron , x::Array{Float64,2} ; eq_tol=1e-4)
    dx =  broadcast( - , p.h , p.G * x )
    if( ~isempty(p.A) )
        eq_viol = broadcast( - , p.A * x , p.b )
        print("Max Eq. Viol:   $(maximum(eq_viol[:]))")
        print("Min Ineq. Dist: $(minimum(dx[:]))")
        return all( dx[:] .> 0 ) && all(eq_viol .< eq_tol)
    else
        print("Min Ineq. Dist: $(minimum(dx[:]))")
        return all( dx[:] .> 0 )
    end

end

function bounding_box_old(lbs::Vector{Float64},ubs::Vector{Float64})
    n = length(lbs)
    idx_lb_not_inf = find( x -> ( !isinf(x) )  , lbs )
    idx_ub_not_inf = find( x -> ( !isinf(x) )  , ubs )
    G::Array{Float64,2} = Array{Float64}(0,length(lbs))
    h::Array{Float64,2} = Array{Float64}(0,1)
    for zi=1:length(idx_lb_not_inf)
        zii = idx_lb_not_inf[zi]
        gi = zeros(1,n); gi[zii] = -1;
        G = [G ;  gi ]
        h = [h ; -lbs[zii] ]
    end
    for zi=1:length(idx_ub_not_inf)
        zii = idx_ub_not_inf[zi]
        gi = zeros(1,n); gi[zii] = 1;
        G = [G ;  gi ]
        h = [h ; ubs[zii] ]
    end
    return Polyhedron( G,h )
end

function bounding_box(lbs::Vector{Float64},ubs::Vector{Float64})
    n = length(lbs)
    idx_lb_not_inf = find( x -> ( !isinf(x) )  , lbs )
    idx_ub_not_inf = find( x -> ( !isinf(x) )  , ubs )
    G_s = spzeros(0,length(lbs))
    h::Array{Float64,2} = Array{Float64}(0,1)
    for zi=1:length(idx_lb_not_inf)
        zii = idx_lb_not_inf[zi]
        gi = spzeros(1,n); gi[zii] = -1;
        G_s = [G_s ;  gi ]
        h = [h ; -lbs[zii] ]
    end
    for zi=1:length(idx_ub_not_inf)
        zii = idx_ub_not_inf[zi]
        gi = spzeros(1,n); gi[zii] = 1;
        G_s = [G_s ;  gi ]
        h = [h ; ubs[zii] ]
    end
    return Polyhedron( full(G_s),h )
end

#function find_minimal_affine_subspace(p::Polyhedron ; min_diam=1e-6 )
function find_minimal_affine_subspace( p::Polyhedron )

    conf_min_diam = 1e-6
    output_level = 2
    nc = size(p.G,1)
    d  = size(p.G,2)

    # sort out "zero constraints.."
    idx_zero_constr = []
    for zi in 1:nc
        if(norm(p.G[zi,:]) < 1e-6 )
            push!(idx_zero_constr,zi)
        end
    end
    if(!isempty(idx_zero_constr))
        print(" WARNING! $(length(idx_zero_constr)) infinite inequality constraint found.. ")
        p.G = p.G[ setdiff( (1:nc),(idx_zero_constr)) , : ]
        p.h = p.h[ setdiff( (1:nc),(idx_zero_constr)) , : ]
        nc = size(p.G,1)
    end

    ## Part 1: Check all inequality constraints
    fmins = Vector{Float64}(nc)
    fmaxs = Vector{Float64}(nc)

    x_boundary_min = Array{Float64,2}(d,nc)
    x_boundary_max = Array{Float64,2}(d,nc)

    #print( size(p.G) , size(p.h) , size(p.A) , size(p.b) )
    if(output_level>=1)
        print("Probe $(d) inequality constraints: ")
    end
    for zi in 1:nc
        if(output_level==1); if( mod( ceil(nc/8) )==0 ); print("."); end; end
        if(output_level==2); print("\n$(@sprintf("%3d",zi)) : "); end

        #(lp_a_f,lp_b_f,lp_a_x,lp_b_x) = solve_lp_minmax(p,p.G[zi,:])
        #print("  " ,norm(p.G[zi,:]), "  ")
        (lp_a_f,lp_b_f,lp_a_x,lp_b_x) = solve_lp_minmax(p, p.G[zi,:] / norm(p.G[zi,:]) )
        fmins[zi] = lp_a_f; fmaxs[zi,1] = lp_b_f
        x_boundary_min[:,zi] = lp_a_x
        x_boundary_max[:,zi] = lp_b_x

        if(output_level==2); print("diam= $(@sprintf("%10.6f",lp_b_f-lp_a_f))"); end
    end
    if(output_level>=1); print("\n"); end

    x_boundary = [x_boundary_min  x_boundary_max]
    f_diams = (fmaxs-fmins)
    c_ineq = f_diams .>= conf_min_diam
    c_ineq_to_eq = .~ c_ineq

    if(output_level>2)
        #println("\nDistances: ")
        #println(f_diams)
    end

    if(output_level>0)
        println( "\nMin considered Diameter: $(minimum(f_diams[c_ineq] ))\n" )
        if(any(c_ineq_to_eq))
            println(   "Max non-considered Diameter: $(maximum(f_diams[c_ineq_to_eq] ))\n" )
        else
            println(   "All dimensions considered, polyhedron is full-dimensional")
        end
    end

    ## Part 2: Create the new polytope
    G = p.G[c_ineq,:]; h = p.h[c_ineq,1]
    A = [ p.A ; p.G[c_ineq_to_eq,:] ]
    b = [ p.b ; p.h[c_ineq_to_eq,:] ]

    ## PART 3: Compute nullspace and reparametrize:
    x0_ip     = sum(x_boundary,2) / (2*nc) # just take the mean of all computed points..
    null_A    = nullspace(A)

    # reparametrize:
    G_1 = G * null_A
    h_1 = h - ( G * x0_ip )
    # That should be it.. (?)
    p_new = Polyhedron(G_1,h_1)

#    return p_new

    at::AffineMap = AffineMap(null_A,x0_ip)
    return p_new , at
    # x0 = zeros(size(G_1,2),1) # origin is in reparametrized polytope.
end

function round_polytope(p::Polyhedron)
    if(!isempty(p.A)); error("Rounding works only for polytopes in ineq. representation!"); end

    (f,E) = mve_run(p.G,p.h[:])

    G_r = p.G * inv(E)
    h_r = p.h - p.G * f
    p_r = Polyhedron(G_r,h_r)

    at::AffineMap = AffineMap(inv(E),f)
    return p_r , at
end

# p = Polyhedron( ones(1,4) , ones(4,1) , [eye(4);-eye(4)] , 2*ones(4,1) )

function solve_lp_minmax(p::Polyhedron,v::Vector{Float64})
    m = Model(solver=ClpSolver())
    @variable(m, x[1:length(v)])
    #@objective(m, Max, AffExpr( x , v' , 0 ) )
    @objective(m, Min, v'*x )
    #@constraint( m , sparse(p.A) * x .<= p.b )
    sG = sparse(p.G)
    sA = sparse(p.A)
    @constraint( m , sG * x .<= p.h )
    @constraint( m , sA * x .== p.b )
    # Solve with any available LP solver
    solve(m)
    lp_a_x = getvalue(x)
    lp_a_f = getobjectivevalue(m)

    @objective( m , Max , v'*x )
    solve(m)
    lp_b_x = getvalue(x)
    lp_b_f = getobjectivevalue(m)

    return (lp_a_f,lp_b_f,lp_a_x,lp_b_x)
end

function solve_lp_max_v(p::Polyhedron,v::Vector{Float64})
    m = Model(solver=ClpSolver())
    @variable(m, x[1:length(v)])
    #@objective(m, Max, AffExpr( x , v' , 0 ) )
    @objective(m, Max, v'*x )
    @constraint( m , p.A * x .<= p.b )
    @constraint( m , p.G * x .<= p.h )
    @constraint( m , p.A * x .== p.b )
    # Solve with any available LP solver
    solve(m)
    lp_a_x = getvalue(x)
    lp_a_f = getobjectivevalue(m)
    return (lp_a_f,lp_a_x)
end

function remove_redundant_constraints(p::Polyhedron)
    red = map( xi -> test_redundancy(p,xi) , 1:size(p.G,1) )
    p2 = Polyhedron( copy( p.G[ .~ red , : ] ) , copy(p.h[ .~ red , : ])  , copy(p.A) , copy(p.b) )
    return p2
end

function test_redundancy(p::Polyhedron,i::Int64)
    p2 = copy(p)
    p_g = p2.G[i,:]
    p_h = p2.h[i,1]
    p2.G = p2.G[ (1:size(p.G,1)).!=i,:]
    p2.h = p2.h[ (1:size(p.G,1)).!=i,:]

    (x_f , x_x) = solve_lp_max_v(p2,p_g)

    if( isnan(x_f) )
        # whatever happend here.. play it safe and say we need it.
        # currently, the solve_lp_max_v function returns NaN for unbounded, so this is important..
        return false
    end

    # SAFETY CUSHION!!
    # to be on the safe side wrt. removing redundant constraints..
    CONF_redundance_test_cushion = 0.001;

    if( x_f > p_h-CONF_redundance_test_cushion)
        return false
    else
        return true
    end
 end

"""
heuristic way of identifying the constraints which make a polyhedron empty.
works by relaxing all constraints, until polyhedron is non-empty
"""
function heuristically_find_infeasible_constraints( p::Polyhedron ; relaxation_step = 0.1 , max_tries=500 )

    G_A_p = [p.A;-p.A]
    h_b_p = [p.b;-p.b]

    P_original = Polyhedron( [p.G;G_A_p] , [p.h;h_b_p] )
    Pi         = copy(P_original)

    test_polytope = (p::Polyhedron) -> !isnan( solve_lp_max_v(p,ones(size(p.G,2)))[1] )

    cnt = 0
    print("Relax:\n")
    polytope_ok = test_polytope(Pi)
    while(cnt<max_tries && ! polytope_ok)
        cnt += 1
        print("."); if(mod(cnt,60)==0);print("\n");end;
        Pi.h[:] = Pi.h + relaxation_step
        polytope_ok = test_polytope(Pi)
    end
    print("\nSucces! Polytope feasible.. Now shrink back:\n")
    # now check for every constraint, if the shrunk back version is feasible / infeasible..
    idx_problems = Vector{Int64}()

    for zi=1:length(Pi.h)
        P2 = copy(Pi)
        P2.h[zi] = P_original.h[zi]
        if(!test_polytope(P2))
            push!(idx_problems,zi)
            print("x"); if(mod(zi,60)==0);print("\n");end;
        else
            print("."); if(mod(zi,60)==0);print("\n");end;
        end
    end

    # and remap indeces:
    problems_eq   = filter( x->x >  size(p.G,1) , idx_problems) .- size(p.G,1)
    problems_ineq = filter( x->x <= size(p.G,1) , idx_problems)

    return (problems_ineq,problems_eq)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libsvm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libsvm")
JLLWrappers.@generate_main_file("libsvm", UUID("08558c22-525a-5d2a-acf6-0ac6658ffce4"))
end  # module libsvm_jll


text:
# Postcompose with scaling and sum

export Postcompose

"""
**Postcomposition with an affine transformation**

    Postcompose(f, a=1, b=0)

Returns the function
```math
g(x) = a\\cdot f(x) + b.
```
"""
struct Postcompose{T <: ProximableFunction, R <: Real, S <: Real} <: ProximableFunction
    f::T
    a::R
    b::S
    function Postcompose{T,R,S}(f::T, a::R, b::S) where {T <: ProximableFunction, R <: Real, S <: Real}
        if a <= 0
            error("parameter `a` must be positive")
        else
            new(f, a, b)
        end
    end
end

is_prox_accurate(f::Postcompose) = is_prox_accurate(f.f)
is_separable(f::Postcompose) = is_separable(f.f)
is_convex(f::Postcompose) = is_convex(f.f)
is_set(f::Postcompose) = is_set(f.f)
is_singleton(f::Postcompose) = is_singleton(f.f)
is_cone(f::Postcompose) = is_cone(f.f)
is_affine(f::Postcompose) = is_affine(f.f)
is_smooth(f::Postcompose) = is_smooth(f.f)
is_quadratic(f::Postcompose) = is_quadratic(f.f)
is_generalized_quadratic(f::Postcompose) = is_generalized_quadratic(f.f)
is_strongly_convex(f::Postcompose) = is_strongly_convex(f.f)

Postcompose(f::T, a::R=1, b::S=0) where {T <: ProximableFunction, R <: Real, S <: Real} = Postcompose{T, R, S}(f, a, b)

Postcompose(f::Postcompose{T, R, S}, a::R=1, b::S=0) where {T <: ProximableFunction, R <: Real, S <: Real} = Postcompose{T, R, S}(f.f, a * f.a, b + a * f.b)

function (g::Postcompose)(x::AbstractArray{T}) where {R <: Real, T <: RealOrComplex{R}}
    return g.a * g.f(x) + g.b
end

function gradient!(y::AbstractArray{T}, g::Postcompose, x::AbstractArray{T}) where {
    R <: Real, T <: RealOrComplex{R}
}
    v = gradient!(y, g.f, x)
    y .*= g.a
    return g.a * v + g.b
end

function prox!(y::AbstractArray{T}, g::Postcompose, x::AbstractArray{T}, gamma=R(1)) where {
    R <: Real, T <: RealOrComplex{R}
}
    v = prox!(y, g.f, x, g.a * gamma)
    return g.a * v + g.b
end

function prox_naive(g::Postcompose, x::AbstractArray{T}, gamma=R(1)) where {
    R <: Real, T <: RealOrComplex{R}
}
    y, v = prox_naive(g.f, x, g.a * gamma)
    return y, g.a * v + g.b
end

fun_name(f::Postcompose) = string("Postcomposition of ", fun_name(f.f))
fun_dom(f::Postcompose) = fun_dom(f.f)
fun_expr(f::Postcompose) = "x ↦ a*f(x)+b"
fun_params(f::Postcompose) = string("f(x) = ", fun_expr(f.f), ", a = $(f.a), b = $(f.b)")


text:
using EodData

const USERNAME = "string"
const PASSWORD = "string"

#=
    EodData Web Service Calls
=#

# login()
# -------
# Call login. This will assign you a token which is needed to
# make EodData web service calls.
# The second call directly assigns a token, which is merely
# an ::ASCIIString.
resp = login(USERNAME, PASSWORD)
println(resp.message)
println(resp.token)
println(typeof(resp))

token = login(USERNAME, PASSWORD).token
println(token)
println(typeof(token))

# country_list()
# --------------
# Call and assign countries
countries = country_list(resp.token)
println(countries)

# data_client_latest_version()
# ----------------------------
# Call and assign the latest version for EodData's data client
version = data_client_latest_version(resp.token)
println(version)

# data_formats()
# --------------
# Call and assign the formats available, then assign the
# Standard CSV format to work with at your leisure
formats = data_formats(resp.token)
println(formats)

csv = formats["CSV"]
println(csv.name)
println(csv.format_header)
for column=values(csv.columns)
    println("$(column.column_name) | $(column.column_header)")
end

# exchange_get()
# --------------
# Call and assign a single exchange, in this case the NASDAQ
# and assign the NASDAQ to a variable to work with at your leisure
nasdaq = exchange_get(resp.token, "NASDAQ")
println(nasdaq.name)
println(nasdaq.advances)
println(nasdaq.declines)
println("Advance/Decline Ratio \= $(nasdaq.advances / nasdaq.declines)")

# exchange_list()
# ---------------
# Call and assign the exchanges available (these can be iterated over
# if you wish), then assign the New York Stock Exchange to work with
# at your leisure
exchanges = exchange_list(resp.token)
println(exchanges)

nyse = exchanges["NYSE"]
println(nyse.name)
println(nyse.advances)
println(nyse.declines)
println("Advance/Decline Ratio \= $(nyse.advances / nyse.declines)")

intraday_exchanges = filter((k,v) -> v.is_intraday == true, exchanges)

# exchange_months()
# -----------------
# Call and assign the number of months history available to download
# for a given exchange.
months = exchange_months(resp.token,"NYSE")
println(months)

# fundamental_list()
# ------------------
# Call and assign the fundamentals of all the listings on a given exchange.
# Here we look at the New Zealand Exchange (NZX) market capitalisation which
# also has bonds listed on it, and from time-to-time options as well.
nzx_fundamentals = fundamental_list(resp.token,"NZX")

nzx_market_cap = 0.0
for listing = values(nzx_fundamentals)
    nzx_market_cap += listing.market_cap
end
println(nzx_market_cap)
@sprintf "%.2f" nzx_market_cap

# login_2()
# ---------
# Call and assign the response.
# Not really necessary for end-users from what I can tell
resp = login_2(USERNAME,PASSWORD,"0.1")

# membership()
# ------------
# Call and assign the users membership level/account type with EodData.com
membership = membership(resp.token)
println(membership)

# quote_get()
# -----------
# Call and assign the end-of-day quote for a given instrument
# Here we get the quote for JP Morgan
jpm = quote_get(resp.token, "NYSE", "JPM")
println(jpm)
println(jpm.close)
println(jpm.previous)
println(jpm.change)
println(jpm.simple_return)

# quote_list()
# ------------
# Call and assign the end-of-day quotes for an exchange.
# The collection can be iterated over if you wish
nyse_quotes = quote_list(resp.token, "NYSE")

# quote_list_2()
# --------------
# Call and assign the end-of-day quotes for a custom group
# on a particular exchange.
# The collection can be iterated over if you wish
my_quotes = quote_list_2(resp.token, "NYSE", "C,MS,JPM,BAC,DB")
println(my_quotes)

# quote_list_by_date()
# --------------------
# Call and assign the end-of-day quotes for a custom date
# on a particular exchange.
# The collection can be iterated over if you wish
nyse_20140605 = quote_list_by_date(resp.token, "NYSE", "20140605")
for qt = values(nyse_20140605)
    println("$(qt.name) | $(qt.close)")
end

# quote_list_by_date_2()
# ----------------------
# Call and assign end-of-day quotes, with a smaller type ::Quote_2
# for a custom date on a particular exchange.
# The collection can be iterated over if you wish
asx_20131203 = quote_list_by_date_2(resp.token, "ASX", "20131203")
for qt = values(asx_20131203)
    println("$(qt.ticker_code) | $(qt.close)")
end

# quote_list_by_date_period()
# ---------------------------
# Call and assign quotes for a custom date, and a custom period
# on a particular exchange. If you choose "h" this will return
# hourly data for the exchange
# The collection can be iterated over if you wish
cme_20141008_h = quote_list_by_date_period(resp.token, "CME", "20141008", "h")

for k = keys(cme_20141008_h)
    println(k)
end

for qt = values(cme_20141008_h)
    println("$(qt.ticker_code)\t|\tDate Time: $(qt.date_time)\t|\tClose: $(qt.close)\t|\tVolume: $(qt.volume)")
end

cme_h = collect(cme_20141008_h)

# quote_list_by_date_period_2()
# -----------------------------
# Call and assign quotes, with a smaller type ::Quote_2, for a custom date, and a custom period
# on a particular exchange. If you choose "h" this will return
# hourly data for the exchange
# The keys, values, and collection can be iterated over if you wish
cme_20141008_h_2 = quote_list_by_date_period_2(resp.token, "CME", "20141008", "h")

for k = keys(cme_20141008_h_2)
    println(k)
end

for qt = values(cme_20141008_h_2)
    println("$(qt.ticker_code)\t|\tDate Time: $(qt.date_time)\t|\tClose: $(qt.close)\t|\tVolume: $(qt.volume)")
end

cme_h_2 = collect(cme_20141008_h_2)

# split_list_by_exchange()
# ------------------------
# Call and assign the most recent splits for a given exchange
nyse_splits = split_list_by_exchange(resp.token, "NYSE")
splits = collect(nyse_splits)
println(splits)

# split_list_by_symbol()
# ----------------------
# Call and assign the most recent splits for a given symbol on a particular exchange
nct_splits = split_list_by_symbol(resp.token, "NYSE", "NCT")
for sp = values(nct_splits)
    println("$(sp.exchange_code)\t|\t$(sp.ticker_code)\t|\tDate Time: $(sp.date_time)\t|\tRatio: $(sp.ratio)\t|\tPrice Multiplier: $(sp.price_multiplier)\t|\tReverse Split: $(sp.is_reverse_split)")
end

# symbol_changes_by_exchange()
# ----------------------------
# Call and assign the most recent changes to stock symbols,
# and changes to exchanges
amex_changes = symbol_changes_by_exchange(resp.token, "AMEX")
for sc = values(amex_changes)
    println(sc)
end

# symbol_chart()
# --------------
# Call and assign the url for a chart of the symbol's
# price history.
# Note: There is no "safety" on this web service call, if you give it
# incorrect exchanges, symbols, or incorrect combinations you will receive
# no error.
url = symbol_chart(resp.token, "NYSE", "A")
println(url)

# symbol_get()
# ------------
# Call and assign the detail for a ticker
fb = symbol_get(resp.token, "NASDAQ", "FB")
println(fb)

# symbol_history()
# ----------------
# Call and assign quotes for a ticker from a start date until "today"
# Due to the web service not returning 100% data, the following fields of the quote type will
# be 0, or NaN:
# * open_interest
# * previous
# * change
# * simple_return
# * bid
# * ask
# * previous_close
# * next_open
# * modified
c_20140601_today = symbol_history(resp.token, "NYSE", "C", "20140601")
println(c_20140601_today)

# symbol_history_period()
# -----------------------
# Call and assign quotes for a ticker, for a date, and a custom period
# Due to the web service not returning 100% data, the following fields of the quote type will
# be 0, or NaN:
# * open_interest
# * previous
# * change
# * simple_return
# * bid
# * ask
# * previous_close
# * next_open
# * modified
pg_2014102_h = symbol_history_period(resp.token, "NYSE", "PG", "20141002", "h")
println(pg_2014102_h)

# symbol_history_period_by_date_range()
# -------------------------------------
# Call and assign quotes for a ticker, between a start date and end date, and a custom period
# Due to the web service not returning 100% data, the following fields of the quote type will
# be "", 0, or NaN:
# * description
# * name
# * open_interest
# * previous
# * change
# * simple_return
# * bid
# * ask
# * previous_close
# * next_open
# * modified
amzn_20141020_20141024_30 =
    symbol_history_period_by_date_range(resp.token, "NASDAQ", "AMZN", "20141020", "20141024", "30")
println(amzn_20141020_20141024_30)

# symbol_list()
# -------------
# Call and assign the tickers for a given exchange
nyse_tickers = symbol_list(resp.token, "NYSE")
println(nyse_tickers)

# symbol_list_2()
# ---------------
# Call and assign the tickers for a given exchange
# This is a "smaller" version of the ticker object with only the
# ticker code and ticker name
nyse_tickers_2 = symbol_list_2(resp.token, "NYSE")
println(nyse_tickers_2)

# technical_list()
# ----------------
# Call and assign the technical indicator values for each ticker on a given exchange
nyse_technicals = technical_list(resp.token, "NYSE")
println(nyse_technicals)

# top_10_gains()
# --------------
# Call and assign the quotes for the top 10 gains for the NZX,
# also collect the ticker codes into an array, as the ticker codes
# are the dictionary keys.
nzx_top_10_gains_dict = top_10_gains(resp.token, "NZX")
nzx_top_10_gains_tickers = collect(keys(nzx_top_10_gains_dict))
println(nzx_top_10_gains_tickers)

# top_10_losses()
# ---------------
# Call and assign the quotes for the top 10 losses for the NZX,
# also collect the ticker codes into an array, as the ticker codes
# are the dictionary keys.
nzx_top_10_losses_dict = top_10_losses(resp.token, "NZX")
nzx_top_10_losses_tickers = collect(keys(nzx_top_10_losses_dict))
println(nzx_top_10_losses_tickers)

# validate_access()
# -----------------
# The first call passes plain strings in.
# The second call uses a utility function `set_date_string()` to
# build a ::String argument based on today's date.
validate_access(resp.token, "NZX", "20141001", "h")
validate_access(resp.token, "NYSE", set_date_string(today()), "h")

#=
    EodData external utility functions
=#
# set_date_string()
# -----------------
# The first call passes in a ::DateTime
# The secondcall passes in a ::Date
# Calls 3-6 show how further dynamic strings can be set
set_date_string(now())
set_date_string(today())
quarter_1 = set_date_string(Date(year(today()), 1, 1))
quarter_2 = set_date_string(Date(year(today()), 4, 1))
quarter_3 = set_date_string(Date(year(today()), 7, 1))
quarter_4 = set_date_string(Date(year(today()), 9, 1))


text:
@testset "Test components" begin
    data = IS.SystemData()

    name = "component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(data, component)

    component = IS.get_component(IS.TestComponent, data, name)
    @test component.name == name

    components = IS.get_components(IS.TestComponent, data)
    @test length(components) == 1

    i = 0
    for component in IS.iterate_components(data)
        i += 1
    end
    @test i == 1

    IS.remove_component!(data, collect(components)[1])
    components = IS.get_components(IS.TestComponent, data)
    @test length(components) == 0

    IS.add_component!(data, component)
    components = IS.get_components_by_name(IS.InfrastructureSystemsType, data, name)
    @test length(components) == 1
    @test components[1].name == name

    IS.remove_components!(IS.TestComponent, data)
    components = IS.get_components(IS.TestComponent, data)
    @test length(components) == 0

    @test_throws ArgumentError IS.remove_components!(IS.TestComponent, data)
end

@testset "Test forecasts" begin
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 8954681e-273b-11ec-344d-9326e736c69f
begin
	using Pkg
	Pkg.develop( url="https://github.com/grahamstark/ScottishTaxBenefitModel.jl")
	using ScottishTaxBenefitModel
	using Plots
	using PlutoUI
	Pkg.add( "DataFrames" )
	using DataFrames
end

# ╔═╡ dd1d1115-41b3-4f3f-8cf9-6d97d944b9fa
begin
	using .ExampleHouseholdGetter
	using .STBParameters
	using .BCCalcs
	using .ModelHousehold
	using .Utils
	using .Definitions
	using .SingleHouseholdCalculations
	using .RunSettings
end

# ╔═╡ 6b4b41a4-2475-4d5d-8ea8-aef6d409c1e6
begin
	ExampleHouseholdGetter.initialise()
	hh = get_household("example_hh1")
	# MODEL_PARAMS_DIR

	# head.age
	spouse = get_spouse(hh)
	head = get_head(hh)
	empty!( head.income )
	if spouse !== nothing
		empty!( spouse.income )
	end
	hh
end

# ╔═╡ a68d4508-8749-41fa-ab49-157a18bc04fa
begin
	#	sys = TaxBenefitSystem{Float64}()
	sys = load_file( joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021_22.jl" ))
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021-uplift-removed.jl"))
	# println( "weeklyise start wpm=$PWPM wpy=52")
	weeklyise!( sys )

end

# ╔═╡ f32a80f9-37df-4a36-90f5-8f83f6e3551e
begin
	settings1 = Settings()
	settings2 = Settings()
	settings1.means_tested_routing = lmt_full
	lbc = BCCalcs.makebc(hh, sys, settings1 )
	settings2.means_tested_routing = uc_full
	ubc = BCCalcs.makebc(hh, sys, settings2 )

end

# ╔═╡ 49f2e5eb-71c1-48f0-909b-d85e2670e5f5
begin
	default(fontfamily="Gill Sans",
		titlefont = (12,:grey),
		legendfont = (9),
		guidefont = (10),
		tickfont = (9),
		annotationfontsize=(8),
		annotationcolor=:blue
	  )
	p1=plot( lbc[:,:gross], lbc[:,:net], label="legacy", widen=false, title="UC vs Legacy", ylims = (0,1200), xlims=(0,1200))
	plot!(p1, ubc[:,:gross], ubc[:,:net], label="UC")
end

# ╔═╡ 2108b252-cb89-434b-811a-f222be16bbf3
begin

	println( lbc[16:17, [:gross,:label]])
	println( lbc[16,:].label)
	println( lbc[17,:].label)
	head.income[wages]=484.797
	hres = do_one_calc( hh, sys, settings1 )
	hres.bus[1].pers[head.pid].it
	# head.income[wages]=484.797

end

# ╔═╡ fb9b81a5-91c5-4f9a-b26a-930ae74b07fc
begin

function bcplot( lbc:: DataFrame, ubc :: DataFrame )
	bl = PlotlyJS.scatter(
           lbc, x=:gross, y=:net,
           mode="line", name="Legacy"
       )
	bu = PlotlyJS.scatter(
		ubc, x=:gross, y=:net,
		mode="line", name="UC"
	)
	layout = PlotlyJS.Layout(title="BC leg vs uc",
        xaxis_title="Gross",
        yaxis_title="Net")
	PlotlyJS.plot( [bl, bu], layout)
end

bcplot( lbc, ubc )

end

# ╔═╡ Cell order:
# ╠═8954681e-273b-11ec-344d-9326e736c69f
# ╠═dd1d1115-41b3-4f3f-8cf9-6d97d944b9fa
# ╠═6b4b41a4-2475-4d5d-8ea8-aef6d409c1e6
# ╠═a68d4508-8749-41fa-ab49-157a18bc04fa
# ╠═f32a80f9-37df-4a36-90f5-8f83f6e3551e
# ╠═49f2e5eb-71c1-48f0-909b-d85e2670e5f5
# ╠═2108b252-cb89-434b-811a-f222be16bbf3
# ╠═fb9b81a5-91c5-4f9a-b26a-930ae74b07fc


text:
using VoronoiCells
using Test

@testset "Sort points" begin
    @testset "Average point" begin
        points = [
            GeometryBasics.Point2(1.0, 1.0),
            GeometryBasics.Point2(3.0, 3.0)
        ]

        avg_point = VoronoiCells.mean(points)

        @test VoronoiCells.mean(points) == GeometryBasics.Point2(2.0, 2.0)
    end

    @testset "Sorting points is independent of translation" begin
        radians = range(0, stop = 2π, length = 10)
        points = [GeometryBasics.Point2(cos(θ), sin(θ)) for θ in radians]
        translated_points = [GeometryBasics.Point2(1, 1) + point for point in points]

        @test sortperm(points) == sortperm(translated_points)
    end
end


text:
#=

struct DataLoader <: MapDataset
    ds::MapDataset
    bs::Int
    shuffle::Bool
    num_workers::Int
end

DataLoader(ds::MapDataset, bs=16, shuffle=true, num_workers=0) = DataLoader(ds,bs,shuffle,num_workers)

Base.getindex(dl::DataLoader,idx::Int) = getindex(dl.ds,idx)

Base.length(dl::DataLoader) = length(dl.ds)

one_batch(dl::DataLoader) = dl.bs >= length(dl.ds) ? dl.ds[1:end] : dl.ds[1:dl.bs]

Base.iterate(dl::DataLoader, offset=1) = offset >= length(dl.ds) ? nothing : offset+dl.bs >= length(dl.ds) ? (dl.ds[offset:end],offset+dl.bs) : (dl.ds[offset:offset+dl.bs-1],offset+dl.bs)
=#


text:
module TestDynamics

import Plots
import VegaLite
import PartiallyObservedInverseGames.DynamicsModelInterface
import PartiallyObservedInverseGames.TrajectoryVisualization
import JuMP
using JuMP: @variable, @constraint, @NLconstraint

include("unicycle.jl")
include("product_system.jl")

end


text:
function lgwt(N0::Integer,a::Real=-1,b::Real=1)
    # This script is for computing definite integrals using Legendre-Gauss
    # Quadrature. Computes the Legendre-Gauss nodes and weights on an interval
    # [a,b] with truncation order N
    #
    # Suppose you have a continuous function f(x) which is defined on [a,b]
    # which you can evaluate at any x in [a,b]. Simply evaluate it at all of
    # the values contained in the x vector to obtain a vector f. Then compute
    # the definite integral using sum(f.*w);
    #
    # Written by Greg von Winckel - 02/25/2004
    # First adapted to Julia by Tyler Ransom on 08-04-2015
    # Updated on 08-17-2020

    N  = N0-1
    N1 = N+1
    N2 = N+2

    xu = range(-1,stop=1,length=N1)

    # Initial guess
    y = cos.((2*(0:N) .+ 1)*pi/(2*N .+ 2))  .+  ( 0.27/N1 ) .* sin.( pi .* xu .* N/N2 )

    # Legendre-Gauss Vandermonde Matrix
    L  = zeros(N1,N2)

    # Derivative of LGVM
    Lp = zeros(N1,N2)

    # Compute the zeros of the N+1 Legendre Polynomial
    # using the recursion relation and the Newton-Raphson method

    y0 = 2

    vareps = 2e-52
    # Iterate until new points are uniformly within epsilon of old points

    i = 0
    tracker=0
    it_max=10
    while (norm(y.-y0,Inf)>vareps && tracker<=it_max)
        d=norm(y.-y0,Inf)

        L[:,1]  .= 1
        Lp[:,1] .= 0

        L[:,2] .= y

        for k=2:N1
            L[:,k+1] = ( (2*k-1)*y .* L[:,k] .- (k-1)*L[:,k-1] )/k
        end

        Lp = (N2)*( L[:,N1] .- y .* L[:,N2] )./(1 .- y.^2)
        y0 = y
        y  = y0 - L[:,N2]./Lp
        if norm(y.-y0,Inf)==d
            tracker+=1
        end
        i+=1

    end

    # Linear map from[-1,1] to [a,b]
    x = (a.*(1 .- y) .+ b .* (1 .+ y))./2

    # Compute the weights
    w=(b-a)./((1 .- y.^2).*Lp.^2)*(N2/N1)^2

    return x,w
end

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

# Dr. Ransom Function for mlogit from PS3 solution
function mlogit_with_Z(theta, X, Z, y)

        alpha = theta[1:end-1]
        gamma = theta[end]
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        T = promote_type(eltype(X),eltype(theta))
        num   = zeros(T,N,J)
        dem   = zeros(T,N)
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end


text:
using SpinDoctor
using LinearAlgebra

## Choose setup script
# include("setups/1axon_analytical.jl")
# include("setups/1sphere_analytical.jl")
include("setups/cylinders.jl")
# include("setups/spheres.jl")
# include("setups/neuron.jl")

## Prepare experiments
model = create_model(setup)
volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)

# Sizes
ncompartment = length(model.mesh.points)
namplitude = length(experiment.gradient.values)
nsequence = length(experiment.gradient.sequences)
ndirection = size(experiment.gradient.directions, 2)

qvalues, bvalues = get_values(experiment.gradient)

## Assemble finite element matrices
@time matrices = assemble_matrices(model);

## Use manual time stepping scheme (theta rule)
if !isnothing(experiment.btpde_midpoint)
    @time btpde = solve_btpde_midpoint(model, matrices, experiment)
end

## Solve BTPDE
if !isnothing(experiment.btpde)
    @time btpde = solve_btpde(model, matrices, experiment)

    refinement_str = ""
    if !isnothing(setup.refinement)
        refinement_str = "_refinement$(setup.refinement)"
    end
    output_dir = "output/" * setup.name * refinement_str
    isdir(output_dir) || mkpath(output_dir)

    if experiment.btpde.nsave == 1
        savefield(
            model.mesh,
            btpde.magnetization[:, end, 1, 1],
            "$output_dir/magnetization_btpde",
        )
    else
        save_btpde_results(model.mesh, experiment, btpde, "$output_dir/magnetization_btpde")
    end

    adc = [
        fit_adc(
            bvalues[:, iseq],
            real(btpde.signal[icmpt, :, iseq, idir]) / (model.ρ' * volumes),
        ) for idir = 1:ndirection, iseq = 1:nsequence, icmpt = 1:ncompartment
    ]
end

## Solve MF
if !isnothing(experiment.mf)
    λ_max = length2eig(experiment.mf.length_scale, D_avg)
    lap_eig = compute_laplace_eig(model, matrices, λ_max, experiment.mf.neig_max)
    length_scales = eig2length.(lap_eig.values, D_avg)
    mf = solve_mf(model, matrices, lap_eig, experiment)

    output_dir = "output/$(setup.name)"
    isdir(output_dir) || mkpath(output_dir)
    savefield(model.mesh, mf.magnetization[:, 1, 1, 1], "$output_dir/magnetization_mf")

    npoint_cmpts = size.(model.mesh.points, 2)
    bounds = cumsum([0; npoint_cmpts])
    ϕ_cmpts = [lap_eig.funcs[bounds[i]+1:bounds[i+1], :] for i = 1:ncompartment]
    savefield(model.mesh, ϕ_cmpts, "$output_dir/laplace_eig", "Laplace eigenfunction")
end

## Solve HADC
if !isnothing(experiment.hadc)
    hadc = solve_hadc(model, experiment)
end

## Solve Karger model
if !isnothing(experiment.karger)
    # Fit difftensors (assuming HADC has been solved)
    difftensors = fit_tensors(experiment.gradient.directions, hadc.adc)

    # Solv Karger
    karger = solve_karger(model, experiment, difftensors)
end

## Solve analytical model
if !isnothing(experiment.analytical)
    analytical_signal = solve_analytical(setup, experiment, volumes)
end


text:
using IcanProj
using NetCDF
using NveData
using PyPlot
using DataFrames

function resolution_ind(res, nrow, ncol)
    @assert 50%res==0
    ind = fill(-9999, nrow, ncol)
    nrow = nrow-nrow%res-res+1
    ncol = ncol-ncol%res-res+1
    counter = 1
    for irow = 1:res:nrow, icol = 1:res:ncol
        ind[irow:irow+res-1, icol:icol+res-1] = counter
        counter += 1
    end
    return ind
end

# Settings

path = "/data02/Ican/vic_sim/fsm_past_1km"

# Load soil parameters

soil_param = read_soil_params(path)

# Get information about senorge extent

ind_senorge, xcoord, ycoord = senorge_info()

# Get julia indices

ind_julia = reshape(collect(1:length(ind_senorge[:])), size(ind_senorge))

# Generate indices for coarser grid resolutions

ind_5km = resolution_ind(5, size(ind_senorge,1), size(ind_senorge,2))

ind_10km = resolution_ind(10, size(ind_senorge,1), size(ind_senorge,2))

ind_25km = resolution_ind(25, size(ind_senorge,1), size(ind_senorge,2))

ind_50km = resolution_ind(50, size(ind_senorge,1), size(ind_senorge,2))

ind_5km = convert(Array{Any}, ind_5km)

ind_10km = convert(Array{Any}, ind_10km)

ind_25km = convert(Array{Any}, ind_25km)

ind_50km = convert(Array{Any}, ind_50km)

ind_5km[ind_5km .< 0] = NA

ind_10km[ind_10km .< 0] = NA

ind_25km[ind_25km .< 0] = NA

ind_50km[ind_50km .< 0] = NA

# Data frame from soil parameters

df_left = DataFrame(ind_senorge = soil_param[:gridcel],
                    lon = soil_param[:lon],
                    lat = soil_param[:lat],
                    elev =soil_param[:elev])

df_right = DataFrame(ind_senorge = ind_senorge[:],
                     ind_julia = ind_julia[:],
                     ind_5km = ind_5km[:],
                     ind_10km = ind_10km[:],
                     ind_25km = ind_25km[:],
                     ind_50km = ind_50km[:])

# Join dataframes

df_links = join(df_left, df_right, on=:ind_senorge)

# Remove missing values

completecases!(df_links)

# Convert columns back to integers

df_links[:ind_senorge] = convert(Array{Int64}, df_links[:ind_senorge])

df_links[:ind_5km] = convert(Array{Int64}, df_links[:ind_5km])

df_links[:ind_10km] = convert(Array{Int64}, df_links[:ind_10km])

df_links[:ind_25km] = convert(Array{Int64}, df_links[:ind_25km])

df_links[:ind_50km] = convert(Array{Int64}, df_links[:ind_50km])

# Save to file

file = Pkg.dir("IcanProj", "data", "df_links.csv")

writetable(file, df_links)

#=

# Do some checks

# Some histograms

grid_ids = df_links[:ind_10km]

ngrids = [length(find(grid_ids .== i)) for i in unique(grid_ids)]

plt[:hist](ngrids,10)

# Some maps

imap = df_links[:ind_julia]

tmp = fill(NaN, size(ind_senorge))

tmp[imap] = df_links[:ind_25km]

PyPlot.imshow(tmp)

=#