text:
#=
IMPORTANT: Update the lists below whenever a 1 or 2 qubit gate is added to this file.
           Update QCO._get_angle_gates_idx in src/data.jl if gates with angle parameters are added to this file.
=#

const ONE_QUBIT_GATES_ANGLE_PARAMETERS     = ["U3", "U2", "U1", "RX", "RY", "RZ", "Phase"]

const ONE_QUBIT_GATES_CONSTANTS            = ["Identity", "I", "H", "X", "Y", "Z", "S",
                                              "Sdagger", "T", "Tdagger", "SX", "SXdagger"]

const TWO_QUBIT_GATES_ANGLE_PARAMETERS     = ["CRX", "CRXRev", "CRY", "CRYRev", "CRZ", "CRZRev",
                                              "CU3", "CU3Rev"]

# Gates invariant to qubit flip
const TWO_QUBIT_GATES_CONSTANTS_SYMMETRIC  = ["Swap", "iSwap", "Sycamore", "DCX", "W", "M",
                                              "QFT2", "HCoin", "GroverDiffusion"]

# Gates non-invariant to qubit flip
const TWO_QUBIT_GATES_CONSTANTS_ASYMMETRIC = ["CNot", "CNotRev", "CX", "CXRev", "CY", "CYRev", "CZ",
                                              "CZRev", "CH", "CHRev", "CV", "CVRev", "CVdagger",
                                              "CVdaggerRev", "CSX", "CSXRev"]

const ONE_QUBIT_GATES           = union(QCO.ONE_QUBIT_GATES_CONSTANTS,
                                        QCO.ONE_QUBIT_GATES_ANGLE_PARAMETERS)

const TWO_QUBIT_GATES_CONSTANTS = union(QCO.TWO_QUBIT_GATES_CONSTANTS_SYMMETRIC,
                                        QCO.TWO_QUBIT_GATES_CONSTANTS_ASYMMETRIC)

const TWO_QUBIT_GATES           = union(QCO.TWO_QUBIT_GATES_CONSTANTS,
                                        QCO.TWO_QUBIT_GATES_ANGLE_PARAMETERS)

#----------------------------------------#
#            Single-qubit gates          #
#----------------------------------------#

@doc raw"""
    IGate(num_qubits::Int64)

Identity matrix for an input number of qubits.

**Matrix Representation (num_qubits = 1)**

```math
I = \begin{pmatrix}
        1 & 0 \\
        0 & 1
    \end{pmatrix}
```
"""
function IGate(num_qubits::Int64)

    return Array{Complex{Float64},2}(Matrix(LA.I, 2^num_qubits, 2^num_qubits))

end

@doc raw"""
    U3Gate(θ::Number, ϕ::Number, λ::Number)

Universal single-qubit rotation gate with three Euler angles, ``\theta``, ``\phi`` and ``\lambda``.

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

U3(\theta, \phi, \lambda) =
    \begin{pmatrix}
        \cos(\th)          & -e^{i\lambda}\sin(\th) \\
        e^{i\phi}\sin(\th) & e^{i(\phi+\lambda)}\cos(\th)
    \end{pmatrix}
```
"""
function U3Gate(θ::Number, ϕ::Number, λ::Number)

    QCO._verify_θ_bounds(θ)
    QCO._verify_ϕ_bounds(ϕ)
    QCO._verify_λ_bounds(λ)

    U3 = Array{Complex{Float64},2}([           cos(θ/2)               -(cos(λ) + (sin(λ))im)*sin(θ/2)
                                    (cos(ϕ) + (sin(ϕ))im)*sin(θ/2)  (cos(λ+ϕ) + (sin(λ+ϕ))im)*cos(θ/2)])

    return QCO.round_complex_values(U3)
end

@doc raw"""
    U2Gate(ϕ::Number, λ::Number)

Universal single-qubit rotation gate with two Euler angles, ``\phi`` and ``\lambda``. U2Gate is the special case of
[U3Gate](@ref).

**Matrix Representation**

```math
U2(\phi, \lambda) = \frac{1}{\sqrt{2}}
\begin{pmatrix}
    1          & -e^{i\lambda} \\
    e^{i\phi} & e^{i(\phi+\lambda)}
\end{pmatrix}
```
"""
function U2Gate(ϕ::Number, λ::Number)

    θ = π/2

    U2 = QCO.U3Gate(θ, ϕ, λ)
    return U2
end

@doc raw"""
    U1Gate(λ::Number)

Universal single-qubit rotation gate with one Euler angle, ``\lambda``. U1Gate represents rotation about the Z axis and
is the special case of [U3Gate](@ref), which also known as the [PhaseGate](@ref). Also note that ``U1(\pi) = ``[ZGate](@ref), ``U1(\pi/2) = ``[SGate](@ref) and
``U1(\pi/4) = ``[TGate](@ref).

**Matrix Representation**

```math
U1(\lambda) =
\begin{pmatrix}
    1 & 0 \\
    0 & e^{i\lambda}
\end{pmatrix}
```
"""
function U1Gate(λ::Number)

    θ = 0
    ϕ = 0

    U1 = QCO.U3Gate(θ, ϕ, λ)
    return U1
end

@doc raw"""
    RXGate(θ::Number)

A single-qubit Pauli gate which represents rotation about the X axis.

**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}

RX(\theta) = exp(-i \th X) =
    \begin{pmatrix}
        \cos{\th}   & -i\sin{\th} \\
        -i\sin{\th} & \cos{\th}
    \end{pmatrix}
```
"""
function RXGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    RX = Array{Complex{Float64},2}([cos(θ/2) -(sin(θ/2))im; -(sin(θ/2))im cos(θ/2)])

    return QCO.round_complex_values(RX)
end

@doc raw"""
    RYGate(θ::Number)

A single-qubit Pauli gate which represents rotation about the Y axis.

**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}

RY(\theta) = exp(-i \th Y) =
    \begin{pmatrix}
        \cos{\th} & -\sin{\th} \\
        \sin{\th} & \cos{\th}
    \end{pmatrix}
```
"""
function RYGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    RY = Array{Complex{Float64},2}([cos(θ/2) -(sin(θ/2)); (sin(θ/2)) cos(θ/2)])

    return QCO.round_complex_values(RY)
end

@doc raw"""
    RZGate(θ::Number)

A single-qubit Pauli gate which represents rotation about the Z axis. This gate is also equivalent to [U1Gate](@ref) up to a phase factor,
that is, ``RZ(\theta) = e^{-i{\theta}/2}U1(\theta)``.

**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}

RZ(\theta) = exp(-i\th Z) =
\begin{pmatrix}
    e^{-i\th} & 0 \\
    0 & e^{i\th}
\end{pmatrix}
```
"""
function RZGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    RZ = Array{Complex{Float64},2}([(cos(θ/2) - (sin(θ/2))im) 0; 0 (cos(θ/2) + (sin(θ/2))im)])

    return QCO.round_complex_values(RZ)
end

@doc raw"""
    HGate(num_qubits::Int64)

Single-qubit Hadamard gate, which is a ``\pi`` rotation about the X+Z axis, thus equivalent to [U3Gate](@ref)(``\frac{\pi}{2},0,\pi``)

**Matrix Representation**

```math
H = \frac{1}{\sqrt{2}}
        \begin{pmatrix}
            1 & 1 \\
            1 & -1
        \end{pmatrix}
```
"""
function HGate()

    return Array{Complex{Float64},2}(1/sqrt(2)*[1 1; 1 -1])

end

@doc raw"""
    XGate()

Single-qubit Pauli-X gate (``\sigma_x``), equivalent to [U3Gate](@ref)(``\pi,0,\pi``)

**Matrix Representation**

```math
X = \begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}
```
"""
function XGate()

    return Array{Complex{Float64},2}([0 1; 1 0])

end

@doc raw"""
    YGate()

Single-qubit Pauli-Y gate (``\sigma_y``), equivalent to [U3Gate](@ref)(``\pi,\frac{\pi}{2},\frac{\pi}{2}``)

**Matrix Representation**

```math
Y = \begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix}
```
"""
function YGate()

    return Array{Complex{Float64},2}([0 -im; im 0])

end

@doc raw"""
    ZGate()

Single-qubit Pauli-Z gate (``\sigma_z``), equivalent to [U3Gate](@ref)(``0,0,\pi``)

**Matrix Representation**

```math
Z = \begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}
```
"""
function ZGate()

    return Array{Complex{Float64},2}([1 0; 0 -1])

end

@doc raw"""
    SGate()

Single-qubit S gate, equivalent to [U3Gate](@ref)(``0,0,\frac{\pi}{2}``). This
gate is also referred to as a Clifford gate, P gate or a square-root of Pauli-[ZGate](@ref). Historically, this is also
called as the phase gate (denoted by P), since it shifts the phase of the one state relative to the zero state.

**Matrix Representation**

```math
S = \begin{pmatrix}
1 & 0 \\
0 & i
\end{pmatrix}
```
"""
function SGate()

    return Array{Complex{Float64},2}([1 0; 0 im])

end

@doc raw"""
    SdaggerGate()

Single-qubit, hermitian conjugate of the [SGate](@ref). This is also an alternative square root of
the [ZGate](@ref).

**Matrix Representation**

```math
S = \begin{pmatrix}
1 & 0 \\
0 & -i
\end{pmatrix}
```
"""
function SdaggerGate()

    return Array{Complex{Float64},2}([1 0; 0 -im])

end

