text:
/*
 *  Authors: Lana Mineh and John Scott
 *  Copyright 2021 Phasecraft Ltd. and John Scott
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 */

/**
 * \file gates.cpp
 * \brief Contains tests for the gate member functions
 */

#include <gtest/gtest.h>
#include <qsl/qubits.hpp>
#include <qsl/concepts.hpp>
#include <armadillo>
#include <complex>
#include <list>
#include "test-utils.hpp"

template<typename T>
struct SimWrapper
{
    using Sim = T;
};

/**
 * \brief Typed test suite for one-qubit gates
 */
template <typename T>
class Gates : public testing::Test {
public:
    using List = std::list<T>;
    static T shared_;
    T value_;
};

template <typename T>
class NPGates : public testing::Test {
public:
    using List = std::list<T>;
    static T shared_;
    T value_;
};

/// List of simulator types to check
using Sim1 = SimWrapper<qsl::Qubits<qsl::Type::Default, float>>;
using Sim2 = SimWrapper<qsl::Qubits<qsl::Type::Default, double>>;
using Sim3 = SimWrapper<qsl::Qubits<qsl::Type::Omp, float>>;
using Sim4 = SimWrapper<qsl::Qubits<qsl::Type::Omp, double>>;
using Sim5 = SimWrapper<qsl::Qubits<qsl::Type::Resize, float>>;
using Sim6 = SimWrapper<qsl::Qubits<qsl::Type::Resize, double>>;
using SimTypes = ::testing::Types<Sim1, Sim2, Sim3, Sim4, Sim5, Sim6>;

using Sim7 = SimWrapper<qsl::Qubits<qsl::Type::NP, float>>;
using Sim8 = SimWrapper<qsl::Qubits<qsl::Type::NP, double>>;
using Sim9 = SimWrapper<qsl::Qubits<qsl::Type::OmpNP, float>>;
using Sim10 = SimWrapper<qsl::Qubits<qsl::Type::OmpNP, double>>;

using NPSimTypes = ::testing::Types<Sim7,Sim8,Sim9,Sim10>;

TYPED_TEST_SUITE(Gates, SimTypes);
TYPED_TEST_SUITE(NPGates, NPSimTypes);

TEST(GateTests, MakeMatrixTestsPauliX)
{
    // Pauli gate, targ = 4
    const unsigned num_qubits{ 7 };
    const unsigned targ{ 4 };
    arma::Mat<std::complex<double>> mat(2, 2, arma::fill::zeros);
    mat(0b0,0b1) = 1;
    mat(0b1,0b0) = 1;

    // Make identity gate
    // Create gate in armadillo
    // Sizes of idenity matrix padding
    std::size_t pre = 1 << targ;
    std::size_t post = 1 << (num_qubits - targ - 1);

    // Tensor to make the gate. Have fun with these lines...
    arma::SpMat<std::complex<double>> gate =
	arma::speye<arma::SpMat<std::complex<double>>>(pre, pre);
    gate = arma::kron(arma::conv_to<arma::SpMat<std::complex<double>>>::from(mat),
		      gate);
    gate = arma::kron(arma::speye<arma::SpMat<std::complex<double>>>(post, post),
		      gate);

    // Compute the same matrix using the makeMatrix function
    const auto gate2{ makeMatrix(mat, num_qubits, {targ}) };
    EXPECT_TRUE(arma::approx_equal(gate2, gate, "absdiff", 1e-8));

}

TEST(GateTests,MakeMatrixTestsCnot)
{

    // CNOT gate, qubit 1 is ctrl, qubit 0 is targ
    arma::Mat<std::complex<double>> gate(4, 4, arma::fill::zeros);
    gate(0b00,0b00) = 1;
    gate(0b01,0b01) = 1;
    gate(0b10,0b11) = 1;
    gate(0b11,0b10) = 1;

    // Apply CNOT to qubits ctrl = 1, targ = 0
    unsigned ctrl = 1;
    unsigned targ = 0;
    const auto m0 = makeMatrix(gate, 3, {targ, ctrl});
    //std::cout << m0 << std::endl;
    // Check all the basis states
    const std::complex<double> one{1,0};
    EXPECT_EQ(m0(0b000,0b000), one); // |000) -> |000)
    EXPECT_EQ(m0(0b001,0b001), one); // |001) -> |001)
    EXPECT_EQ(m0(0b010,0b011), one); // |010) -> |011)
    EXPECT_EQ(m0(0b011,0b010), one); // |011) -> |010)
    EXPECT_EQ(m0(0b100,0b100), one); // |100) -> |100)
    EXPECT_EQ(m0(0b101,0b101), one); // |101) -> |101)
    EXPECT_EQ(m0(0b110,0b111), one); // |110) -> |111)
    EXPECT_EQ(m0(0b111,0b110), one); // |111) -> |110)

    // Apply CNOT to qubits ctrl = 0, targ = 2
    ctrl = 0;
    targ = 2;
    const auto m1 = makeMatrix(gate, 3, {targ, ctrl});
    //std::cout << m1 << std::endl;
    // Check all the basis states
    EXPECT_EQ(m1(0b000,0b000), one);
    EXPECT_EQ(m1(0b001,0b101), one);
    EXPECT_EQ(m1(0b010,0b010), one);
    EXPECT_EQ(m1(0b011,0b111), one);
    EXPECT_EQ(m1(0b100,0b100), one);
    EXPECT_EQ(m1(0b101,0b001), one);
    EXPECT_EQ(m1(0b110,0b110), one);
    EXPECT_EQ(m1(0b111,0b011), one);

    // Check the degnerate cases
    ctrl = 0;
    targ = 1;
    const auto m2 = makeMatrix(gate, 2, {targ, ctrl});
    //std::cout << m2 << std::endl;
    // Check all the basis states
    EXPECT_EQ(m2(0b00,0b00), one);
    EXPECT_EQ(m2(0b01,0b11), one);
    EXPECT_EQ(m2(0b10,0b10), one);
    EXPECT_EQ(m2(0b11,0b01), one);

    // Check the degnerate case the other way round
    ctrl = 1;
    targ = 0;
    const auto m3 = makeMatrix(gate, 2, {targ, ctrl});
    //std::cout << m3 << std::endl;
    // Check all the basis states
    EXPECT_EQ(m3(0b00,0b00), one);
    EXPECT_EQ(m3(0b01,0b01), one);
    EXPECT_EQ(m3(0b10,0b11), one);
    EXPECT_EQ(m3(0b11,0b10), one);
}

TEST(GateTests,MakeMatrixTestsSwap)
{

    // SWAP gate, qubit 1 is ctrl, qubit 0 is targ
    arma::Mat<std::complex<double>> gate(4, 4, arma::fill::zeros);
    gate(0b00,0b00) = 1;
    gate(0b01,0b10) = 1;
    gate(0b10,0b01) = 1;
    gate(0b11,0b11) = 1;

    // Apply SWAP to qubits ctrl = 1, targ = 0
    unsigned ctrl = 1;
    unsigned targ = 0;
    const auto m0 = makeMatrix(gate, 3, {targ, ctrl});
    //std::cout << m0 << std::endl;
    // Check all the basis states
    const std::complex<double> one{1,0};
    EXPECT_EQ(m0(0b000,0b000), one); // |000) -> |000)
    EXPECT_EQ(m0(0b001,0b010), one); // |001) -> |010)
    EXPECT_EQ(m0(0b010,0b001), one); // |010) -> |001)
    EXPECT_EQ(m0(0b011,0b011), one); // |011) -> |011)
    EXPECT_EQ(m0(0b100,0b100), one); // |100) -> |100)
    EXPECT_EQ(m0(0b101,0b110), one); // |101) -> |110)
    EXPECT_EQ(m0(0b110,0b101), one); // |110) -> |101)
    EXPECT_EQ(m0(0b111,0b111), one); // |111) -> |111)

    // Apply SWAP to qubits ctrl = 0, targ = 2
    ctrl = 0;
    targ = 2;
    const auto m1 = makeMatrix(gate, 3, {targ, ctrl});
    //std::cout << m1 << std::endl;
    // Check all the basis states
    EXPECT_EQ(m1(0b000,0b000), one);
    EXPECT_EQ(m1(0b001,0b100), one);
    EXPECT_EQ(m1(0b010,0b010), one);
    EXPECT_EQ(m1(0b011,0b110), one);
    EXPECT_EQ(m1(0b100,0b001), one);
    EXPECT_EQ(m1(0b101,0b101), one);
    EXPECT_EQ(m1(0b110,0b011), one);
    EXPECT_EQ(m1(0b111,0b111), one);

}

