sparsematrix.cc

一种聚类算法,名字是cocluster

  /* compute the squared L-2 norms for each vec in the matrix
     first check if array 'norm' has been given memory space
   */
{
  if (L1_norm == NULL){
     L1_norm = new double[numCol];
     memoryUsed += numCol * sizeof(double);
  }
  for (int i = 0; i < numCol; i++){
    L1_norm[i] =0.0;
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      L1_norm[i] += value[j] ;
  }
}


void SparseMatrix::computeNorm_KL(int laplace)
  // the norm[i] is in unit of nats NOT bits
{
  double row_inv_alpha = smoothingFactor / numRow;
  if (norm == NULL){
    norm = new double[numCol];
    memoryUsed += numCol * sizeof(double);
  }
  Norm_sum = 0;
  switch (laplace){
    case NO_SMOOTHING:
    case UNIFORM_SMOOTHING:
      for (int i = 0; i < numCol; i++){
        norm[i] = 0.0;
        for (int j = colPtr[i]; j < colPtr[i+1]; j++){
	  double tempValue = value[j];
          norm[i] += tempValue * log(tempValue);
	}
        Norm_sum += norm[i] * priors[i];
      }
      break;
    case MAXIMUM_ENTROPY_SMOOTHING:
      for (int i = 0; i < numCol; i++){
        norm[i] = 0.0;
        for (int j = colPtr[i]; j < colPtr[i+1]; j++){
	  double tempValue = value[j];
          norm[i] += tempValue * log(tempValue);
	}
        Norm_sum +=norm[i] * priors[i];
      }
      if (p_x == NULL){
        p_x = new double[numRow];
	memoryUsed += numRow * sizeof(double);
      }
      for (int i = 0; i < numRow; i++)
        p_x[i] = 0;
      for (int i = 0; i < numCol; i++)
        for (int j = colPtr[i]; j < colPtr[i+1]; j++)
          p_x[rowIdx[j]] += priors[i] * value[j];
      break;
    case LAPLACE_SMOOTHING:
      for (int i = 0; i < numCol; i++){
        norm[i] =0.0;
        for (int j = colPtr[i]; j < colPtr[i+1]; j++){
	  double tempValue = value[j];
          norm[i] += (tempValue + row_inv_alpha) * log(tempValue + row_inv_alpha);
	}
        norm[i] += (numRow - (colPtr[i+1] - colPtr[i])) * row_inv_alpha * log(row_inv_alpha) ;
        norm[i] = norm[i] / (1+smoothingFactor) + log(1+smoothingFactor);
        Norm_sum += norm[i] * priors[i];
      }
    default:
      break;
  }
  Norm_sum /= log(2.0);
}


void SparseMatrix::normalize_mat_L2()
  /* compute the L_2 norms for each vec in the matrix and L_2-normalize them
     first check if array 'norm' has been given memory space
   */
{
  for (int i = 0; i < numCol; i++){
    double norm = 0.0;
    for (int j = colPtr[i]; j < colPtr[i+1]; j++){
      double tempValue = value[j];
      norm += tempValue * tempValue;
    }
    if( norm > 0.0){
      norm = sqrt(norm);
      for (int j = colPtr[i]; j < colPtr[i+1]; j++)
        value[j] /= norm;
    }
  }
}


void SparseMatrix::normalize_mat_L1()
  /* compute the L_1 norms for each vec in the matrix and L_1-normalize them
     first check if array 'L1_norm' has been given memory space
   */
{
  for (int i = 0; i < numCol; i++){
    double norm = 0.0;
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      norm += fabs(value[j]);
    if(norm > 0)
      for (int j = colPtr[i]; j < colPtr[i+1]; j++)
         value[j] /= norm;
  }
}


void SparseMatrix::ith_add_CV(int i, double *CV)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    CV[rowIdx[j]] += value[j];
}


void SparseMatrix::ith_add_CV_prior(int i, double *CV)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    CV[rowIdx[j]] += priors[i] * value[j];
}


