// $Id: pattern.C,v 1.13 2001/03/06 22:32:42 tgkolda Exp $ 
// $Source: /usr/local/cvsroot/appspack/apps/src/pattern.C,v $ 

// ASYNCHRONOUS PARALLEL PATTERN SEARCH (APPS)
// COPYRIGHT (2000) Sandia Corporation.
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// LICENSE & WARRANTY INFORMATION in README.txt and LICENSE.txt.
// CONTACT Tamara G. Kolda (tgkolda@sandia.gov).

// ---------------------------------------------------------------------
// See pattern.H for description.
// ---------------------------------------------------------------------

#include <iomanip>		// for setw
#include <fstream>		// file i/o

#ifndef BRAINDEAD_NAMESPACE
using namespace std;
#endif

#ifndef BRAINDEAD_HEADERS
#include <cstdlib>		// for srand
#include <ctime>		// for clock()
#include <cmath>		// for sqrt(double)
#include <cstring>		// for strlen
#else
#include <stdlib.h>	
#include <time.h>	
#include <math.h>	
#include <string.h>
#endif

#include "pattern.H"
#include "apps.H"		// for idstr, precision, debug
#include "gci.H"		// for pack and unpack

// Initialize pattern types
const int COORD_MAX = 0;
const int SPHERICAL_MIN = 1;
const int COORD_MIN = 2;
const int SPHERICAL_MAX = 3;
const int USER_DEF = 4;

Pattern::Pattern()
{
  ndim = nsearch = 0;		// placeholder
  type = -1;			// placeholder
  minnsearch = -1;		// placeholder
  maxnsearch = -1;		// default
  rnormtol = 0.0;		// placeholder
  isfile = false;		// default
  isset = false;		// default
  acopy = NULL;
  bcopy = NULL;
  index = NULL;			// placeholder
  solution = NULL;		// placeholder
  work1 = NULL;			// placeholder
  work2 = NULL;			// placeholder
}

Pattern::~Pattern()
{
  delete [] acopy;
  delete [] bcopy;
  delete [] index;
  delete [] solution;
  delete [] work1;
  delete [] work2;
}

bool Pattern::isSet() const
{
  return (isset);
}

const vector<double>& Pattern::operator[](int i) const
{
  // Get the ith column of the pattern
  if (!p.isSet())
    throw "Must call Pattern::setPattern() before Pattern[]";
  if ((i < 0) || (i >= nsearch))
    throw "Invalid index in Pattern[]";
  return p.column(i);
}

void Pattern::setNdim(int n)
{
  if (ndim != 0)
    throw "Cannot call Pattern::setNdim() more than once";
  if (n <= 0)
    throw "Invalid parameter for Pattern::setNdim()";
  ndim = n;
}

void Pattern::setFile(const char* fname)
{
  if (ndim <= 0)
    throw "Must call Pattern::setNdim() before Pattern::setFileName()";

  if (type != -1)
    throw "Must call Pattern::setFileName() before Pattern::setType()";

  if (strlen(fname) == 0)
    return;

  // Set isfile to true (default is false)
  isfile = true;

  // Open file
  ifstream input;
  input.open(fname);
  if (!input)
    throw "Unable to open pattern file";

  // First entry should be maxnsearch, i.e., the number of search
  // directions defined in the file
  input >> maxnsearch;

  // Second entry should be minnsearch, i.e., the number of search
  // directions that must converge to guarantee convergence.
  input >> minnsearch;

  // Error check these two values
  if ((maxnsearch <= 0) || (minnsearch <= 0) || (minnsearch > maxnsearch))
    throw "Error specifying maxnsearch or minnsearch in pattern file";

  // Set pattern size
  p.setMatrix(ndim, maxnsearch);

  // Read pattern, one COLUMN per row
  for (int j = 0; j < maxnsearch; j ++)
    for (int i = 0; i < ndim; i ++)
      input >> p.entry(i, j);

  // Close input file
  input.close();

  // Delete memory for filename
  delete[] fname;
  fname = NULL;
}