TEST(GateTests,MakeMatrixTestsCHadamard)
{
    // Controlled Hadamard gate, qubit 1 is ctrl, qubit 0 is targ
    arma::Mat<std::complex<double>> gate(4, 4, arma::fill::zeros);
    double one_sqrt2{1/std::sqrt(2)};
    gate(0b00,0b00) = 1;
    gate(0b01,0b01) = 1;
    gate(0b10,0b10) = one_sqrt2;
    gate(0b10,0b11) = one_sqrt2;
    gate(0b11,0b10) = one_sqrt2;
    gate(0b11,0b11) = -one_sqrt2;

    // Apply controlled Hadamard to qubits ctrl = 1, targ = 0
    unsigned ctrl = 1;
    unsigned targ = 0;
    const auto m0 = makeMatrix(gate, 3, {targ, ctrl});
    //std::cout << m0 << std::endl;
    // Check all the basis states
    const std::complex<double> one{1,0};
    const std::complex<double> val{one_sqrt2,0};
    EXPECT_EQ(m0(0b000,0b000), one); // |000) -> |000)
    EXPECT_EQ(m0(0b001,0b001), one); // |001) -> |010)
    EXPECT_EQ(m0(0b010,0b010), val); // |010) -> [ |010) + |011) ] / sqrt2
    EXPECT_EQ(m0(0b010,0b011), val);
    EXPECT_EQ(m0(0b011,0b010), val); // |011) -> [ |010) - |011) ] / sqrt2
    EXPECT_EQ(m0(0b011,0b011), -val);
    EXPECT_EQ(m0(0b100,0b100), one); // |100) -> |100)
    EXPECT_EQ(m0(0b101,0b101), one); // |101) -> |110)
    EXPECT_EQ(m0(0b110,0b110), val); // |110) -> [ |110) + |111) ] / sqrt2
    EXPECT_EQ(m0(0b110,0b111), val);
    EXPECT_EQ(m0(0b111,0b110), val); // |111) -> [ |110) - |111) ] / sqrt2
    EXPECT_EQ(m0(0b111,0b111), -val);

    // Apply controlled Hadamard to qubits ctrl = 1, targ = 0
    ctrl = 0;
    targ = 2;
    const auto m1 = makeMatrix(gate, 3, {targ, ctrl});
    //std::cout << m1 << std::endl;

    // Check all the basis states
    // Identity
    EXPECT_EQ(m1(0b000,0b000), one);
    EXPECT_EQ(m1(0b010,0b010), one);
    EXPECT_EQ(m1(0b100,0b100), one);
    EXPECT_EQ(m1(0b110,0b110), one);

    // To the plus state
    EXPECT_EQ(m1(0b001,0b001), val);
    EXPECT_EQ(m1(0b001,0b101), val);
    EXPECT_EQ(m1(0b011,0b011), val);
    EXPECT_EQ(m1(0b011,0b111), val);

    // To the minus state
    EXPECT_EQ(m1(0b101,0b001), val);
    EXPECT_EQ(m1(0b101,0b101), -val);
    EXPECT_EQ(m1(0b111,0b011), val);
    EXPECT_EQ(m1(0b111,0b111), -val);

}

/// Test the makeMatrix function against tensor product matrix
TEST(GateTests,MakeMatrixTestsTensor)
{
    const unsigned num_qubits{ 8 };
    unsigned ctrl = 6;
    unsigned targ = 5;

    // Controlled Hadamard gate, qubit 1 is ctrl, qubit 0 is targ
    arma::Mat<std::complex<double>> mat(4, 4, arma::fill::zeros);
    double one_sqrt2{1/std::sqrt(2)};
    mat(0b00,0b00) = 1;
    mat(0b01,0b01) = 1;
    mat(0b10,0b10) = one_sqrt2;
    mat(0b10,0b11) = one_sqrt2;
    mat(0b11,0b10) = one_sqrt2;
    mat(0b11,0b11) = -one_sqrt2;

    // Make identity gate
    // Create gate in armadillo
    // Sizes of idenity matrix padding
    std::size_t pre = 1 << targ;
    std::size_t post = 1 << (num_qubits - ctrl - 1);

    // Tensor to make the gate. Have fun with these lines...
    arma::SpMat<std::complex<double>> gate =
	arma::speye<arma::SpMat<std::complex<double>>>(pre, pre);
    gate = arma::kron(arma::conv_to<arma::SpMat<std::complex<double>>>::from(mat),
		      gate);
    gate = arma::kron(arma::speye<arma::SpMat<std::complex<double>>>(post, post),
		      gate);

    // Compute the same matrix using the makeMatrix function
    const auto gate2{ makeMatrix(mat, num_qubits, {targ, ctrl}) };
    EXPECT_TRUE(arma::approx_equal(gate2, gate, "absdiff", 1e-8));
}

TYPED_TEST(Gates, OneQubitGate)
{
    const unsigned num_qubits = 8;
    const unsigned targ = 4;
    using Sim = TypeParam::Sim;
    using Fp = TypeParam::Sim::Fp_type;
    // Generate random angle to test parameterised gates
    Fp angle = arma::datum::pi * arma::randn<Fp>();

    // Create list of gates mapped to matrices
    std::vector<std::pair<std::function<void(Sim &, unsigned)>,
			  arma::Mat<std::complex<Fp>>>> gates;

    // PauliX
    arma::Mat<std::complex<Fp>> pauliX(2, 2, arma::fill::zeros);
    pauliX(0, 1) = 1;
    pauliX(1, 0) = 1;
    auto fn_pauliX = [](Sim & sim, unsigned targ) {
			 sim.pauliX(targ);
		     };
    gates.push_back({fn_pauliX, pauliX});

    // PauliY
    arma::Mat<std::complex<Fp>> pauliY(2, 2, arma::fill::zeros);
    pauliY(0, 1) = std::complex<Fp>{0, -1};
    pauliY(1, 0) = std::complex<Fp>{0, 1};
    auto fn_pauliY = [](Sim & sim, unsigned targ) {
    			 sim.pauliY(targ);
    		     };
    gates.push_back({fn_pauliY, pauliY});

    // PauliZ
    arma::Mat<std::complex<Fp>> pauliZ(2, 2, arma::fill::zeros);
    pauliZ(0, 0) = 1;
    pauliZ(1, 1) = -1;
    auto fn_pauliZ = [](Sim & sim, unsigned targ) {
			 sim.pauliZ(targ);
		     };
    gates.push_back({fn_pauliZ, pauliZ});

    // Hadamard
    arma::Mat<std::complex<Fp>> hadamard(2, 2, arma::fill::zeros);
    Fp sqrt2 = 1/std::sqrt(2);
    hadamard(0, 0) = sqrt2;
    hadamard(0, 1) = sqrt2;
    hadamard(1, 0) = sqrt2;
    hadamard(1, 1) = -sqrt2;
    auto fn_hadamard = [](Sim & sim, unsigned targ) {
			   sim.hadamard(targ);
		       };
    gates.push_back({fn_hadamard, hadamard});

    // phase shift
    arma::Mat<std::complex<Fp>> phase(2, 2, arma::fill::zeros);
    phase(0, 0) = 1;
    phase(1, 1) = std::complex<Fp>{std::cos(angle), std::sin(angle)};
    auto fn_phase = [=](Sim & sim, unsigned targ) {
			sim.phase(targ, angle);
		       };
    gates.push_back({fn_phase, phase});

    // rotateX
    arma::Mat<std::complex<Fp>> rotateX(2, 2, arma::fill::zeros);
    rotateX(0, 0) = std::cos(angle/2);
    rotateX(0, 1) = std::complex<Fp>{0,-std::sin(angle/2)};
    rotateX(1, 0) = std::complex<Fp>{0,-std::sin(angle/2)};
    rotateX(1, 1) = std::cos(angle/2);
    auto fn_rotateX = [=](Sim & sim, unsigned targ) {
			sim.rotateX(targ, angle);
		       };
    gates.push_back({fn_rotateX, rotateX});

    // rotateY
    arma::Mat<std::complex<Fp>> rotateY(2, 2, arma::fill::zeros);
    rotateY(0, 0) = std::cos(angle/2);
    rotateY(0, 1) = -std::sin(angle/2);
    rotateY(1, 0) = std::sin(angle/2);
    rotateY(1, 1) = std::cos(angle/2);
    auto fn_rotateY = [=](Sim & sim, unsigned targ) {
			sim.rotateY(targ, angle);
		       };
    gates.push_back({fn_rotateY, rotateY});

    // rotateZ
    arma::Mat<std::complex<Fp>> rotateZ(2, 2, arma::fill::zeros);
    rotateZ(0, 0) = std::complex<Fp>{std::cos(angle/2), std::sin(-angle/2)};
    rotateZ(1, 1) = std::complex<Fp>{std::cos(angle/2), std::sin(angle/2)};
    auto fn_rotateZ = [=](Sim & sim, unsigned targ) {
			sim.rotateZ(targ, angle);
		       };
    gates.push_back({fn_rotateZ, rotateZ});

    for (const auto & [fn, mat] : gates) {

	// Make a random state
	Sim q{num_qubits};
	const std::vector<qsl::complex<Fp>> state
	    = qsl::makeRandomState<Fp>(num_qubits);
	q.setState(state);

	// Set an armadillo vector to the same state
	std::size_t dim = 1 << num_qubits;
	arma::Col<std::complex<Fp>> v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    v(i) = std::complex<Fp>{state[i].real, state[i].imag};
	}

	// Apply gate to qubits
	std::invoke(fn, q, targ);

	// Apply gate in armadillo
	v = makeMatrix(mat, num_qubits, {targ}) * v;

	// Read qubit state into armadillo
	std::vector<qsl::complex<Fp>> res = q.getState();
	arma::Col<std::complex<Fp>> qubit_v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    qubit_v(i) = std::complex<Fp>{res[i].real, res[i].imag};
	}

    	EXPECT_TRUE(arma::approx_equal(v, qubit_v, "both", 1e-6, 1e-8));
    }
}

