src/common/tree_decomp.cc

/*

This tree decomposition based library was written by 
Kyle Ellrott (kellrott@csbl.bmb.uga.edu)
and is distributed under the GPL v2.0 license

*/


#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "tree_decomp.h"

td_graph * td_graph_init( long vertex_count, long *state_count, char **matrix ) {
        td_graph *graph = (td_graph *)malloc(sizeof(td_graph));
        
        long edge_count = 0;
        for (int i = 0; i < vertex_count; i++) {
                for (int j = i+1; j < vertex_count; j++) {
                        if ( matrix[i][j] )
                                edge_count++;
                }
        }
        
        
        graph->vertex_count = vertex_count;
        graph->edge_count = edge_count;
        
        if ( state_count !=  NULL ) {
                long cur_loc = 0;
                graph->state_count = (long *)malloc(sizeof(long) * vertex_count);
                for (int i = 0; i < vertex_count; i++) {
                        graph->state_count[i] = state_count[i];
                        cur_loc += state_count[i];
                }
                graph->vertex_score = (double **)malloc(sizeof(double *) * vertex_count );
                graph->vertex_score[0] = (double *)malloc(sizeof(double) * cur_loc );
                cur_loc = 0;
                for (int i = 0; i < vertex_count; i++) {
                        graph->vertex_score[i] = graph->vertex_score[0] + cur_loc;
                        cur_loc += graph->state_count[i];
                }
                graph->edge_score = (double ***)malloc(sizeof(double **) * vertex_count * vertex_count );
                for (int i = 0; i < vertex_count; i++) {
                        for (int j = 0; j < vertex_count; j++) {
                                long offset = td_edge_offset( graph, i, j ); 
                                if (i < j) {
                                        if ( matrix[i][j] ) {
                                                graph->edge_score[ offset ] = (double **)malloc(sizeof(double *) * state_count[i] );
                                                graph->edge_score[ offset ][0] = (double *)malloc(sizeof(double) * state_count[i] * state_count[j]);
                                                for (int k = 0; k < state_count[i]; k++) {
                                                        graph->edge_score[ offset ][k] = graph->edge_score[ offset ][0] + k * state_count[j];
                                                }
                                        } else {
                                                graph->edge_score[ offset ] = NULL;
                                        }
                                } else {
                                        graph->edge_score[ offset ] = NULL;
                                }
                        }
                }
        }
        graph->matrix = matrix;
        return graph;
}


char **td_matrix_init( long vertex_count ) {
        char **matrix = (char **)malloc(sizeof(char *) * vertex_count );
        matrix[0] = (char *)malloc(sizeof(char) * vertex_count * vertex_count);
        for (int i = 0; i < vertex_count; i++) {
                matrix[i] = matrix[0] + vertex_count * i;
                for (int j = 0; j < vertex_count; j++) {
                        matrix[i][j] = 0;
                }
        }       
        return matrix;  
}



void td_matrix_free( char ** graph) {
        free( graph[0] );
        free( graph );
}


void  td_graph_free( td_graph *graph) {
        free( graph->vertex_score[0] );
        free( graph->vertex_score );
        free( graph->state_count );
        for (int i = 0; i < graph->vertex_count; i++) {
                for (int j = 0; j < graph->vertex_count; j++) {
                        if ( graph->edge_score[i][j] ) {
                                free( graph->edge_score[ td_edge_offset(graph, i,j) ][0] );
                                free( graph->edge_score[ td_edge_offset(graph, i,j) ] );
                        }
                }       
        }       
}



long td_dist_count(char **matrix, long matrix_size, long cur_node, long source) {
        int i;
        long dist = 0, new_dist;
        for (i = 0; i < matrix_size; i++) {
                if (cur_node != i && i != source && matrix[cur_node][i]) {
                        new_dist = td_dist_count(matrix, matrix_size, i, cur_node);
                        if (new_dist > dist)
                                dist = new_dist;
                }
        }
        return dist + 1;        
}

void td_tree_clean_direction( td_decomp *decomp, int cur_bag, int parent) {
        for (int i=0; i < decomp->bag_count; i++) {
                if ( decomp->bag_matrix[cur_bag][i] ) {
                        if ( i == cur_bag || i == parent ) {
                                decomp->bag_matrix[cur_bag][i] = 0;
                        } else {
                                td_tree_clean_direction(decomp, i, cur_bag);
                        }
                }
        }       
}




