Daniel@0: /* triangulate.c written by Ilya Shpitser  */
Daniel@0: 
Daniel@0: #include <stdlib.h>
Daniel@0: 
Daniel@0: #ifdef UNIX
Daniel@0: #include "matlab.h"
Daniel@0: #endif
Daniel@0: 
Daniel@0: #include "matrix.h"
Daniel@0: #include "mex.h"
Daniel@0: 
Daniel@0: #include "elim.h"
Daniel@0: #include "map.h"
Daniel@0: #include "misc.h"
Daniel@0: 
Daniel@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
Daniel@0: 
Daniel@0:   int dims [2];
Daniel@0:   int i, j, k, m, n;
Daniel@0:   long index;
Daniel@0:   double * G_pr;
Daniel@0:   double * stage_pr;
Daniel@0:   double * answer_G_pr, * fill_ins_pr;
Daniel@0:   double * matlab_clique_pr;
Daniel@0:   mxArray * matlab_clique;
Daniel@0:   Elimination e;
Daniel@0:   float ** adj_mat;
Daniel@0:   int ** order = (int **) NULL;
Daniel@0:   Iterator iter, iter2;
Daniel@0:   word w, w2;
Daniel@0:   int ** fill_ins;
Daniel@0:   Map cliques;
Daniel@0:   Map clique;
Daniel@0:   mxArray * fill_ins_mat;
Daniel@0:   int * nodes;
Daniel@0:   mxArray * full;
Daniel@0: 
Daniel@0: // (original)  full = mlfFull((mxArray *) prhs[0]);
Daniel@0:   full = (mxArray *) mlfFull((mxArray *) prhs[0]);  // added typecasting
Daniel@0:   /* Obtain graph matrix information. */
Daniel@0:   m = mxGetM(full);
Daniel@0:   n = mxGetN(full);
Daniel@0:   G_pr = mxGetPr(full);
Daniel@0: 
Daniel@0:   if(n < 1 || m < 1){
Daniel@0:     return;
Daniel@0:   }
Daniel@0: 
Daniel@0:   /* Allocate and populate the log weight adjacency matrix corresponding
Daniel@0:      to the input graph. */
Daniel@0:   adj_mat = (float **) malloc(sizeof(float *) * m);
Daniel@0:   adj_mat[0] = (float *) malloc(sizeof(float) * m * n);
Daniel@0:   for(i = 1; i < m; i++){
Daniel@0:     adj_mat[i] = adj_mat[i - 1] + n;
Daniel@0:   }
Daniel@0:   /* We no longer have log weight info, but we have a (total) ordering on
Daniel@0:      the nodes already, so we do not need this information. */
Daniel@0:   for(i = 0; i < m; i++){
Daniel@0:     for(j = 0; j < n; j++){
Daniel@0:       index = j * m + i;
Daniel@0:       if(G_pr[index] > 0){
Daniel@0:         adj_mat[i][j] = 1;
Daniel@0:       } else {
Daniel@0:         adj_mat[i][j] = 0;
Daniel@0:       }
Daniel@0:     }
Daniel@0:   }
Daniel@0: 
Daniel@0:   /* Convert the total elimination ordering into a partial order argument
Daniel@0:      for the elimination routine.  The elimination routine's purpose in this
Daniel@0:      mode of operation is to return cliques and fill-in edges. */
Daniel@0:   if(nrhs > 1){
Daniel@0:     order = (int **) malloc(sizeof(int *) * m);
Daniel@0:     order[0] = (int *) malloc(sizeof(int) * m * n);
Daniel@0:     for(i = 1; i < m; i++){
Daniel@0:       order[i] = order[i - 1] + n;
Daniel@0:     }
Daniel@0:     for(i = 0; i < m; i++){
Daniel@0:       for(j = 0; j < n; j++){
Daniel@0:         order[i][j] = 0;
Daniel@0:       }
Daniel@0:     }
Daniel@0:     stage_pr = mxGetPr(prhs[1]);
Daniel@0:     for(i = 0; i < mxGetN(prhs[1]) - 1; i++){
Daniel@0:       order[(int) stage_pr[i] - 1][(int) stage_pr[i + 1] - 1] = 1;
Daniel@0:     }
Daniel@0:   }
Daniel@0: 
Daniel@0:   /* Find the elimination ordering. */
Daniel@0:   e = find_elim(n, adj_mat, order, -1);
Daniel@0: 
Daniel@0:   /* Allocate memory for the answer, and set the answer. */
Daniel@0:   plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
Daniel@0:   answer_G_pr = mxGetPr(plhs[0]);
Daniel@0:   cliques = get_cliques(e);
Daniel@0: /* 
Daniel@0:   dims[0] = 1;
Daniel@0:   dims[1] = get_size_Map(cliques);
Daniel@0:   plhs[1] = mxCreateCellArray(2, (const int *) dims);*/
Daniel@0:   plhs[1] = mxCreateCellMatrix(get_size_Map(cliques), 1);
Daniel@0:   fill_ins = get_fill_ins(e);
Daniel@0:   fill_ins_mat = mxCreateDoubleMatrix(m, n, mxREAL);
Daniel@0:   fill_ins_pr = mxGetPr(fill_ins_mat);
Daniel@0: 
Daniel@0:   for(i = 0; i < n; i++){
Daniel@0:     for(j = 0; j < m; j++){
Daniel@0:       index = j * m + i;
Daniel@0:       answer_G_pr[index] = G_pr[index];
Daniel@0:       if(fill_ins[i][j] > 0){
Daniel@0:         answer_G_pr[index] = 1;
Daniel@0:         fill_ins_pr[index] = 1;
Daniel@0:       }
Daniel@0:     }
Daniel@0:   }
Daniel@0:   mxDestroyArray(full);
Daniel@0: // (original)  plhs[2] = mlfSparse(fill_ins_mat, NULL, NULL, NULL, NULL, NULL);
Daniel@0:   plhs[2] = (mxArray *) mlfSparse(fill_ins_mat, NULL, NULL, NULL, NULL, NULL); // added typecasting
Daniel@0:   mxDestroyArray(fill_ins_mat);
Daniel@0:   nodes = (int *) malloc(sizeof(int) * n);
Daniel@0:   k = 0;
Daniel@0:   iter = get_Iterator(cliques);
Daniel@0:   while(!is_empty(iter)){
Daniel@0:     w = next_key(iter);
Daniel@0:     clique = (Map) w.v;
Daniel@0:     matlab_clique = mxCreateDoubleMatrix(1, get_size_Map(clique), mxREAL);
Daniel@0:     matlab_clique_pr = mxGetPr(matlab_clique);
Daniel@0:     for(i = 0; i < n; i++){
Daniel@0:       nodes[i] = 0;
Daniel@0:     }
Daniel@0:     iter2 = get_Iterator(clique);
Daniel@0:     while(!is_empty(iter2)){
Daniel@0:       w2 = next_key(iter2);
Daniel@0:       nodes[w2.i] = w2.i + 1;
Daniel@0:     }
Daniel@0:     j = 0;
Daniel@0:     for(i = 0; i < n; i++){
Daniel@0:       if(nodes[i] > 0){
Daniel@0:         matlab_clique_pr[j++] = nodes[i];
Daniel@0:       }
Daniel@0:     }
Daniel@0:     mxSetCell(plhs[1], k++, matlab_clique);
Daniel@0:   }
Daniel@0:   free(nodes);
Daniel@0: 
Daniel@0:   /* Finally, free the allocated memory. */
Daniel@0:   destroy_Elimination(e);
Daniel@0:   if(adj_mat){
Daniel@0:     if(adj_mat[0]){
Daniel@0:       free(adj_mat[0]);
Daniel@0:     }
Daniel@0:     free(adj_mat);
Daniel@0:   }
Daniel@0:   if(order){
Daniel@0:     if(order[0]){
Daniel@0:       free(order[0]);
Daniel@0:     }
Daniel@0:     free(order);
Daniel@0:   }
Daniel@0: }