/* This file is part of the Gudhi Library. The Gudhi library
* (Geometric Understanding in Higher Dimensions) is a generic C++
* library for computational topology.
*
* Author(s): Pawel Dlotko
*
* Copyright (C) 2015 INRIA (France)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 this program. If not, see .
*/
#ifndef FILL_IN_MISSING_DATA_H
#define FILL_IN_MISSING_DATA_H
namespace Gudhi
{
namespace Gudhi_stat
{
/**
* Quite often in biological sciences we are facing a problem of missing data. We may have for instance a number of sequences of observations made in between times A and B in a discrete
* collection of times A = t1, t2,...,tn = B. But quite typically some of the observations may be missing. Then quite often it is hard to estimate the values in the missing times.
* The procedure below assumes that we compute some topological descriptor of the observations we are given. Most typically this will be some type of persistence homology representation.
* The the values in the missing points are filled in by linear interpolation and extrapolation. The procedure below have minimal requirements: it assumes that we have two elements that
* filled in. The rest will be filled in by using them.
* The fill-in process is done based on the idea of linear approximation. Let us assume that we have two positions A and B which are filled in with proper objects of a type Representation_of_topology.
* Any intermediate time step can be interpolated by taking a linear combination t*A + (1-t)*B, where t \in [0,1].
* If the missing data are located at the beginning of at the end of vector of Representation_of_topology, then we use the following extrapolation scheme:
* First we make sure that there are no missing data except at the very beginning and at the very end of a vector of data. Then, we pick two constitutive closest filled-in data A and B and we extrapolate
* by using a formula: t*A-(t-1)*B, where t is a natural number > 1.
* Note that both vector of data and the vector is_the_position_filled_in are modified by this procedure. Upon successful termination of the procedure, the vector is_the_position_filled_in has only 'true' entries,
* and the vector data do not have missing data.
**/
template < typename Representation_of_topology >
void fill_in_missing_data( std::vector< Representation_of_topology* >& data , std::vector< bool >& is_the_position_filled_in )
{
bool dbg = false;
//first check if at least two positions are filled in:
size_t number_of_positions_that_are_filled_in = 0;
for ( size_t i = 0 ; i != is_the_position_filled_in.size() ; ++i )
{
if ( is_the_position_filled_in[i] )++number_of_positions_that_are_filled_in;
}
if ( number_of_positions_that_are_filled_in < 2 )
{
std::cerr << "There are too few positions filled in to do extrapolation / interpolation. The program will now terminate.\n";
throw "There are too few positions filled in to do extrapolation / interpolation. The program will now terminate.\n";
}
for ( size_t i = 0 ; i != data.size() ; ++i )
{
if ( !is_the_position_filled_in[i] )
{
//This position is not filled in. Find the next position which is nonzero.
size_t j = 1;
while ( (is_the_position_filled_in[i+j] == false) && ( i+j != data.size() ) )++j;
if ( dbg )
{
std::cout << "The position number : " << i << " is not filled in. The next filled-in position is : " << i+j << std::endl;
}
if ( i != 0 )
{
//this is not the first position of the data:
if ( i + j != data.size() )
{
//this is not the last position of the data either
for ( size_t k = 0 ; k != j ; ++k )
{
double weight1 = double(j-k)/(double)(j+1);
double weight2 = double(k+1)/(double)(j+1);
data[i+k] = new Representation_of_topology(weight1*(*data[i-1]) + weight2*(*data[i+j]));
is_the_position_filled_in[i+k] = true;
if ( dbg )
{
std::cerr << "We fill in a position : " << i+k << " with: position " << i-1 << " with weight " << weight1 << " and position " << i+j << " with weight " << weight2 << std::endl;
}
}
}
}
else
{
//this is the first position of the data, i.e. i == 0.
while ( is_the_position_filled_in[i] == 0 )++i;
}
}
}
if ( is_the_position_filled_in[0] == false )
{
//find the first nonzero (then, we know that the second one will be nonzero too, since we filled it in above):
size_t i = 0;
while ( is_the_position_filled_in[i] == false )++i;
//the data at position i is declared. Since, we made sure that all other positions, except maybe a few first and a few last, are declared too. Therefore, if we find a
//first declared, then the next one will be declared too. So, we can do a telescopic declaration backward and forward.
for ( size_t j = i ; j != 0 ; --j )
{
data[j-1] = new Representation_of_topology(2*(*data[j]) + (-1)*(*data[j+1]));
is_the_position_filled_in[j-1] = true;
if ( dbg )
{
std::cerr << "Filling in a position : " << j-1 << " by using: " << j << " and " << j+1 <