@doc raw"""
    TGate()

Single-qubit T gate, equivalent to [U3Gate](@ref)(``0,0,\frac{\pi}{4}``). This
gate is also referred to as a ``\frac{\pi}{8}`` gate or as a fourth-root of Pauli-[ZGate](@ref).

**Matrix Representation**

```math
T = \begin{pmatrix}
1 & 0 \\
0 & e^{i\pi/4}
\end{pmatrix}
```
"""
function TGate()

    return Array{Complex{Float64},2}([1 0; 0 (1/sqrt(2)) + (1/sqrt(2))im])

end

@doc raw"""
    TdaggerGate()

Single-qubit, hermitian conjugate of the [TGate](@ref). This gate is equivalent to [U3Gate](@ref)(``0,0,-\frac{\pi}{4}``). This
gate is also referred to as the fourth-root of Pauli-[ZGate](@ref).

**Matrix Representation**

```math
T^{\dagger} = \begin{pmatrix}
1 & 0 \\
0 & e^{-i\pi/4}
\end{pmatrix}
```
"""
function TdaggerGate()

    return Array{Complex{Float64},2}([1 0; 0 (1/sqrt(2)) - (1/sqrt(2))im])

end

@doc raw"""
    SXGate()

Single-qubit square root of pauli-[XGate](@ref).

**Matrix Representation**

```math
\sqrt{X} = \frac{1}{2} \begin{pmatrix}
1 + i & 1 - i \\
1 - i & 1 + i
\end{pmatrix}
```
"""
function SXGate()

    return Array{Complex{Float64},2}(1/2*[1+im 1-im; 1-im 1+im])

end

@doc raw"""
    SXdaggerGate()

Single-qubit hermitian conjugate of the square root of pauli-[XGate](@ref), or the [SXGate](@ref).

**Matrix Representation**

```math
\sqrt{X}^{\dagger} = \frac{1}{2} \begin{pmatrix}
1 - i & 1 + i \\
1 + i & 1 - i
\end{pmatrix}
```
"""
function SXdaggerGate()

    return Array{Complex{Float64},2}(1/2*[1-im 1+im; 1+im 1-im])

end

@doc raw"""
    PhaseGate()

Single-qubit rotation gate about the Z axis. This is also equivalent to [U3Gate](@ref)(``0,0,\lambda``). This
gate is also referred to as the [U1Gate](@ref).

**Matrix Representation**

```math
P(\lambda) = \begin{pmatrix}
    1 & 0 \\
    0 & e^{i\lambda}
\end{pmatrix}
```
"""
function PhaseGate(λ::Number)
    return QCO.U1Gate(λ)
end

#-------------------------------------#
#            Two-qubit gates          #
#-------------------------------------#
@doc raw"""
    CNotGate()

Two-qubit controlled NOT gate with control and target on first and second qubits, respectively. This is also
called the controlled X gate ([CXGate](@ref)).

**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ X ├
     └───┘
```

**Matrix Representation**

```math
CNot = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0
    \end{pmatrix}
```
"""
function CNotGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0])

end

@doc raw"""
    CNotRevGate()

Two-qubit reverse controlled NOT gate, with target and control on first and second qubits, respectively.

**Circuit Representation**
```
     ┌───┐
q_0: ┤ X ├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CNotRev = \begin{pmatrix}
            1 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 \\
            0 & 0 & 1 & 0 \\
            0 & 1 & 0 & 0
            \end{pmatrix}
```
"""
function CNotRevGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0])

end

@doc raw"""
    DCXGate()

Two-qubit double controlled NOT gate consisting of two back-to-back [CNotGate](@ref)s with alternate controls.

**Circuit Representation**
```
          ┌───┐
q_0: ──■──┤ X ├
     ┌─┴─┐└─┬─┘
q_1: ┤ X ├──■──
     └───┘
```

**Matrix Representation**

```math
DCX = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 1 & 0
    \end{pmatrix}
```
"""
function DCXGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 0 0 1; 0 1 0 0; 0 0 1 0])

end

@doc raw"""
    CXGate()

Two-qubit controlled [XGate](@ref), which is also the same as [CNotGate](@ref).

**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ X ├
     └───┘
```

**Matrix Representation**

```math
CX = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0
    \end{pmatrix}
```
"""
function CXGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0])

end

@doc raw"""
    CXRevGate()

Two-qubit reverse controlled-X gate, with target and control on first and second qubits, respectively.
This is also the same as [CNotRevGate](@ref).

**Circuit Representation**
```
     ┌───┐
q_0: ┤ X ├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CXRev = I \otimes |0 \rangle\langle 0| + X \otimes |1 \rangle\langle 1| = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 0 & 0 & 1 \\
        0 & 0 & 1 & 0 \\
        0 & 1 & 0 & 0
    \end{pmatrix}
```
"""
function CXRevGate()

    # I ⊗ |0⟩⟨0|
    control_0 = kron(QCO.IGate(1), Array{Complex{Float64},2}([1 0; 0 0]))
    # X ⊗ |1⟩⟨1|
    control_1 = kron(QCO.XGate(), Array{Complex{Float64},2}([0 0; 0 1]))

    return control_0 + control_1
end

@doc raw"""
    CYGate()

Two-qubit controlled [YGate](@ref).

**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ Y ├
     └───┘
```

**Matrix Representation**

```math
CY = |0 \rangle\langle 0| \otimes I + |1 \rangle\langle 1| \otimes Y = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & -i \\
    0 & 0 & i & 0
    \end{pmatrix}
```
"""
function CYGate()

    # |0⟩⟨0| ⊗ I
    control_0 = kron(Array{Complex{Float64},2}([1 0; 0 0]) , QCO.IGate(1))
    # |1⟩⟨1| ⊗ Y
    control_1 = kron(Array{Complex{Float64},2}([0 0; 0 1]) , QCO.YGate())

    return control_0 + control_1
end

@doc raw"""
    CYRevGate()

Two-qubit reverse controlled-Y gate, with target and control on first and second qubits, respectively.

**Circuit Representation**
```
     ┌───┐
q_0: ┤ Y ├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CYRev = I \otimes |0 \rangle\langle 0| + Y \otimes |1 \rangle\langle 1| = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 0 & 0 & -i \\
        0 & 0 & 1 & 0 \\
        0 & i & 0 & 0
    \end{pmatrix}
```
"""
function CYRevGate()

    # I ⊗ |0⟩⟨0|
    control_0 = kron(QCO.IGate(1), Array{Complex{Float64},2}([1 0; 0 0]))
    # Y ⊗ |1⟩⟨1|
    control_1 = kron(QCO.YGate(), Array{Complex{Float64},2}([0 0; 0 1]))

    return control_0 + control_1
end

@doc raw"""
    CZGate()

Two-qubit, symmetric, controlled [ZGate](@ref).

**Circuit Representation**
```
q_0: ──■──     ─■─
     ┌─┴─┐  ≡   │
q_1: ┤ Z ├     ─■─
     └───┘
```

**Matrix Representation**

```math
CZ = |0 \rangle\langle 0| \otimes I + |1 \rangle\langle 1| \otimes Z = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 1 & 0 \\
    0 & 0 & 0 & -1
    \end{pmatrix}
```
"""
function CZGate()

    # |0⟩⟨0| ⊗ I
    control_0 = kron(Array{Complex{Float64},2}([1 0; 0 0]) , QCO.IGate(1))
    # |1⟩⟨1| ⊗ Z
    control_1 = kron(Array{Complex{Float64},2}([0 0; 0 1]) , QCO.ZGate())

    return control_0 + control_1
end

@doc raw"""
    CZRevGate()

Two-qubit reverse controlled-Z gate, with target and control on first and second qubits, respectively.

**Circuit Representation**
```
     ┌───┐
q_0: ┤ Z ├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CZRev = I \otimes |0\rangle\langle0| + Z \otimes |1\rangle\langle1| = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 1 & 0 \\
    0 & 0 & 0 & -1
    \end{pmatrix}
```
"""
function CZRevGate()

    # I ⊗ |0⟩⟨0|
    control_0 = kron(QCO.IGate(1), Array{Complex{Float64},2}([1 0; 0 0]))
    # Z ⊗ |1⟩⟨1|
    control_1 = kron(QCO.ZGate(), Array{Complex{Float64},2}([0 0; 0 1]))

    return control_0 + control_1
end

@doc raw"""
    CHGate()

Two-qubit, symmetric, controlled Hadamard gate ([HGate](@ref)).

**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ H ├
     └───┘
```

**Matrix Representation**

```math
CH = |0\rangle\langle 0| \otimes I + |1\rangle\langle 1| \otimes H = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\
0 & 0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}
\end{pmatrix}
```
"""
function CHGate()

    # |0⟩⟨0| ⊗ I
    control_0 = kron(Array{Complex{Float64},2}([1 0; 0 0]) , QCO.IGate(1))
    # |1⟩⟨1| ⊗ H
    control_1 = kron(Array{Complex{Float64},2}([0 0; 0 1]) , QCO.HGate())

    return control_0 + control_1
end

