gbppreprocessor.cpp

The package includes 3 Matlab-interfaces to the c-code: 1. inference.m An interface to the full

  }
  
  delete[] pp_doubleCount;
  pp_doubleCount = 0;
  delete[] pp_parents;
  pp_parents = 0;
}

void GBPPreProcessor::calcNumOfStates() {
  int const* V = pp_mrf->V;

  for (int i=0; i<pp_numAllRegs; i++) {
    Region* region = pp_allRegionsPtr[i];

    // calc i's number of possible states
    int i_numStates = 1;
    
    nodes_iterator reg_niter = region->begin();
    while (reg_niter != region->end()) {
      int& node = (*reg_niter);
      i_numStates *= V[node];

      reg_niter++;
    }

    pp_regions_mrf->V[i] = i_numStates;

  }
}

void GBPPreProcessor::calcAssignTable() {

  int const& N = pp_mrf->N;
  int const* V = pp_mrf->V;
  
  pp_assignTable = new int**[pp_numAllRegs];
  for (int i=0; i<pp_numAllRegs; i++) {
    pp_assignTable[i] = new int*[pp_regions_mrf->neighbNum(i)];

    Region* i_region = pp_allRegionsPtr[i];
    int i_numStates = pp_regions_mrf->V[i];
    
    Assignment assignment;
    assignment.clear();
    assignment.resize(N);

    for (int n=0; n<pp_regions_mrf->neighbNum(i); n++) {
      
      int j = (*pp_adjMat)[i][n];
      
      pp_assignTable[i][n] = 0;
      
      if (i<j) {

	Region* j_region = pp_allRegionsPtr[j];
	
	pp_assignTable[i][n] = new int[i_numStates];
	for (int i_assignInd=0; i_assignInd<i_numStates; i_assignInd++) {
	  i_region->indToAssign(i_assignInd, assignment, V);
	  int j_assignInd = j_region->assignToInd(assignment, V);
	  pp_assignTable[i][n][i_assignInd] = j_assignInd;
	}
      }
    }
    
  }
}

