densematrix.cc

一种聚类算法,名字是cocluster

    double norm = 0.0;
    for (int j = 0; j < numRow; j++){
      double tempValue = value[j][i];
      norm += tempValue * tempValue;
    }
    if (norm > 0.0){
      norm = sqrt(norm);
      for (int j = 0; j < numRow; j++)
        value[j][i] /= norm;
    }
  }
}


void DenseMatrix::normalize_mat_L1()
  /* L1 normalize each vector in the dense matrix to have L1 norm 1
   */
{
  for (int i = 0; i < numCol; i++){
    double norm = 0.0;
    for (int j = 0; j < numRow; j++)
      norm += fabs(value[j][i]);
    if (norm > 0)
      for (int j = 0; j < numRow; j++)
        value[j][i] /= norm;
  }
}


void DenseMatrix::ith_add_CV(int i, double *CV)
{
  for (int j = 0; j < numRow; j++)
    CV[j] += value[j][i];
}


void DenseMatrix::ith_add_CV_prior(int i, double *CV)
{
  for (int j = 0; j < numRow; j++)
    CV[j] += priors[i] * value[j][i];
}


void DenseMatrix::CV_sub_ith(int i, double *CV)
{
  for (int j = 0; j < numRow; j++)
    CV[j] -= value[j][i];
}


void DenseMatrix::CV_sub_ith_prior(int i, double *CV)
{
  for (int j = 0; j < numRow; j++)
    CV[j] -= priors[i] * value[j][i];
}


double DenseMatrix::computeMutualInfo()
{
  double *rowSum = new double[numRow], MI = 0.0;
  for (int i = 0; i < numRow; i++){
    rowSum[i] = 0.0;
    for (int j = 0; j < numCol; j++)
      rowSum[i] += value[i][j] * priors[i];
  }
  for (int i = 0; i < numRow; i++){
    double temp = 0;
    for (int j = 0; j < numCol; j++){
      double tempValue = value[i][j];
      if (tempValue > 0.0)
        temp += tempValue * log(tempValue / rowSum[i]);
    }
    MI += temp * priors[i];
  }
  delete [] rowSum;
  return MI / log(2.0);
}


double DenseMatrix::exponential_kernel(double *x, int i, double norm_x, double sigma_squared)
  // this function computes the exponential kernel distance between i_th data with the centroid v
{
  double result = 0.0;
  for (int j = 0; j < numRow; j++)
    result += x[j] * value[j][i];
  result *= -2.0;
  result += norm[i] + norm_x;
  return exp(result * 0.5 / sigma_squared);
}


void DenseMatrix::exponential_kernel(double *x, double norm_x, double *result, double sigma_squared)
  //this function computes the exponential kernel distance between all data with the centroid x
{
  for (int i = 0; i < numCol; i++)
    result[i] = exponential_kernel(x, i, norm_x, sigma_squared);
}


double DenseMatrix::i_j_dot_product(int i, int j)
//this function computes  dot product between vectors i and j
{
  double result = 0;
  if (i == j)
    for (int l = 0; l < numRow; l++){
      double tempValue = value[l][i];
      result += tempValue * tempValue;
    }
  else
    for (int l = 0; l < numRow; l++)
      result += value[l][i] * value[l][j];
  return result;
}


double DenseMatrix::squared_i_j_euc_dis(int i, int j)
  //computes squared Euc_dis between vec i and vec j
{
  return (norm[i] - 2.0 * i_j_dot_product(i,j) + norm[j]);
}


void  DenseMatrix::pearson_normalize()
//this function shift each vector so that it has 0 mean of all its entries
  //then the vector gets L2 normalized.
{
  for (int i = 0; i < numCol; i++){
    double mean = 0;
    for (int j = 0; j < numRow; j++)
      mean += value[j][i];
    mean /= numRow;
    for (int j = 0; j < numRow; j++)
      value[j][i] -= mean;
  }
/*
  for (int i = 0; i < numCol; i++) {
    for (int j = 0; j < numRow; j++)
      cout << value[j][i] << " ";
    cout << endl;
  }
*/
  normalize_mat_L2();
/*
  for (int i = 0; i < numCol; i++){
    for (j = 0; j < numRow; j++)
      cout << value[j][i] << " ";
    cout << endl;
  }
*/
}


bool DenseMatrix::isHavingNegative()
{
  bool havingNegative = false;
  for (int r = 0; r < numRow && !havingNegative; r++)
    for (int c = 0; c < numCol && !havingNegative; c++)
      if (value[r][c] < 0)
        havingNegative = true;
  return havingNegative;
}


