source: trunk/c++_tools/utility/kNNI.cc @ 616

// $Id: kNNI.cc 616 2006-08-31 08:52:02Z jari $

/*
  Copyright (C) 2004 Jari Häkkinen
  Copyright (C) 2005 Peter Johansson
  Copyright (C) 2006 Jari Häkkinen

  This file is part of the thep c++ tools library,
                                http://lev.thep.lu.se/trac/c++_tools

  The c++ tools library 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 2 of
  the License, or (at your option) any later version.

  The c++ tools library 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, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include <c++_tools/utility/kNNI.h>

#include <algorithm>
#include <cmath>
#include <fstream>
#include <vector>

#include <c++_tools/utility/stl_utility.h>

namespace theplu {
namespace utility {

  kNNI::kNNI(const utility::matrix& matrix,const utility::matrix& flag,
             const u_int neighbours)
    : NNI(matrix,flag,neighbours)
  {
    for (unsigned int i=0; i<weight_.rows(); i++)
      for (unsigned int j=0; j<weight_.columns(); j++)
        if (!weight_(i,j)) {
          mv_rows_.push_back(i);
          break;
        }
    //estimate();
  }



  // \hat{x_{ij}}=\frac{ \sum_{k=1,N} \frac{x_{kj}}{d_{ki}} }
  //                   { \sum_{k=1,N} \frac{1     }{d_{ki}} },
  // where N is defined in the paper cited in the NNI class definition
  // documentation.
  u_int kNNI::estimate(void)
  {
    for (unsigned int i=0; i<mv_rows_.size(); i++) {
      // Jari, avoid copying in next line
      std::vector<std::pair<u_int,double> > distance=
        calculate_distances(mv_rows_[i]);
      std::sort(distance.begin(),distance.end(),
                pair_value_compare<u_int,double>());
      for (unsigned int j=0; j<data_.columns(); j++)
        if (!weight_(mv_rows_[i],j)) {
          std::vector<u_int> knn=nearest_neighbours(j,distance);
          double new_value=0.0;
          double norm=0.0;
          for (std::vector<u_int>::const_iterator k=knn.begin(); k!=knn.end();
               ++k) {
            // Jari, a small number needed here, use something standardized.
            // Avoid division with zero (perfect match vectors)
            double d=(distance[*k].second ? distance[*k].second : 1e-10);
            new_value+=data_(distance[*k].first,j)/d;
            norm+=1.0/d;
          }
          // No impute if no contributions from neighbours.
          if (norm)
            imputed_data_(mv_rows_[i],j)=new_value/norm;
          else {
            not_imputed_.push_back(i);
            // if norm is zero for one column it is zero for all columns
            // having zero weight
            break;
          }
        }
    }
    return not_imputed_.size();
  }


}} // of namespace utility and namespace theplu