void Pattern::setType(int i, bool bounds)
{
  if (ndim == 0)
    throw "Cannot call Pattern::setType() before Pattern::setNdim()";

  if (type != -1)
    throw "Cannot call Pattern::setType() more than once";

  // If a user-specified pattern was used, everything is already
  // set. No checking.
  if (isfile) {
    type = USER_DEF;
    return;
  }

  // If bounds are being used, then we must use the +/- unit vectors.
  if (bounds && (i != COORD_MAX)) {
    cout << "Usage: only pattern type 0 is valid with bound constraints" << endl;
    exit(1);
  }

  // Check that we have a valid type.
  if ((i != COORD_MIN) && (i != SPHERICAL_MIN) &&
      (i != COORD_MAX) && (i != SPHERICAL_MAX)) {
    cout << "Usage: pattern type must be 0, 1, 2 or 3." << endl;
    exit(1);
  }

  // Set the type.
  type = i;

  // Set minnsearch.
  if ((type == COORD_MIN) || (type == SPHERICAL_MIN))
    minnsearch = ndim + 1;
  else	// ((type == COORD_MAX) || (type == SPHERICAL_MAX))
    minnsearch = 2 * ndim;
}

void Pattern::setNsearch(int n) 
{
  if (nsearch != 0)
    throw "Cannot call Pattern::setNsearch() more than once";
  if (n <= 0)
    throw "Invalid parameter for Pattern::setNsearch()";
  nsearch = n;
}

int Pattern::getNsearch() const
{
  if (nsearch == 0)
    throw "Must call Pattern::getNsearch() before Pattern::getNsearch()";
  return nsearch;
}

int Pattern::getMinNsearch() const
{
  if (type == -1)
    throw "Must call Pattern::setNdim() and Pattern::setType() before Pattern::getMinSize()";
  return minnsearch;
}

const Matrix& Pattern::getPattern() const
{
  if (!isSet())
    throw "Cannot call getPattern() before setPattern()";
  return p;
}

void Pattern::setRnormtol(double d)
{
  if (d < 0) {
    cout << "Usage: NNLS tolerance (" << d << ") must be positive" << endl;
    exit(1);
  }

  rnormtol = d;
}

void Pattern::setPattern()
{
  if (rnormtol == 0)
    throw "Must call Pattern::setRnormtol() before Pattern::setPattern()";
  if (type == -1)
    throw "Must call Pattern::setType() before Pattern::setPattern()";
  if (ndim == 0)
    throw "Must call Pattern::setNdim() before Pattern::setPattern()";
  if (nsearch == 0)
    throw "Must call Pattern::setNsearch() before Pattern::setPattern()";
  if (isset)
    throw "Cannot call Pattern::setPattern() more than once";

  // Counters
  int i, j, k;

  // Number of currently filled columns
  int nfill = 0;
  
  if (isfile) {

    // The pattern was read from a file.

    // Check that nsearch is large enough. If not, reset it. (Recall
    // that minnsearch is read from the pattern file.
    if (nsearch < minnsearch) {
      cerr << APPS::idstr << "WARNING! ";
      cerr << "Resetting nsearch from " << nsearch << " to " << minnsearch << endl;
      nsearch = minnsearch;
    }
    
    // Enlarge/Shrink the pattern 
    if (nsearch != maxnsearch) { 

      // Copy p into temporary matrix
      Matrix tmp = p;

      // Copy the first min(nsearch,maxnsearch) columns of the old p
      // into the new p
      vector<bool> b(maxnsearch);
      for (i = 0; i < maxnsearch; i ++)
	b[i] = (i < nsearch) ? true : false;
      p.copy(tmp, b);

    }      

    // Set the number of filled columns to min(nsearch,maxnsearch)
    nfill = (nsearch < maxnsearch) ? nsearch : maxnsearch;

  }
  else {

    // The pattern was not read from a file, so we must create it.

    // Allocate space and zero out all entries
    p.setMatrix(ndim, nsearch);

    if ((type == COORD_MIN) || (type == COORD_MAX)) {
      
      // Pattern based on coordinate vectors

      if (type == COORD_MIN) 
	for (i = 0; i < ndim; i ++) {
	  p.entry(i, i) = 1;
	  p.entry(i, ndim) = -1;
	}
      else
	for (i = 0; i < ndim; i ++) {
	  p.entry(i, i) = 1;
	  p.entry(i, i + ndim) = -1;
	}
    }
    else {
      
      // Pattern based on regular simplex
      
      double centroid = ndim / sqrt(ndim + 1.0e0);
      double pbar = (sqrt(ndim + 1.0e0) - 1.0e0 + ndim) - centroid;
      double qbar = (sqrt(ndim + 1.0e0) - 1.0e0) - centroid;
      double divisor = ndim * sqrt(1.0e0 * ndim) / sqrt(ndim + 1.0e0);
      
      for (i = 0; i < ndim; i ++) 
	for (j = 0; j < ndim; j ++)
	  if (i == j)
	    p.entry(i, j) = pbar / divisor;
	  else
	    p.entry(i, j) = qbar / divisor;
      
      if (type == SPHERICAL_MIN) 
	for (i = 0; i < ndim; i++) 
	  p.entry(i, ndim) = - centroid / divisor;
      else 
	for (i = 0; i < ndim; i ++) 
	  for (j = 0; j < ndim; j ++)
	    p.entry(i, j + ndim) = - p.entry(i, j);
    }

    // Number of columns filled nsearch based on ndim and the type
    nfill = minnsearch;
  }

  // Fill the remaining vectors (between n and nsearch) with random
  // binary combinations of the previous search directions

  srand(clock());
  j = nfill;
  int cnt;
  while (j < nsearch) {

    for (i = 0; i < ndim; i ++) 
      p.entry(i, j) = 0;
    
    for (k = 0; k < nfill; k ++) 
      if (rand() < (RAND_MAX / 2)) 
	for (i = 0; i < ndim; i ++) 
	  p.entry(i,j) += p.entry(i,k);
    
    cnt = 0;
    for (i = 0; i < ndim; i ++) 
      if (p.entry(i, j) != 0) 
	cnt ++;
    
    if (cnt >= 2)
      j ++;
  }

  // Now the pattern is set
  isset = true;
}