void GBPPreProcessor::calcPotentials(Potentials* bigRegsPot) {

  int* regCard = pp_regions_mrf->V;

  pp_regions_mrf->initLocalPotentials();
  Potentials* regPot = pp_regions_mrf->localMat;
  if (pp_mrf->logspace) {
    for (int i=0; i<pp_numAllRegs; i++) {
      for (int xi=0; xi<regCard[i]; xi++) {
	regPot[i][xi] = 0.0;
      }
    }
  }
  else {
    for (int i=0; i<pp_numAllRegs; i++) {
      for (int xi=0; xi<regCard[i]; xi++) {
	regPot[i][xi] = 1.0;
      }
    }
  }
  int const& N = pp_mrf->N;
  int const* V = pp_mrf->V;
  
  pp_extractSingle = new int[N];
  pp_extractPairs = new int*[N];
  for (int i=0; i<N; i++) {
    pp_extractSingle[i] = -1;
    pp_extractPairs[i] = new int[pp_mrf->neighbNum(i)];
    for (int n=0; n<pp_mrf->neighbNum(i); n++) {
      pp_extractPairs[i][n] = -1;
    }
  }

  int numLocalUsed = 0;
  bool* recievedPot = new bool[pp_numBigRegs];
  for (int i=0; i<pp_numBigRegs; i++) {
    recievedPot[i] = false;
  }  

  Assignment assignment;
  assignment.clear();
  assignment.resize(N);

  // deal the local potentials of the nodes between the big regions  

  for (int i=0; i<pp_numBigRegs; i++) {
    
    Region* region = pp_allRegionsPtr[i];
    const_nodes_iterator niter = region->begin();
    while (niter != region->end()) {
      int n = *niter;
      if (pp_extractSingle[n]==-1) {
	if (pp_mrf->isLocal()) {
	  for (int i_assignInd=0; i_assignInd<regCard[i]; i_assignInd++) {
	    region->indToAssign(i_assignInd, assignment, V);
	    int xn = assignment[n];
	    if (pp_mrf->logspace) {
	      // 	    cout << "i = " << i << " n = " << n << " exp(-pp_mrf->localMat[n][xn]) = " << exp(-pp_mrf->localMat[n][xn]) << endl;
	      regPot[i][i_assignInd] += pp_mrf->localMat[n][xn];
	    }
	    else {
	      // 	    cout << "i = " << i << " n = " << n << " pp_mrf->localMat[n][xn] = " << pp_mrf->localMat[n][xn] << endl;
	      regPot[i][i_assignInd] *= pp_mrf->localMat[n][xn];
	    }
	  }
	}
	numLocalUsed++;

	recievedPot[i] = true;
	pp_extractSingle[n] = i;
      }
      niter++;
    }

    if (numLocalUsed == N) {
      break;
    }
  }
  // deal the pairwise potentials of the nodes between the big regions
  for (int i=0; i<N; i++) {
    for (int n=0; n<pp_mrf->neighbNum(i); n++) {
      int j = pp_mrf->adjMat[i][n];
      if (i<j) {
	for (int reg=0; reg<pp_numBigRegs; reg++) {
	  Region* region = pp_allRegionsPtr[reg];

	  if (region->contains(i) && region->contains(j)) {
	    if (pp_mrf->isPairwise()) {
	      for (int reg_assignInd=0; reg_assignInd<regCard[reg]; reg_assignInd++) {
		region->indToAssign(reg_assignInd, assignment, V);
		int xi = assignment[i];
		int xj = assignment[j];
		if (pp_mrf->logspace) {
		  regPot[reg][reg_assignInd] += pp_mrf->pairPotential(i,n,xi,xj);
		}
		else {
		  regPot[reg][reg_assignInd] *= pp_mrf->pairPotential(i,n,xi,xj);
		}
	      }
	    }
	    
	    recievedPot[reg] = true;
	    pp_extractPairs[i][n] = reg;
	    break;
	  }
	}
      }
    }
  }
  // deal the regions potentials, if given
  if (bigRegsPot != 0) {
    for (int reg=0; reg<pp_numBigRegs; reg++) {
      for (int reg_assignInd=0; reg_assignInd<regCard[reg]; reg_assignInd++) {
	if (pp_mrf->logspace) {
	  regPot[reg][reg_assignInd] += bigRegsPot[reg][reg_assignInd];
	}
	else {
	  regPot[reg][reg_assignInd] *= bigRegsPot[reg][reg_assignInd];
	}
	recievedPot[reg] = true;
      }
    }
  }
  
  // check wether exists a big region which is redundant, and issue an
  // error message

  for (int i=0; i<pp_numBigRegs; i++) {
    if (!recievedPot[i]) {
      //      cerr << "region number " << i << " is redundant: " << (*pp_allRegionsPtr[i]);
      mexPrintf("Warning: region number %d is redundant: ",i);
      Region* region = pp_allRegionsPtr[i];
      const_nodes_iterator niter = region->begin();
      while (niter != region->end()) {
	mexPrintf(" %d",*niter);
	niter++;
      }
      mexPrintf("\n");
    }
  }
  
  delete[] recievedPot;
  recievedPot = 0;

}

void GBPPreProcessor::extractSingle(double** regsBeliefs,
				    double** singleBeliefs,
				    SumOrMax smFlag) const {

  if (pp_mrf->logspace) {
    extractSingleLog(regsBeliefs,singleBeliefs,smFlag);
  }
  else {
    int const& N = pp_mrf->N;
    int const* V = pp_mrf->V;
    int* regCard = pp_regions_mrf->V;
  
    Assignment assignment;
    assignment.clear();
    assignment.resize(N);

    for (int n=0; n<N; n++) {
      int reg = pp_extractSingle[n];
      if (reg<0) {
	for (int xn=0; xn<V[n]; xn++) {
	  singleBeliefs[n][xn] = 1.0 / V[n];
	}
      }
      else {
	Region* region = pp_allRegionsPtr[reg];
	switch (smFlag) {
	  case SUM:
	    for (int xn=0; xn<V[n]; xn++) {
	      singleBeliefs[n][xn] = 0.0;
	    }
	    break;
	  case MAX:
	    for (int xn=0; xn<V[n]; xn++) {
	      singleBeliefs[n][xn] = -1.0;
	    }
	    break;
	  default:
	    break;
	}

	for (int reg_assignInd=0; reg_assignInd<regCard[reg]; reg_assignInd++) {
	  region->indToAssign(reg_assignInd, assignment, V);
	  int xn = assignment[n];
	  switch (smFlag) {
	    case SUM:
	      singleBeliefs[n][xn] += regsBeliefs[reg][reg_assignInd];
	      break;
	    case MAX:
	      if (regsBeliefs[reg][reg_assignInd] > singleBeliefs[n][xn]) {
		singleBeliefs[n][xn] = regsBeliefs[reg][reg_assignInd];
	      }
	      break;
	    default:
	      break;
	  }
	}
	// normalize for max-gbp
	if (smFlag == MAX) {
	  double sum_bel_n = 0.0;
	  for (int xn=0; xn<V[n]; xn++) {
	    sum_bel_n += singleBeliefs[n][xn];      
	  }
	  for (int xn=0; xn<V[n]; xn++) {
	    if (sum_bel_n > 0.0) {
	      singleBeliefs[n][xn] /= sum_bel_n;
	    }
	  }    
	}
      }
    }
  }
}