TYPED_TEST(Gates, TwoQubitGate)
{
    const unsigned num_qubits = 8;
    const unsigned ctrl = 2;
    const unsigned targ = 4;
    using Sim = TypeParam::Sim;
    using Fp = TypeParam::Sim::Fp_type;
    // Generate random angle to test parameterised gates
    Fp angle = arma::datum::pi * arma::randn<Fp>();

    // Create list of gates mapped to matrices
    std::vector<std::pair<std::function<void(Sim &, unsigned, unsigned)>,
			  arma::Mat<std::complex<Fp>>>> gates;

    // controlNot
    arma::Mat<std::complex<Fp>> cnot(4, 4, arma::fill::zeros);
    cnot(0b00,0b00) = 1;
    cnot(0b01,0b01) = 1;
    cnot(0b10,0b11) = 1;
    cnot(0b11,0b10) = 1;

    auto fn_cnot = [](Sim & sim, unsigned ctrl, unsigned targ) {
			 sim.controlNot(ctrl, targ);
		     };
    gates.push_back({fn_cnot, cnot});

    // controlY
    arma::Mat<std::complex<Fp>> cy(4, 4, arma::fill::zeros);
    cy(0b00,0b00) = 1;
    cy(0b01,0b01) = 1;
    cy(0b10,0b11) = std::complex<Fp>{0, -1};
    cy(0b11,0b10) = std::complex<Fp>{0, 1};;
    auto fn_cy = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		     sim.controlY(ctrl, targ);
		   };
    gates.push_back({fn_cy, cy});

    // controlZ
    arma::Mat<std::complex<Fp>> cz(4, 4, arma::fill::zeros);
    cz(0b00,0b00) = 1;
    cz(0b01,0b01) = 1;
    cz(0b10,0b10) = 1;
    cz(0b11,0b11) = -1;
    auto fn_cz = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		     sim.controlZ(ctrl, targ);
		   };
    gates.push_back({fn_cz, cz});

    // controlRotateX
    arma::Mat<std::complex<Fp>> crx(4, 4, arma::fill::zeros);
    crx(0b00,0b00) = 1;
    crx(0b01,0b01) = 1;
    crx(0b10,0b10) = std::complex<Fp>{std::cos(angle/2), 0};
    crx(0b10,0b11) = std::complex<Fp>{0, -std::sin(angle/2)};
    crx(0b11,0b10) = std::complex<Fp>{0, -std::sin(angle/2)};
    crx(0b11,0b11) = std::complex<Fp>{std::cos(angle/2), 0};

    auto fn_crx = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.controlRotateX(ctrl, targ, angle);
		   };
    gates.push_back({fn_crx, crx});

    // controlRotateY
    arma::Mat<std::complex<Fp>> cry(4, 4, arma::fill::zeros);
    cry(0b00,0b00) = 1;
    cry(0b01,0b01) = 1;
    cry(0b10,0b10) = std::cos(angle/2);
    cry(0b10,0b11) = -std::sin(angle/2);
    cry(0b11,0b10) = std::sin(angle/2);
    cry(0b11,0b11) = std::cos(angle/2);

    auto fn_cry = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.controlRotateY(ctrl, targ, angle);
		   };
    gates.push_back({fn_cry, cry});

    // controlRotateZ
    arma::Mat<std::complex<Fp>> crz(4, 4, arma::fill::zeros);
    crz(0b00,0b00) = 1;
    crz(0b01,0b01) = 1;
    crz(0b10,0b10) = std::complex<Fp>{std::cos(angle/2), -std::sin(angle/2)};
    crz(0b11,0b11) = std::complex<Fp>{std::cos(angle/2), std::sin(angle/2)};
    auto fn_crz = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.controlRotateZ(ctrl, targ, angle);
		   };
    gates.push_back({fn_crz, crz});

    // controlPhase
    arma::Mat<std::complex<Fp>> cphase(4, 4, arma::fill::zeros);
    cphase(0b00,0b00) = 1;
    cphase(0b01,0b01) = 1;
    cphase(0b10,0b10) = 1;
    cphase(0b11,0b11) = std::complex<Fp>{std::cos(angle), std::sin(angle)};
    auto fn_cphase = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		       sim.controlPhase(ctrl, targ, angle);
		   };
    gates.push_back({fn_cphase, cphase});

    // controlHadamard
    Fp sqrt2 = 1/std::sqrt(2);
    arma::Mat<std::complex<Fp>> ch(4, 4, arma::fill::zeros);
    ch(0b00,0b00) = 1;
    ch(0b01,0b01) = 1;
    ch(0b10,0b10) = sqrt2;
    ch(0b10,0b11) = sqrt2;
    ch(0b11,0b10) = sqrt2;
    ch(0b11,0b11) = -sqrt2;

    auto fn_ch = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.controlHadamard(ctrl, targ);
		   };
    gates.push_back({fn_ch, ch});

    // swap
    arma::Mat<std::complex<Fp>> swap(4, 4, arma::fill::zeros);
    swap(0b00,0b00) = 1;
    swap(0b01,0b10) = 1;
    swap(0b10,0b01) = 1;
    swap(0b11,0b11) = 1;
    auto fn_swap = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		       sim.swap(ctrl, targ);
		   };
    gates.push_back({fn_swap, swap});

    // fswap
    arma::Mat<std::complex<Fp>> fswap(4, 4, arma::fill::zeros);
    fswap(0b00,0b00) = 1;
    fswap(0b01,0b10) = 1;
    fswap(0b10,0b01) = 1;
    fswap(0b11,0b11) = -1;
    auto fn_fswap = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.fswap(ctrl, targ);
		   };
    gates.push_back({fn_fswap, fswap});

    // npRotateX
    arma::Mat<std::complex<Fp>> nprx(4, 4, arma::fill::zeros);
    nprx(0b00,0b00) = 1;
    nprx(0b01,0b01) = std::complex<Fp>{std::cos(angle/2), 0};
    nprx(0b01,0b10) = std::complex<Fp>{0, -std::sin(angle/2)};
    nprx(0b10,0b01) = std::complex<Fp>{0, -std::sin(angle/2)};
    nprx(0b10,0b10) = std::complex<Fp>{std::cos(angle/2), 0};
    nprx(0b11,0b11) = 1;

    auto fn_nprx = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.npRotateX(targ, ctrl, angle);
		   };
    gates.push_back({fn_nprx, nprx});

    // npRotateY
    arma::Mat<std::complex<Fp>> npry(4, 4, arma::fill::zeros);
    npry(0b00,0b00) = 1;
    npry(0b01,0b01) = std::cos(angle/2);
    npry(0b01,0b10) = -std::sin(angle/2);
    npry(0b10,0b01) = std::sin(angle/2);
    npry(0b10,0b10) = std::cos(angle/2);
    npry(0b11,0b11) = 1;

    auto fn_npry = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.npRotateY(targ, ctrl, angle);
		   };
    gates.push_back({fn_npry, npry});

    // npHadamard
    arma::Mat<std::complex<Fp>> nph(4, 4, arma::fill::zeros);
    nph(0b00,0b00) = 1;
    nph(0b01,0b01) = sqrt2;
    nph(0b01,0b10) = sqrt2;
    nph(0b10,0b01) = sqrt2;
    nph(0b10,0b10) = -sqrt2;
    nph(0b11,0b11) = 1;

    auto fn_nph = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.npHadamard(targ, ctrl);
		   };
    gates.push_back({fn_nph, nph});

    for (const auto & [fn, mat] : gates) {

	// Make a random state
	Sim q{num_qubits};
	const std::vector<qsl::complex<Fp>> state
	    = qsl::makeRandomState<Fp>(num_qubits);
	q.setState(state);

	// Set an armadillo vector to the same state
	std::size_t dim = 1 << num_qubits;
	arma::Col<std::complex<Fp>> v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    v(i) = std::complex<Fp>{state[i].real, state[i].imag};
	}

	// Apply gate to qubits
	std::invoke(fn, q, ctrl, targ);

	// Apply gate in armadillo
	v = makeMatrix(mat, num_qubits, {targ, ctrl}) * v;

	// Read qubit state into armadillo
	std::vector<qsl::complex<Fp>> res = q.getState();
	arma::Col<std::complex<Fp>> qubit_v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    qubit_v(i) = std::complex<Fp>{res[i].real, res[i].imag};
	}

    	EXPECT_TRUE(arma::approx_equal(v, qubit_v, "both", 1e-6, 1e-8));
    }
}

