// $Id$
//
//  Copyright (C) 2004-2006 Rational Discovery LLC
//
//   @@ All Rights Reserved @@
//  This file is part of the RDKit.
//  The contents are covered by the terms of the BSD license
//  which is included in the file license.txt, found at the root
//  of the RDKit source tree.
//
#include "BoundsMatrix.h"
#include "TriangleSmooth.h"

namespace DistGeom {
bool triangleSmoothBounds(BoundsMatPtr boundsMat, double tol) {
  return triangleSmoothBounds(boundsMat.get(), tol);
}
bool triangleSmoothBounds(BoundsMatrix *boundsMat, double tol) {
  int npt = boundsMat->numRows();
  int i, j, k;
  double Uik, Lik, Ukj, sumUikUkj, diffLikUjk, diffLjkUik;

  for (k = 0; k < npt; k++) {
    for (i = 0; i < npt - 1; i++) {
      if (i == k) {
        continue;
      }
      Uik = boundsMat->getUpperBound(i, k);
      Lik = boundsMat->getLowerBound(i, k);
      for (j = i + 1; j < npt; j++) {
        if (j == k) {
          continue;
        }
        Ukj = boundsMat->getUpperBound(k, j);
        sumUikUkj = Uik + Ukj;
        if (boundsMat->getUpperBound(i, j) > sumUikUkj) {
          boundsMat->setUpperBound(i, j, sumUikUkj);
        }

        diffLikUjk = Lik - Ukj;
        diffLjkUik = boundsMat->getLowerBound(j, k) - Uik;
        if (boundsMat->getLowerBound(i, j) < diffLikUjk) {
          boundsMat->setLowerBound(i, j, diffLikUjk);
        } else if (boundsMat->getLowerBound(i, j) < diffLjkUik) {
          boundsMat->setLowerBound(i, j, diffLjkUik);
        }
        double lBound = boundsMat->getLowerBound(i, j);
        double uBound = boundsMat->getUpperBound(i, j);
        // std::cerr << "     smooth: " << i << "-" << j << "(" << k
        //           << "): " << boundsMat->getLowerBound(i, j) << " "
        //           << boundsMat->getUpperBound(i, j) << " lik: " << Lik
        //           << " uik: " << Uik
        //           << " lkj: " << boundsMat->getLowerBound(j, k)
        //           << " ukj: " << Ukj << "\n";

        if (tol > 0. && (lBound - uBound) / lBound > 0. &&
            (lBound - uBound) / lBound < tol) {
          boundsMat->setUpperBound(i, j, lBound);
          uBound = lBound;
        }
        if (lBound - uBound > 0.) {
          // std::cerr << std::endl;
          // for (unsigned int ii = 0; ii < npt; ++ii) {
          //   for (unsigned int jj = 0; jj < npt; ++jj) {
          //     std::cerr << " " << std::setprecision(3)
          //               << boundsMat->getVal(ii, jj);
          //   }
          //   std::cerr << std::endl;
          // }
          // std::cerr << std::endl;
          // std::cerr << " Fail: " << i << "-" << j << ": "
          //           << boundsMat->getLowerBound(i, j) << " "
          //           << boundsMat->getUpperBound(i, j) << "\n";
          return false;
        }
      }
    }
  }
  return true;
}
}