double DenseMatrix::getPlogQ(double **pxhatyhat, int *rowCL, int *colCL, double *pXhat, double *pYhat)
{
  double PlogQ = 0;
  for (int i = 0; i < numRow; i++)
    for (int j = 0; j < numCol; j++){
      double tempValue = value[i][j];
      if (tempValue > 0.0){
        int tempRowCLI = rowCL[i];
	int tempColCLJ = colCL[j];
        PlogQ += tempValue * log(pxhatyhat[tempRowCLI][tempColCLJ] * pX[i] * pY[j] / (pXhat[tempRowCLI] * pYhat[tempColCLJ]));
      }
    }
  return PlogQ / log(2.0);
}


void DenseMatrix::preprocess()
{
  double totalSum = 0;
  PlogP = mutualInfo = 0;
  pX = new double[numRow];
  pY = new double[numCol];
  for (int i = 0; i < numRow; i++)
    pX[i] = 0;
  for (int j = 0; j < numCol; j++)
    pY[j] = 0;
  for (int i = 0; i < numRow; i++)
    for (int j = 0; j < numCol; j++){
      double tempValue = value[i][j];
      totalSum += tempValue;
      pX[i] += tempValue;
      pY[j] += tempValue;
    }
  for (int i = 0; i < numRow; i++)
    for (int j = 0; j < numCol; j++)
      value[i][j] /= totalSum;
  for (int i = 0; i < numRow; i++)
    pX[i] /= totalSum;
  for (int j = 0; j < numCol; j++)
    pY[j] /= totalSum;
  xnorm = new double[numRow];
  ynorm = new double[numCol];
  for (int i = 0; i < numRow; i++)
    xnorm[i] = 0;
  for (int j = 0; j < numCol; j++)
    ynorm[j] = 0;
  for (int i = 0; i < numRow; i++){
    for (int j = 0; j < numCol; j++){
      double tempValue = value[i][j];
      if(tempValue > 0){
        double temp = tempValue * log(tempValue);
        xnorm[i] += temp;
        ynorm[j] += temp;
        mutualInfo += tempValue * log(pX[i] * pY[j]);
      }
    }
    PlogP += xnorm[i];
  }
  mutualInfo = PlogP - mutualInfo;
}


void DenseMatrix::condenseMatrix(int *rowCL, int *colCL, int numRowCluster, int numColCluster, double **cM)
{
  for (int i = 0; i < numRowCluster; i++)
    for (int j = 0; j < numColCluster; j++)
      cM[i][j] = 0;
  for (int i = 0; i < numRow; i++)
    for (int j = 0; j < numCol; j++)
      cM[rowCL[i]][colCL[j]] += value[i][j];
}


void DenseMatrix::condenseMatrix(int *rowCL, int *colCL, int numRowCluster, int numColCluster, double **cM, bool *isReversed)
{
  int rowId = 0;
  for (int i = 0; i < numRowCluster; i++)
    for (int j = 0; j < numColCluster; j++)
      cM[i][j] = 0;
  for (int i = 0; i < numRow; i++){
    rowId = rowCL[i];
    if (isReversed[i])
      for (int j = 0; j < numCol; j++)
        cM[rowId][colCL[j]] += (0 - value[i][j]);
    else
      for (int j = 0; j < numCol; j++)
        cM[rowId][colCL[j]] += value[i][j];
  }    
}


double DenseMatrix::Kullback_leibler(double *x, int i, int priorType, int clusterDimension)
{
  double result = 0.0, uni_alpha;
  if (clusterDimension == COL_DIMENSION){
    switch (priorType){
      case NO_SMOOTHING:
        for (int j = 0; j < numRow; j++){
          if(x[j] > 0)
            result += value[j][i] * log(x[j]);
          else if (value[j][i] >0)
            return MY_DBL_MAX; // if KL(vec[i], x) is inf. give it a huge number 1.0e20
        }
        result = (ynorm[i] - result)/ pY[i] - log(pY[i]);
	break;
      case UNIFORM_SMOOTHING:
        uni_alpha = smoothingFactor / numRow;
        for (int j = 0; j < numRow; j++)
          result += value[j][i] * log(x[j] + uni_alpha) ;
        result = (ynorm[i] - result) / pY[i] - log(pY[i]) + log(1 + smoothingFactor);
        break;
      case MAXIMUM_ENTROPY_SMOOTHING:
        for (int j = 0; j < numRow; j++)
          result += value[j][i] * log(x[j] + pX[j] * smoothingFactor);
        result = (ynorm[i] - result)/ pY[i] - log(pY[i]) + log(1 + smoothingFactor);
        break;
      default:
        break;
    }
  } else {
    switch (priorType){
      case NO_SMOOTHING:
        for (int j = 0; j < numCol; j++)
          if (value[i][j] > 0){
            if(x[j] > 0)
              result += value[i][j] * log(x[j]);
            else 
              return MY_DBL_MAX; // if KL(vec[i], x) is inf. give it DBL_MAX.
	  }
        result = (xnorm[i] - result) / pX[i] - log(pX[i]);
        break;
      case UNIFORM_SMOOTHING:
        uni_alpha = smoothingFactor / numCol;
        for (int j = 0; j < numCol; j++)
          result += value[i][j] * log(x[j] + uni_alpha) ;
	result = (xnorm[i] - result) / pX[i] - log(pX[i]) + log(1 + smoothingFactor);
        break;
      case MAXIMUM_ENTROPY_SMOOTHING:
        for (int j = 0; j < numCol; j++)
          result += value[i][j] * log(x[j] + pY[j] * smoothingFactor) ;
        result = (xnorm[i] - result) / pX[i] - log(pX[i]) + log(1 + smoothingFactor);
	break;
      default:
        break;
    }
  }
//cout << result << endl;
  return result;
}