bool Pattern::isPosBasis(const Matrix& q, const Matrix& rhs)
{
  // Return true if the search directions indicated by isIncluded (a
  // boolean array of length nsearch) form a positive spanning set.

  if (!isSet())
    throw "Must call Pattern::setPattern() before Pattern::isPosBasis()";

  // Must have at least one column in q to proceed.
  if (q.getNcols() < 1) 
    return false;

  if (APPS::debug >= 9) 
    q.print(cout, "Projected Pos Basis Test Col ");

  // Solve a NNLS problem for the n+1 rhs's. The residual should be
  // zero (i.e., less than RNORMTOL) for every problem.

  for (int j = 0; j < rhs.getNcols(); j ++) {

    double rnorm = Pattern::nnls(q, rhs.column(j));

    if (APPS::debug >= 9) {
      cout.precision(0);
      cout << APPS::idstr << "Solution of NNLS problem " << j + 1 << " of " << rhs.getNcols() 
	   << " with rhs = [ ";
      for (int i = 0; i < rhs.getNrows(); i ++)
	cout << setw(2) << rhs.entry(i,j) << " ";
      cout.setf(ios::scientific);
      cout.precision(APPS::precision);
      cout << "] is rnorm = " << setw(APPS::precision + 6) << rnorm << endl;
      cout.unsetf(ios::scientific);
    }

    if (rnorm > rnormtol)
      return false;
  }

  return true;
}

void Pattern::pack() const
{
  if (!isSet())
    throw "Must call Pattern::setPattern() before Pattern::pack()";

  GCI::pack(ndim);
  GCI::pack(nsearch);
  GCI::pack(type);
  GCI::pack(rnormtol);
  p.pack();
}

void Pattern::unpack()
{
  GCI::unpack(ndim);
  GCI::unpack(nsearch);
  GCI::unpack(type);
  GCI::unpack(rnormtol);
  p.unpack();
}

void Pattern::print(ostream& stream) const
{
  cout << APPS::idstr << "Pattern ndim=" << ndim << " nsearch=" << nsearch
       << " type=" << type << " rnormtol=" << rnormtol << "\n";
  p.print(stream, "Pattern column ");
}


double Pattern::nnls(const Matrix& A, const vector<double>& b)
{
  // Solve the NNLS problem min { |Ax - b| s.t. x >= 0 } and return
  // the residual.

  int nrows = A.getNrows();
  int ncols = A.getNcols();

  if (b.size() != nrows)
    throw "Size mismatch in Pattern::nnls()";

  // Create work vectors, if they have not already been created.
  if (acopy == NULL) {
    int maxncols = p.getNcols();
    acopy = new double[nrows * maxncols];
    bcopy = new double[nrows];
    index = new int[maxncols];
    solution = new double[maxncols];
    work1 = new double[maxncols];
    work2 = new double[maxncols];
  }

  double rnorm;
  int mode;

  int i,j;

  // Copy A into acopy
  for (i = 0; i < nrows; i ++)
    for (j = 0; j < ncols; j ++) 
      acopy[i + (j * nrows)] = A.entry(i,j);

  // Copy b into bcopy
  for (i = 0; i < nrows; i ++)
    bcopy[i] = b[i];

  // Solve the NNLS problem
  nnls_c(acopy, &nrows, &nrows, &ncols, bcopy, solution, &rnorm, work1, work2, index, &mode);
  
  if (mode > 1)
    cout << APPS::idstr << "Warning: NNLS may have given inaccurate result." << endl;

  return rnorm;
}