int bag_intersect( td_decomp *decomp, int bag_1, int bag_2, char *answer ) {
        char mask[ decomp->graph->vertex_count ];
        int count = 0;  
        memset(mask, 0x00, decomp->graph->vertex_count );       
        if (answer != NULL) {
                memset(answer, 0x00, decomp->graph->vertex_count );     
        }
        if ( bag_1 != -1) {
                for (int i = 0; i < decomp->bag_size[bag_1]; i++) {
                        mask[ decomp->bag[bag_1][i] ] |= 0x01;
                }
        }
        if ( bag_2 != -1 ) {
                for (int i = 0; i < decomp->bag_size[bag_2]; i++) {
                        mask[ decomp->bag[bag_2][i] ] |= 0x02;
                }       
        }
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                if (mask[i] == 0x03) {
                        if (answer != NULL)
                                answer[i] = 0x01;
                        count++;
                }
        }       
        return count;
}


char bag_contains_intersect(  td_decomp *decomp, int bag_mid, int bag_1, int bag_2 ) {  
        char mask[ decomp->graph->vertex_count ];
        memset(mask, 0x00, decomp->graph->vertex_count );       
        for (int i = 0; i < decomp->bag_size[bag_1]; i++) {
                mask[ decomp->bag[bag_1][i] ] |= 0x01;
        }
        for (int i = 0; i < decomp->bag_size[bag_2]; i++) {
                mask[ decomp->bag[bag_2][i] ] |= 0x02;
        }       
        for (int i = 0; i < decomp->bag_size[bag_mid]; i++) {
                mask[ decomp->bag[bag_2][i] ] |= 0x04;
        }               
        
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                if (mask[i] == 0x03) {
                        return 0;
                }               
        }
        return 1;
}

