/*
This file is part of the Kernel Quantum Probability library (KQP).

KQP is free software: you can redistribute it and/or modify
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

KQP is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with KQP.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <kqp/feature_matrix/dense.hpp>
#include <kqp/kernel_evd/incremental.hpp>
#include <kqp/probabilities.hpp>

int main(int, const char**) {

// --- Compute a density at random

// Definitions
using namespace kqp;
typedef Eigen::Matrix<double,Dynamic,Dynamic> Matrix;
int dim = 10;

// Creating an incremental builder
IncrementalKernelEVD<DenseMatrix<double>> kevd;

// Add 10 vectors with $\alpha_i=1$
for(int i = 0; i < 10; i++) {
// Adds a random $\varphi_i$
Matrix m = Matrix::Random(dim, 1);
}

// Get the result $\rho \approx X Y D Y^\dagger X^\dagger$
Decomposition<DenseMatrix<double>> d = kevd.getDecomposition();

// --- Compute a kEVD for a subspace

IncrementalKernelEVD<DenseMatrix<double>> kevd_event;
for(int i = 0; i < 3; i++) {
// Adds a random $\varphi_i$
Matrix m = Matrix::Random(dim, 1);
}

// --- Compute some probabilities

// Setup densities and events
Density<DenseMatrix<double>> rho(kevd);
rho.normalize();
Event<DenseMatrix<double>> event(kevd_event);

// Compute the probability
std::cout << "Probability = " << rho.probability(event) << std::endl;

// Conditional probability
Density<DenseMatrix<double>> rho_cond = event.project(kevd).normalize();
std::cout << "Entropy rho/E = " << rho_cond.entropy() << std::endl;

// Conditional probability (orthogonal event)
Density<DenseMatrix<double>> rho_cond_orth = event.project(kevd, true).normalize();
std::cout << "Entropy rho/not E = " << rho_cond_orth.entropy() << std::endl;

return 0;
}