TYPED_TEST(NPGates, TwoQubitGate)
{
    const unsigned num_qubits = 8;
    const unsigned num_ones = 4;
    const unsigned ctrl = 2;
    const unsigned targ = 4;
    using Sim = TypeParam::Sim;
    using Fp = TypeParam::Sim::Fp_type;
    // Generate random angle to test parameterised gates
    Fp angle = arma::datum::pi * arma::randn<Fp>();

    // Create list of gates mapped to matrices
    std::vector<std::pair<std::function<void(Sim &, unsigned, unsigned)>,
			  arma::Mat<std::complex<Fp>>>> gates;

    // controlZ
    arma::Mat<std::complex<Fp>> cz(4, 4, arma::fill::zeros);
    cz(0b00,0b00) = 1;
    cz(0b01,0b01) = 1;
    cz(0b10,0b10) = 1;
    cz(0b11,0b11) = -1;
    auto fn_cz = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		     sim.controlZ(ctrl, targ);
		   };
    gates.push_back({fn_cz, cz});

    // controlRotateZ
    arma::Mat<std::complex<Fp>> crz(4, 4, arma::fill::zeros);
    crz(0b00,0b00) = 1;
    crz(0b01,0b01) = 1;
    crz(0b10,0b10) = std::complex<Fp>{std::cos(angle/2), -std::sin(angle/2)};
    crz(0b11,0b11) = std::complex<Fp>{std::cos(angle/2), std::sin(angle/2)};
    auto fn_crz = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.controlRotateZ(ctrl, targ, angle);
		   };
    gates.push_back({fn_crz, crz});

    // controlPhase
    arma::Mat<std::complex<Fp>> cphase(4, 4, arma::fill::zeros);
    cphase(0b00,0b00) = 1;
    cphase(0b01,0b01) = 1;
    cphase(0b10,0b10) = 1;
    cphase(0b11,0b11) = std::complex<Fp>{std::cos(angle), std::sin(angle)};
    auto fn_cphase = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		       sim.controlPhase(ctrl, targ, angle);
		   };
    gates.push_back({fn_cphase, cphase});

    // swap
    arma::Mat<std::complex<Fp>> swap(4, 4, arma::fill::zeros);
    swap(0b00,0b00) = 1;
    swap(0b01,0b10) = 1;
    swap(0b10,0b01) = 1;
    swap(0b11,0b11) = 1;
    auto fn_swap = [=](Sim & sim, unsigned ctrl, unsigned targ) {
		       sim.swap(ctrl, targ);
		   };
    gates.push_back({fn_swap, swap});

    // fswap
    arma::Mat<std::complex<Fp>> fswap(4, 4, arma::fill::zeros);
    fswap(0b00,0b00) = 1;
    fswap(0b01,0b10) = 1;
    fswap(0b10,0b01) = 1;
    fswap(0b11,0b11) = -1;
    auto fn_fswap = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.fswap(ctrl, targ);
		   };
    gates.push_back({fn_fswap, fswap});

    // npRotateX
    arma::Mat<std::complex<Fp>> nprx(4, 4, arma::fill::zeros);
    nprx(0b00,0b00) = 1;
    nprx(0b01,0b01) = std::complex<Fp>{std::cos(angle/2), 0};
    nprx(0b01,0b10) = std::complex<Fp>{0, -std::sin(angle/2)};
    nprx(0b10,0b01) = std::complex<Fp>{0, -std::sin(angle/2)};
    nprx(0b10,0b10) = std::complex<Fp>{std::cos(angle/2), 0};
    nprx(0b11,0b11) = 1;

    auto fn_nprx = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.npRotateX(targ, ctrl, angle);
		   };
    gates.push_back({fn_nprx, nprx});

    // npRotateY
    arma::Mat<std::complex<Fp>> npry(4, 4, arma::fill::zeros);
    npry(0b00,0b00) = 1;
    npry(0b01,0b01) = std::cos(angle/2);
    npry(0b01,0b10) = -std::sin(angle/2);
    npry(0b10,0b01) = std::sin(angle/2);
    npry(0b10,0b10) = std::cos(angle/2);
    npry(0b11,0b11) = 1;

    auto fn_npry = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.npRotateY(targ, ctrl, angle);
		   };
    gates.push_back({fn_npry, npry});

    // npHadamard
    arma::Mat<std::complex<Fp>> nph(4, 4, arma::fill::zeros);
    Fp sqrt2 = 1/std::sqrt(2);
    nph(0b00,0b00) = 1;
    nph(0b01,0b01) = sqrt2;
    nph(0b01,0b10) = sqrt2;
    nph(0b10,0b01) = sqrt2;
    nph(0b10,0b10) = -sqrt2;
    nph(0b11,0b11) = 1;

    auto fn_nph = [=](Sim & sim, unsigned ctrl, unsigned targ) {
	sim.npHadamard(targ, ctrl);
		   };
    gates.push_back({fn_nph, nph});

    for (const auto & [fn, mat] : gates) {

	// Make a random state
	Sim q{num_qubits};
	const std::vector<qsl::complex<Fp>> state
	    = qsl::makeRandomNPState<Fp>(num_qubits, num_ones);
	q.setState(state);

	// Set an armadillo vector to the same state
	std::size_t dim = 1 << num_qubits;
	arma::Col<std::complex<Fp>> v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    v(i) = std::complex<Fp>{state[i].real, state[i].imag};
	}

	// Apply gate to qubits
	std::invoke(fn, q, ctrl, targ);

	// Apply gate in armadillo
	v = makeMatrix(mat, num_qubits, {targ, ctrl}) * v;

	// Read qubit state into armadillo
	std::vector<qsl::complex<Fp>> res = q.getState();
	arma::Col<std::complex<Fp>> qubit_v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    qubit_v(i) = std::complex<Fp>{res[i].real, res[i].imag};
	}

    	EXPECT_TRUE(arma::approx_equal(v, qubit_v, "both", 1e-6, 1e-8));
    }
}