td_decomp * td_graph_decomp_tri( td_graph *graph ) {
        td_decomp *decomp = (td_decomp *)malloc(sizeof(td_decomp));
        
        
        char **connect_matrix = td_matrix_init( graph->vertex_count );  
        for (int i = 0; i < graph->vertex_count; i++) {
                for (int j = 0; j < graph->vertex_count; j++) {
                        connect_matrix[i][j] = graph->matrix[i][j];
                }
        }
        
        char *active_vertex = (char *)malloc(sizeof(char) * graph->vertex_count);
        for (int i = 0; i < graph->vertex_count; i++) {
                active_vertex[i] = 1;
        }
        
        long **bag_list = (long **)malloc(sizeof(long *) * graph->vertex_count);
        long *bag_size = (long *)malloc(sizeof(long) * graph->vertex_count);
        bag_list[0]  =(long *)malloc(sizeof(long) * graph->vertex_count * graph->vertex_count );
        for (int i = 0; i < graph->vertex_count; i++) {
                bag_size[i] = 0;
                bag_list[i] = bag_list[0] + i * (graph->vertex_count);
        }
        
        long active_vertex_count = graph->vertex_count;
        long bag_count = 0;
        long least_edge, least_count;
        do {
                //find the edge with the least number of connections
                least_edge = -1;
                least_count = 0;
                for (int i = 0; i < graph->vertex_count; i++) {
                        if (active_vertex[i]) {
                                int connect_count = 0;
                                for (int j = 0; j < graph->vertex_count; j++) {
                                        if (active_vertex[j] && i != j && connect_matrix[i][j]) {
                                                connect_count++;
                                        }
                                }
                                if (least_edge == -1 || least_count > connect_count) {
                                        least_edge = i;
                                        least_count = connect_count;
                                }
                        }
                }
                
                if (least_edge != -1 && active_vertex_count > least_count) {
                        active_vertex[least_edge] = 0;
                        int l = 0;
                        //begin removing vertex
                        //fprintf(stderr, "%d: %d", bag_count, least_edge);
                        bag_list[ bag_count ][0] = least_edge;
                        active_vertex_count--;
                        l++;
                        //list connected edges, and connect them to each other
                        for (int i = 0; i < graph->vertex_count; i++) {
                                if (active_vertex[i] && connect_matrix[i][least_edge]) {
                                        //fprintf(stderr, " %d", i);
                                        bag_list[ bag_count ][ l  ] = i;
                                        l++;
                                        for (int j = 0; j < graph->vertex_count; j++) {
                                                if (active_vertex[j] && i != j && connect_matrix[j][least_edge]) {
                                                        connect_matrix[i][j] = 1;
                                                        connect_matrix[j][i] = 1;
                                                }
                                        }
                                }
                        }
                        //fprintf(stderr, "\n");
                        bag_size[bag_count] = l;
                        bag_count++;
                }               
        } while (least_edge != -1 && active_vertex_count > least_count);
        free(active_vertex);
        char **bag_matrix = td_matrix_init( bag_count );
        
        decomp->bag_matrix = bag_matrix;
        decomp->bag = bag_list;
        decomp->bag_size = bag_size;
        decomp->bag_count = bag_count;  
        decomp->graph = graph;
        

        char *active_bag = (char *)malloc(sizeof(char) * bag_count);
        memset(active_bag, 0, sizeof(char) * bag_count);
        
        int active_bag_count = 1;
        active_bag[0] = 1;
        while ( active_bag_count < bag_count ) {
                //find the bag that has the maximal cardinality of the seperator set with one of the bags already in the active set             
                int max_src = -1, max_dst = -1, max_count = 0;
                for (int i = 0; i < bag_count; i++) {
                        if ( active_bag[i] ) {
                                for (int j = 0; j < bag_count; j++) {
                                        if ( !active_bag[j] ) {
                                                int set_size = bag_intersect( decomp, i, j, NULL );
                                                if ( set_size > max_count ) {
                                                        max_src = i;
                                                        max_dst = j;
                                                        max_count = set_size;
                                                }
                                        }
                                }
                        }                       
                }
                if (max_count > 0) {
                        active_bag[ max_dst ] = 1;
                        decomp->bag_matrix[max_src][max_dst] = 1;
                        decomp->bag_matrix[max_dst][max_src] = 1;
                        active_bag_count++;
                } else {
                        if ( active_bag_count < bag_count ) {
                                //we should only reach this point if there are disconnected bags on the graph
                                //we're going to attach one of them at an arbitrary point
                                for (int i = 0; i < bag_count; i++) {
                                        if ( !active_bag[i] ) {
                                                for (int j = 0; j < bag_count; j++) {
                                                        if (active_bag[j]) {
                                                                active_bag[i] = 1;
                                                                decomp->bag_matrix[i][j] = 1;
                                                                decomp->bag_matrix[j][i] = 1;
                                                                active_bag_count++;
                                                                i = bag_count;
                                                                j = bag_count;
                                                        }
                                                }                                               
                                        }
                                }                               
                        }
                }
        }
        free(active_bag);
                
        //pick the root

        long dist = -1;
        for (int i = 0; i < bag_count; i++) {
                int k = td_dist_count(bag_matrix, bag_count, i, -1);
                if (dist == -1 || k < dist) {
                        decomp->root = i;
                        dist = k;
                }
        }

        decomp->root = 0;
        
        //make the connection matrix respect the directions
        td_tree_clean_direction( decomp, decomp->root, -1);     
        
        return decomp;
}


long mask_enum( td_decomp *decomp, char *mask, long *set_pos) {
        long out = 0;
        long row_size = 1;
        int mask_count = 0;
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                if (mask[i]) {
                        out += set_pos[i] * row_size;
                        row_size *= decomp->graph->state_count[i];
                        mask_count++;
                }
        }       
        if (mask_count == 0) {
                return -1;
        }
        return out;
}

void  mask_de_enum( td_decomp *decomp, char *mask, long pos, long *set_pos) {
        long in_pos = pos;
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                if (mask[i]) {
                        set_pos[i] = in_pos % decomp->graph->state_count[i];
                        in_pos = in_pos / decomp->graph->state_count[i];
                } else {
                        set_pos[i] = 0;
                }
        }       
}


typedef struct td_row {
        int chosen;
        double  val;
} td_row;



