Statistics
| Branch: | Tag: | Revision:

sssimulator / Matrix / matrix.c @ d6fa79fc

History | View | Annotate | Download (10.2 KB)

1
/*
2
 * This is SSSim: the Simple & Stupid Simulator
3
 *
4
 *  Copyright (c) 2015 Luca Baldesi
5
 *
6
 *  This is free software; see gpl-3.0.txt
7
 */
8

    
9
#include <string.h>
10
#include <stdio.h>
11
#include <malloc.h>
12
#include <stdlib.h>
13
#include <math.h>
14

    
15
#include "matrix.h"
16
#include "string_indexer.h"
17
#include "tokens.h"
18

    
19
#define MAX_EDGE_RECORD_STRING_LEN 255
20
#define PRECISION 0.0001
21

    
22
#define MAT_GET(m,i,j) \
23
        (m->values[(i)+(j)*(m->rows)])
24
#define MAT_SET(m,i,j,v) \
25
        (m->values[(i)+(j)*(m->rows)] = (v))
26
#define MAX(a,b) ((a) > (b) ? (a) : (b))
27
#define MIN(a,b) ((a) < (b) ? (a) : (b))
28

    
29
struct matrix {
30
        matsize rows;
31
        matsize cols;
32
        matvalue * values;
33
};
34

    
35
matvalue matrix_element_get(const struct matrix * m, matsize i, matsize j)
36
{
37
        if (m && i < m->rows && j < m->cols)
38
                return MAT_GET(m, i, j);
39
        return 0;
40
}
41

    
42
int matrix_element_set(const struct matrix * m, matsize i, matsize j, matvalue v)
43
{
44
        int res = 0;
45

    
46
        if (m && i < m->rows && j < m->cols)
47
                MAT_SET(m, i, j, v);
48
        else
49
                res = -1;
50
        return res;
51
}
52

    
53
matsize matrix_num_rows(const struct matrix * m)
54
{
55
        if (m)
56
                return m->rows;
57
        return 0;
58
}
59

    
60
matsize matrix_num_cols(const struct matrix * m)
61
{
62
        if (m)
63
                return m->cols;
64
        return 0;
65
}
66

    
67
int matrix_init(struct matrix *m, matsize nrow, matsize ncol)
68
{
69
        if (m)
70
        {
71
                m->rows = nrow;
72
                m->cols = ncol;
73
                m->values = (matvalue *) malloc (nrow*ncol*sizeof(matvalue));
74
                memset(m->values, 0, nrow*ncol*sizeof(matvalue));
75
        }
76
        return 0;
77
}
78

    
79
int matrix_make_bidir(struct matrix * m)
80
{
81
        int res = 0;
82
        matsize i, j;
83
        matvalue v;
84

    
85
        if (m)
86
        {
87
                for(i = 0; i < m->rows; i++)
88
                        for(j = 0; j < m->cols; j++)
89
                        {
90
                                v = MIN(MAT_GET(m, i, j), MAT_GET(m, j, i));
91
                                if (v == 0)
92
                                {
93
                                        v = MAX(MAT_GET(m, i, j), MAT_GET(m, j, i));
94
                                        MAT_SET(m, i, j, v);
95
                                        MAT_SET(m, j, i, v);
96
                                }
97
                        }
98
        } else
99
                res = -1;
100

    
101
        return res;
102
}
103

    
104
int matrix_add_edge_record(struct matrix * m, struct string_indexer * si, char * rec)
105
        //input: 0 4 0.56