@doc raw"""
    CHRevGate()

Two-qubit reverse controlled-H gate, with target and control on first and second qubits, respectively.

**Circuit Representation**
```
     ┌───┐
q_0: ┤ H ├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CHRev = I \otimes |0\rangle\langle 0| + H \otimes |1\rangle\langle 1| = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & \frac{1}{\sqrt{2}} & 0 & \frac{1}{\sqrt{2}} \\
        0 & 0 & 1 & 0 \\
        0 & \frac{1}{\sqrt{2}} & 0 & -\frac{1}{\sqrt{2}}
    \end{pmatrix}
```
"""
function CHRevGate()

    # I ⊗ |0⟩⟨0|
    control_0 = kron(QCO.IGate(1), Array{Complex{Float64},2}([1 0; 0 0]))
    # H ⊗ |1⟩⟨1|
    control_1 = kron(QCO.HGate(), Array{Complex{Float64},2}([0 0; 0 1]))

    return control_0 + control_1
end

@doc raw"""
    CVGate()

Two-qubit, controlled-V gate, which is also the same as Controlled square-root of X gate ([CSXGate](@ref)).

**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ V ├
     └───┘
```

**Matrix Representation**

```math
CV = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & 0.5+0.5i & 0.5-0.5i \\
        0 & 0 & 0.5-0.5i & 0.5+0.5i
    \end{pmatrix}
```
"""
function CVGate()

    return QCO.CSXGate()

end

@doc raw"""
    CVRevGate()

Two-qubit reverse controlled-V gate, with target and control on first and second qubits, respectively.

**Circuit Representation**
```
     ┌───┐
q_0: ┤ V ├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CVRev = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 0.5+0.5i & 0 & 0.5-0.5i \\
        0 & 0 & 1 & 0 \\
        0 & 0.5-0.5i & 0 & 0.5+0.5i
    \end{pmatrix}
```
"""
function CVRevGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 0.5+0.5im 0 0.5-0.5im; 0 0 1 0; 0 0.5-0.5im 0 0.5+0.5im])

end

@doc raw"""
    CVdaggerGate()

Two-qubit hermitian conjugate of controlled-V gate, which is also the same as hermitian conjugate Controlled square-root of X gate ([CSXGate](@ref)).

**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ V'├
     └───┘
```

**Matrix Representation**

```math
CVdagger = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & 0.5-0.5i & 0.5+0.5i \\
        0 & 0 & 0.5+0.5i & 0.5-0.5i
    \end{pmatrix}
```
"""
function CVdaggerGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 1 0 0; 0 0 0.5-0.5im 0.5+0.5im; 0 0 0.5+0.5im 0.5-0.5im])

end

@doc raw"""
    CVdaggerRevGate()

Two-qubit hermitian conjugate of reverse controlled-V gate, with target and control on first and second qubits, respectively.

**Circuit Representation**
```
     ┌───┐
q_0: ┤ V'├
     └─┬─┘
q_1: ──■──
```

**Matrix Representation**

```math
CVdaggerRev = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 0.5-0.5i & 0 & 0.5+0.5i \\
        0 & 0 & 1 & 0 \\
        0 & 0.5+0.5i & 0 & 0.5-0.5i
    \end{pmatrix}
```
"""
function CVdaggerRevGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 0.5-0.5im 0 0.5+0.5im; 0 0 1 0; 0 0.5+0.5im 0 0.5-0.5im])

end

@doc raw"""
    WGate()

Two-qubit, W hermitian gate, typically useful to diagonlize the ([SwapGate](@ref)).

**Matrix Representation**

```math
W = \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 \\
        0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} & 0 \\
        0 & 0 & 0 & 1
    \end{pmatrix}
```
"""
function WGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 1/sqrt(2) 1/sqrt(2) 0; 0 1/sqrt(2) -1/sqrt(2) 0; 0 0 0 1])

end

@doc raw"""
    CRXGate(θ::Number)

Two-qubit controlled [RXGate](@ref).

**Circuit Representation**
```
q_0: ────■────
     ┌───┴───┐
q_1: ┤ RX(ϴ) ├
     └───────┘
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CRX(\theta)\ q_1, q_0 =
|0\rangle\langle0| \otimes I + |1\rangle\langle1| \otimes RX(\theta) =
    \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & \cos{\th}   & -i\sin{\th} \\
        0 & 0 & -i\sin{\th} & \cos{\th}
    \end{pmatrix}
```
"""
function CRXGate(θ::Number)

    CRX = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 1 0 0
                                      0 0  cos(θ/2)  -(sin(θ/2))im
                                      0 0  -(sin(θ/2))im  cos(θ/2)])

    return QCO.round_complex_values(CRX)
end

@doc raw"""
    CRXRevGate(θ::Number)

Two-qubit controlled reverse [RXGate](@ref).

**Circuit Representation**
```
     ┌───────┐
q_1: ┤ RX(ϴ) ├
     └───┬───┘
q_0: ────■────
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CRXRev(\theta)\ q_1, q_0 =
|0\rangle\langle0| \otimes I + |1\rangle\langle1| \otimes RX(\theta) =
    \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & \cos{\th} & 0 & -i\sin{\th} \\
        0 & 0 & 1 & 0\\
        0 & -i\sin{\th} & 0 & \cos{\th}
    \end{pmatrix}
```
"""
function CRXRevGate(θ::Number)

    CRXRev = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 cos(θ/2) 0 -(sin(θ/2))im
                                      0 0 1 0
                                      0 -(sin(θ/2))im 0 cos(θ/2)])

    return QCO.round_complex_values(CRXRev)
end

@doc raw"""
    CRYGate(θ::Number)

Two-qubit controlled [RYGate](@ref).

**Circuit Representation**
```
q_0: ────■────
     ┌───┴───┐
q_1: ┤ RY(ϴ) ├
     └───────┘
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CRY(\theta)\ q_1, q_0 =
|0\rangle\langle0| \otimes I + |1\rangle\langle1| \otimes RY(\theta) =
    \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & \cos{\th}   & -\sin{\th} \\
        0 & 0 & \sin{\th} & \cos{\th}
    \end{pmatrix}
```
"""
function CRYGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    CRY = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 1 0 0
                                      0 0  cos(θ/2)  -(sin(θ/2))
                                      0 0  (sin(θ/2))  cos(θ/2)])

    return QCO.round_complex_values(CRY)
end

@doc raw"""
    CRYRevGate(θ::Number)

Two-qubit controlled reverse [RYGate](@ref).

**Circuit Representation**
```
     ┌───────┐
q_1: ┤ RY(ϴ) ├
     └───┬───┘
q_0: ────■────
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CRYRev(\theta)\ q_1, q_0 =
|0\rangle\langle0| \otimes I + |1\rangle\langle1| \otimes RY(\theta) =
    \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & \cos{\th} & 0 & -\sin{\th} \\
        0 & 0 & 1 & 0 \\
        0 & \sin{\th} & 0 & \cos{\th}
    \end{pmatrix}
```
"""
function CRYRevGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    CRYRev = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 cos(θ/2) 0 -(sin(θ/2))
                                      0 0 1 0
                                      0 (sin(θ/2)) 0 cos(θ/2)])

    return QCO.round_complex_values(CRYRev)
end

@doc raw"""
    CRZGate(θ::Number)

Two-qubit controlled [RZGate](@ref).

**Circuit Representation**
```
q_0: ────■────
     ┌───┴───┐
q_1: ┤ RZ(ϴ) ├
     └───────┘
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CRZ(\theta)\ q_1, q_0 =
|0\rangle\langle0| \otimes I + |1\rangle\langle1| \otimes RZ(\theta) =
    \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & e^{-i\th}  &  0 \\
        0 & 0 & 0 & e^{i\th}
    \end{pmatrix}
```
"""
function CRZGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    CRZ = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 1 0 0
                                      0 0  (cos(θ/2) - (sin(θ/2))im)  0
                                      0 0  0  (cos(θ/2) + (sin(θ/2))im)])

    return QCO.round_complex_values(CRZ)
end

@doc raw"""
    CRZRevGate(θ::Number)

Two-qubit controlled reverse [RZGate](@ref).

**Circuit Representation**
```
     ┌───────┐
q_1: ┤ RZ(ϴ) ├
     └───┬───┘
q_0: ────■────
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CRZRev(\theta)\ q_1, q_0 =
|0\rangle\langle0| \otimes I + |1\rangle\langle1| \otimes RZ(\theta) =
    \begin{pmatrix}
        1 & 0 & 0 & 0 \\
        0 & e^{-i\th} & 0 & 0 \\
        0 & 0 & 1 &  0 \\
        0 & 0 & 0 & e^{i\th}
    \end{pmatrix}
```
"""
function CRZRevGate(θ::Number)

    QCO._verify_θ_bounds(θ)

    CRZRev = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 (cos(θ/2) - (sin(θ/2))im) 0 0
                                      0 0 1 0
                                      0 0  0  (cos(θ/2) + (sin(θ/2))im)])

    return QCO.round_complex_values(CRZRev)
end

@doc raw"""
    CU3Gate(θ::Number, ϕ::Number, λ::Number)

Two-qubit, controlled version of the universal rotation gate with three Euler angles ([U3Gate](@ref)).

**Circuit Representation**
```
q_0: ──────■──────
     ┌─────┴─────┐
q_1: ┤ U3(ϴ,φ,λ) ├
     └───────────┘
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CU3(\theta, \phi, \lambda)\ q_1, q_0 =
                |0\rangle\langle 0| \otimes I +
                |1\rangle\langle 1| \otimes U3(\theta,\phi,\lambda) =
                \begin{pmatrix}
                    1 & 0   & 0                  & 0 \\
                    0 & 1   & 0                  & 0 \\
                    0 & 0   & \cos(\th)          & -e^{i\lambda}\sin(\th) \\
                    0 & 0   & e^{i\phi}\sin(\th) & e^{i(\phi+\lambda)}\cos(\th)
                \end{pmatrix}
```
"""
function CU3Gate(θ::Number, ϕ::Number, λ::Number)

    QCO._verify_θ_bounds(θ)
    QCO._verify_ϕ_bounds(ϕ)
    QCO._verify_λ_bounds(λ)

    CU3 = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 1 0 0
                                      0 0  cos(θ/2)  -(cos(λ)+(sin(λ))im)*sin(θ/2)
                                      0 0 (cos(ϕ)+(sin(ϕ))im)*sin(θ/2)  (cos(λ+ϕ)+(sin(λ+ϕ))im)*cos(θ/2)])

    return QCO.round_complex_values(CU3)
