//! Computing the interpolative decomposition of a matrix.
//!
//! This example computes the rank k interpolative decomposition of an $m\times n$ matrix.
//!

use rusty_compression::*;

pub fn main() {

    // We initialize a random number generator.
    let mut rng = rand::thread_rng();

    // The dimension of the matrix for which we want to
    // compute a low-rank approximation.
    let dimension = (500, 100);

    // The compression rank.
    let k = 20;

    // Generate a random matrix with singular values logarithmically
    // distributed between 1 and 1E-10.
    let mat = f64::random_approximate_low_rank_matrix(dimension, 1.0, 1E-10, &mut rng);

    // Compute the pivoted QR decomposition of the matrix.
    let qr = QR::<f64>::compute_from(mat.view()).expect("Could not compress the matrix.");

    // We now compress the pivoted QR decomposition to only include the k most significant
    // basis vectors of the range.
    let qr_compressed = qr.compress(CompressionType::RANK(k)).unwrap();

    // From the compressed representation we now compute the column interpolative decomposition
    let col_int_decomp = qr_compressed.column_id().unwrap();

    // We can also compute a two sided interpolative decomposition.
    let two_sided_int_decomp = col_int_decomp.two_sided_id().unwrap();

    // Let us compare the relative difference of the compressed operator to the original matrix.

    // First we multiply the factors of the two-sided decomposition to get back a matrix.
    // This should only be done for debugging purposes and non-probabilistic relative error
    // computations.
    let mat_approx = two_sided_int_decomp.to_mat();

    // We now compue the relative difference.
    let rel_diff = f64::rel_diff_fro(mat.view(), mat_approx.view());

    // Print the result.
    println!("The relative difference of the compressed and original matrix is {:1.2E}", rel_diff);


}