int td_graph_decomp_extract( td_decomp *decomp, int cur_bag, int parent, td_row **tables, long *sol ) {
        if ( parent != -1 ) {           
                char isec_mask[ decomp->graph->vertex_count ];
                bag_intersect( decomp, cur_bag, parent, isec_mask );
                long chosen = mask_enum( decomp, isec_mask, sol ); 
                //printf("bag %d - %d\n", cur_bag, chosen);
                char uniq_mask[ decomp->graph->vertex_count ];
                memcpy( uniq_mask, isec_mask, decomp->graph->vertex_count );
                for (int i = 0; i < decomp->bag_size[ cur_bag ]; i++) {
                        if ( ! uniq_mask[ decomp->bag[cur_bag][i] ] ) {
                                uniq_mask[ decomp->bag[cur_bag][i] ] = 1;
                        } else {
                                uniq_mask[ decomp->bag[cur_bag][i] ] = 0;
                        }
                }
                for (int i = 0; i < decomp->bag_size[ parent ]; i++) {
                        uniq_mask[ decomp->bag[parent][i] ] = 0;
                }
                long sel_pos[ decomp->graph->vertex_count ];
                mask_de_enum( decomp, uniq_mask, tables[ cur_bag ][ chosen ].chosen, sel_pos);
                for (int i = 0; i < decomp->graph->vertex_count; i++) {
                        if ( uniq_mask[i] )
                                sol[i] = sel_pos[i];
                }
                //now determine for children
                for (int i = 0; i < decomp->bag_count; i++) {
                        if ( decomp->bag_matrix[cur_bag][i] ) {
                                td_graph_decomp_extract( decomp, i, cur_bag, tables, sol);
                        }
                }
        } else {
                char uniq_mask[ decomp->graph->vertex_count ];
                memset(uniq_mask, 0x00, decomp->graph->vertex_count );
                for (int i = 0; i < decomp->bag_size[ cur_bag ]; i++) {
                        uniq_mask[ decomp->bag[cur_bag][i] ] = 1;
                }
                long cur_pos[ decomp->graph->vertex_count ];
                mask_de_enum( decomp, uniq_mask, tables[ cur_bag ][0].chosen, cur_pos ); 
                /*
                printf("::bag %d - ", cur_bag);
                for (int i = 0; i < decomp->graph->vertex_count; i++) {
                        if ( uniq_mask[i] )
                                printf("%d(%d) ", i, cur_pos[i]);
                }
                printf("\n");
                 */
                for (int i = 0; i < decomp->graph->vertex_count; i++) {
                        if (uniq_mask[i])
                                sol[i] = cur_pos[i];
                }                         
                for (int i = 0; i < decomp->bag_count; i++) {
                        if ( decomp->bag_matrix[cur_bag][i] ) {
                                td_graph_decomp_extract( decomp, i, cur_bag, tables, sol);
                        }
                }
        }
}