106
{
107
        char ** tokens;
108
        uint32_t ntok;
109
        int res = 0;
110

    
111
        if (m && rec)
112
        {
113
                tokens = tokens_create(rec,' ',&ntok);
114
                if (tokens[0][0] != '#' && ntok >= 3)
115
                        MAT_SET(m, string_indexer_id(si, (tokens[0])),\
116
                                string_indexer_id(si, (tokens[1])),\
117
                                atof(tokens[2]));
118
                else
119
                        res = 1;
120

    
121
                tokens_destroy(&tokens, ntok);
122
        } else
123
                res = -1;
124
        return res;
125
}
126

    
127
matsize matrix_parse_edfile_names(struct string_indexer * si, FILE *fp)
128
{
129
        char edge_record[MAX_EDGE_RECORD_STRING_LEN];
130
        char ** tokens;
131
        uint32_t ntok;
132

    
133
        if (fp && si)
134
        {
135
                while (fgets(edge_record, MAX_EDGE_RECORD_STRING_LEN, fp) != NULL)
136
                {
137
                        tokens = tokens_create(edge_record,' ',&ntok);
138
                        if (tokens[0][0] != '#' && ntok >= 3)
139
                        {
140
                                string_indexer_id(si, tokens[0]);
141
                                string_indexer_id(si, tokens[1]);
142
                        }
143
                        tokens_destroy(&tokens, ntok);
144
                }
145
        }
146

    
147
        return (matsize) string_indexer_size(si);
148
}
149

    
150
void matrix_print(const struct matrix *m)
151
{
152
        matsize i, j;
153

    
154
        if(m && m->values)
155
        {
156
                fprintf(stderr,"----------------------------------\n");
157
                for(i = 0; i < m->rows; i++)
158
                {
159
                        for(j = 0; j < m->cols; j++)
160
                                fprintf(stderr,"%f ", MAT_GET(m, i, j));
161
                        fprintf(stderr,"\n");
162
                }
163
                fprintf(stderr,"----------------------------------\n");
164
        }
165
}
166

    
167
struct matrix * matrix_from_edgefile(const char * edgefile)
168
{
169
        struct matrix * m = NULL;
170
        FILE *fp;
171
        matsize num_nodes;
172
        char edge_record[MAX_EDGE_RECORD_STRING_LEN];
173
        struct string_indexer * si;
174
        
175
        if (edgefile)
176
        {
177
                //fprintf(stderr, "[DEBUG] file: %s\n", edgefile);
178
                fp = fopen(edgefile, "r");
179
                if (fp != NULL)
180
                {
181
                        m = (struct matrix *) malloc(sizeof(struct matrix));
182
                        m->values = NULL; // pseudo init to prevent unexpected free
183
                        si = string_indexer_new(100);
184

    
185
                        num_nodes = matrix_parse_edfile_names(si, fp);
186
                        //fprintf(stderr, "[DEBUG] nodes: %d\n", num_nodes);
187
                        if(num_nodes > 0)
188
                        {
189
                                matrix_init(m, num_nodes, num_nodes);
190

    
191
                                rewind(fp);
192
                                while (fgets(edge_record, MAX_EDGE_RECORD_STRING_LEN, fp) != NULL)
193
                                        matrix_add_edge_record(m, si, edge_record);
194
                                //matrix_print(m);
195
                        } else
196
                                matrix_destroy(&m);
197
                        string_indexer_destroy(&si);
198
                        fclose(fp);
199
                }
200
        }
201
        return m;
202
}
203

    
204
void matrix_destroy(struct matrix ** m)
205
{
206
        if (m && *m)
207
        {
208
                if ((*m)->values)
209
                        free((*m)->values);
210
                free(*m);
211
        }
212
        *m = NULL;
213
}
214

    
215
struct matrix * matrix_single_value(matsize rows, matsize cols, matvalue v)
216
{
217
        struct matrix * m = NULL;
218
        matsize i,j;
219

    
220
        if (cols > 0 && rows > 0)
221
        {
222
                m = (struct matrix *) malloc(sizeof(struct matrix));
223
                matrix_init(m, rows, cols);
224
                for(i = 0; i < m->rows; i++)
225
                        for(j = 0; j < m->cols; j++)
226
                                MAT_SET(m, i, j, v);
227
        }
228
        return m;
229
}
230

    
231
struct matrix * matrix_ones(matsize rows, matsize cols)
232
{
233
        return matrix_single_value(rows, cols, 1);
234
}
235

    
236
struct matrix * matrix_multiply(const struct matrix * m1, const struct matrix * m2)
237
{
238
        struct matrix * res = NULL;
239
        matsize i, j, k;
240
        matvalue v;
241

    
242
        if (m1 && m2 && m1->cols == m2->rows)
243
        {
244
                res = (struct matrix *) malloc(sizeof(struct matrix));
245
                matrix_init(res, m1->rows, m2->cols);
246

    
247
                for(i = 0; i < m1->rows; i++)
248
                        for(j = 0; j < m2->cols; j++)
249
                        {
250
                                v = 0;
251
                                for(k = 0; k < m1->cols; k++)
252
                                        v += MAT_GET(m1, i, k) * MAT_GET(m2, k, j);
253
                                MAT_SET(res, i, j, v);
254
                        }
255
        }
256
        return res;
257
}
258

    
259
matvalue matrix_sum_values(const struct matrix *m)
260
{
261
        matsize i, j;
262
        matvalue sum = 0;
263

    
264
        if (m)
265
                for(i = 0; i < m->rows; i++)
266
                        for(j = 0; j < m->cols; j++)
267
                                sum += MAT_GET(m, i, j);
268
        return sum;
269
}
270

    
271
void matrix_divide(struct matrix *m, matvalue v)
272
{
273
        matsize i, j;
274

    
275
        if (m)
276
                for(i = 0; i < m->rows; i++)
277
                        for(j = 0; j < m->cols; j++)
278
                                MAT_SET(m, i, j, MAT_GET(m, i, j)/v);
279
}
280

    
281
struct matrix * matrix_difference(const struct matrix * m1, const struct matrix * m2)
282
{
283
        struct matrix * res = NULL;
284
        matsize i, j;
285

    
286
        if (m1 && m2 && m1->rows == m2->rows && m1->cols == m2->cols)
287
        {
288
                res = malloc(sizeof(struct matrix));
289
                matrix_init(res, m1->rows, m2->cols);
290
                for (i = 0; i < res->rows; i++)
291
                        for (j = 0; j < res->cols; j++)
292
                                MAT_SET(res, i, j, \
293
                                                MAT_GET(m1, i, j) -\
294
                                                MAT_GET(m2, i, j));
295
        }
296
        return res;
297
}
298

    
299
matvalue matrix_distance(const struct matrix * m1, const struct matrix * m2)
300
{
301
        struct matrix * diff;
302
        matsize i, j;
303
        matvalue v = 0;
304

    
305
        if (m1 && m2)
306
        {
307
                diff = matrix_difference(m1, m2);
308
                if (diff)
309
                {
310
                        for(i = 0; i < diff->rows; i++)
311
                                for(j = 0; j < diff->cols; j++)
312
                                        v += MAT_GET(diff, i, j) * MAT_GET(diff, i, j);
313
                        v = sqrt(v);
314
                }
315
                matrix_destroy(&diff);
316
        }
317
        return v;
318
}
319
                