TYPED_TEST(NPGates, OneQubitGate)
{
    const unsigned num_qubits = 8;
    const unsigned num_ones = 4;
    const unsigned targ = 4;
    using Sim = TypeParam::Sim;
    using Fp = TypeParam::Sim::Fp_type;
    // Generate random angle to test parameterised gates
    Fp angle = arma::datum::pi * arma::randn<Fp>();

    // Create list of gates mapped to matrices
    std::vector<std::pair<std::function<void(Sim &, unsigned)>,
			  arma::SpMat<std::complex<Fp>>>> gates;

    // PauliZ
    arma::SpMat<std::complex<Fp>> pauliZ(2, 2);
    pauliZ(0, 0) = 1;
    pauliZ(1, 1) = -1;
    auto fn_pauliZ = [](Sim & sim, unsigned targ) {
			 sim.pauliZ(targ);
		     };
    gates.push_back({fn_pauliZ, pauliZ});

    // phase shift
    arma::SpMat<std::complex<Fp>> phase(2, 2);
    phase(0, 0) = 1;
    phase(1, 1) = std::complex<Fp>{std::cos(angle), std::sin(angle)};
    auto fn_phase = [=](Sim & sim, unsigned targ) {
			sim.phase(targ, angle);
		       };
    gates.push_back({fn_phase, phase});

    // rotateZ
    arma::SpMat<std::complex<Fp>> rotateZ(2, 2);
    rotateZ(0, 0) = std::complex<Fp>{std::cos(angle/2), std::sin(-angle/2)};
    rotateZ(1, 1) = std::complex<Fp>{std::cos(angle/2), std::sin(angle/2)};
    auto fn_rotateZ = [=](Sim & sim, unsigned targ) {
			sim.rotateZ(targ, angle);
		       };
    gates.push_back({fn_rotateZ, rotateZ});

    for (const auto & [fn, mat] : gates) {

	// Make a random state
	Sim q{num_qubits};
	const std::vector<qsl::complex<Fp>> state
	    = qsl::makeRandomNPState<Fp>(num_qubits, num_ones);
	q.setState(state);

	// Set an armadillo vector to the same state
	std::size_t dim = 1 << num_qubits;
	arma::Col<std::complex<Fp>> v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    v(i) = std::complex<Fp>{state[i].real, state[i].imag};
	}

	// Apply gate to qubits
	std::invoke(fn, q, targ);

	// Create gate in armadillo
	// Sizes of idenity matrix padding
	std::size_t pre = 1 << targ;
	std::size_t post = 1 << (num_qubits - targ - 1);

	// Tensor to make the gate
	arma::SpMat<std::complex<Fp>> gate =
	    arma::speye<arma::SpMat<std::complex<Fp>>>(pre, pre);
	gate = arma::kron(mat, gate);
	gate = arma::kron(arma::speye<arma::SpMat<std::complex<Fp>>>(post, post),
			  gate);

	// Apply gate in armadillo
	v = gate * v;

	// Read qubit state into armadillo
	std::vector<qsl::complex<Fp>> res = q.getState();
	arma::Col<std::complex<Fp>> qubit_v(dim);
	for (std::size_t i = 0; i < dim; i++) {
	    qubit_v(i) = std::complex<Fp>{res[i].real, res[i].imag};
	}

    	EXPECT_TRUE(arma::approx_equal(v, qubit_v, "both", 1e-6, 1e-8));
    }
}

max_stars_repo_name: lanamineh/qsl

max_issues_repo_name: lanamineh/qsl

max_forks_repo_name: lanamineh/qsl


text: #include <boost/archive/binary_iarchive_impl.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#define BOOST_TEST_MODULE tests

#include "lib.h"
#include <boost/test/unit_test.hpp>

BOOST_AUTO_TEST_SUITE(test_version)

BOOST_AUTO_TEST_CASE(test_valid_version)
{
    auto [major, minor, patch] = version();
    BOOST_CHECK(major == 0);
    BOOST_CHECK(minor == 0);
    BOOST_CHECK(patch > 0);
}

BOOST_AUTO_TEST_SUITE_END()

BOOST_AUTO_TEST_SUITE(test_message)

BOOST_AUTO_TEST_CASE(test_valid_message)
{
    BOOST_CHECK(message() == "Hello world!\n");
}

BOOST_AUTO_TEST_SUITE_END()

ext: cpp


text:
/* $Id: SPStream.cpp,v 1.99 2016/08/09 05:35:56 severin Exp $ */

/***

NAME - EEDB::SPStream

SYNOPSIS

DESCRIPTION

Abstract superclass for all stream signal-processing modules.
SPStream is short hand for Signal-Process-Stream

CONTACT

Jessica Severin <severin@gsc.riken.jp>

LICENSE

 * Software License Agreement (BSD License)
 * ZENBU [EEDB] system
 * copyright (c) 2007-2013 Jessica Severin RIKEN OSC
 * All rights reserved.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Neither the name of Jessica Severin RIKEN OSC nor the
 *       names of its contributors may be used to endorse or promote products
 *       derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

APPENDIX

The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _

***/

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string>
#include <stdarg.h>
#include <rapidxml.hpp>  //rapidxml must be include before boost
#include <boost/algorithm/string.hpp>
#include <MQDB/MappedQuery.h>
#include <EEDB/SPStream.h>
#include <EEDB/Experiment.h>

//for dynamic stream creation from XML
#include <EEDB/SPStreams/Dummy.h>
#include <EEDB/SPStreams/TemplateCluster.h>
#include <EEDB/SPStreams/MergeStreams.h>
#include <EEDB/SPStreams/FeatureEmitter.h>
#include <EEDB/SPStreams/MultiMergeStream.h>
#include <EEDB/SPStreams/FederatedSourceStream.h>
#include <EEDB/SPStreams/ObjectCount.h>
#include <EEDB/SPStreams/Proxy.h>
#include <EEDB/SPStreams/TemplateFilter.h>
#include <EEDB/SPStreams/CutoffFilter.h>
#include <EEDB/SPStreams/CalcFeatureSignificance.h>
#include <EEDB/SPStreams/IDFilter.h>
#include <EEDB/SPStreams/NormalizePerMillion.h>
#include <EEDB/SPStreams/CalcInterSubfeatures.h>
#include <EEDB/SPStreams/StreamSubfeatures.h>
#include <EEDB/SPStreams/UniqueFeature.h>
#include <EEDB/SPStreams/ResizeFeatures.h>
#include <EEDB/SPStreams/FeatureLengthFilter.h>
#include <EEDB/SPStreams/MakeStrandless.h>
#include <EEDB/SPStreams/TopHits.h>
#include <EEDB/SPStreams/ExpressionDatatypeFilter.h>
#include <EEDB/SPStreams/Paraclu.h>
#include <EEDB/SPStreams/NeighborCutoff.h>
#include <EEDB/SPStreams/NormalizeRPKM.h>
#include <EEDB/SPStreams/DevNull.h>
#include <EEDB/SPStreams/StreamTKs.h>
#include <EEDB/SPStreams/NormalizeByFactor.h>
#include <EEDB/SPStreams/FeatureRename.h>
#include <EEDB/SPStreams/RenameExperiments.h>
#include <EEDB/SPStreams/RescalePseudoLog.h>
#include <EEDB/SPStreams/FilterSubfeatures.h>
#include <EEDB/SPStreams/OverlapAnnotate.h>
#include <EEDB/SPStreams/OverlapMerge.h>
#include <EEDB/SPStreams/MetadataFilter.h>
#include <EEDB/SPStreams/MetadataManipulate.h>
#include <EEDB/SPStreams/StrandFilter.h>
#include <EEDB/SPStreams/FlipStrand.h>
#include <EEDB/SPStreams/MannWhitneyRanksum.h>
#include <EEDB/SPStreams/CachePoint.h>