void GBPPreProcessor::extractSingleLog(double** regsBeliefs,
				       double** singleBeliefs,
				       SumOrMax smFlag) const {

  int const& N = pp_mrf->N;
  int const* V = pp_mrf->V;
  int* regCard = pp_regions_mrf->V;
  
  Assignment assignment;
  assignment.clear();
  assignment.resize(N);


  for (int n=0; n<N; n++) {
    int reg = pp_extractSingle[n];
    if (reg<0) {
      for (int xn=0; xn<V[n]; xn++) {
	singleBeliefs[n][xn] = 0.0;
      }
    }
    else{
      Region* region = pp_allRegionsPtr[reg];
      switch (smFlag) {
	case SUM:
	  for (int xn=0; xn<V[n]; xn++) {
	    singleBeliefs[n][xn] = -HUGE_VAL;
	  }
	  break;
	case MAX:
	  for (int xn=0; xn<V[n]; xn++) {
	    singleBeliefs[n][xn] = HUGE_VAL;
	  }
	  break;
	default:
	  break;
      }

      for (int reg_assignInd=0; reg_assignInd<regCard[reg]; reg_assignInd++) {
	region->indToAssign(reg_assignInd, assignment, V);
	int xn = assignment[n];
	switch (smFlag) {
	  case SUM:
	    singleBeliefs[n][xn] = AddLogFactor(singleBeliefs[n][xn],regsBeliefs[reg][reg_assignInd],-pp_mrf->getTemperature());
	    break;
	  case MAX:
	    if (regsBeliefs[reg][reg_assignInd] < singleBeliefs[n][xn]) {
	      singleBeliefs[n][xn] = regsBeliefs[reg][reg_assignInd];
	    }
	    break;
	  default:
	    break;
	}
      }

      // normalize
      double norm_bel_n = HUGE_VAL;
      double sum_bel_n = 0.0;
      for (int xn=0; xn<V[n]; xn++) {
	if (singleBeliefs[n][xn] < norm_bel_n) {
	  norm_bel_n = singleBeliefs[n][xn];
	}
      }
      for (int xn=0; xn<V[n]; xn++) {
	singleBeliefs[n][xn] -= norm_bel_n;
      }
    }
  }
}

void GBPPreProcessor::extractPairs(double** regsBeliefs,
				   double**** pairBeliefs,
				   SumOrMax smFlag) const {
  if (pp_mrf->logspace) {
    extractPairsLog(regsBeliefs,pairBeliefs,smFlag);
  }
  else {
    int const& N = pp_mrf->N;
    int const* V = pp_mrf->V;
    int* regCard = pp_regions_mrf->V;
  
    Assignment assignment;
    assignment.clear();
    assignment.resize(N);


    for (int i=0; i<N; i++) {
      for (int n=0; n<pp_mrf->neighbNum(i); n++) {
	int j = pp_mrf->adjMat[i][n];
	if (i<j) {
	  int reg = pp_extractPairs[i][n];
	  if (reg<0) {
	    for (int xi=0; xi<V[i]; xi++) {
	      for (int xj=0; xj<V[j]; xj++) {
		pairBeliefs[i][n][xi][xj] = 1.0 / (V[i] * V[j]);
	      }
	    }
	  }
	  else {
	    Region* region = pp_allRegionsPtr[reg];
	
	    for (int xi=0; xi<V[i]; xi++) {
	      for (int xj=0; xj<V[j]; xj++) {
	    
		switch (smFlag) {
		  case SUM:
		    pairBeliefs[i][n][xi][xj] = 0.0;
		    break;
		  case MAX:
		    pairBeliefs[i][n][xi][xj] = -1.0;
		    break;
		  default:
		    break;
		}
	      }
	    }
	    for (int reg_assignInd=0; reg_assignInd<regCard[reg]; reg_assignInd++) {
	      region->indToAssign(reg_assignInd, assignment, V);
	      int xi = assignment[i];
	      int xj = assignment[j];
	      switch (smFlag) {
		case SUM:
		  pairBeliefs[i][n][xi][xj] += regsBeliefs[reg][reg_assignInd];
		  break;
		case MAX:
		  if (regsBeliefs[reg][reg_assignInd] > pairBeliefs[i][n][xi][xj]) {
		    pairBeliefs[i][n][xi][xj] = regsBeliefs[reg][reg_assignInd];
		  }
		  break;
		default:
		  break;
	      }
	    }
      
	    // normalize for max-gbp
	    if (smFlag == MAX) {
	      double sum_bel_ij = 0.0;
	      for (int xi=0; xi<V[i]; xi++) {
		for (int xj=0; xj<V[j]; xj++) {
		  sum_bel_ij += pairBeliefs[i][n][xi][xj];
		}
	      }
	      for (int xi=0; xi<V[i]; xi++) {
		for (int xj=0; xj<V[j]; xj++) {
		  pairBeliefs[i][n][xi][xj] /= sum_bel_ij;
		}
	      }
	    }
	  }
	}
      }
    }
  }
}

