Statistics
| Branch: | Tag: | Revision:

## sssimulator / Matrix / matrix.c @ d6fa79fc

 1 ```/* ``` ``` * This is SSSim: the Simple & Stupid Simulator ``` ``` * ``` ``` * Copyright (c) 2015 Luca Baldesi ``` ``` * ``` ``` * This is free software; see gpl-3.0.txt ``` ``` */ ``` ```#include ``` ```#include ``` ```#include ``` ```#include ``` ```#include ``` ```#include "matrix.h" ``` ```#include "string_indexer.h" ``` ```#include "tokens.h" ``` ```#define MAX_EDGE_RECORD_STRING_LEN 255 ``` ```#define PRECISION 0.0001 ``` ```#define MAT_GET(m,i,j) \ ``` ``` (m->values[(i)+(j)*(m->rows)]) ``` ```#define MAT_SET(m,i,j,v) \ ``` ``` (m->values[(i)+(j)*(m->rows)] = (v)) ``` ```#define MAX(a,b) ((a) > (b) ? (a) : (b)) ``` ```#define MIN(a,b) ((a) < (b) ? (a) : (b)) ``` ```struct matrix { ``` ``` matsize rows; ``` ``` matsize cols; ``` ``` matvalue * values; ``` ```}; ``` ```matvalue matrix_element_get(const struct matrix * m, matsize i, matsize j) ``` ```{ ``` ``` if (m && i < m->rows && j < m->cols) ``` ``` return MAT_GET(m, i, j); ``` ``` return 0; ``` ```} ``` ```int matrix_element_set(const struct matrix * m, matsize i, matsize j, matvalue v) ``` ```{ ``` ``` int res = 0; ``` ``` if (m && i < m->rows && j < m->cols) ``` ``` MAT_SET(m, i, j, v); ``` ``` else ``` ``` res = -1; ``` ``` return res; ``` ```} ``` ```matsize matrix_num_rows(const struct matrix * m) ``` ```{ ``` ``` if (m) ``` ``` return m->rows; ``` ``` return 0; ``` ```} ``` ```matsize matrix_num_cols(const struct matrix * m) ``` ```{ ``` ``` if (m) ``` ``` return m->cols; ``` ``` return 0; ``` ```} ``` ```int matrix_init(struct matrix *m, matsize nrow, matsize ncol) ``` ```{ ``` ``` if (m) ``` ``` { ``` ``` m->rows = nrow; ``` ``` m->cols = ncol; ``` ``` m->values = (matvalue *) malloc (nrow*ncol*sizeof(matvalue)); ``` ``` memset(m->values, 0, nrow*ncol*sizeof(matvalue)); ``` ``` } ``` ``` return 0; ``` ```} ``` ```int matrix_make_bidir(struct matrix * m) ``` ```{ ``` ``` int res = 0; ``` ``` matsize i, j; ``` ``` matvalue v; ``` ``` if (m) ``` ``` { ``` ``` for(i = 0; i < m->rows; i++) ``` ``` for(j = 0; j < m->cols; j++) ``` ``` { ``` ``` v = MIN(MAT_GET(m, i, j), MAT_GET(m, j, i)); ``` ``` if (v == 0) ``` ``` { ``` ``` v = MAX(MAT_GET(m, i, j), MAT_GET(m, j, i)); ``` ``` MAT_SET(m, i, j, v); ``` ``` MAT_SET(m, j, i, v); ``` ``` } ``` ``` } ``` ``` } else ``` ``` res = -1; ``` ``` return res; ``` ```} ``` ```int matrix_add_edge_record(struct matrix * m, struct string_indexer * si, char * rec) ``` ``` //input: 0 4 0.56 ``` ```{ ``` ``` char ** tokens; ``` ``` uint32_t ntok; ``` ``` int res = 0; ``` ``` if (m && rec) ``` ``` { ``` ``` tokens = tokens_create(rec,' ',&ntok); ``` ``` if (tokens[0][0] != '#' && ntok >= 3) ``` ``` MAT_SET(m, string_indexer_id(si, (tokens[0])),\ ``` ``` string_indexer_id(si, (tokens[1])),\ ``` ``` atof(tokens[2])); ``` ``` else ``` ``` res = 1; ``` ``` tokens_destroy(&tokens, ntok); ``` ``` } else ``` ``` res = -1; ``` ``` return res; ``` ```} ``` ```matsize matrix_parse_edfile_names(struct string_indexer * si, FILE *fp) ``` ```{ ``` ``` char edge_record[MAX_EDGE_RECORD_STRING_LEN]; ``` ``` char ** tokens; ``` ``` uint32_t ntok; ``` ``` if (fp && si) ``` ``` { ``` ``` while (fgets(edge_record, MAX_EDGE_RECORD_STRING_LEN, fp) != NULL) ``` ``` { ``` ``` tokens = tokens_create(edge_record,' ',&ntok); ``` ``` if (tokens[0][0] != '#' && ntok >= 3) ``` ``` { ``` ``` string_indexer_id(si, tokens[0]); ``` ``` string_indexer_id(si, tokens[1]); ``` ``` } ``` ``` tokens_destroy(&tokens, ntok); ``` ``` } ``` ``` } ``` ``` return (matsize) string_indexer_size(si); ``` ```} ``` ```void matrix_print(const struct matrix *m) ``` ```{ ``` ``` matsize i, j; ``` ``` if(m && m->values) ``` ``` { ``` ``` fprintf(stderr,"----------------------------------\n"); ``` ``` for(i = 0; i < m->rows; i++) ``` ``` { ``` ``` for(j = 0; j < m->cols; j++) ``` ``` fprintf(stderr,"%f ", MAT_GET(m, i, j)); ``` ``` fprintf(stderr,"\n"); ``` ``` } ``` ``` fprintf(stderr,"----------------------------------\n"); ``` ``` } ``` ```} ``` ```struct matrix * matrix_from_edgefile(const char * edgefile) ``` ```{ ``` ``` struct matrix * m = NULL; ``` ``` FILE *fp; ``` ``` matsize num_nodes; ``` ``` char edge_record[MAX_EDGE_RECORD_STRING_LEN]; ``` ``` struct string_indexer * si; ``` ``` ``` ``` if (edgefile) ``` ``` { ``` ``` //fprintf(stderr, "[DEBUG] file: %s\n", edgefile); ``` ``` fp = fopen(edgefile, "r"); ``` ``` if (fp != NULL) ``` ``` { ``` ``` m = (struct matrix *) malloc(sizeof(struct matrix)); ``` ``` m->values = NULL; // pseudo init to prevent unexpected free ``` ``` si = string_indexer_new(100); ``` ``` num_nodes = matrix_parse_edfile_names(si, fp); ``` ``` //fprintf(stderr, "[DEBUG] nodes: %d\n", num_nodes); ``` ``` if(num_nodes > 0) ``` ``` { ``` ``` matrix_init(m, num_nodes, num_nodes); ``` ``` rewind(fp); ``` ``` while (fgets(edge_record, MAX_EDGE_RECORD_STRING_LEN, fp) != NULL) ``` ``` matrix_add_edge_record(m, si, edge_record); ``` ``` //matrix_print(m); ``` ``` } else ``` ``` matrix_destroy(&m); ``` ``` string_indexer_destroy(&si); ``` ``` fclose(fp); ``` ``` } ``` ``` } ``` ``` return m; ``` ```} ``` ```void matrix_destroy(struct matrix ** m) ``` ```{ ``` ``` if (m && *m) ``` ``` { ``` ``` if ((*m)->values) ``` ``` free((*m)->values); ``` ``` free(*m); ``` ``` } ``` ``` *m = NULL; ``` ```} ``` ```struct matrix * matrix_single_value(matsize rows, matsize cols, matvalue v) ``` ```{ ``` ``` struct matrix * m = NULL; ``` ``` matsize i,j; ``` ``` if (cols > 0 && rows > 0) ``` ``` { ``` ``` m = (struct matrix *) malloc(sizeof(struct matrix)); ``` ``` matrix_init(m, rows, cols); ``` ``` for(i = 0; i < m->rows; i++) ``` ``` for(j = 0; j < m->cols; j++) ``` ``` MAT_SET(m, i, j, v); ``` ``` } ``` ``` return m; ``` ```} ``` ```struct matrix * matrix_ones(matsize rows, matsize cols) ``` ```{ ``` ``` return matrix_single_value(rows, cols, 1); ``` ```} ``` ```struct matrix * matrix_multiply(const struct matrix * m1, const struct matrix * m2) ``` ```{ ``` ``` struct matrix * res = NULL; ``` ``` matsize i, j, k; ``` ``` matvalue v; ``` ``` if (m1 && m2 && m1->cols == m2->rows) ``` ``` { ``` ``` res = (struct matrix *) malloc(sizeof(struct matrix)); ``` ``` matrix_init(res, m1->rows, m2->cols); ``` ``` for(i = 0; i < m1->rows; i++) ``` ``` for(j = 0; j < m2->cols; j++) ``` ``` { ``` ``` v = 0; ``` ``` for(k = 0; k < m1->cols; k++) ``` ``` v += MAT_GET(m1, i, k) * MAT_GET(m2, k, j); ``` ``` MAT_SET(res, i, j, v); ``` ``` } ``` ``` } ``` ``` return res; ``` ```} ``` ```matvalue matrix_sum_values(const struct matrix *m) ``` ```{ ``` ``` matsize i, j; ``` ``` matvalue sum = 0; ``` ``` if (m) ``` ``` for(i = 0; i < m->rows; i++) ``` ``` for(j = 0; j < m->cols; j++) ``` ``` sum += MAT_GET(m, i, j); ``` ``` return sum; ``` ```} ``` ```void matrix_divide(struct matrix *m, matvalue v) ``` ```{ ``` ``` matsize i, j; ``` ``` if (m) ``` ``` for(i = 0; i < m->rows; i++) ``` ``` for(j = 0; j < m->cols; j++) ``` ``` MAT_SET(m, i, j, MAT_GET(m, i, j)/v); ``` ```} ``` ```struct matrix * matrix_difference(const struct matrix * m1, const struct matrix * m2) ``` ```{ ``` ``` struct matrix * res = NULL; ``` ``` matsize i, j; ``` ``` if (m1 && m2 && m1->rows == m2->rows && m1->cols == m2->cols) ``` ``` { ``` ``` res = malloc(sizeof(struct matrix)); ``` ``` matrix_init(res, m1->rows, m2->cols); ``` ``` for (i = 0; i < res->rows; i++) ``` ``` for (j = 0; j < res->cols; j++) ``` ``` MAT_SET(res, i, j, \ ``` ``` MAT_GET(m1, i, j) -\ ``` ``` MAT_GET(m2, i, j)); ``` ``` } ``` ``` return res; ``` ```} ``` ```matvalue matrix_distance(const struct matrix * m1, const struct matrix * m2) ``` ```{ ``` ``` struct matrix * diff; ``` ``` matsize i, j; ``` ``` matvalue v = 0; ``` ``` if (m1 && m2) ``` ``` { ``` ``` diff = matrix_difference(m1, m2); ``` ``` if (diff) ``` ``` { ``` ``` for(i = 0; i < diff->rows; i++) ``` ``` for(j = 0; j < diff->cols; j++) ``` ``` v += MAT_GET(diff, i, j) * MAT_GET(diff, i, j); ``` ``` v = sqrt(v); ``` ``` } ``` ``` matrix_destroy(&diff); ``` ``` } ``` ``` return v; ``` ```} ``` ``` ``` ```int matrix_sum(struct matrix * dst, const struct matrix * op1) ``` ```{ ``` ``` int res = 0; ``` ``` matsize i, j; ``` ``` if (dst && op1 && dst->cols == op1->cols && dst->rows == op1->rows) ``` ``` { ``` ``` for(i = 0; i < dst->rows; i++) ``` ``` for(j = 0; j < dst->cols; j++) ``` ``` MAT_SET(dst, i, j, MAT_GET(dst, i, j)\ ``` ``` + MAT_GET(op1, i, j)); ``` ``` } else ``` ``` res = -1; ``` ``` return res; ``` ```} ``` ```struct matrix * matrix_clone(const struct matrix * m) ``` ```{ ``` ``` struct matrix * c = NULL; ``` ``` if(m) ``` ``` { ``` ``` c = malloc(sizeof(struct matrix)); ``` ``` matrix_init(c, m->rows, m->cols); ``` ``` memmove(c->values, m->values, m->rows*m->cols*sizeof(matvalue)); ``` ``` } ``` ``` return c; ``` ```} ``` ``` ``` ```struct matrix * matrix_eigencentrality(const struct matrix * A, matvalue * dominant_eigenvalue) ``` ```{ ``` ``` struct matrix *x, *y = NULL, *S, *M = NULL; ``` ``` float m = 0.15; // decided by Google guys ``` ``` matvalue l; ``` ``` if (A && A->cols == A->rows) ``` ``` { ``` ``` //matrix_print(A); ``` ``` S = matrix_ones(A->rows, A->cols); ``` ``` matrix_divide(S, A->rows/m); ``` ``` M = matrix_clone(A); ``` ``` matrix_divide(M, 1/(1-m)); ``` ``` matrix_sum(M, S); ``` ``` x = matrix_single_value(A->rows, 1, 0); ``` ``` MAT_SET(x, 0, 0, 1); ``` ``` //matrix_print(M); ``` ``` y = matrix_multiply(M, x); ``` ``` l = matrix_sum_values(y); ``` ``` matrix_divide(y, l); ``` ``` while (matrix_distance(x, y) > PRECISION) ``` ``` { ``` ``` matrix_destroy(&x); ``` ``` x = y; ``` ``` y = matrix_multiply(M, x); ``` ``` l = matrix_sum_values(y); ``` ``` matrix_divide(y, l); ``` ``` } ``` ``` matrix_destroy(&x); ``` ``` matrix_destroy(&S); ``` ``` matrix_destroy(&M); ``` ``` //matrix_print(y); ``` ``` //fprintf(stderr,"[DEBUG] dominant eigenvalue: %f\n", l); ``` ``` if (dominant_eigenvalue) ``` ``` *dominant_eigenvalue = l; ``` ``` } ``` ``` return y; ``` ```} ``` ```struct matrix * matrix_diag(const struct matrix * m) ``` ```{ ``` ``` matsize i; ``` ``` struct matrix * diag = NULL; ``` ``` if (m && m->rows == 1) ``` ``` { ``` ``` diag = matrix_single_value(m->cols, m->cols, 0); ``` ``` for(i = 0; i < diag->cols; i++) ``` ``` { ``` ``` MAT_SET(diag, i, i, MAT_GET(m, 0, i)); ``` ``` } ``` ``` } else if (m && m->cols == 1) ``` ``` { ``` ``` diag = matrix_single_value(m->rows, m->rows, 0); ``` ``` for(i = 0; i < diag->rows; i++) ``` ``` MAT_SET(diag, i, i, MAT_GET(m, i, 0)); ``` ``` } ``` ``` return diag; ``` ```} ``` ```struct matrix * matrix_invdiag(const struct matrix * m) ``` ```{ ``` ``` matsize i; ``` ``` struct matrix * diag = NULL; ``` ``` if (m && m->rows == 1) ``` ``` { ``` ``` diag = matrix_single_value(m->cols, m->cols, 0); ``` ``` for(i = 0; i < diag->cols; i++) ``` ``` MAT_SET(diag, i, i, 1/MAT_GET(m, 0, i)); ``` ``` } else if (m && m->cols == 1) ``` ``` { ``` ``` diag = matrix_single_value(m->rows, m->rows, 0); ``` ``` for(i = 0; i < diag->rows; i++) ``` ``` MAT_SET(diag, i, i, 1/MAT_GET(m, i, 0)); ``` ``` } ``` ``` return diag; ``` ```} ``` ```int matrix_stochastify(struct matrix *m) ``` ``` // makes the matrix COLUMN-stochastic ``` ```{ ``` ``` int res = 0; ``` ``` struct matrix * u, *sums; ``` ``` matsize i, j; ``` ``` if(m) ``` ``` { ``` ``` u = matrix_ones(1, m->rows); ``` ``` sums = matrix_multiply(u, m); ``` ``` for(i = 0; i < m->rows; i++) ``` ``` for(j = 0; j < m->cols; j++) ``` ``` MAT_SET(m, i, j, MAT_GET(m, i, j)/MAT_GET(sums, 0, j)); ``` ``` matrix_destroy(&u); ``` ``` matrix_destroy(&sums); ``` ``` } else ``` ``` res = -1; ``` ``` return res; ``` ```} ``` ```void matrix_shrink(struct matrix **ptr, matsize pos) ``` ```{ ``` ``` struct matrix *neo, *m; ``` ``` matsize i, j; ``` ``` if(ptr && *ptr) ``` ``` { ``` ``` m = *ptr; ``` ``` neo = malloc(sizeof(struct matrix)); ``` ``` if (m->rows > 1 && m->cols > 1) // matrix case ``` ``` { ``` ``` matrix_init(neo, pos < m->rows ? m->rows - 1 : m->rows, pos < m->cols ? m->cols - 1 : m->cols); ``` ``` for(i = 0; i < m->rows; i++) ``` ``` for(j = 0; j < m->cols; j++) ``` ``` { ``` ``` if(i < pos) ``` ``` { ``` ``` if(j < pos) MAT_SET(neo, i, j, MAT_GET(m, i, j)); ``` ``` if(j > pos) MAT_SET(neo, i, j-1, MAT_GET(m, i, j)); ``` ``` } ``` ``` if(i > pos) ``` ``` { ``` ``` if(j < pos) MAT_SET(neo, i - 1, j, MAT_GET(m, i, j)); ``` ``` if(j > pos) MAT_SET(neo, i - 1, j-1, MAT_GET(m, i, j)); ``` ``` } ``` ``` } ``` ``` } else // vector case ``` ``` { ``` ``` if(m->rows == 1) // column vector ``` ``` { ``` ``` matrix_init(neo, m->rows, m->cols - 1); ``` ``` for(i = 0; i < m->cols; i++) ``` ``` { ``` ``` if(i < pos) MAT_SET(neo, 0, i, MAT_GET(m, 0, i)); ``` ``` if(i > pos) MAT_SET(neo, 0, i-1, MAT_GET(m, 0, i)); ``` ``` } ``` ``` } ``` ``` if(m->cols == 1) // row vector ``` ``` { ``` ``` matrix_init(neo, m->rows - 1, m->cols); ``` ``` for(i = 0; i < m->rows; i++) ``` ``` { ``` ``` if(i < pos) MAT_SET(neo, i, 0, MAT_GET(m, i, 0)); ``` ``` if(i > pos) MAT_SET(neo, i-1, 0, MAT_GET(m, i, 0)); ``` ``` } ``` ``` } ``` ``` } ``` ``` matrix_destroy(ptr); ``` ``` *ptr = neo; ``` ``` } ``` ```} ```