Mauch MIREX 2010

/* triangulate.c written by Ilya Shpitser  */

#include <stdlib.h>

#ifdef UNIX

#include "matlab.h"

#endif

#include "matrix.h"

#include "mex.h"

#include "elim.h"

#include "map.h"

#include "misc.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){

  int dims [2];

  int i, j, k, m, n;

  long index;

  double * G_pr;

  double * stage_pr;

  double * answer_G_pr, * fill_ins_pr;

  double * matlab_clique_pr;

  mxArray * matlab_clique;

  Elimination e;

  float ** adj_mat;

  int ** order = (int **) NULL;

  Iterator iter, iter2;

  word w, w2;

  int ** fill_ins;

  Map cliques;

  Map clique;

  mxArray * fill_ins_mat;

  int * nodes;

  mxArray * full;

// (original)  full = mlfFull((mxArray *) prhs[0]);

  full = (mxArray *) mlfFull((mxArray *) prhs[0]);  // added typecasting

  /* Obtain graph matrix information. */

  m = mxGetM(full);

  n = mxGetN(full);

  G_pr = mxGetPr(full);

  if(n < 1 || m < 1){

    return;

  }

  /* Allocate and populate the log weight adjacency matrix corresponding

     to the input graph. */

  adj_mat = (float **) malloc(sizeof(float *) * m);

  adj_mat[0] = (float *) malloc(sizeof(float) * m * n);

  for(i = 1; i < m; i++){

    adj_mat[i] = adj_mat[i - 1] + n;

  }

  /* We no longer have log weight info, but we have a (total) ordering on

     the nodes already, so we do not need this information. */

  for(i = 0; i < m; i++){

    for(j = 0; j < n; j++){

      index = j * m + i;

      if(G_pr[index] > 0){

        adj_mat[i][j] = 1;

      } else {

        adj_mat[i][j] = 0;

      }

    }

  }

  /* Convert the total elimination ordering into a partial order argument

     for the elimination routine.  The elimination routine's purpose in this

     mode of operation is to return cliques and fill-in edges. */

  if(nrhs > 1){

    order = (int **) malloc(sizeof(int *) * m);

    order[0] = (int *) malloc(sizeof(int) * m * n);

    for(i = 1; i < m; i++){

      order[i] = order[i - 1] + n;

    }

    for(i = 0; i < m; i++){

      for(j = 0; j < n; j++){

        order[i][j] = 0;

      }

    }

    stage_pr = mxGetPr(prhs[1]);

    for(i = 0; i < mxGetN(prhs[1]) - 1; i++){

      order[(int) stage_pr[i] - 1][(int) stage_pr[i + 1] - 1] = 1;

    }

  }

  /* Find the elimination ordering. */

  e = find_elim(n, adj_mat, order, -1);

  /* Allocate memory for the answer, and set the answer. */

  plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);

  answer_G_pr = mxGetPr(plhs[0]);

  cliques = get_cliques(e);

/* 

  dims[0] = 1;

  dims[1] = get_size_Map(cliques);

  plhs[1] = mxCreateCellArray(2, (const int *) dims);*/

  plhs[1] = mxCreateCellMatrix(get_size_Map(cliques), 1);

  fill_ins = get_fill_ins(e);

  fill_ins_mat = mxCreateDoubleMatrix(m, n, mxREAL);

  fill_ins_pr = mxGetPr(fill_ins_mat);

  for(i = 0; i < n; i++){

    for(j = 0; j < m; j++){

      index = j * m + i;

      answer_G_pr[index] = G_pr[index];

      if(fill_ins[i][j] > 0){

        answer_G_pr[index] = 1;

        fill_ins_pr[index] = 1;

      }

    }

  }

  mxDestroyArray(full);

// (original)  plhs[2] = mlfSparse(fill_ins_mat, NULL, NULL, NULL, NULL, NULL);

  plhs[2] = (mxArray *) mlfSparse(fill_ins_mat, NULL, NULL, NULL, NULL, NULL); // added typecasting

  mxDestroyArray(fill_ins_mat);

  nodes = (int *) malloc(sizeof(int) * n);

  k = 0;

  iter = get_Iterator(cliques);

  while(!is_empty(iter)){

    w = next_key(iter);

    clique = (Map) w.v;

    matlab_clique = mxCreateDoubleMatrix(1, get_size_Map(clique), mxREAL);

    matlab_clique_pr = mxGetPr(matlab_clique);

    for(i = 0; i < n; i++){

      nodes[i] = 0;

    }

    iter2 = get_Iterator(clique);

    while(!is_empty(iter2)){

      w2 = next_key(iter2);

      nodes[w2.i] = w2.i + 1;

    }

    j = 0;

    for(i = 0; i < n; i++){

      if(nodes[i] > 0){

        matlab_clique_pr[j++] = nodes[i];

      }

    }

    mxSetCell(plhs[1], k++, matlab_clique);

  }

  free(nodes);

  /* Finally, free the allocated memory. */

  destroy_Elimination(e);

  if(adj_mat){

    if(adj_mat[0]){

      free(adj_mat[0]);

    }

    free(adj_mat);

  }

  if(order){

    if(order[0]){

      free(order[0]);

    }

    free(order);

  }

}