void GBPPreProcessor::extractPairsLog(double** regsBeliefs,
				      double**** pairBeliefs,
				      SumOrMax smFlag) const {
  int const& N = pp_mrf->N;
  int const* V = pp_mrf->V;
  int* regCard = pp_regions_mrf->V;
  
  Assignment assignment;
  assignment.clear();
  assignment.resize(N);


  for (int i=0; i<N; i++) {
    for (int n=0; n<pp_mrf->neighbNum(i); n++) {
      int j = pp_mrf->adjMat[i][n];
      if (i<j) {
	int reg = pp_extractPairs[i][n];
	if (reg<0) {
	  for (int xi=0; xi<V[i]; xi++) {
	    for (int xj=0; xj<V[j]; xj++) {
	      pairBeliefs[i][n][xi][xj] = 0.0;
	    }
	  }
	}
	else {
	  Region* region = pp_allRegionsPtr[reg];
	
	  for (int xi=0; xi<V[i]; xi++) {
	    for (int xj=0; xj<V[j]; xj++) {
	    
	      switch (smFlag) {
		case SUM:
		  pairBeliefs[i][n][xi][xj] = -HUGE_VAL;
		  break;
		case MAX:
		  pairBeliefs[i][n][xi][xj] = HUGE_VAL;
		  break;
		default:
		  break;
	      }
	    }
	  }
	  for (int reg_assignInd=0; reg_assignInd<regCard[reg]; reg_assignInd++) {
	    region->indToAssign(reg_assignInd, assignment, V);
	    int xi = assignment[i];
	    int xj = assignment[j];
	    switch (smFlag) {
	      case SUM:
		pairBeliefs[i][n][xi][xj] = AddLogFactor(pairBeliefs[i][n][xi][xj],regsBeliefs[reg][reg_assignInd],
							 -pp_mrf->getTemperature());
		break;
	      case MAX:
		if (regsBeliefs[reg][reg_assignInd] < pairBeliefs[i][n][xi][xj]) {
		  pairBeliefs[i][n][xi][xj] = regsBeliefs[reg][reg_assignInd];
		}
		break;
	      default:
		break;
	    }
	  }
      
	  // normalize for max-gbp
	  double norm_bel_ij = HUGE_VAL;
	  for (int xi=0; xi<V[i]; xi++) {
	    for (int xj=0; xj<V[j]; xj++) {
	      if (pairBeliefs[i][n][xi][xj] < norm_bel_ij) {
		norm_bel_ij = pairBeliefs[i][n][xi][xj];
	      }
	    }
	  }
	  for (int xi=0; xi<V[i]; xi++) {
	    for (int xj=0; xj<V[j]; xj++) {
	      pairBeliefs[i][n][xi][xj] -= norm_bel_ij;
	    }
	  }
	}
      }
    }
  }
}

void GBPPreProcessor::printAllRegions() const {

  for (int i=0; i<pp_numAllRegs; i++) {
    cout << i << '\t' << (*pp_allRegionsPtr[i]);
  }
  cout << endl;
}

void GBPPreProcessor::printAdj() const {

  for (int i=0; i<pp_numAllRegs; i++) {
    Nodes& reg = (*pp_adjMat)[i];
    cout << i << '\t';
    for (unsigned int j=0; j<reg.size(); j++) {
      cout << reg[j] << ' ';
    }
    cout << endl;
  }
  cout << endl;
}

void GBPPreProcessor::printBethe() const {
  if (pp_bethe != 0) {
    for (int i=0; i<pp_numAllRegs; i++) {
      cout << i << " bethe: " << pp_bethe[i] << endl;
    }
    cout << endl;
  }
}