end

@doc raw"""
    CU3RevGate(θ::Number, ϕ::Number, λ::Number)

Two-qubit, reverse controlled version of the universal rotation gate with three Euler angles ([U3Gate](@ref)).

**Circuit Representation**
```
     ┌────────────┐
q_1: ┤  U3(ϴ,φ,λ) ├
     └──────┬─────┘
q_0: ───────■──────
```

**Matrix Representation**

```math
\newcommand{\th}{\frac{\theta}{2}}

CU3(\theta, \phi, \lambda)\ q_1, q_0 =
                |0\rangle\langle 0| \otimes I +
                |1\rangle\langle 1| \otimes U3(\theta,\phi,\lambda) =
                \begin{pmatrix}
                    1 & 0   & 0  & 0 \\
                    0 & \cos(\th)   & 0 & -e^{i\lambda}\sin(\th) \\
                    0 & 0   &  1 & 0 \\
                    0 & e^{i\phi}\sin(\th)  & 0  & e^{i(\phi+\lambda)}\cos(\th)
                \end{pmatrix}
```
"""
function CU3RevGate(θ::Number, ϕ::Number, λ::Number)

    QCO._verify_θ_bounds(θ)
    QCO._verify_ϕ_bounds(ϕ)
    QCO._verify_λ_bounds(λ)

    CU3Rev = Array{Complex{Float64},2}([ 1 0 0 0
                                      0 cos(θ/2) 0 -(cos(λ)+(sin(λ))im)*sin(θ/2)
                                      0 0 1 0
                                      0 (cos(ϕ)+(sin(ϕ))im)*sin(θ/2) 0 (cos(λ+ϕ)+(sin(λ+ϕ))im)*cos(θ/2)])

    return QCO.round_complex_values(CU3Rev)
end

@doc raw"""
    SwapGate()

Two-qubit, symmetric, SWAP gate.

**Circuit Representation**
```
q_0: ─X─
      │
q_1: ─X─
```

**Matrix Representation**

```math
SWAP = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
```
"""
function SwapGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1])

end

@doc raw"""
    iSwapGate()

Two-qubit, symmetric and clifford, iSWAP gate. This is an entangling swapping gate where the qubits
obtain a phase of ``i`` if the state of the qubits is swapped.

**Circuit Representation**
```
q_0: ─⨂─
      │
q_1: ─⨂─
```
Minimum depth representation
```
      ┌───┐     ┌───┐ ┌───┐
q_0: ─┤ X ├──■──┤ S ├─┤ X ├─
      └─┬─┘┌─┴─┐└───┘ └─┬─┘
q_1: ───■──┤ X ├────────■──
           └───┘
```

**Matrix Representation**

```math
iSWAP = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & i & 0 \\
0 & i & 0 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
```
"""
function iSwapGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1])

end

@doc raw"""
    CSXGate()

Two-qubit controlled version of ([SXGate](@ref)).

**Circuit Representation**
```
q_0: ─────■─────
     ┌────┴────┐
q_1: ┤ sqrt(X) ├
     └─────────┘
```

**Matrix Representation**

```math
CSXGate = |0 \rangle\langle 0| \otimes I + |1 \rangle\langle 1| \otimes SX = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0.5+0.5i & 0.5-0.5i \\
0 & 0 & 0.5-0.5i & 0.5+0.5i
    \end{pmatrix}
```
"""
function CSXGate()

    # |0⟩⟨0| ⊗ I
    control_0 = kron(Array{Complex{Float64},2}([1 0; 0 0]) , QCO.IGate(1))
    # |1⟩⟨1| ⊗ SX
    control_1 = kron(Array{Complex{Float64},2}([0 0; 0 1]) , QCO.SXGate())

    return control_0 + control_1
end

@doc raw"""
    CSXRevGate()

Two-qubit controlled version of the reverse ([SXGate](@ref)).

**Circuit Representation**
```
     ┌─────────┐
q_1: ┤ sqrt(X) ├
     └────┬────┘
q_0: ─────■────
```

**Matrix Representation**

```math
CSXRevGate = I \otimes |0\rangle\langle 0| + SX \otimes |1\rangle\langle 1| = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 0.5+0.5i & 0 & 0.5-0.5i \\
0 & 0 & 1 & 0 \\
0 & 0.5-0.5i & 0 & 0.5+0.5i
    \end{pmatrix}
```
"""
function CSXRevGate()

    # I ⊗ |0⟩⟨0|
    control_0 = kron(QCO.IGate(1), Array{Complex{Float64},2}([1 0; 0 0]))
    # SX ⊗ |1⟩⟨1|
    control_1 = kron(QCO.SXGate(), Array{Complex{Float64},2}([0 0; 0 1]))

    return control_0 + control_1
end

@doc raw"""
    MGate()

Two-qubit Magic gate, also known as the Ising coupling or the XX gate.

Reference: [https://doi.org/10.1103/PhysRevA.69.032315](https://doi.org/10.1103/PhysRevA.69.032315)

**Circuit Representation**
```
      ┌───┐        ┌───┐
q_0: ─┤ X ├────────┤ S ├
      └─┬─┘        └─┬─┘
        │   ┌───┐  ┌─┴─┐
q_1: ───■───┤ H ├──┤ S ├
            └───┘  └───┘
```

**Matrix Representation**

```math
M = \frac{1}{\sqrt{2}} \begin{pmatrix}
1 & i & 0 & 0 \\
0 & 0 & i & 1 \\
0 & 0 & i & -1 \\
1 & -i & 0 & 0
\end{pmatrix}
```
"""
function MGate()

    return Array{Complex{Float64},2}(1/sqrt(2)*[1 im 0 0; 0 0 im 1; 0 0 im -1; 1 -im 0 0])

end

@doc raw"""
    QFT2Gate()

Two-qubit Quantum Fourier Transform (QFT) gate, where the QFT operation on n-qubits is given by:
```math
|j\rangle \mapsto \frac{1}{2^{n/2}} \sum_{k=0}^{2^n - 1} e^{2\pi ijk / 2^n} |k\rangle
```

**Circuit Representation**
```
     ┌──────┐
q_0: ┤      ├
     │ QFT2 │
q_1: ┤      ├
     └──────┘
```

**Matrix Representation**

```math
M = \frac{1}{2} \begin{pmatrix}
1 & 1 & 1 & 1 \\
1 & i & -1 & -i \\
1 & -1 & 1 & -1 \\
1 & -i & -1 & i
\end{pmatrix}
```
"""
function QFT2Gate()

    return Array{Complex{Float64},2}(0.5*[1 1 1 1; 1 im -1 -im; 1 -1 1 -1; 1 -im -1 im])

end

@doc raw"""
    HCoinGate()

Two-qubit, Hadamard Coin gate when implemented in tune with the quantum cellular automata.
Reference: [https://doi.org/10.1007/s11128-018-1983-x](https://doi.org/10.1007/s11128-018-1983-x), [https://arxiv.org/pdf/2106.03115.pdf](https://arxiv.org/pdf/2106.03115.pdf)

**Matrix Representation**

```math
HCoinGate = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 \\
0 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} & 0 \\
0 & 0 & 0 & 1
\end{pmatrix}
```
"""
function HCoinGate()

    return Array{Complex{Float64},2}([1 0 0 0; 0 1/sqrt(2) 1/sqrt(2) 0; 0 1/sqrt(2) -1/sqrt(2) 0; 0 0 0 1])

end

@doc raw"""
    GroverDiffusionGate()

Two-qubit, Grover's diffusion operator, a key building block of the Glover's algorithm used to find a specific
item (with probability > 0.5) within a randomly ordered database of N items in O(sqrt(N)) operations.
Reference: [https://arxiv.org/pdf/1804.03719.pdf](https://arxiv.org/pdf/1804.03719.pdf)

**Matrix Representation**

```math
GroverDiffusionGate = \frac{1}{2}\begin{pmatrix}
1 & -1 & -1 & -1 \\
-1 & 1 & -1 & -1 \\
-1 & -1 & 1 & -1 \\
-1 & -1 & -1 & 1
\end{pmatrix}
```
"""
function GroverDiffusionGate()

    return Array{Complex{Float64},2}(0.5*[1 -1 -1 -1; -1 1 -1 -1; -1 -1 1 -1; -1 -1 -1 1])

end

@doc raw"""
    SycamoreGate()

Two-qubit Sycamore Gate, native to Google's universal quantum processor.
Reference: [quantumai.google/cirq/google/devices](https://quantumai.google/cirq/google/devices)

**Circuit Representation**
```
     ┌──────┐
q_0: ┤      ├
     │ SYC  │
q_1: ┤      ├
     └──────┘
```

**Matrix Representation**

```math

SycamoreGate() = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & -i & 0 \\
0 & -i & 0 & 0 \\
0 & 0 & 0 & e^{-i \frac{\pi}{6}}
\end{pmatrix}

```
"""
function SycamoreGate()

    return  Array{Complex{Float64},2}([ 1 0 0 0
                                      0 0 -im 0
                                      0 -im 0 0
                                      0 0 0 cos(pi/6)-(sin(pi/6))im])