using namespace std;
using namespace MQDB;

const char*               EEDB::SPStream::class_name = "EEDB::SPStream";

//function prototypes
//these are functions so need to pass obj as parameter
//this is also outside of class scope so only public attributes/methods available
MQDB::DBObject* _spstream_default_next_in_stream_func(EEDB::SPStream* node) {
  return ((EEDB::SPStream*)node)->_next_in_stream();
}
MQDB::DBObject* _spstream_default_fetch_object_by_id_func(EEDB::SPStream* node, string fid) {
  return ((EEDB::SPStream*)node)->_fetch_object_by_id(fid);
}
void _spstream_default_stream_clear_func(EEDB::SPStream* node) {
  ((EEDB::SPStream*)node)->_stream_clear();
}
bool _spstream_default_stream_by_named_region_func(EEDB::SPStream* node, string assembly_name, string chrom_name, long int start, long int end) {
  return ((EEDB::SPStream*)node)->_stream_by_named_region(assembly_name, chrom_name, start, end);
}
void _spstream_default_stream_features_by_metadata_search_func(EEDB::SPStream* node, string search_logic) {
  ((EEDB::SPStream*)node)->_stream_features_by_metadata_search(search_logic);
}
void _spstream_default_stream_data_sources_func(EEDB::SPStream* node, string classname, string filter_logic) {
  ((EEDB::SPStream*)node)->_stream_data_sources(classname, filter_logic);
}
void _spstream_default_stream_chromosomes_func(EEDB::SPStream* node, string assembly_name, string chrom_name) {
  ((EEDB::SPStream*)node)->_stream_chromosomes(assembly_name, chrom_name);
}
void _spstream_default_stream_peers_func(EEDB::SPStream* node) {
  ((EEDB::SPStream*)node)->_stream_peers();
}
void _spstream_default_reload_stream_data_sources_func(EEDB::SPStream* node) {
  ((EEDB::SPStream*)node)->_reload_stream_data_sources();
}
void _spstream_default_disconnect_func(EEDB::SPStream* node) {
  ((EEDB::SPStream*)node)->_disconnect();
}
void _spstream_default_get_proxies_by_name(EEDB::SPStream* node, string proxy_name, vector<EEDB::SPStream*> &proxies) {
  ((EEDB::SPStream*)node)->_get_proxies_by_name(proxy_name, proxies);
}
void _spstream_default_get_dependent_datasource_ids(EEDB::SPStream* node, map<string, bool> &source_ids) {
  ((EEDB::SPStream*)node)->_get_dependent_datasource_ids(source_ids);
}
void _spstream_default_reset_stream_node_func(EEDB::SPStream* node) {
  //this is not passed down the stream, but can be replaced by subclasses
  //with a function which resets pointers,counters or other internal
  //variables to bring the stream back to a clean starting point
  ((EEDB::SPStream*)node)->_reset_stream_node();
}

void _eedb_spstream_delete_func(MQDB::DBObject *obj) {
  delete (EEDB::SPStream*)obj;
}
void _eedb_spstream_xml_func(MQDB::DBObject *obj, string &xml_buffer) {
  ((EEDB::SPStream*)obj)->_xml(xml_buffer);
}
void _eedb_spstream_simple_xml_func(MQDB::DBObject *obj, string &xml_buffer) {
  ((EEDB::SPStream*)obj)->_xml(xml_buffer);
}
string _eedb_spstream_display_desc_func(MQDB::DBObject *obj) {
  return ((EEDB::SPStream*)obj)->_display_desc();
}

//
////////////////////////////////////////////////////////////////////////////////////////
//

EEDB::SPStream::SPStream() {
  init();
}

EEDB::SPStream::~SPStream() {
}

void EEDB::SPStream::init() {
  MQDB::DBObject::init();
  _classname                 = EEDB::SPStream::class_name;
  _funcptr_delete            = _eedb_spstream_delete_func;
  _funcptr_xml               = _eedb_spstream_xml_func;
  _funcptr_simple_xml        = _eedb_spstream_simple_xml_func;
  _funcptr_display_desc      = _eedb_spstream_display_desc_func;

  _source_stream  = NULL;
  _side_stream    = NULL;
  _region_start   = -1;
  _region_end     = -1;

  //function pointer code
  _funcptr_next_in_stream                     = _spstream_default_next_in_stream_func;
  _funcptr_fetch_object_by_id                 = _spstream_default_fetch_object_by_id_func;
  _funcptr_disconnect                         = _spstream_default_disconnect_func;
  _funcptr_stream_clear                       = _spstream_default_stream_clear_func;
  _funcptr_stream_by_named_region             = _spstream_default_stream_by_named_region_func;
  _funcptr_stream_features_by_metadata_search = _spstream_default_stream_features_by_metadata_search_func;
  _funcptr_stream_data_sources                = _spstream_default_stream_data_sources_func;
  _funcptr_stream_chromosomes                 = _spstream_default_stream_chromosomes_func;
  _funcptr_stream_peers                       = _spstream_default_stream_peers_func;
  _funcptr_reload_stream_data_sources         = _spstream_default_reload_stream_data_sources_func;
  _funcptr_reset_stream_node                  = _spstream_default_reset_stream_node_func;
  _funcptr_get_proxies_by_name                = _spstream_default_get_proxies_by_name;
  _funcptr_get_dependent_datasource_ids       = _spstream_default_get_dependent_datasource_ids;

}

void EEDB::SPStream::display_info() {
  printf("%s\n", display_desc().c_str());
}

string EEDB::SPStream::_display_desc() {
  string xml = "SPStream[";
  xml += _classname;
  xml += "]";
  return xml;
}

string EEDB::SPStream::display_contents() {
  return display_desc();
}

void EEDB::SPStream::_xml_start(string &xml_buffer) {
  //depth first search to reverse order
  if(_source_stream) {
    _source_stream->xml(xml_buffer);
  }

  xml_buffer.append("<spstream module=\"");
  xml_buffer.append(_module_name);
  xml_buffer.append("\">");
}

void EEDB::SPStream::_xml_end(string &xml_buffer) {
  if(_side_stream) {
    xml_buffer.append("<side_stream>");
    _side_stream->xml(xml_buffer);
    xml_buffer.append("</side_stream>");
  }
  xml_buffer.append("</spstream>");
}

void EEDB::SPStream::_xml(string &xml_buffer) {
  _xml_start(xml_buffer);
  _xml_end(xml_buffer);
}

////////////////////////////////////////////////////////////////////////////

void EEDB::SPStream::unparse_eedb_id(string fid, string &uuid, long int &objID, string &objClass) {
  uuid.clear();
  objClass="Feature";
  objID = -1;

  size_t   p1;
  if((p1 = fid.find("::")) != string::npos) {
    uuid = fid.substr(0, p1);
    string idclass = fid.substr(p1+2);
    if((p1 = idclass.find(":::")) != string::npos) {
      objClass = idclass.substr(p1+3);
      string t_id = idclass.substr(0, p1);
      objID = atoi(t_id.c_str());
    } else {
      objID = atoi(idclass.c_str());
    }
  }
}

////////////////////////////////////////////////////////////////////////////
// Instance methods
////////////////////////////////////////////////////////////////////////////

/***** source_stream

  Description: set the input or source stream feeding objects into this level of the stream stack.
  Returntype : either an MQdb::DBStream or subclass of EEDB::SPStream object or NULL if not set
  Exceptions : none

*****/

EEDB::SPStream*  EEDB::SPStream::source_stream() { return (EEDB::SPStream*)_source_stream; }

