Ticket #43269: sparceReduced.cc

File sparceReduced.cc, 15.3 KB (added by py@…, 10 years ago)

g++ source

Line 
1
2#include "sparse.h"
3#include "PrismSparseGlob.h"
4#include <new>
5
6//------------------------------------------------------------------------------
7
8// local function prototypes
9
10static void split_mdp_rec(DdManager *ddman, DdNode *dd, DdNode **ndvars, int num_ndvars, int level, DdNode **matrices);
11static void split_mdp_and_sub_mdp_rec(DdManager *ddman, DdNode *dd, DdNode *subdd, DdNode **ndvars, int num_ndvars, int level, DdNode **matrices, DdNode **submatrices);
12static void traverse_mtbdd_matr_rec(DdManager *ddman, DdNode *dd, DdNode **rvars, DdNode **cvars, int num_vars, int level, ODDNode *row, ODDNode *col, int r, int c, int code, bool transpose);
13static void traverse_mtbdd_vect_rec(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, int level, ODDNode *odd, int i, int code);
14
15// global variables (used by local functions)
16static int count;
17static int *starts, *starts2;
18static int *actions;
19static RMSparseMatrix *rmsm;
20static CMSparseMatrix *cmsm;
21static RCSparseMatrix *rcsm;
22static CMSRSparseMatrix *cmsrsm;
23static CMSCSparseMatrix *cmscsm;
24static NDSparseMatrix *ndsm;
25
26
27//-----------------------------------------------------------------------------------
28
29// build nondeterministic (mdp) sparse matrix
30// throws std::bad_alloc on out-of-memory
31
32NDSparseMatrix *build_nd_sparse_matrix(DdManager *ddman, DdNode *mdp, DdNode **rvars, DdNode **cvars, int num_vars, DdNode **ndvars, int num_ndvars, ODDNode *odd)
33{
34        int i, n, nm, nc, nnz, max, max2;
35        DdNode *tmp = NULL, **matrices = NULL, **matrices_bdds = NULL;
36       
37        // try/catch for memory allocation/deallocation
38        try {
39       
40        // create new data structure
41        ndsm = NULL; ndsm = new NDSparseMatrix();
42       
43        // get number of states from odd
44        n = ndsm->n = odd->eoff+odd->toff;
45        // get num of choices (prob. distributions)
46        Cudd_Ref(mdp);
47        tmp = DD_ThereExists(ddman, DD_Not(ddman, DD_Equals(ddman, mdp, 0)), cvars, num_vars);
48        nc = ndsm->nc = (int)DD_GetNumMinterms(ddman, tmp, num_vars+num_ndvars);
49        // get num of transitions
50        nnz = ndsm->nnz = (int)DD_GetNumMinterms(ddman, mdp, num_vars*2+num_ndvars);
51        // break the mdp mtbdd into several (nm) mtbdds
52        tmp = DD_ThereExists(ddman, tmp, rvars, num_vars);
53        nm = (int)DD_GetNumMinterms(ddman, tmp, num_ndvars);
54        Cudd_RecursiveDeref(ddman, tmp);
55        matrices = new DdNode*[nm];
56        count = 0;
57        split_mdp_rec(ddman, mdp, ndvars, num_ndvars, 0, matrices);
58        // and for each one create a bdd storing which rows/choices are non-empty
59        matrices_bdds = new DdNode*[nm];
60        for (i = 0; i < nm; i++) {
61                Cudd_Ref(matrices[i]);
62                matrices_bdds[i] = DD_ThereExists(ddman, DD_Not(ddman, DD_Equals(ddman, matrices[i], 0)), cvars, num_vars);
63        }
64       
65        // create arrays
66        ndsm->non_zeros = new double[nnz];
67        ndsm->cols = new unsigned int[nnz];
68        starts = NULL; starts = new int[n+1];
69        starts2 = NULL; starts2 = new int[nc+1];
70       
71        // first traverse mtbdds to compute how many choices are in each row
72        for (i = 0; i < n+1; i++) starts[i] = 0;
73        for (i = 0; i < nm; i++) {
74                traverse_mtbdd_vect_rec(ddman, matrices_bdds[i], rvars, num_vars, 0, odd, 0, 1);
75        }
76        // and use this to compute the starts information
77        // (and at same time, compute max num choices in a state)
78        max = 0;
79        for (i = 1 ; i < n+1; i++) {
80                if (starts[i] > max) max = starts[i];
81                starts[i] += starts[i-1];
82        }
83        ndsm->k = max;
84       
85        // now traverse mtbdds to compute how many transitions in each choice
86        for (i = 0; i < nc+1; i++) starts2[i] = 0;
87        for (i = 0; i < nm; i++) {
88                traverse_mtbdd_matr_rec(ddman, matrices[i], rvars, cvars, num_vars, 0, odd, odd, 0, 0, 10, false);
89                traverse_mtbdd_vect_rec(ddman, matrices_bdds[i], rvars, num_vars, 0, odd, 0, 2);
90        }
91        // and use this to compute the starts2 information
92        // (and at same time, compute max num transitions in a choice)
93        max2 = 0;
94        for (i = 1; i < nc+1; i++) {
95                if (starts2[i] > max2) max2 = starts2[i];
96                starts2[i] += starts2[i-1];
97        }
98        // recompute starts (because we altered them during last traversal)
99        for (i = n; i > 0; i--) {
100                starts[i] = starts[i-1];
101        }
102        starts[0] = 0;
103       
104        // max num choices/transitions determines whether we store counts or starts:
105        ndsm->use_counts = (max < (unsigned int)(1 << (8*sizeof(unsigned char))));
106        ndsm->use_counts &= (max2 < (unsigned int)(1 << (8*sizeof(unsigned char))));
107       
108        // now traverse the mtbdd again to get the actual matrix entries
109        for (i = 0; i < nm; i++) {
110                traverse_mtbdd_matr_rec(ddman, matrices[i], rvars, cvars, num_vars, 0, odd, odd, 0, 0, 11, false);
111                traverse_mtbdd_vect_rec(ddman, matrices_bdds[i], rvars, num_vars, 0, odd, 0, 2);
112        }
113        // recompute starts (because we altered them during last traversal)
114        for (i = n; i > 0; i--) {
115                starts[i] = starts[i-1];
116        }
117        starts[0] = 0;
118        // recompute starts2 (likewise)
119        for (i = nc; i > 0; i--) {
120                starts2[i] = starts2[i-1];
121        }
122        starts2[0] = 0;
123       
124        // if it's safe to do so, replace starts/starts2 with (smaller) arrays of counts
125        if (ndsm->use_counts) {
126                ndsm->row_counts = new unsigned char[n];
127                for (i = 0; i < n; i++) ndsm->row_counts[i] = (unsigned char)(starts[i+1] - starts[i]);
128                delete[] starts; starts = NULL;
129                ndsm->choice_counts = new unsigned char[nc];
130                for (i = 0; i < nc; i++) ndsm->choice_counts[i] = (unsigned char)(starts2[i+1] - starts2[i]);
131                delete[] starts2; starts2 = NULL;
132                ndsm->mem = (nnz * (sizeof(double) + sizeof(unsigned int)) + (n+nc) * sizeof(unsigned char)) / 1024.0;
133        } else {
134                ndsm->row_counts = (unsigned char*)starts;
135                ndsm->choice_counts = (unsigned char*)starts2;
136                ndsm->mem = (nnz * (sizeof(double) + sizeof(unsigned int)) + (n+nc) * sizeof(int)) / 1024.0;
137        }
138       
139        // try/catch for memory allocation/deallocation
140        } catch(std::bad_alloc e) {
141                if (ndsm) delete ndsm;
142                if (matrices) delete[] matrices;
143                if (matrices_bdds) {
144                        for (i = 0; i < nm; i++) Cudd_RecursiveDeref(ddman, matrices_bdds[i]);
145                        delete[] matrices_bdds;
146                }
147                if (starts) delete[] starts;
148                if (starts2) delete[] starts2;
149                throw e;
150        }
151       
152        // clear up memory
153        for (i = 0; i < nm; i++) {
154                Cudd_RecursiveDeref(ddman, matrices_bdds[i]);
155                // nb: don't deref matrices array because that was just pointers, not new copies
156        }
157        delete[] matrices;
158        delete[] matrices_bdds;
159       
160        return ndsm;
161}
162
163//-----------------------------------------------------------------------------------
164
165// Build nondeterministic (MDP) sparse matrix for "sub-MDP".
166// This function basically exists to construct a sparse matrix representing the transition rewards for an MDP.
167// The complication is that we need to use the nondeterministic choice indexing of the main
168// MDP matrix, not the rewards matrix, otherwise we can't tell which reward is on which transition.
169// throws std::bad_alloc on out-of-memory
170
171NDSparseMatrix *build_sub_nd_sparse_matrix(DdManager *ddman, DdNode *mdp, DdNode *submdp, DdNode **rvars, DdNode **cvars, int num_vars, DdNode **ndvars, int num_ndvars, ODDNode *odd)
172{
173        int i, n, nm, nc, nnz, max, max2;
174        DdNode *tmp = NULL, **matrices = NULL, **submatrices = NULL, **matrices_bdds = NULL;
175       
176        // try/catch for memory allocation/deallocation
177        try {
178       
179        // create new data structure
180        ndsm = NULL; ndsm = new NDSparseMatrix();
181       
182        // get number of states from odd
183        n = ndsm->n = odd->eoff+odd->toff;
184        // get num of choices (prob. distributions) (USING MDP)
185        Cudd_Ref(mdp);
186        tmp = DD_ThereExists(ddman, DD_Not(ddman, DD_Equals(ddman, mdp, 0)), cvars, num_vars);
187        nc = ndsm->nc = (int)DD_GetNumMinterms(ddman, tmp, num_vars+num_ndvars);
188        // get num of transitions (USING SUB-MDP)
189        nnz = ndsm->nnz = (int)DD_GetNumMinterms(ddman, submdp, num_vars*2+num_ndvars);
190        // break the two mtbdds (MDP AND SUB-MDP) into several (nm) mtbdds
191        tmp = DD_ThereExists(ddman, tmp, rvars, num_vars);
192        nm = (int)DD_GetNumMinterms(ddman, tmp, num_ndvars);
193        Cudd_RecursiveDeref(ddman, tmp);
194        matrices = new DdNode*[nm];
195        submatrices = new DdNode*[nm];
196        count = 0;
197        split_mdp_and_sub_mdp_rec(ddman, mdp, submdp, ndvars, num_ndvars, 0, matrices, submatrices);
198        // and for each one create a bdd storing which rows/choices are non-empty (USING MDP)
199        matrices_bdds = new DdNode*[nm];
200        for (i = 0; i < nm; i++) {
201                Cudd_Ref(matrices[i]);
202                matrices_bdds[i] = DD_ThereExists(ddman, DD_Not(ddman, DD_Equals(ddman, matrices[i], 0)), cvars, num_vars);
203        }
204       
205        // create arrays
206        ndsm->non_zeros = new double[nnz];
207        ndsm->cols = new unsigned int[nnz];
208        starts = NULL; starts = new int[n+1];
209        starts2 = NULL; starts2 = new int[nc+1];
210       
211        // first traverse mtbdds to compute how many choices are in each row (USING MDP)
212        for (i = 0; i < n+1; i++) starts[i] = 0;
213        for (i = 0; i < nm; i++) {
214                traverse_mtbdd_vect_rec(ddman, matrices_bdds[i], rvars, num_vars, 0, odd, 0, 1);
215        }
216        // and use this to compute the starts information
217        // (and at same time, compute max num choices in a state)
218        max = 0;
219        for (i = 1 ; i < n+1; i++) {
220                if (starts[i] > max) max = starts[i];
221                starts[i] += starts[i-1];
222        }
223        ndsm->k = max;
224       
225        // now traverse mtbdds to compute how many transitions in each choice (USING SUB-MDP)
226        for (i = 0; i < nc+1; i++) starts2[i] = 0;
227        for (i = 0; i < nm; i++) {
228                traverse_mtbdd_matr_rec(ddman, submatrices[i], rvars, cvars, num_vars, 0, odd, odd, 0, 0, 10, false);
229                traverse_mtbdd_vect_rec(ddman, matrices_bdds[i], rvars, num_vars, 0, odd, 0, 2);
230        }
231        // and use this to compute the starts2 information
232        // (and at same time, compute max num transitions in a choice)
233        max2 = 0;
234        for (i = 1; i < nc+1; i++) {
235                if (starts2[i] > max2) max2 = starts2[i];
236                starts2[i] += starts2[i-1];
237        }
238        // recompute starts (because we altered them during last traversal)
239        for (i = n; i > 0; i--) {
240                starts[i] = starts[i-1];
241        }
242        starts[0] = 0;
243       
244        // max num choices/transitions determines whether we store counts or starts:
245        ndsm->use_counts = (max < (unsigned int)(1 << (8*sizeof(unsigned char))));
246        ndsm->use_counts &= (max2 < (unsigned int)(1 << (8*sizeof(unsigned char))));
247       
248        // now traverse the mtbdd again to get the actual matrix entries (USING SUB-MDP)
249        for (i = 0; i < nm; i++) {
250                traverse_mtbdd_matr_rec(ddman, submatrices[i], rvars, cvars, num_vars, 0, odd, odd, 0, 0, 11, false);
251                traverse_mtbdd_vect_rec(ddman, matrices_bdds[i], rvars, num_vars, 0, odd, 0, 2);
252        }
253        // recompute starts (because we altered them during last traversal)
254        for (i = n; i > 0; i--) {
255                starts[i] = starts[i-1];
256        }
257        starts[0] = 0;
258        // recompute starts2 (likewise)
259        for (i = nc; i > 0; i--) {
260                starts2[i] = starts2[i-1];
261        }
262        starts2[0] = 0;
263       
264        // if it's safe to do so, replace starts/starts2 with (smaller) arrays of counts
265        if (ndsm->use_counts) {
266                ndsm->row_counts = new unsigned char[n];
267                for (i = 0; i < n; i++) ndsm->row_counts[i] = (unsigned char)(starts[i+1] - starts[i]);
268                delete[] starts; starts = NULL;
269                ndsm->choice_counts = new unsigned char[nc];
270                for (i = 0; i < nc; i++) ndsm->choice_counts[i] = (unsigned char)(starts2[i+1] - starts2[i]);
271                delete[] starts2; starts2 = NULL;
272                ndsm->mem = (nnz * (sizeof(double) + sizeof(unsigned int)) + (n+nc) * sizeof(unsigned char)) / 1024.0;
273        } else {
274                ndsm->row_counts = (unsigned char*)starts;
275                ndsm->choice_counts = (unsigned char*)starts2;
276                ndsm->mem = (nnz * (sizeof(double) + sizeof(unsigned int)) + (n+nc) * sizeof(int)) / 1024.0;
277        }
278       
279        // try/catch for memory allocation/deallocation
280        } catch(std::bad_alloc e) {
281                if (ndsm) delete ndsm;
282                if (matrices) delete[] matrices;
283                if (submatrices) delete[] submatrices;
284                if (matrices_bdds) {
285                        for (i = 0; i < nm; i++) Cudd_RecursiveDeref(ddman, matrices_bdds[i]);
286                        delete[] matrices_bdds;
287                }
288                if (starts) delete[] starts;
289                if (starts2) delete[] starts2;
290                throw e;
291        }
292       
293        // clear up memory
294        for (i = 0; i < nm; i++) {
295                Cudd_RecursiveDeref(ddman, matrices_bdds[i]);
296                // nb: don't deref matrices/submatrices array because that was just pointers, not new copies
297        }
298        delete[] matrices;
299        delete[] submatrices;
300        delete[] matrices_bdds;
301       
302        return ndsm;
303}
304
305//------------------------------------------------------------------------------
306
307// traverses the mtbdd and gets all the MATRIX entries out
308// does different things depending on the value of 'code'
309// if tranpose flag is true, actually extract from tranpose of matrix
310
311void traverse_mtbdd_matr_rec(DdManager *ddman, DdNode *dd, DdNode **rvars, DdNode **cvars, int num_vars, int level, ODDNode *row, ODDNode *col, int r, int c, int code, bool transpose)
312{
313        DdNode *e, *t, *ee, *et, *te, *tt;
314        int i, dist_num;
315        double *dist, d;
316       
317        // base case - zero terminal
318        if (dd == Cudd_ReadZero(ddman)) return;
319       
320        // base case - non zero terminal
321        if (level == num_vars) {
322                switch (code) {
323               
324                // row major - first pass
325                case 1:
326                        starts[(transpose?c:r)+1]++;
327                        break;
328                       
329                // row major - second pass
330                case 2:
331                        rmsm->non_zeros[starts[(transpose?c:r)]] = Cudd_V(dd);
332                        rmsm->cols[starts[(transpose?c:r)]] = (transpose?r:c);
333                        starts[(transpose?c:r)]++;
334                        break;
335                       
336                // column major - first pass
337                case 3:
338                        starts[(transpose?r:c)+1]++;
339                        break;
340                       
341                // column major - second pass
342                case 4:
343                        cmsm->non_zeros[starts[(transpose?r:c)]] = Cudd_V(dd);
344                        cmsm->rows[starts[(transpose?r:c)]] = (transpose?c:r);
345                        starts[(transpose?r:c)]++;
346                        break;
347                       
348                // row/column - only pass
349                case 5:
350                        rcsm->non_zeros[count] = Cudd_V(dd);
351                        rcsm->rows[count] = (transpose?c:r);
352                        rcsm->cols[count] = (transpose?r:c);
353                        count++;
354                        break;
355                       
356                // compact modified sparse row - first pass
357                case 6:
358                        starts[(transpose?c:r)+1]++;
359                        break;
360                       
361                // compact modified sparse row - second pass
362                case 7:
363                        // try and find value
364                        dist = cmsrsm->dist;
365                        dist_num = cmsrsm->dist_num;
366                        d = Cudd_V(dd);
367                        for (i = 0; i < dist_num; i++) if (dist[i] == d) break;
368                        // if it's not there, add it
369                        if (i == dist_num) {
370                                dist[i] = d;
371                                cmsrsm->dist_num++;
372                        }
373                        // store info
374                        cmsrsm->cols[starts[(transpose?c:r)]] = (unsigned int)(((unsigned int)(transpose?r:c) << cmsrsm->dist_shift) + (unsigned int)i);
375                        starts[(transpose?c:r)]++;
376                        break;
377                       
378                // compact modified sparse column - first pass
379                case 8:
380                        starts[(transpose?r:c)+1]++;
381                        break;
382                       
383                // compact modified sparse column - second pass
384                case 9:
385                        // try and find value
386                        dist = cmscsm->dist;
387                        dist_num = cmscsm->dist_num;
388                        d = Cudd_V(dd);
389                        for (i = 0; i < dist_num; i++) if (dist[i] == d) break;
390                        // if it's not there, add it
391                        if (i == dist_num) {
392                                dist[i] = d;
393                                cmscsm->dist_num++;
394                        }
395                        // store info
396                        cmscsm->rows[starts[(transpose?r:c)]] = (unsigned int)(((unsigned int)(transpose?c:r) << cmscsm->dist_shift) + (unsigned int)i);
397                        starts[(transpose?r:c)]++;
398                        break;
399                       
400                // mdp - first pass
401                case 10:
402                        starts2[starts[(transpose?c:r)]+1]++;
403                        break;
404                       
405                // mdp - second pass
406                case 11:
407                        ndsm->non_zeros[starts2[starts[(transpose?c:r)]]] = Cudd_V(dd);
408                        ndsm->cols[starts2[starts[(transpose?c:r)]]] = (transpose?r:c);
409                        starts2[starts[(transpose?c:r)]]++;
410                        break;
411                }
412                return;
413        }
414       
415        // recurse
416        if (dd->index > cvars[level]->index) {
417                ee = et = te = tt = dd;
418        }
419        else if (dd->index > rvars[level]->index) {
420                ee = te = Cudd_E(dd);
421                et = tt = Cudd_T(dd);
422        }
423        else {
424                e = Cudd_E(dd);
425                if (e->index > cvars[level]->index) {
426                        ee = et = e;
427                }
428                else {
429                        ee = Cudd_E(e);
430                        et = Cudd_T(e);
431                }
432                t = Cudd_T(dd);
433                if (t->index > cvars[level]->index) {
434                        te = tt = t;
435                }
436                else {
437                        te = Cudd_E(t);
438                        tt = Cudd_T(t);
439                }
440        }
441
442        traverse_mtbdd_matr_rec(ddman, ee, rvars, cvars, num_vars, level+1, row->e, col->e, r, c, code, transpose);
443        traverse_mtbdd_matr_rec(ddman, et, rvars, cvars, num_vars, level+1, row->e, col->t, r, c+col->eoff, code, transpose);
444        traverse_mtbdd_matr_rec(ddman, te, rvars, cvars, num_vars, level+1, row->t, col->e, r+row->eoff, c, code, transpose);
445        traverse_mtbdd_matr_rec(ddman, tt, rvars, cvars, num_vars, level+1, row->t, col->t, r+row->eoff, c+col->eoff, code, transpose);
446}
447