end

#---------------------------------------#
#            Three-qubit gates          #
#---------------------------------------#

@doc raw"""
    ToffoliGate()

Three-qubit Toffoli gate, also known as the CCX (controlled-controlled-NOT) gate.

**Circuit Representation**
```
q_0: ──■──
       │
q_1: ──■──
     ┌─┴─┐
q_2: ┤ X ├
     └───┘
```

**Matrix Representation**

```math
Toffoli     =
            |0 \rangle \langle 0| \otimes I \otimes I + |1 \rangle \langle 1| \otimes CXGate =
            \begin{pmatrix}
                1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\
                0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\
                0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\
                0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
                0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\
                0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\
                0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\
                0 & 0 & 0 & 0 & 0 & 0 & 1 & 0
            \end{pmatrix}
```
"""
function ToffoliGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  1  0  0  0
                                      0  0  0  0  0  1  0  0
                                      0  0  0  0  0  0  0  1
                                      0  0  0  0  0  0  1  0])

end

@doc raw"""
    CSwapGate()

Three-qubit, controlled [SwapGate](@ref), also known as the [Fredkin gate](https://en.wikipedia.org/wiki/Fredkin_gate).

**Circuit Representation**
```
q_0: ─■─
      │
q_1: ─X─
      │
q_2: ─X─
```

**Matrix Representation**

```math
CSwapGate =
            \begin{pmatrix}
            1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
            0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
        \end{pmatrix}
```
"""
function CSwapGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  1  0  0  0
                                      0  0  0  0  0  0  1  0
                                      0  0  0  0  0  1  0  0
                                      0  0  0  0  0  0  0  1])

end

@doc raw"""
    CCZGate()

Three-qubit controlled-controlled Z gate.

**Circuit Representation**
```
q_0: ─■─
      │
q_1: ─■─
      │
q_2: ─■─
```

**Matrix Representation**

```math
CCZGate =
            \begin{pmatrix}
            1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 \\
        \end{pmatrix}
```
"""
function CCZGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  1  0  0  0
                                      0  0  0  0  0  1  0  0
                                      0  0  0  0  0  0  1  0
                                      0  0  0  0  0  0  0  -1])

end

@doc raw"""
    PeresGate()

Three-qubit Peres gate. This gate is equivalent to [ToffoliGate](@ref) followed by the [CNotGate](@ref) in 3 qubits.
Reference: [https://doi.org/10.1103/PhysRevA.32.3266](https://doi.org/10.1103/PhysRevA.32.3266)

**Circuit Representation**
```
q_0: ──■─────■──
       │   ┌─┴─┐
q_1: ──■───┤ X ├
     ┌─┴─┐ └───┘
q_2: ┤ X ├──────
     └───┘
```

**Matrix Representation**

```math
PeresGate =
            \begin{pmatrix}
            1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
            0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
            0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
        \end{pmatrix}
```
"""
function PeresGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  0  0  1  0
                                      0  0  0  0  0  0  0  1
                                      0  0  0  0  0  1  0  0
                                      0  0  0  0  1  0  0  0])

end

@doc raw"""
    RCCXGate()

Three-qubit relative (or simplified) Toffoli gate, or the CCX gate. This gate is equivalent to [ToffoliGate](@ref) upto relative phases.
The advantage of this gate is that it's implementation requires only three CNot (or CX) gates.
Reference: [https://arxiv.org/pdf/1508.03273.pdf](https://arxiv.org/pdf/1508.03273.pdf)

**Matrix Representation**

```math
RCCXGate =
            \begin{pmatrix}
            1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 0 & -i \\
            0 & 0 & 0 & 0 & 0 & 0 & i & 0 \\
        \end{pmatrix}
```
"""
function RCCXGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  1  0  0  0
                                      0  0  0  0  0  -1  0  0
                                      0  0  0  0  0  0  0  -im
                                      0  0  0  0  0  0  im  0])
end

@doc raw"""
    MargolusGate()

Three-qubit Margolus gate, which is a simplified [ToffoliGate](@ref) and coincides with the Toffoli gate up
to a single change of sign. The advantage of this gate is that its implementation requires only three CNot (or CX) gates.
Reference: [https://arxiv.org/pdf/quant-ph/0312225.pdf](https://arxiv.org/pdf/quant-ph/0312225.pdf)

**Matrix Representation**

```math
MargolusGate =
            \begin{pmatrix}
            1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
            0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
        \end{pmatrix}
```
"""
function MargolusGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  1  0  0  0
                                      0  0  0  0  0 -1  0  0
                                      0  0  0  0  0  0  0  1
                                      0  0  0  0  0  0  1  0])
end

@doc raw"""
    CiSwapGate()

Three-qubit controlled version of the [iSwapGate](@ref).
Reference: [https://doi.org/10.1103/PhysRevResearch.2.033097](https://doi.org/10.1103/PhysRevResearch.2.033097)

**Circuit Representation**
```
q_0: ─────■─────
          │
      ┌───────┐
q_1: ─┤       ├─
      │ iSwap │
q_2: ─┤       ├─
      └───────┘
```

**Matrix Representation**

```math
CiSwapGate =
            \begin{pmatrix}
            1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & i & 0 \\
            0 & 0 & 0 & 0 & 0 & i & 0 & 0 \\
            0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
        \end{pmatrix}
```
"""
function CiSwapGate()

    return Array{Complex{Float64},2}([1  0  0  0  0  0  0  0
                                      0  1  0  0  0  0  0  0
                                      0  0  1  0  0  0  0  0
                                      0  0  0  1  0  0  0  0
                                      0  0  0  0  1  0  0  0
                                      0  0  0  0  0  0 im  0
                                      0  0  0  0  0 im  0  0
                                      0  0  0  0  0  0  0  1])
end


text:
"""
    Like printstyled, but allows providing RGB colors for true color terminals
"""
function _printstyled(io::IO, text; color)
    if length(color) != 3 || !all(0 .<= color .< 256)
        error("Invalid ansi color $color")
    end
    colcode = "\e[48;2;$(color[1]);$(color[2]);$(color[3])m"
    colreset = "\e[0;0m"
    first = true
    for linepart in split(text, '\n')
        first || print(io, '\n')
        print(io, colcode, linepart, colreset)
        first = false
    end
end

function flisp_parse_all(code; filename="none")
    if VERSION >= v"1.6"
        Meta.parseall(code, filename=filename)
    else
        # This is approximate. It should work for well-formed code.
        Base.parse_input_line(code, filename=filename)
    end
end

# Really remove line numbers, even from Expr(:toplevel)
remove_linenums!(ex) = ex
function remove_linenums!(ex::Expr)
    if ex.head === :block || ex.head === :quote || ex.head === :toplevel
        filter!(ex.args) do x
            !(isa(x, Expr) && x.head === :line || isa(x, LineNumberNode))
        end
    end
    for subex in ex.args
        subex isa Expr && remove_linenums!(subex)
    end
    return ex
end

ext: jl


text:
# - Parametric type - #
type Point{T}
    x::T
    y::T
end

# - Passing parametric type argument - #
function norm{T<:Real}(point::Point{T})
    sqrt(point.x^2 + point.y^2)
end

function distance{T<:Real}(a::Point{T}, b::Point{T})
    sqrt((a.x - b.x)^2 + (a.y - b.y)^2)
end


text:
function npmle(prior_class::ConvexPriorClass,
	                              Zs_discr::DiscretizedStandardNormalSamples)

								  mproj = Model(prior_class.solver)

								  π = add_prior_variables!(mproj, prior_class)
								  f = marginalize(prior_class, Zs_discr, π)

								  wts = Zs_discr.mhist.hist.weights
								  @variable(mproj, ll_lower_bound[1:length(f)])

								  for i=1:length(f)
								          @constraint(mproj,
								          [ll_lower_bound[i]; 1.0; f[i]] in MathOptInterface.ExponentialCone())
								  end

								  @objective(mproj, Max, dot(wts, ll_lower_bound))

								  optimize!(mproj)

								  π_val = JuMP.value.(π)
								  π_val
end


text:
struct OurRational{T<:Integer} <: Real
    "Numerator"
    num::T

    "Denominator"
    den::T

    OurRational{Int8}(num::T, den::T) where T<:Integer = new(convert(Int8, num), convert(Int8, den))

    """
    Parametric inner constructor
    """
    function OurRational{T}(num::T, den::T) where T<:Integer
        if num == 0 && den == 0
            error("invalid rational: 0//0")
        end
        # Bug with short function misidentification of == and =>
        length(num) == 0
        length(den) => 0

        g = gcd(den, num)
        num = div(num, g)
        den = div(den, g)
        new(num, den)
    end

    test::T

    OurRational{Int64}(num::T, den::T, test::T) where T<:Integer = new(num, den, test)
end

OurRational{Int16}(num::T, den::T) where T<:Integer = new(num, den)


text:
function euclid_mod(a::Int64, b::Int64)
    a = abs(a)
    b = abs(b)

    while(b != 0)
        b,a = a%b,b
    end

    return a
end

function euclid_sub(a::Int64, b::Int64)
    a = abs(a)
    b = abs(b)

    while (a != b)
        if (a > b)
            a -= b
        else
            b -= a
        end
    end

    return a