int td_graph_decomp_solve_explore( td_decomp *decomp, int cur_bag, int parent, td_row **tables ) {
        //first solve the children
        for (int i = 0; i < decomp->bag_count; i++) {
                if ( decomp->bag_matrix[cur_bag][i] ) {
                        td_graph_decomp_solve_explore( decomp, i, cur_bag, tables );
                }
        }       
        
        //start by find intersections with the parent, if we don't have any,
        //or we don't have a parent, special steps will be taken
        char isec_mask[ decomp->graph->vertex_count ];
        int isec_count = bag_intersect( decomp, cur_bag, parent, isec_mask );
        long row_count = 1;
        //if these is no intersection with the parent, we only have one row, which is the
        //basis of the answer
        if ( isec_count > 0 ) {
                //figure out how many rows it will take to enumerate all the combinations
                //of the vertices that intersect with out parent
                row_count = 1;
                for (int i = 0; i < decomp->graph->vertex_count; i++) {
                        if ( isec_mask[i] ) {
                                row_count *= decomp->graph->state_count[ i ];
                        }
                }
        } 
                
        tables[ cur_bag ] = (td_row *)malloc( sizeof(td_row) * row_count );
        //figure out what vertices this bag has that it parent doesnt
        //we will have to find the best combination of these for each of the  
        //rows in the table of shared vertices combinations
        char uniq_mask[ decomp->graph->vertex_count ];
        memcpy( uniq_mask, isec_mask, decomp->graph->vertex_count );
        for (int i = 0; i < decomp->bag_size[ cur_bag ]; i++) {
                if ( ! uniq_mask[ decomp->bag[cur_bag][i] ] ) {
                        uniq_mask[ decomp->bag[cur_bag][i] ] = 1;
                } else {
                        uniq_mask[ decomp->bag[cur_bag][i] ] = 0;
                }
        }
        if ( parent != -1) {
                for (int i = 0; i < decomp->bag_size[ parent ]; i++) {
                        uniq_mask[ decomp->bag[parent][i] ] = 0;
                }
        }
        int uniq_count = 0; 
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                if (uniq_mask[i])
                        uniq_count++;
        }
        //figure out the unique row count, for ever row in the table
        //this is how many possibilities we will have to search
        long uniq_row_count = 1;
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                if ( uniq_mask[i] ) {
                        uniq_row_count *= decomp->graph->state_count[ i  ];
                }
        }
        
        //printf("Bag %d : (rows: %d)  (uniq_rows: %d) \n", cur_bag, row_count, uniq_row_count );
        
        //setup the children
        //we will be looking at them when we try to figure out the best combination of 
        //vertices that don't occur in the parent
        //remember, each bag has a table that enumerate all possible combinations
        //of vertices that occur in the parent annd the child, in this case, the parent would be the 
        //current bag
        char **child_isec = (char **)malloc(sizeof(char *) * decomp->bag_count );               
        int child_count = 0;
        for (int i = 0; i < decomp->bag_count; i++) {
                if ( decomp->bag_matrix[cur_bag][i] ) {
                        child_isec[i] = (char *)malloc(sizeof(char) * decomp->graph->vertex_count);
                        bag_intersect( decomp, cur_bag, i, child_isec[ i ] );
                        child_count++;
                } else {
                        child_isec[i] = NULL;
                }
        }
        
        //fill in the rows 
        //remember, for each row, find the best combination of vertices that don't occur in the parent
        //each row is a different combination of vertices that do occur in the parent
        for (int i = 0; i < row_count; i++) {
                long cur_pos[ decomp->graph->vertex_count ];
                if ( isec_count > 0 )
                        mask_de_enum( decomp, isec_mask, i, cur_pos );
                else 
                        mask_de_enum( decomp, uniq_mask, i, cur_pos );
                
                int min = -1;
                double min_val = 1e20; //make the assumption that the min val will be less then 1e20
                for (int j = 0; j < uniq_row_count; j++) {
                        //figure out what enumeration the current row is eqivelant to
                        //and add those scores into the value
                        long row_pos[ decomp->graph->vertex_count ];
                        mask_de_enum( decomp, uniq_mask, j, row_pos );
                        double val = 0;
                        for (int k = 0; k < decomp->graph->vertex_count; k++) {
                                if (uniq_mask[k]) {
                                        val = decomp->graph->vertex_score[k][ row_pos[k] ];
                                        for (int l = k+1; l < decomp->graph->vertex_count; l++) {
                                                if ( uniq_mask[l] ) {
                                                        if ( decomp->graph->matrix[k][l] ) {
                                                                int offset = td_edge_offset( decomp->graph, k, l);
                                                                val += decomp->graph->edge_score[ offset ][ row_pos[k] ][ row_pos[l] ];
                                                                //val += td_edge_get( decomp->graph, k, l, row_pos[k], row_pos[l] );
                                                        }
                                                }
                                        }
                                }
                        }
                        //for each of the children, look at the child row that corresponds to the current row
                        for (int k = 0; k < decomp->bag_count; k++) {
                                if ( decomp->bag_matrix[cur_bag][k] ) {
                                        int child_row_num = mask_enum( decomp, child_isec[k], cur_pos);
                                        if ( child_row_num != -1 )
                                                val += tables[ k ][ child_row_num ].val;                                        
                                }
                        }
                        if ( val < min_val ) {
                                min = j;
                                min_val = val;
                        }                               
                }
                tables[ cur_bag ][ i ].chosen = min;
                tables[ cur_bag ][ i ].val = min_val;
        }
        
        for (int i = 0; i < decomp->bag_count; i++) {
                if (child_isec[i] != NULL) {
                        free(child_isec[i]);
                }
        }
        free(child_isec);
}


long *td_graph_decomp_solve( td_decomp *decomp  ) {
        if ( decomp->bag_count == 0)
                return NULL;
        
        td_row **tables = (td_row **)malloc(sizeof(td_row *) * decomp->bag_count );
        for (int i = 0; i < decomp->bag_count; i++) {
                tables[i] = NULL;
        }       
        td_graph_decomp_solve_explore( decomp, decomp->root, -1, tables );
        long *solution = (long *)malloc(sizeof(long) *decomp->graph->vertex_count );
        for (int i = 0; i <  decomp->graph->vertex_count ; i++) {
                solution[i] = -1;
        }
        td_graph_decomp_extract( decomp, decomp->root, -1, tables, solution );  
        /*
        printf("solution: ");
        for (int i = 0; i < decomp->graph->vertex_count; i++) {
                printf("%d(%d) ", i, solution[i]); 
        }
        printf("\n");   
         */
        for (int i = 0; i < decomp->bag_count; i++) {
                if ( tables[i] != NULL) {
                        free( tables[i] );
                }
        }       
        return solution;
}

---