Statistics
| Branch: | Tag: | Revision:

sssimulator / wuli.c @ 743e659d

History | View | Annotate | Download (2.64 KB)

1
#include"wuli.h"
2
#include<lp_lib.h>
3

    
4
#define x(n, i, j) (i+j*n + 1)
5
#define b(i, t, s_id) ((i == s_id) ? 1 : (i == t) ? -1 : 0)
6

    
7
int f(matsize n, matsize i, matsize j, int t, int source_id)
8
{
9
        if(t < source_id)
10
                return (n*n + i + j*n + t*n*n + 1);
11
        if(t > source_id)
12
                return (n*n + i + j*n + (t-1)*n*n + 1);
13
        return -1;
14
}
15

    
16
REAL node_output(const struct matrix * A, matsize i)
17
{
18
        matsize n;
19
        REAL res;
20
        struct matrix *u, *product;
21
        n = matrix_num_rows(A);  // number of nodes
22

    
23
        u = matrix_zeros(1, n);
24
        matrix_element_set(u, 0, i, 1);
25

    
26
        product = matrix_multiply(u, A);
27
        res = matrix_sum_values(product);
28

    
29
        matrix_destroy(&u);
30
        matrix_destroy(&product);
31
        return res;
32
}
33

    
34
void overlay_matrix_wuli(struct matrix ** A, struct matrix ** x, int source_id)
35
{
36
        lprec *lp;
37
        matsize n, N, i, j;
38
        int *var_index, k, t;
39
        REAL *coeff, *sol, sum;
40

    
41
        if(!A || !x)
42
                return;
43

    
44
        n = matrix_num_rows(*A);  // number of nodes
45
        N = n*n + n*n*(n-1);  // number of variables
46
        coeff = malloc(sizeof(REAL) * (N+1) );
47
        var_index = malloc(sizeof(int) * (N+1) );
48

    
49
        for(i=1; i < (n*n)+1; i++)
50
                coeff[i] = 1;
51
        for(i=n*n+1; i < N+1; i++)
52
                coeff[i] = 0;
53

    
54
        lp = make_lp(0, N);
55
        set_obj_fn(lp, coeff);
56

    
57
        k=0;
58
        for(i=0; i<n; i++)
59
        {
60
                for(t=0; t<n; t++)
61
                {
62
                        if (t != source_id)
63
                        {
64
                                // flux constraint
65
                                for(k=0, j=0; j<n; j++)
66
                                {
67
                                        if(matrix_element_get(*A, i, j))
68
                                        {
69
                                                var_index[k] = f(n, i, j, t, source_id);
70
                                                coeff[k++] = 1;
71
                                        }
72
                                        if(matrix_element_get(*A, j, i))
73
                                        {
74
                                                var_index[k] = f(n, j, i, t, source_id);
75
                                                coeff[k++] = -1;
76
                                        }
77
                                }
78
                                if (k)
79
                                        add_constraintex(lp, k, coeff, var_index, EQ, b(i, t, source_id));
80
                                // f_ij <= x_ij constraint
81
                                for(j=0; j<n; j++)
82
                                {
83
                                        k = 0;
84
                                        var_index[k] = f(n, i, j, t, source_id);
85
                                        coeff[k++] = 1;
86
                                        var_index[k] = x(n, i, j);
87
                                        coeff[k++] = -1;
88
                                        if (k)
89
                                                add_constraintex(lp, k, coeff, var_index, LE, 0);
90
                                }
91
                        }
92
                }
93
                // node output constraint
94
                for(k=0, j=0; j<n; j++)
95
                {
96
                        if(matrix_element_get(*A, i, j))
97
                        {
98
                                var_index[k] = x(n, i, j);
99
                                coeff[k++] = 1;
100
                        }
101
                }
102
                if (k)
103
                        add_constraintex(lp, k, coeff, var_index, LE, node_output(*A, i));
104
                for(j=0; j<n; j++)
105
                        if(!matrix_element_get(*A, i, j))
106
                                set_upbo(lp, x(n, i, j), 0);
107
        }
108

    
109
        k = solve(lp);
110
        get_ptr_variables(lp, &sol);
111

    
112
        if (k == 0)
113
        {
114
                for(i=0; i < n; i++)
115
                {
116
                        for(sum=0, j=0; j < n; j++)
117
                                sum += sol[x(n, i, j)-1];
118
                        matrix_element_set(*x, i, 0, sum);
119

    
120
                        for(j=0; j < n; j++)
121
                                if (sum)
122
                                        matrix_element_set(*A, i, j, sol[x(n, i, j)-1]/sum);
123
                                else
124
                                        matrix_element_set(*A, i, j, 0);
125
                }
126
        }
127
        else
128
                for(i=0; i < n; i++)
129
                        matrix_element_set(*x, i, 0, -1);
130

    
131
        free(coeff);
132
        free(var_index);
133
        delete_lp(lp);
134
}