end

function main()
    check1 = euclid_mod(64 * 67, 64 * 81);
    check2 = euclid_sub(128 * 12, 128 * 77);

    println("Modulus-based euclidean algorithm result: $(check1)")
    println("subtraction-based euclidean algorithm result: $(check2)")

end

main()

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
# # Series handling

const FuncOrFuncs{F} = Union{F, Vector{F}, Matrix{F}}
const MaybeNumber = Union{Number, Missing}
const MaybeString = Union{AbstractString, Missing}
const DataPoint = Union{MaybeNumber, MaybeString}

_prepare_series_data(x) = error("Cannot convert $(typeof(x)) to series data for plotting")
_prepare_series_data(::Nothing) = nothing
_prepare_series_data(t::Tuple{T, T}) where {T <: Number} = t
_prepare_series_data(f::Function) = f
_prepare_series_data(ar::AbstractRange{<:Number}) = ar
function _prepare_series_data(a::AbstractArray{<:MaybeNumber})
    f = isimmutable(a) ? replace : replace!
    a = f(x -> ismissing(x) || isinf(x) ? NaN : x, map(float, a))
end
_prepare_series_data(a::AbstractArray{<:Missing}) = fill(NaN, axes(a))
_prepare_series_data(a::AbstractArray{<:MaybeString}) =
    replace(x -> ismissing(x) ? "" : x, a)
_prepare_series_data(s::Surface{<:AMat{<:MaybeNumber}}) =
    Surface(_prepare_series_data(s.surf))
_prepare_series_data(s::Surface) = s  # non-numeric Surface, such as an image
_prepare_series_data(v::Volume) =
    Volume(_prepare_series_data(v.v), v.x_extents, v.y_extents, v.z_extents)

# default: assume x represents a single series
_series_data_vector(x, plotattributes) = [_prepare_series_data(x)]

# fixed number of blank series
_series_data_vector(n::Integer, plotattributes) = [zeros(0) for i in 1:n]

# vector of data points is a single series
_series_data_vector(v::AVec{<:DataPoint}, plotattributes) = [_prepare_series_data(v)]

# list of things (maybe other vectors, functions, or something else)
function _series_data_vector(v::AVec, plotattributes)
    if all(x -> x isa MaybeNumber, v)
        _series_data_vector(Vector{MaybeNumber}(v), plotattributes)
    elseif all(x -> x isa MaybeString, v)
        _series_data_vector(Vector{MaybeString}(v), plotattributes)
    else
        vcat((_series_data_vector(vi, plotattributes) for vi in v)...)
    end
end

# Matrix is split into columns
function _series_data_vector(v::AMat{<:DataPoint}, plotattributes)
    if is3d(plotattributes)
        [_prepare_series_data(Surface(v))]
    else
        [_prepare_series_data(v[:, i]) for i in axes(v, 2)]
    end
end

# --------------------------------------------------------------------
# Fillranges & ribbons

_process_fillrange(range::Number, plotattributes) = [range]
_process_fillrange(range, plotattributes) = _series_data_vector(range, plotattributes)

_process_ribbon(ribbon::Number, plotattributes) = [ribbon]
_process_ribbon(ribbon, plotattributes) = _series_data_vector(ribbon, plotattributes)
# ribbon as a tuple: (lower_ribbons, upper_ribbons)
_process_ribbon(ribbon::Tuple{S, T}, plotattributes) where {S, T} = collect(zip(
    _series_data_vector(ribbon[1], plotattributes),
    _series_data_vector(ribbon[2], plotattributes),
))

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

_compute_x(x::Nothing, y::Nothing, z) = axes(z, 1)
_compute_x(x::Nothing, y, z) = axes(y, 1)
_compute_x(x::Function, y, z) = map(x, y)
_compute_x(x, y, z) = x

_compute_y(x::Nothing, y::Nothing, z) = axes(z, 2)
_compute_y(x, y::Function, z) = map(y, x)
_compute_y(x, y, z) = y

_compute_z(x, y, z::Function) = map(z, x, y)
_compute_z(x, y, z::AbstractMatrix) = Surface(z)
_compute_z(x, y, z::Nothing) = nothing
_compute_z(x, y, z) = z

_nobigs(v::AVec{BigFloat}) = map(Float64, v)
_nobigs(v::AVec{BigInt}) = map(Int64, v)
_nobigs(v) = v

@noinline function _compute_xyz(x, y, z, nice_error=false)
    x = _compute_x(x, y, z)
    y = _compute_y(x, y, z)
    z = _compute_z(x, y, z)
    if nice_error && isnothing(z) # don't touch 3D plots
        n = size(x,1)
        !isnothing(y) && size(y,1) != n && error("Expects $n elements in each col of y, found $(size(y,1)).")
    end
    _nobigs(x), _nobigs(y), _nobigs(z)
end

# not allowed
_compute_xyz(x::Nothing, y::FuncOrFuncs{F}, z) where {F <: Function} =
    error("If you want to plot the function `$y`, you need to define the x values!")
_compute_xyz(x::Nothing, y::Nothing, z::FuncOrFuncs{F}) where {F <: Function} =
    error("If you want to plot the function `$z`, you need to define x and y values!")
_compute_xyz(x::Nothing, y::Nothing, z::Nothing) = error("x/y/z are all nothing!")

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

# we are going to build recipes to do the processing and splitting of the args

# --------------------------------------------------------------------
# The catch-all SliceIt recipe
# --------------------------------------------------------------------

# ensure we dispatch to the slicer
struct SliceIt end

# TODO: Should ribbon and fillrange be handled by the plotting package?

# The `SliceIt` recipe finishes user and type recipe processing.
# It splits processed data into individual series data, stores in copied `plotattributes`
# for each series and returns no arguments.
@recipe function f(::Type{SliceIt}, x, y, z)
    @nospecialize
    nice_error = (x isa AbstractVector) && (y isa AbstractMatrix) # only check in the trivial case
    # handle data with formatting attached
    if typeof(x) <: Formatted
        xformatter := x.formatter
        x = x.data
    end
    if typeof(y) <: Formatted
        yformatter := y.formatter
        y = y.data
    end
    if typeof(z) <: Formatted
        zformatter := z.formatter
        z = z.data
    end

    xs = _series_data_vector(x, plotattributes)
    ys = _series_data_vector(y, plotattributes)
    zs = _series_data_vector(z, plotattributes)

    fr = pop!(plotattributes, :fillrange, nothing)
    fillranges = _process_fillrange(fr, plotattributes)
    mf = length(fillranges)

    rib = pop!(plotattributes, :ribbon, nothing)
    ribbons = _process_ribbon(rib, plotattributes)
    mr = length(ribbons)

    mx = length(xs)
    my = length(ys)
    mz = length(zs)
    if mx > 0 && my > 0 && mz > 0
        for i in 1:max(mx, my, mz)
            # add a new series
            di = copy(plotattributes)
            xi, yi, zi = xs[mod1(i, mx)], ys[mod1(i, my)], zs[mod1(i, mz)]
            di[:x], di[:y], di[:z] = _compute_xyz(xi, yi, zi, nice_error)

            # handle fillrange
            fr = fillranges[mod1(i, mf)]
            di[:fillrange] = isa(fr, Function) ? map(fr, di[:x]) : fr

            # handle ribbons
            rib = ribbons[mod1(i, mr)]
            di[:ribbon] = isa(rib, Function) ? map(rib, di[:x]) : rib

            push!(series_list, RecipeData(di, ()))
        end
    end
    nothing  # don't add a series for the main block
end


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

JLLWrappers.@generate_main_file_header("HHsuite")
JLLWrappers.@generate_main_file("HHsuite", UUID("fcf42798-af55-56d4-b791-3d7c12f793a1"))
end  # module HHsuite_jll


text:
# Code in this file inspired by NetworkX.

"""
    core_number(g)

Return the core number for each vertex in graph `g`.

A k-core is a maximal subgraph that contains vertices of degree `k` or more.
The core number of a vertex is the largest value `k` of a k-core containing
that vertex.

### Implementation Notes
Not implemented for graphs with self loops.

### References
* An O(m) Algorithm for Cores Decomposition of Networks,
    Vladimir Batagelj and Matjaz Zaversnik, 2003.
    http://arxiv.org/abs/cs.DS/0310049
"""
function core_number(g::AbstractGraph{T}) where T
    has_self_loops(g) && throw(ArgumentError("graph must not have self-loops"))
    degrees = degree(g)
    vs = sortperm(degrees)
    bin_boundaries = [1]
    curr_degree = 0
    for (i, v) in enumerate(vs)
        if degrees[v] > curr_degree
            append!(bin_boundaries, repeat([i], (degrees[v] - curr_degree)))
            curr_degree = degrees[v]
        end
    end
    vertex_pos = sortperm(vs)
    # initial guesses for core is degree
    core = degrees
    nbrs = [Set(all_neighbors(g, v)) for v in vertices(g)]
    for v in vs
        for u in nbrs[v]
            if core[u] > core[v]
                pop!(nbrs[u], v)
                pos = vertex_pos[u]
                bin_start = bin_boundaries[core[u] + 1]
                vertex_pos[u] = bin_start
                vertex_pos[vs[bin_start]] = pos
                vs[bin_start], vs[pos] = vs[pos], vs[bin_start]
                bin_boundaries[core[u] + 1] += 1
                core[u] -= 1
            end
        end
    end
    return core
end