void SparseMatrix::CV_sub_ith(int i, double *CV)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    CV[rowIdx[j]] -= value[j];
}


void SparseMatrix::CV_sub_ith_prior(int i, double *CV)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    CV[rowIdx[j]] -= priors[i]*value[j];
}


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


double SparseMatrix::exponential_kernel(double *v, int i, double norm_v, 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 = colPtr[i]; j < colPtr[i+1]; j++)
    result += value[j] * v[rowIdx[j]];
  result *= -2.0;
  result += norm[i] + norm_v;
  return exp(result * 0.5 / sigma_squared);
}


void SparseMatrix::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 SparseMatrix::i_j_dot_product(int i, int j)
//this function computes  dot product between vectors i and j
{
  double result = 0.0;
  if (i == j)
    for (int l = colPtr[i]; l < colPtr[i+1]; l++){
      double tempValue = value[l];
      result += tempValue * tempValue;
    }
  else
    for (int l = colPtr[i]; l < colPtr[i+1]; l++)
      for (int k = colPtr[j]; k < colPtr[j+1]; k++)
	if(rowIdx[l] == rowIdx[k])
	  result += value[l] * value[k];
  return result;
}


double SparseMatrix::squared_i_j_euc_dis(int i, int j)
{
  return (norm[i] - 2.0 * i_j_dot_product(i,j) + norm[j]);
}


void SparseMatrix::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.0;
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      mean += value[j];
    mean /= numRow;
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      value[j] -= mean;
  }
  normalize_mat_L2();
}


bool SparseMatrix::isHavingNegative()
{
  bool havingNegative = false;
  for (int c = 0; c < numCol && !havingNegative; c++)
    for (int k = colPtr[c]; k < colPtr[c+1] && !havingNegative; k++)
      if (value[k] < 0)
        havingNegative = true;
  return havingNegative;
}


double SparseMatrix::getPlogQ(double **pxhatyhat, int *rowCL, int *colCL, double *pXhat, double *pYhat)
{

  double PlogQ = 0;
  for (int i = 0; i < numCol; i++){
    int colID = colCL[i];
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      PlogQ += value[j] * log(pxhatyhat[rowCL[rowIdx[j]]][colID] * pX[rowIdx[j]]* pY[i] / (pXhat[rowCL[rowIdx[j]]] * pYhat[colID]));
  }
  return PlogQ / log(2.0);
}


void SparseMatrix::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;
  //cout<<numRow<<" "<<numCol<<endl;
  for (int i = 0; i < numCol; i++)
    for (int j = colPtr[i]; j < colPtr[i+1]; j++){
      totalSum += value[j];
      pY[i] += value[j];
      pX[rowIdx[j]] += value[j];
    }
  //cout<<"ts="<<totalSum<<endl;
  for (int i = 0; i < numCol; i++)
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      value[j] /= totalSum;
  for (int i = 0; i < numRow; i++)
    pX[i] /= totalSum;
  for (int j = 0; j < numCol; j++)
    pY[j] /= totalSum;
  norm = new double[numCol];
  for (int i = 0; i < numCol; i++){
    norm[i] = 0;
    for (int j = colPtr[i]; j < colPtr[i+1]; j++){
      double tempValue = value[j];
      norm[i] += tempValue*log(tempValue);
      mutualInfo += tempValue*log(pX[rowIdx[j]] * pY[i]);
    }
    PlogP += norm[i];
  }
  mutualInfo = PlogP - mutualInfo;
}
	  	  

void SparseMatrix::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 < numCol; i++)
    for (int j = colPtr[i]; j < colPtr[i+1]; j++)
      cM[rowCL[rowIdx[j]]][colCL[i]] += value[j]; 
}


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