320
int matrix_sum(struct matrix  * dst, const struct matrix * op1)
321
{
322
        int res = 0;
323
        matsize i, j;
324

    
325
        if (dst && op1 && dst->cols == op1->cols && dst->rows == op1->rows)
326
        {
327
                for(i = 0; i < dst->rows; i++)
328
                        for(j = 0; j < dst->cols; j++)
329
                                MAT_SET(dst, i, j, MAT_GET(dst, i, j)\
330
                                                + MAT_GET(op1, i, j));
331
        } else
332
                res = -1;
333
        return res;
334
}
335

    
336
struct matrix * matrix_clone(const struct matrix * m)
337
{
338
        struct matrix * c = NULL;
339

    
340
        if(m)
341
        {
342
                c = malloc(sizeof(struct matrix));
343
                matrix_init(c, m->rows, m->cols);
344
                memmove(c->values, m->values, m->rows*m->cols*sizeof(matvalue));
345
        }
346
        return c;
347
}
348
        
349
struct matrix * matrix_eigencentrality(const struct matrix * A, matvalue * dominant_eigenvalue)
350
{
351
        struct matrix *x, *y = NULL, *S, *M = NULL;
352
        float m = 0.15; // decided by Google guys
353
        matvalue l;
354

    
355
        if (A && A->cols == A->rows)
356
        {
357
                //matrix_print(A);
358
                S = matrix_ones(A->rows, A->cols);
359
                matrix_divide(S, A->rows/m);
360
                M = matrix_clone(A);
361
                matrix_divide(M, 1/(1-m));
362
                matrix_sum(M, S);
363
                x = matrix_single_value(A->rows, 1, 0);
364
                MAT_SET(x, 0, 0, 1);
365
                //matrix_print(M);
366

    
367
                y = matrix_multiply(M, x);
368
                l = matrix_sum_values(y);
369
                matrix_divide(y, l);
370
                while (matrix_distance(x, y) > PRECISION)
371
                {
372
                        matrix_destroy(&x);
373
                        x = y;
374
                        y = matrix_multiply(M, x);
375
                        l = matrix_sum_values(y);
376
                        matrix_divide(y, l);
377
                }
378
                matrix_destroy(&x);
379
                matrix_destroy(&S);
380
                matrix_destroy(&M);
381
                //matrix_print(y);
382
                //fprintf(stderr,"[DEBUG] dominant eigenvalue: %f\n", l);
383
                if (dominant_eigenvalue)
384
                        *dominant_eigenvalue = l;
385
        }
386
        return y;
387
}
388

    
389
struct matrix * matrix_diag(const struct matrix * m)
390
{
391
        matsize i;
392
        struct matrix * diag = NULL;
393

    
394
        if (m && m->rows == 1)
395
        {
396
                diag = matrix_single_value(m->cols, m->cols, 0);
397
                for(i = 0; i < diag->cols; i++)
398
                {
399
                        MAT_SET(diag, i, i, MAT_GET(m, 0, i));
400
                }
401
        } else if (m && m->cols == 1)
402
        {
403
                diag = matrix_single_value(m->rows, m->rows, 0);
404
                for(i = 0; i < diag->rows; i++)
405
                        MAT_SET(diag, i, i, MAT_GET(m, i, 0));
406
        }
407
        return diag;
408
}
409

    
410
struct matrix * matrix_invdiag(const struct matrix * m)
411
{
412
        matsize i;
413
        struct matrix * diag = NULL;
414

    
415
        if (m && m->rows == 1)
416
        {
417
                diag = matrix_single_value(m->cols, m->cols, 0);
418
                for(i = 0; i < diag->cols; i++)
419
                        MAT_SET(diag, i, i, 1/MAT_GET(m, 0, i));
420
        } else if (m && m->cols == 1)
421
        {
422
                diag = matrix_single_value(m->rows, m->rows, 0);
423
                for(i = 0; i < diag->rows; i++)
424
                        MAT_SET(diag, i, i, 1/MAT_GET(m, i, 0));
425
        }
426
        return diag;
427
}
428

    
429
int matrix_stochastify(struct matrix *m)
430
        // makes the matrix COLUMN-stochastic