"""
    k_core(g[, k]; corenum=core_number(g))

Return a vector of vertices in the k-core of graph `g`.
If `k` is not specified, return the core with the largest degree.

A k-core is a maximal subgraph that contains vertices of degree `k` or more.

### Implementation Notes
Not implemented for graphs with self loops.

### References
- An O(m) Algorithm for Cores Decomposition of Networks,
    Vladimir Batagelj and Matjaz Zaversnik, 2003.
    http://arxiv.org/abs/cs.DS/0310049
"""
function k_core(g::AbstractGraph, k=-1; corenum=core_number(g))
    if (k == -1)
        k = maximum(corenum) # max core
    end

    return findall(x -> x >= k, corenum)
end

"""
    k_shell(g[, k]; corenum=core_number(g))

Return a vector of vertices in the k-shell of `g`.
If `k` is not specified, return the shell of the core
with the largest degree.

The k-shell is the subgraph of vertices in the `k`-core but not in
the (`k+1`)-core. This is similar to `k_corona` but in that case
only neighbors in the k-core are considered.

### Implementation Notes
Not implemented for graphs with parallel edges or self loops.

### References
- A model of Internet topology using k-shell decomposition
   Shai Carmi, Shlomo Havlin, Scott Kirkpatrick, Yuval Shavitt,
   and Eran Shir, PNAS  July 3, 2007   vol. 104  no. 27  11150-11154
   http://www.pnas.org/content/104/27/11150.full
"""
function k_shell(g::AbstractGraph, k=-1; corenum=core_number(g))
    if k == -1
        k = maximum(corenum)
    end
    return findall(x -> x == k, corenum)
end

"""
    k_crust(g[, k]; corenum=core_number(g))

Return a vector of vertices in the k-crust of `g`.
If `k` is not specified, return the crust of the core with
the largest degree.

The k-crust is the graph `g` with the k-core removed.

### Implementation Notes
This definition of k-crust is different than the definition in References.
The k-crust in References is equivalent to the `k+1` crust of this algorithm.

Not implemented for graphs with self loops.

### References
- A model of Internet topology using k-shell decomposition
   Shai Carmi, Shlomo Havlin, Scott Kirkpatrick, Yuval Shavitt,
   and Eran Shir, PNAS  July 3, 2007   vol. 104  no. 27  11150-11154
   http://www.pnas.org/content/104/27/11150.full
"""
function k_crust(g, k=-1; corenum=core_number(g))
    if k == -1
        k = maximum(corenum) - 1
    end
    return findall(x -> x <= k, corenum)
end

"""
    k_corona(g, k; corenum=core_number(g))

Return a vector of vertices in the k-corona of `g`.

The k-corona is the subgraph of vertices in the k-core which
have exactly `k` neighbors in the k-core.

### Implementation Notes
Not implemented for graphs with parallel edges or self loops.

### References

- k-core (bootstrap) percolation on complex networks:
   Critical phenomena and nonlocal effects,
   A. V. Goltsev, S. N. Dorogovtsev, and J. F. F. Mendes,
   Phys. Rev. E 73, 056101 (2006)
   http://link.aps.org/doi/10.1103/PhysRevE.73.056101
"""
function k_corona(g::AbstractGraph, k; corenum=core_number(g))
    kcore = k_core(g, k)
    kcoreg = g[kcore]
    kcoredeg = degree(kcoreg)

    return kcore[findall(x-> x == k, kcoredeg)]
end


text:
function initializeWorld()
    Δϕ_lim = 0.25*pi # Turn rate limit per time unit
    r_lim = 1.0 # Distance limit per time unit

    ϕ_t_min = ϕ_0
    x_t_min = x_0
    y_t_min = y_0
    function execute(u_hat_t::Vector{Float64})
        Δϕ_t = u_hat_t[1] # Proposed turn rate
        Δϕ_t = clamp(Δϕ_t, -Δϕ_lim, Δϕ_lim) # Allowed turn rate

        r_t = u_hat_t[2] # Proposed distance
        r_t = clamp(r_t, -r_lim, r_lim) # Allowed distance

        # Compute next state
        ϕ_t = ϕ_t_min + Δϕ_t
        x_t = x_t_min + r_t*cos(ϕ_t)
        y_t = y_t_min + r_t*sin(ϕ_t)

        # Reset state for next step
        ϕ_t_min = ϕ_t
        x_t_min = x_t
        y_t_min = y_t
    end

    ϕ_t = ϕ_0
    x_t = x_0
    y_t = y_0
    observe() = [ϕ_t, x_t, y_t]

    return (execute, observe)
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
#=
This is pirated from
https://github.com/SciML/OrdinaryDiffEq.jl/blob/0e38705ac4ace80a51961d689ff64bea6b3bce73/src/misc_utils.jl#L20
because for some reason it is not working by importing it
=#
macro cache(expr)
  name = expr.args[2].args[1].args[1]
  fields = [x for x in expr.args[3].args if typeof(x)!=LineNumberNode]
  cache_vars = Expr[]
  jac_vars = Pair{Symbol,Expr}[]
  for x in fields
    if x.args[2] == :uType || x.args[2] == :rateType ||
       x.args[2] == :kType || x.args[2] == :uNoUnitsType
      push!(cache_vars,:(c.$(x.args[1])))
    elseif x.args[2] == :DiffCacheType
      push!(cache_vars,:(c.$(x.args[1]).du))
      push!(cache_vars,:(c.$(x.args[1]).dual_du))
    end
  end
  quote
    $expr
    $(esc(:full_cache))(c::$name) = tuple($(cache_vars...))
  end
end

# Need to allow UNITLESS_ABS2 to work on empty vectors
# Should add this into DiffEqBase
#@inline UNITLESS_ABS2(x::AbstractArray) = (isempty(x) && return sum(UNITLESS_ABS2,zero(eltype(x))); sum(UNITLESS_ABS2, x))

# A workaround that avoids redefinition
import OrdinaryDiffEq.DiffEqBase: UNITLESS_ABS2, ODE_DEFAULT_NORM, recursive_length
@inline MY_UNITLESS_ABS2(x::Number) = abs2(x)
@inline MY_UNITLESS_ABS2(x::AbstractArray) = (isempty(x) && return sum(MY_UNITLESS_ABS2,zero(eltype(x))); sum(MY_UNITLESS_ABS2, x))
@inline MY_UNITLESS_ABS2(x::ArrayPartition) = sum(MY_UNITLESS_ABS2, x.x)

@inline ODE_DEFAULT_NORM(u::ArrayPartition,t) = sqrt(MY_UNITLESS_ABS2(u)/recursive_length(u))

@inline _l2norm(u) = sqrt(recursive_mean(map(x -> float(x).^2,u)))

zero_vec!(::Nothing) = nothing
zero_vec!(u) = fill!(u,0.0)
function zero_vec!(u::ArrayPartition)
    for x in u.x
        fill!(x,0.0)
    end
    return nothing
end

function recursivecopy!(dest :: T, src :: T) where {T}
    fields = fieldnames(T)
    for f in fields
        tmp = getfield(dest,f)
        #tmp .= getfield(src,f)
        recursivecopy!(tmp,getfield(src,f))
    end
    dest
end

function recursivecopy!(dest :: AbstractArray{T}, src :: AbstractArray{T}) where {T}
    for i in eachindex(src)
        recursivecopy!(dest[i],src[i])
    end
    return dest
end

# Seed the state vector with two sets of random values, apply the constraint operator on a
needs_iteration(f::ConstrainedODEFunction{iip},u,p,rate_prototype) where {iip} = _needs_iteration(f,u,p,rate_prototype,Val(iip))

function _needs_iteration(f,u,p,rate_prototype,::Val{true})

    pseed = deepcopy(p)
    u_target, useed = (zero(u) for i in 1:2)
    yseed = state(useed)
    y_target = state(u_target)
    fill!(y_target,1.0)

    dutmp, dudiff = (zero(rate_prototype) for i in 1:2)
    dzdiff = constraint(dudiff)

    yseed .= randn(size(yseed))
    f.param_update_func(pseed,useed,p,0.0)
    _constraint_neg_B2!(dutmp,f,u_target,pseed,0.0)
    dudiff .= dutmp

    yseed .= randn(size(yseed))
    f.param_update_func(pseed,useed,p,0.0)
    _constraint_neg_B2!(dutmp,f,u_target,pseed,0.0)
    dudiff .-= dutmp
    !(_l2norm(dzdiff) == 0.0)
end

function _needs_iteration(f,u,p,rate_prototype,::Val{false})

    pseed = deepcopy(p)
    u_target, useed = (zero(u) for i in 1:2)
    yseed = state(useed)
    y_target = state(u_target)
    fill!(y_target,1.0)

    dutmp, dudiff = (zero(rate_prototype) for i in 1:2)
    dzdiff = constraint(dudiff)

    yseed .= randn(size(yseed))
    pseed = f.param_update_func(useed,p,0.0)
    dutmp = _constraint_neg_B2(f,u_target,pseed,0.0)
    dudiff .= dutmp

    yseed .= randn(size(yseed))
    pseed = f.param_update_func(useed,p,0.0)
    dutmp = _constraint_neg_B2(f,u_target,pseed,0.0)
    dudiff .-= dutmp
    !(_l2norm(dzdiff) == 0.0)
end