double SparseMatrix::Kullback_leibler(double *x, int i, int priorType, int clusterDimension)
// In fact, clusteringDimension is not used in the body...
{
  double result = 0.0, uni_alpha = smoothingFactor / numRow;
  switch (priorType){
    case NO_SMOOTHING:
      for (int j = colPtr[i]; j < colPtr[i+1]; j++){
        double tempValue = x[rowIdx[j]];
        if (tempValue > 0.0)
          result += value[j] * log(tempValue);
        else
          return MY_DBL_MAX; // if KL(vec[i], x) is inf. give it DBL_MAX
      }
      result = (norm[i] - result) / pY[i] - log(pY[i]);
      break;
    case UNIFORM_SMOOTHING:
      for (int j = colPtr[i]; j < colPtr[i+1]; j++)
        result += value[j] * log(x[rowIdx[j]] + uni_alpha);
      result = (norm[i] - result) / pY[i] - log(pY[i]) + log(1+smoothingFactor);
      break;
    case MAXIMUM_ENTROPY_SMOOTHING:
      for (int j = colPtr[i]; j < colPtr[i+1]; j++){
        int tempRowId = rowIdx[j];
        result += value[j] * log(x[tempRowId] + pX[tempRowId] * smoothingFactor) ;
      }
      result = (norm[i] - result)/ pY[i] - log(pY[i]) + log(1+smoothingFactor);
      break;
    default:
      break;
  }
//cout << "SMOOTHING_FACTOR = " << smoothingFactor << endl;  
  return result;
}


void SparseMatrix::subtractRow(double **x, int row, int i, int *colCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[row][colCL[rowIdx[j]]] -= value[j];
}


void SparseMatrix::subtractRow(double *x, int i, int *colCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[colCL[rowIdx[j]]] -= value[j];
}


void SparseMatrix::addRow(double **x, int row, int i, int *colCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[row][colCL[rowIdx[j]]] += value[j];
}


void SparseMatrix::addRow(double *x, int i, int *colCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[colCL[rowIdx[j]]] += value[j];
}


void SparseMatrix::subtractCol(double **x, int col, int i, int *rowCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[rowCL[rowIdx[j]]][col] -= value[j];
}


void SparseMatrix::subtractCol(double *x, int i, int *rowCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[rowCL[rowIdx[j]]] -= value[j];
}


void SparseMatrix::addCol(double **x, int col, int i, int *rowCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[rowCL[rowIdx[j]]][col] += value[j];
}


void SparseMatrix::addCol(double *x, int i, int *rowCL)
{
  for (int j = colPtr[i]; j < colPtr[i+1]; j++)
    x[rowCL[rowIdx[j]]] += value[j];
}


double SparseMatrix::squaredFNorm()
  /* compute the squared Frobenius norm for the matrix.
   */
{
  double temp = 0;
  for (int c = 0; c < numCol; c++)
    for (int k = colPtr[c]; k < colPtr[c+1]; k++){
      double tempValue = value[k];
      temp += tempValue * tempValue;
    }
  return temp;
}


double SparseMatrix::squaredL2Norm4Row(int r)
  /* compute the squared L-2 norm for row vec i in the matrix.
   */
{
  double temp = 0;
  for (int k = colPtr[r]; k < colPtr[r+1]; k++){
    double tempValue = value[k];
    temp += tempValue * tempValue;
  }
  return temp;
}


double SparseMatrix::squaredL2Norm4Col(int c)
  /* compute the squared L-2 norm for column vec j in the matrix.
   */
{
  double temp = 0;
  for (int k = colPtr[c]; k < colPtr[c+1]; k++){
    double tempValue = value[k];
    temp += tempValue * tempValue;
  }
  return temp;
}


double SparseMatrix::computeRowDistance(int rowId, int rowCluster, int *colCL, double **cM, double rowQuality4Compressed)
{
  double temp = 0;
  for (int k = colPtr[rowId]; k < colPtr[rowId+1]; k++)
    temp += value[k] * cM[rowCluster][colCL[rowIdx[k]]];
  return (-2 * temp + rowQuality4Compressed);
}