void  EEDB::SPStream::source_stream(EEDB::SPStream* value) {
  if(value == NULL) {
    _source_stream = NULL;
  } else {
    if(value == this) {
      fprintf(stderr, "ERROR: can not set source_stream() to itself\n");
      return;
    } else {
      _source_stream = value;
    }

  }
}

/***** side_stream

 Description: set the input or source stream feeding objects into this level of the stream stack.
 Returntype : either an MQdb::DBStream or subclass of EEDB::SPStream object or NULL if not set
 Exceptions : none

 *****/

EEDB::SPStream*  EEDB::SPStream::side_stream() {
  return (EEDB::SPStream*)_side_stream;
}

void  EEDB::SPStream::side_stream(EEDB::SPStream* value) {
  if(value == NULL) {
    _side_stream = NULL;
  } else {
    if(value == this) {
      fprintf(stderr, "ERROR: can not set side_stream() to itself\n");
      return;
    } else {
      _side_stream = value;
    }

  }
}

////////////////////////////////////////////////////////////////////////////
//
// override method for subclasses which will
// do all the work
//
////////////////////////////////////////////////////////////////////////////

/***** next_in_stream
  Description: return the next object in the stream stack
               depending on the configuration this will either be an on-the-fly created object
               or an object passed through filters, or a primary object streamed out of database/peer
  Returntype : either a DBObject or NULL if end of stream
*****/

MQDB::DBObject*  EEDB::SPStream::next_in_stream() {
  return _funcptr_next_in_stream(this);
}

/***** fetch_object_by_id
  Description: fetch single feature object from stream using global federated-ID.
               passes query down the stream until it is satisfied
  Arg (1)    : $fedID (federated ID <uuid>::<id>:::<class>)
*****/

MQDB::DBObject* EEDB::SPStream::fetch_object_by_id(string id) {
  reset_stream_node();
  return _funcptr_fetch_object_by_id(this, id);
}

/***** stream_clear
  Description: re-initialize the stream-stack back to a clear/empty state
*****/

void EEDB::SPStream::stream_clear() {
  _funcptr_stream_clear(this);
}

/***** stream_by_named_region
  Description: configure/initialize the stream-stack to a stream from a named region.
               This will be passed down the stack until it reaches an SPStream::SourceStream instance
               which then knows how to create a new stream from a database.
  Arg (1)    : $assembly_name (string)
  Arg (2)    : $chrom_name (string)
  Arg (3)    : $chrom_start (integer)
  Arg (4)    : $chrom_end (integer)
*****/

bool  EEDB::SPStream::stream_by_named_region(string assembly_name, string chrom_name, long int start, long int end) {
  reset_stream_node();
  return _funcptr_stream_by_named_region(this, assembly_name, chrom_name, start, end);
}

/***** stream_data_sources
  Description: stream all sources(FeatureSource, EdgeSource, and Experiment) out of database
*****/

void EEDB::SPStream::stream_data_sources() {
  reset_stream_node();
  _funcptr_stream_data_sources(this, string(""), string(""));
}

void EEDB::SPStream::stream_data_sources(string classname) {
  reset_stream_node();
  _funcptr_stream_data_sources(this, classname, "");
}

void EEDB::SPStream::stream_data_sources(string classname, string filter_logic) {
  reset_stream_node();
  _funcptr_stream_data_sources(this, classname, filter_logic);
}

void EEDB::SPStream::get_dependent_datasource_ids(map<string,bool> &source_ids) {
  reset_stream_node();
  _funcptr_get_dependent_datasource_ids(this, source_ids);
}

/***** reload_stream_data_sources
  Description: SourceStream will cache sources(FeatureSource, EdgeSource, and Experiment)
               System automatically will do a reload ever hour to check for new sources
               added to previously cached peer/databases. But it may be useful to manually
               force a cache refresh.
*****/

void EEDB::SPStream::reload_stream_data_sources() {
  reset_stream_node();
  _funcptr_reload_stream_data_sources(this);
}

/***** stream_chromosomes
  Description: stream all EEDB::Chrom chromosomes from databases on the stream
*****/

void EEDB::SPStream::stream_chromosomes(string assembly_name, string chrom_name) {
  reset_stream_node();
  _funcptr_stream_chromosomes(this, assembly_name, chrom_name);
}

/***** stream_peers
  Description: stream all known peers from database
*****/

void EEDB::SPStream::stream_peers() {
  reset_stream_node();
  _funcptr_stream_peers(this);
}

/***** stream_features_by_metadata_search
  Description: perform search of features through metadata filter logic
  Arg (1)    : a keyword/logic string which is applied to the metadata (optional)
*****/

void  EEDB::SPStream::stream_features_by_metadata_search(string search_logic) {
  reset_stream_node();
  _funcptr_stream_features_by_metadata_search(this, search_logic);
}

/***** disconnect
  Description: send "disconnect" message to all spstream modules
*****/

void  EEDB::SPStream::disconnect() {
  _funcptr_disconnect(this);
}

/***** reset_stream_node
  Description: can be reimplemented by subclasses.
               This is called at the beginning of every stream_xxxxx() call.
               Useful for reseting/clearing internal caches/variables between
               different streams.
*****/

void EEDB::SPStream::reset_stream_node() {
  _funcptr_reset_stream_node(this);
}

////////////////////////////////////////////////////////////////////////////
//
//  creation from XML section
//
////////////////////////////////////////////////////////////////////////////

void  EEDB::SPStream::create_stream_from_xml(void *xml_node, SPStream* &head, SPStream* &tail) {
  //linear order in XML is order of appending onto 'head'
  rapidxml::xml_node<>  *root_node = (rapidxml::xml_node<>*)xml_node;

  head = NULL;
  tail = NULL;

  if(root_node == NULL) { return; }

  rapidxml::xml_node<> *stream_node = root_node->first_node("spstream");
  while(stream_node) {
    EEDB::SPStream *spstream = EEDB::SPStream::_xmlnode_create_spstream(stream_node);
    if(spstream!=NULL) {
      if(tail==NULL) { tail = spstream; }
      if(head!=NULL) { spstream->source_stream(head); }
      head = spstream;
    }
    stream_node = stream_node->next_sibling("spstream");
  }
}

void  EEDB::SPStream::create_side_stream_from_xml(void *xml_node) {
  //used by stream modules which need to create a side_stream on them selves
  rapidxml::xml_node<>  *node = (rapidxml::xml_node<>*)xml_node;

  EEDB::SPStream *head  = NULL;
  EEDB::SPStream *tail  = NULL;
  EEDB::SPStream::create_stream_from_xml(node, head, tail);
  if(head) { side_stream(head); }
}

EEDB::SPStream*  EEDB::SPStream::reverse_stream(SPStream* node) {
  //recursive algorithm to reverse node and its source_stream()
  // 'node' is the parent of this, so in reverse this->source_stream should point at parent 'node'
  // returns the new head which is the deepest node of the stream

  if(_side_stream) {
    _side_stream = _side_stream->reverse_stream(NULL);
  }

  EEDB::SPStream *new_head = this;
  if(_source_stream!=NULL) {
    new_head = _source_stream->reverse_stream(this);
  }

  _source_stream = node;
  return new_head;
}