void DenseMatrix::subtractRow(double **x, int row, int i, int *colCL)
{
  for (int j = 0; j < numCol; j++)
    x[row][colCL[j]] -= value[i][j];
}


void DenseMatrix::addRow(double **x, int row, int i, int *colCL)
{
  for (int j = 0; j < numCol; j++)
    x[row][colCL[j]] += value[i][j];
}


void DenseMatrix::subtractCol(double **x, int col, int j, int *rowCL)
{
  for (int i = 0; i < numRow; i++)
    x[rowCL[i]][col] -= value[i][j];
}


void DenseMatrix::addCol(double **x, int col, int j, int *rowCL)
{
  for (int i = 0; i < numRow; i++)
    x[rowCL[i]][col] += value[i][j];
}


void DenseMatrix::subtractRow(double *x, int i, int *colCL)
{
  for (int j = 0; j < numCol; j++)
    x[colCL[j]] -= value[i][j];
}


void DenseMatrix::addRow(double *x, int i, int *colCL)
{
  for (int j = 0; j < numCol; j++)
    x[colCL[j]] += value[i][j];
}


void DenseMatrix::subtractCol(double *x, int j, int *rowCL)
{
  for (int i = 0; i < numRow; i++)
    x[rowCL[i]] -= value[i][j];
}


void DenseMatrix::addCol(double *x, int j, int *rowCL)
{
  for (int i = 0; i < numRow; i++)
    x[rowCL[i]] += value[i][j];
}


/*
// Does y = A * x, A is mxn, and y & x are mx1 and nx1 vectors
//   respectively  
void DenseMatrix::dmatvec(int m, int n, double **a, double *x, double *y)
{
  double yi;
  for(int i = 0;i < m; i++){
    yi = 0.0;
    for(int j = 0;j < n; j++)
      yi += a[i][j] * x[j];
    y[i] = yi;
  }
}


// Does y = A' * x, A is mxn, and y & x are nx1 and mx1 vectors
//   respectively  
void DenseMatrix::dmatvecat(int m, int n, double **a, double *x, double *y)
{
  double yi;
  for(int i = 0;i < n; i++){
    yi = 0.0;
    for(int j = 0;j < m; j++)
      yi += a[j][i] * x[j];
    y[i] = yi;
  }
}


// The procedure dqrbasis outputs an orthogonal basis spanned by the rows
//   of matrix a (using the QR Factorization of a ).  
void DenseMatrix::dqrbasis(double **q)
{
  double *work = new double[numRow];
  for(int i = 0; i < numRow;i++){
    dmatvec(i, numCol, q, value[i], work);
    dmatvecat(i, numCol, q, work, q[i]);
    for(int j = 0; j < numCol; j++)
      q[i][j] = value[i][j] - q[i][j];
    dvec_l2normalize(numCol, q[i]);
  }
}


// The function dvec_l2normsq computes the square of the
//  Euclidean length (2-norm) of the double precision vector v 
double DenseMatrix::dvec_l2normsq( int dim, double *v )
{
  double length = 0.0, tmp;
  for(int i = 0; i < dim; i++){
    tmp = *v++;
    length += tmp*tmp;
  }
  return length;
}


// The function dvec_l2normalize normalizes the double precision
//   vector v to have 2-norm equal to 1 
void DenseMatrix::dvec_l2normalize( int dim, double *v )
{
  double nrm;