431
{
432
        int res = 0;
433
        struct matrix * u, *sums;
434
        matsize i, j;
435

    
436
        if(m)
437
        {
438
                u = matrix_ones(1, m->rows);
439
                sums = matrix_multiply(u, m);
440
                for(i = 0; i < m->rows; i++)
441
                        for(j = 0; j < m->cols; j++)
442
                                MAT_SET(m, i, j, MAT_GET(m, i, j)/MAT_GET(sums, 0, j));
443

    
444
                matrix_destroy(&u);
445
                matrix_destroy(&sums);
446
        } else
447
                res = -1;
448
        return res;
449
}
450

    
451
void matrix_shrink(struct matrix **ptr, matsize pos)
452
{
453
        struct matrix *neo, *m;
454
        matsize i, j;
455

    
456
        if(ptr && *ptr)
457
        {
458
                m = *ptr;
459
                neo = malloc(sizeof(struct matrix));
460
                if (m->rows > 1 && m->cols > 1) // matrix case
461
                {
462
                        matrix_init(neo, pos < m->rows ? m->rows - 1 : m->rows, pos < m->cols ? m->cols - 1 : m->cols);
463
                        for(i = 0; i < m->rows; i++)
464
                                for(j = 0; j < m->cols; j++)
465
                                {
466
                                        if(i < pos)
467
                                        {
468
                                                if(j < pos) MAT_SET(neo, i, j, MAT_GET(m, i, j));
469
                                                if(j > pos) MAT_SET(neo, i, j-1, MAT_GET(m, i, j));
470
                                        }
471
                                        if(i > pos)
472
                                        {
473
                                                if(j < pos) MAT_SET(neo, i - 1, j, MAT_GET(m, i, j));
474
                                                if(j > pos) MAT_SET(neo, i - 1, j-1, MAT_GET(m, i, j));
475
                                        }
476
                                }
477
                } else // vector case
478
                {
479
                        if(m->rows == 1) // column vector
480
                        {
481
                                matrix_init(neo, m->rows, m->cols - 1);
482
                                for(i = 0; i < m->cols; i++)
483
                                {
484
                                        if(i < pos) MAT_SET(neo, 0, i, MAT_GET(m, 0, i));
485
                                        if(i > pos) MAT_SET(neo, 0, i-1, MAT_GET(m, 0, i));
486
                                }
487
                        }
488
                        if(m->cols == 1) // row vector
489
                        {
490
                                matrix_init(neo, m->rows - 1, m->cols);
491
                                for(i = 0; i < m->rows; i++)
492
                                {
493
                                        if(i < pos) MAT_SET(neo, i, 0, MAT_GET(m, i, 0));
494
                                        if(i > pos) MAT_SET(neo, i-1, 0, MAT_GET(m, i, 0));
495
                                }
496
                        }
497
                }
498
                matrix_destroy(ptr);
499
                *ptr = neo;
500
        }
501
}