EEDB::SPStream* EEDB::SPStream::_xmlnode_create_spstream(void *xml_node) {
  //master mapper from named modules to actual class binding (not a dynamically typed language)
  if(xml_node==NULL) { return NULL; }

  rapidxml::xml_node<>      *root_node = (rapidxml::xml_node<>*)xml_node;
  rapidxml::xml_attribute<> *attr;
  EEDB::SPStream            *spstream = NULL;

  if(string(root_node->name()) != "spstream") { return NULL; }

  attr = root_node->first_attribute("module");
  if(!attr) { return NULL;}
  string modname = attr->value();

  size_t p1;
  string tstr;
  if((p1 = modname.find("EEDB::SPStream::")) != string::npos)  { tstr=modname; modname = tstr.substr(16); }
  if((p1 = modname.find("EEDB::SPStreams::")) != string::npos) { tstr=modname; modname = tstr.substr(17); }

  //massive if block now....but in the future maybe a directory of runtime linked plugins

  if(modname == "TemplateCluster")   { spstream = new EEDB::SPStreams::TemplateCluster(xml_node); }
  if(modname == "CollateExpression") { spstream = new EEDB::SPStreams::TemplateCluster(xml_node); }

  //OverlapCluster and OverlapMerge do the same thing, changed namecd
  if(modname == "OverlapCluster")    { spstream = new EEDB::SPStreams::OverlapMerge(xml_node); } //OverlapMerge performs denovo clustering
  if(modname == "OverlapMerge")      { spstream = new EEDB::SPStreams::OverlapMerge(xml_node); } //OverlapMerge performs denovo clustering

  //OverlapFilter and TemplateFilter do the same thing
  if(modname == "TemplateFilter")    { spstream = new EEDB::SPStreams::TemplateFilter(xml_node); }
  if(modname == "OverlapFilter")     { spstream = new EEDB::SPStreams::TemplateFilter(xml_node); }

  if(modname == "Dummy")             { spstream = new EEDB::SPStreams::Dummy(xml_node); }
  if(modname == "MergeStreams")      { spstream = new EEDB::SPStreams::MergeStreams(xml_node); }
  if(modname == "FeatureEmitter")    { spstream = new EEDB::SPStreams::FeatureEmitter(xml_node); }
  if(modname == "ObjectCount")       { spstream = new EEDB::SPStreams::ObjectCount(xml_node); }
  if(modname == "Proxy")             { spstream = new EEDB::SPStreams::Proxy(xml_node); }

  if(modname == "FederatedSourceStream")     { spstream = new  EEDB::SPStreams::FederatedSourceStream(xml_node); }
  if(modname == "CachePoint")                { spstream = new  EEDB::SPStreams::CachePoint(xml_node); }

  if(modname == "CutoffFilter")              { spstream = new  EEDB::SPStreams::CutoffFilter(xml_node); }
  if(modname == "CalcFeatureSignificance")   { spstream = new  EEDB::SPStreams::CalcFeatureSignificance(xml_node); }
  if(modname == "IDFilter")                  { spstream = new  EEDB::SPStreams::IDFilter(xml_node); }

  if(modname == "NormalizePerMillion")       { spstream = new  EEDB::SPStreams::NormalizePerMillion(xml_node); }
  if(modname == "CalcInterSubfeatures")      { spstream = new  EEDB::SPStreams::CalcInterSubfeatures(xml_node); }
  if(modname == "StreamSubfeatures")         { spstream = new  EEDB::SPStreams::StreamSubfeatures(xml_node); }
  if(modname == "UniqueFeature")             { spstream = new  EEDB::SPStreams::UniqueFeature(xml_node); }
  if(modname == "ResizeFeatures")            { spstream = new  EEDB::SPStreams::ResizeFeatures(xml_node); }
  if(modname == "ExpressionDatatypeFilter")  { spstream = new  EEDB::SPStreams::ExpressionDatatypeFilter(xml_node); }
  if(modname == "FeatureLengthFilter")       { spstream = new  EEDB::SPStreams::FeatureLengthFilter(xml_node); }
  if(modname == "MakeStrandless")            { spstream = new  EEDB::SPStreams::MakeStrandless(xml_node); }
  if(modname == "TopHits")                   { spstream = new  EEDB::SPStreams::TopHits(xml_node); }
  if(modname == "NormalizeRPKM")             { spstream = new  EEDB::SPStreams::NormalizeRPKM(xml_node); }
  if(modname == "NormalizeByFactor")         { spstream = new  EEDB::SPStreams::NormalizeByFactor(xml_node); }
  if(modname == "DevNull")                   { spstream = new  EEDB::SPStreams::DevNull(xml_node); }
  if(modname == "StreamTKs")                 { spstream = new  EEDB::SPStreams::StreamTKs(xml_node); }
  if(modname == "FeatureRename")             { spstream = new  EEDB::SPStreams::FeatureRename(xml_node); }
  if(modname == "RenameExperiments")         { spstream = new  EEDB::SPStreams::RenameExperiments(xml_node); }
  if(modname == "RescalePseudoLog")          { spstream = new  EEDB::SPStreams::RescalePseudoLog(xml_node); }
  if(modname == "FilterSubfeatures")         { spstream = new  EEDB::SPStreams::FilterSubfeatures(xml_node); }
  if(modname == "StrandFilter")              { spstream = new  EEDB::SPStreams::StrandFilter(xml_node); }
  if(modname == "FlipStrand")                { spstream = new  EEDB::SPStreams::FlipStrand(xml_node); }

  if(modname == "Paraclu")                   { spstream = new  EEDB::SPStreams::Paraclu(xml_node); }
  if(modname == "NeighborCutoff")            { spstream = new  EEDB::SPStreams::NeighborCutoff(xml_node); }

  if(modname == "OverlapAnnotate")           { spstream = new  EEDB::SPStreams::OverlapAnnotate(xml_node); }
  if(modname == "MetadataFilter")            { spstream = new  EEDB::SPStreams::MetadataFilter(xml_node); }
  if(modname == "MetadataManipulate")        { spstream = new  EEDB::SPStreams::MetadataManipulate(xml_node); }

  if(modname == "MannWhitneyRanksum")        { spstream = new  EEDB::SPStreams::MannWhitneyRanksum(xml_node); }

  return spstream;
}

EEDB::SPStream::SPStream(void *xml_node) {
  //constructor using a rapidxml <spstream> description
  // default super class method which does nothing.
  // subclasses need only implement this method
  init();
}

void  EEDB::SPStream::get_proxies_by_name(string proxy_name, vector<EEDB::SPStream*> &proxies) {
  reset_stream_node();
  _funcptr_get_proxies_by_name(this, proxy_name, proxies);
}

//
// default stream API methods
//

void  EEDB::SPStream::_reset_stream_node() {
  _region_start   = -1;
  _region_end     = -1;
}

MQDB::DBObject* EEDB::SPStream::_next_in_stream() {
  if(source_stream() != NULL) { return source_stream()->next_in_stream(); }
  return NULL;
}

MQDB::DBObject* EEDB::SPStream::_fetch_object_by_id(string fid) {
  if(source_stream() != NULL) { return source_stream()->fetch_object_by_id(fid); }
  return NULL;
}

void EEDB::SPStream::_stream_clear() {
  if(source_stream() != NULL) { source_stream()->stream_clear(); }
}

bool EEDB::SPStream::_stream_by_named_region(string assembly_name, string chrom_name, long int start, long int end) {
  _region_start = start;
  _region_end   = end;
  if(source_stream() != NULL) {
    return source_stream()->stream_by_named_region(assembly_name, chrom_name, start, end);
  }
  return false;
}

void EEDB::SPStream::_stream_features_by_metadata_search(string search_logic) {
  if(source_stream() != NULL) { source_stream()->stream_features_by_metadata_search(search_logic); }
}

void EEDB::SPStream::_stream_data_sources(string classname, string filter_logic) {
  if(source_stream() != NULL) { source_stream()->stream_data_sources(classname, filter_logic); }
}

void EEDB::SPStream::_get_dependent_datasource_ids(map<string, bool> &source_ids) {
  if(source_stream() != NULL) { source_stream()->get_dependent_datasource_ids(source_ids); }
}

void EEDB::SPStream::_stream_chromosomes(string assembly_name, string chrom_name) {
  if(source_stream() != NULL) { source_stream()->stream_chromosomes(assembly_name, chrom_name); }
}

void EEDB::SPStream::_stream_peers() {
  if(source_stream() != NULL) { source_stream()->stream_peers(); }
}

void EEDB::SPStream::_reload_stream_data_sources() {
  if(source_stream() != NULL) { source_stream()->reload_stream_data_sources(); }
}

void EEDB::SPStream::_disconnect() {
  if(source_stream() != NULL) { source_stream()->disconnect(); }
}

void EEDB::SPStream::_get_proxies_by_name(string proxy_name, vector<EEDB::SPStream*> &proxies) {
  if(source_stream() != NULL) { source_stream()->get_proxies_by_name(proxy_name, proxies); }
}