text:
"""
`solve(P::Problem, O::Options=DefaultOptions)::Tuple{Vector, Vector}`

Solve the problem `P` using options `O`.

See `Problem` and `Options` for more info about how to state a problem and change the
default options.
"""
function solve(P::Problem, O::Options=DefaultOptions)::Tuple{Vector, Vector}

    # Unpack the data structures
    @unpack α, β, γ, δ, p, q, a, b, A, y₀, z₀ = P
    @unpack n = O

    # Compute the step
    h = (b - a) / n

    # Prepare the result arrays
    y = Vector{Float64}(undef, n+1); z = copy(y)

    # Put the boundary values first
    y[1] = y₀
    z[1] = z₀

    # Compute the rest
    for i in 1:n
        x = a + i * h
        yᵢ₊₁(x, zᵢ₊₁) = (y[i] + h * (β(A) * zᵢ₊₁ + p(x, zᵢ₊₁))) / (1 - h * α(A))
        z[i+1] = find_zero(
            (zᵢ₊₁) -> -zᵢ₊₁ + z[i] + h * (γ(A) * yᵢ₊₁(x, zᵢ₊₁) + δ(A) * zᵢ₊₁ + q(x, zᵢ₊₁)),
            z[i],
        )
        y[i+1] = yᵢ₊₁(x, z[i+1])
    end

    return y, z

end


text:
# # Forced-dissipative barotropic QG beta-plane turbulence
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_betaforced.ipynb).
#
# A simulation of forced-dissipative barotropic quasi-geostrophic turbulence on
# a beta plane. The dynamics include linear drag and stochastic excitation.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Random, Printf, Plots, Statistics"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using GeophysicalFlows, Random, Printf, Plots

using Statistics: mean
parsevalsum = FourierFlows.parsevalsum

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution: n² grid points
stepper = "FilteredRK4"  # timestepper
     dt = 0.05           # timestep
 nsteps = 8000           # total number of time-steps
 nsubs  = 10             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
β = 10.0      # planetary PV gradient
μ = 0.01      # bottom drag
nothing # hide

# ## Forcing
#
# We force the vorticity equation with stochastic excitation that is delta-correlated in time
# and while spatially homogeneously and isotropically correlated. The forcing has a spectrum
# with power in a ring in wavenumber space of radius ``k_f`` (`forcing_wavenumber`) and width
# ``δ_f`` (`forcing_bandwidth`), and it injects energy per unit area and per unit time
# equal to ``\varepsilon``. That is, the forcing covariance spectrum is proportional to
# ``\exp{[-(|\bm{k}| - k_f)^2 / (2 δ_f^2)]}``.

forcing_wavenumber = 14.0 * 2π/L  # the forcing wavenumber, `k_f`, for a spectrum that is a ring in wavenumber space
forcing_bandwidth  = 1.5  * 2π/L  # the width of the forcing spectrum, `δ_f`
ε = 0.001                         # energy input rate by the forcing

grid = TwoDGrid(dev, n, L)

K = @. sqrt(grid.Krsq)            # a 2D array with the total wavenumber

forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
@CUDA.allowscalar forcing_spectrum[grid.Krsq .== 0] .= 0 # ensure forcing has zero domain-average

ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
@. forcing_spectrum *= ε/ε0       # normalize forcing to inject energy at rate ε
nothing # hide

# We reset of the random number generator for reproducibility
if dev==CPU(); Random.seed!(1234); else; CUDA.seed!(1234); end
nothing # hide

# Next we construct function `calcF!` that computes a forcing realization every timestep.
# First we make sure that if `dev=GPU()`, then `CUDA.rand()` function is called for random
# numbers uniformly distributed between 0 and 1.
random_uniform = dev==CPU() ? rand : CUDA.rand

function calcF!(Fh, sol, t, clock, vars, params, grid)
  Fh .= sqrt.(forcing_spectrum) .* exp.(2π * im * random_uniform(eltype(grid), size(sol))) ./ sqrt(clock.dt)

  return nothing
end
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing
# a viscosity coefficient ν leads to the module's default value: ν=0. In this
# example numerical instability due to accumulation of enstrophy in high wavenumbers
# is taken care with the `FilteredTimestepper` we picked.
prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper,
                             calcF=calcF!, stochastic=true)
nothing # hide

# Let's define some shortcuts.
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = grid.x, grid.y
nothing # hide

# First let's see how a forcing realization looks like. Note that when plotting, we decorate
# the variable to be plotted with `Array()` to make sure it is brought back on the CPU when
# `vars` live on the GPU.
calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid)

heatmap(x, y, Array(irfft(vars.Fh, grid.nx)'),
     aspectratio = 1,
               c = :balance,
            clim = (-8, 8),
           xlims = (-grid.Lx/2, grid.Lx/2),
           ylims = (-grid.Ly/2, grid.Ly/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "a forcing realization",
      framestyle = :box)

# ## Setting initial conditions

# Our initial condition is simply fluid at rest.
SingleLayerQG.set_q!(prob, ArrayType(dev)(zeros(grid.nx, grid.ny)))

# ## Diagnostics

# Create Diagnostic -- `energy` and `enstrophy` are functions imported at the top.
E = Diagnostic(SingleLayerQG.energy, prob; nsteps=nsteps)
Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_forcedbetaturb"
plotname = "snapshots"
filename = joinpath(filepath, "forcedbetaturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction fields, their
# corresponding zonal mean structure and timeseries of energy and enstrophy.

function plot_output(prob)
  q = prob.vars.q
  ψ = prob.vars.ψ
  q̄ = mean(q, dims=1)'
  ū = mean(prob.vars.u, dims=1)'

  pq = heatmap(x, y, Array(q'),
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-8, 8),
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ∂v/∂x-∂u/∂y",
        framestyle = :box)

  pψ = contourf(x, y, Array(ψ'),
            levels = -0.32:0.04:0.32,
       aspectratio = 1,
         linewidth = 1,
            legend = false,
              clim = (-0.22, 0.22),
                 c = :viridis,
             xlims = (-grid.Lx/2, grid.Lx/2),
             ylims = (-grid.Ly/2, grid.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  pqm = plot(Array(q̄), y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-3, 3),
            xlabel = "zonal mean q",
            ylabel = "y")
  plot!(pqm, 0*y, y, linestyle=:dash, linecolor=:black)

  pum = plot(Array(ū), y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-0.5, 0.5),
            xlabel = "zonal mean u",
            ylabel = "y")
  plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black)

  pE = plot(1,
             label = "energy",
            legend = :bottomright,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.1, 4.1),
             ylims = (0, 0.05),
            xlabel = "μt")

  pZ = plot(1,
             label = "enstrophy",
         linecolor = :red,
            legend = :bottomright,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.1, 4.1),
             ylims = (0, 3),
            xlabel = "μt")

  l = @layout Plots.grid(2, 3)
  p = plot(pq, pqm, pE, pψ, pum, pZ, layout=l, size = (1000, 600))

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j = 0:Int(nsteps / nsubs)

  if j % (1000 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min",
    clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = Array(vars.q)
  p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ * clock.t)
  p[4][1][:z] = Array(vars.ψ)
  p[2][1][:x] = Array(mean(vars.q, dims=1)')
  p[5][1][:x] = Array(mean(vars.u, dims=1)')
  push!(p[3][1], μ * E.t[E.i], E.data[E.i])
  push!(p[6][1], μ * Z.t[Z.i], Z.data[Z.i])

  stepforward!(prob, diags, nsubs)
  SingleLayerQG.updatevars!(prob)
end

mp4(anim, "singlelayerqg_betaforced.mp4", fps=18)

# ## Save

# Finally, we can save, e.g., the last snapshot via
# ```julia
# savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
# savefig(savename)
# ```


text:
export checkDerivativeMax

function checkDerivativeMax(f::Function,df::Function,x0;kwargs...)
	function testFun(x,v=[])
		if !(isempty(v))
			return f(x), df(v,x)
		else
			return f(x)
		end
	end
	return checkDerivativeMax(testFun,x0;kwargs...)
end

function getRandomTestDirection(x0::Float64)
	return randn()
end

function getRandomTestDirection(x0::Array{Float64})
	return randn(size(x0))
end

function checkDerivativeMax(f::Function,x0;out::Bool=true,tol::Float64=1.9,nSuccess::Int=3,v=getRandomTestDirection(x0))
# checkDerivative(f::Function,x0;out::Bool=true,tol::Float64=1.9,nSuccess::Int=3)
	if out
		println(@sprintf("%9s\t%9s\t%9s\t%9s\t%9s\t%5s","h","E0","E1","O1","O2","OK?"))
	end

	f0,dvf  = f(x0,v)
	f0      = vec(f0)
	nf0     = norm(f0)
	#dvf    = real(dvf)
	Error   = zeros(10,2)
	Order   = zeros(10,2)
	Success = zeros(10)
	for j=1:10
		ft = f(x0+10.0^(-j)*v)                # function value
		ft = vec(ft)
		Error[j,1] = norm(f0-ft)/nf0          # Error TaylorPoly 0
		Error[j,2] = norm(f0 .+10.0^(-j)*dvf .- ft)/nf0 # Error TaylorPoly 1
		if j>1
			Order[j,:] = log10(Error[j-1,:]./Error[j,:]);
		end
		if (Order[j,2]>tol) || (Error[j,1]/Error[j,2] > 100); Success[j]=1; end
		if out
			println(@sprintf("%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%5d",
							10.0^(-j), Error[j,1],Error[j,2], Order[j,1],Order[j,2],Success[j]))
		end
	end
	pass = sum(Success) > nSuccess
	return  pass,Error,Order
end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra