← Back to branch summary

~tapaal-ltl/verifypn/scc-optimise

« back to all changes in this revision

Viewing changes to src/PetriEngine/Simplification/LinearProgram.cpp

  • Committer: srba.jiri at gmail
  • Date: 2020-09-11 14:23:39 UTC
  • mfrom: (213.1.151 interval_tar)
  • Revision ID: srba.jiri@gmail.com-20200911142339-bq9328s1gppw24uj
merged in lp:~verifypn-maintainers/verifypn/interval_tar doing 
- Implements TAR w/o z3, but using a simple integer inference engine for Hoare logic.
 - Replaces LP-Solve with GLPK, reduces computation-time and memory overhead
 - Implements new global properties, translated into CTL formulae.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/PetriEngine/Simplification/LinearProgram.cpp
 
1
#include <cassert>
 
2
#include <cmath>
 
3
#include <glpk.h>
 
4
#include <fstream>
 
5
 
 
6
#include "PetriEngine/Simplification/LinearProgram.h"
 
7
#include "PetriEngine/Simplification/LPCache.h"
 
8
#include "PetriEngine/PQL/Contexts.h"
 
9
 
 
10
namespace PetriEngine {
 
11
    namespace Simplification {
 
12
        using REAL = double;
 
13
        LinearProgram::~LinearProgram() {
 
14
        }
 
15
 
 
16
        LinearProgram::LinearProgram(Vector* vec, int constant, op_t op, LPCache* factory){
 
17
            // TODO fix memory-management here!
 
18
            equation_t c;
 
19
            switch(op)
 
20
            {
 
21
                case Simplification::OP_LT:
 
22
                    c.upper = constant - 1;
 
23
                    break;
 
24
                case Simplification::OP_GT:
 
25
                    c.lower = constant + 1;
 
26
                    break;
 
27
                case Simplification::OP_LE:
 
28
                    c.upper = constant;
 
29
                    break;
 
30
                case Simplification::OP_GE:
 
31
                    c.lower = constant;
 
32
                    break;
 
33
                case Simplification::OP_EQ:
 
34
                    c.lower = constant;
 
35
                    c.upper = constant;
 
36
                    break;
 
37
                default:
 
38
                    // We ignore this operator for now by not adding any equation.
 
39
                    // This is untested, however.
 
40
                    assert(false);
 
41
            }
 
42
            c.row = vec;
 
43
            vec->inc();
 
44
            _equations.push_back(c);
 
45
        }
 
46
 
 
47
        constexpr auto infty = std::numeric_limits<REAL>::infinity();
 
48
        
 
49
        bool LinearProgram::isImpossible(const PQL::SimplificationContext& context, uint32_t solvetime) {
 
50
            bool use_ilp = true;
 
51
            auto net = context.net();
 
52
 
 
53
 
 
54
            if(_result != result_t::UKNOWN)
 
55
            {
 
56
                if(_result == result_t::IMPOSSIBLE)
 
57
                    return _result == result_t::IMPOSSIBLE;
 
58
            }
 
59
 
 
60
            if(_equations.size() == 0 || context.timeout()){
 
61
                return false;
 
62
            }
 
63
 
 
64
            const uint32_t nCol = net->numberOfTransitions();
 
65
            const uint32_t nRow = net->numberOfPlaces() + _equations.size();
 
66
 
 
67
            std::vector<REAL> row = std::vector<REAL>(nCol + 1);
 
68
            std::vector<int32_t> indir(std::max(nCol, nRow) + 1);
 
69
            for(size_t i = 0; i <= nCol; ++i)
 
70
                indir[i] = i;
 
71
 
 
72
            auto lp = context.makeBaseLP();
 
73
            if(lp == nullptr)
 
74
                return false;
 
75
            
 
76
            int rowno = 1 + net->numberOfPlaces();
 
77
            glp_add_rows(lp, _equations.size());
 
78
            for(const auto& eq : _equations){
 
79
                auto l = eq.row->write_indir(row, indir);
 
80
                assert(!(std::isinf(eq.upper) && std::isinf(eq.lower)));                
 
81
                glp_set_mat_row(lp, rowno, l-1, indir.data(), row.data());
 
82
                if(!std::isinf(eq.lower) && !std::isinf(eq.upper))
 
83
                {
 
84
                    if(eq.lower == eq.upper)
 
85
                        glp_set_row_bnds(lp, rowno, GLP_FX, eq.lower, eq.upper);
 
86
                    else
 
87
                    {
 
88
                        if(eq.lower > eq.upper)
 
89
                        {
 
90
                            _result = result_t::IMPOSSIBLE;
 
91
                            glp_delete_prob(lp);
 
92
                            return true;
 
93
                        }
 
94
                        glp_set_row_bnds(lp, rowno, GLP_DB, eq.lower, eq.upper);
 
95
                    }
 
96
                }
 
97
                else if(std::isinf(eq.lower))
 
98
                    glp_set_row_bnds(lp, rowno, GLP_UP, -infty, eq.upper);
 
99
                else
 
100
                    glp_set_row_bnds(lp, rowno, GLP_LO, eq.lower, -infty);
 
101
                ++rowno;
 
102
 
 
103
                if(context.timeout())
 
104
                {
 
105
                    std::cerr << "glpk: construction timeout" << std::endl;
 
106
                    glp_delete_prob(lp);
 
107
                    return false;
 
108
                }
 
109
            }
 
110
            
 
111
            // Set objective, kind and bounds
 
112
            for(size_t i = 1; i <= nCol; i++) {
 
113
                glp_set_obj_coef(lp, i, 0);
 
114
                glp_set_col_kind(lp, i, use_ilp ? GLP_IV :GLP_CV);
 
115
                glp_set_col_bnds(lp, i, GLP_LO, 0, infty);
 
116
            }
 
117
 
 
118
            // Minimize the objective
 
119
            glp_set_obj_dir(lp, GLP_MIN);
 
120
            auto stime = glp_time();
 
121
            glp_smcp settings;
 
122
            glp_init_smcp(&settings);
 
123
            auto timeout = std::min(solvetime, context.getLpTimeout())*1000;
 
124
            settings.tm_lim = timeout;
 
125
            settings.presolve = GLP_OFF;
 
126
            settings.msg_lev = 0;
 
127
            auto result = glp_simplex(lp, &settings);
 
128
            if (result == GLP_ETMLIM)
 
129
            {
 
130
                _result = result_t::UKNOWN;
 
131
                std::cerr << "glpk: timeout" << std::endl;
 
132
            }
 
133
            else if(result == 0)
 
134
            {
 
135
                auto status = glp_get_status(lp);
 
136
                if(status == GLP_OPT) {
 
137
                    glp_iocp iset;
 
138
                    glp_init_iocp(&iset);
 
139
                    iset.msg_lev = 0;
 
140
                    iset.tm_lim = std::max<uint32_t>(timeout-(stime - glp_time()), 1);
 
141
                    iset.presolve = GLP_OFF;
 
142
                    auto ires = glp_intopt(lp, &iset);
 
143
                    if(ires == GLP_ETMLIM)
 
144
                    {
 
145
                        _result = result_t::UKNOWN;
 
146
                        std::cerr << "glpk mip: timeout" << std::endl;
 
147
                    }
 
148
                    else if(ires == 0)
 
149
                    {
 
150
                        auto ist = glp_mip_status(lp);
 
151
                        if(ist == GLP_OPT || ist == GLP_FEAS || ist == GLP_UNBND) {
 
152
                            _result = result_t::POSSIBLE;
 
153
                        }
 
154
                        else
 
155
                        {
 
156
                            _result = result_t::IMPOSSIBLE;
 
157
                        }
 
158
 
 
159
                    }
 
160
                }
 
161
                else if(status == GLP_FEAS || status == GLP_UNBND)
 
162
                {
 
163
                    _result = result_t::POSSIBLE;
 
164
                }
 
165
                else
 
166
                    _result = result_t::IMPOSSIBLE;
 
167
            }
 
168
            else if (result == GLP_ENOPFS || result == GLP_ENODFS || result == GLP_ENOFEAS)
 
169
            {
 
170
                _result = result_t::IMPOSSIBLE;
 
171
            }
 
172
            glp_delete_prob(lp);
 
173
 
 
174
            return _result == result_t::IMPOSSIBLE;
 
175
        }
 
176
 
 
177
        std::vector<std::pair<double,bool>> LinearProgram::bounds(const PQL::SimplificationContext& context, uint32_t solvetime, const std::vector<uint32_t>& places)
 
178
        {
 
179
            std::vector<std::pair<double,bool>> result(places.size() + 1, std::make_pair(std::numeric_limits<double>::infinity(), false));
 
180
            auto net = context.net();
 
181
            auto m0 = context.marking();
 
182
            auto timeout = std::min(solvetime, context.getLpTimeout());
 
183
 
 
184
            const uint32_t nCol = net->numberOfTransitions();
 
185
            std::vector<REAL> row = std::vector<REAL>(nCol + 1);
 
186
 
 
187
 
 
188
            glp_smcp settings;
 
189
            glp_init_smcp(&settings);
 
190
            settings.tm_lim = timeout*1000;
 
191
            settings.presolve = GLP_OFF;
 
192
            settings.msg_lev = 0;
 
193
 
 
194
            for(size_t it = 0; it <= places.size(); ++it)
 
195
            {
 
196
                // we want to start with the overall bound, most important
 
197
                // Spend time on rest after
 
198
                auto stime = glp_time();
 
199
                size_t pi;
 
200
                if(it == 0)
 
201
                    pi = places.size();
 
202
                else
 
203
                    pi = it - 1;
 
204
 
 
205
                if(context.timeout())
 
206
                {
 
207
                    return result;
 
208
                }
 
209
                // Create objective
 
210
                memset(row.data(), 0, sizeof (REAL) * net->numberOfTransitions() + 1);
 
211
                double p0 = 0;
 
212
                bool all_le_zero = true;
 
213
                bool all_zero = true;
 
214
                if(pi < places.size())
 
215
                {
 
216
                    auto tp = places[pi];
 
217
                    p0 = m0[tp];
 
218
                    for (size_t t = 0; t < net->numberOfTransitions(); ++t)
 
219
                    {
 
220
                        row[1 + t] = net->outArc(t, tp) - net->inArc(tp, t);
 
221
                        all_le_zero &= row[1 + t] <= 0;
 
222
                        all_zero &= row[1 + t] == 0;
 
223
                    }
 
224
                }
 
225
                else
 
226
                {
 
227
                    for (size_t t = 0; t < net->numberOfTransitions(); ++t)
 
228
                    {
 
229
                        double cnt = 0;
 
230
                        for(auto tp : places)
 
231
                            cnt += net->outArc(t, tp) - net->inArc(tp, t);
 
232
                        row[1 + t] = cnt;
 
233
                        all_le_zero &= row[1 + t] <= 0;
 
234
                        all_zero &= row[1 + t] == 0;
 
235
                    }
 
236
                    for(auto tp : places)
 
237
                        p0 += m0[tp];
 
238
                }
 
239
 
 
240
                if(all_le_zero)
 
241
                {
 
242
                    result[pi].first = p0;
 
243
                    result[pi].second = all_zero;
 
244
                    if(pi == places.size())
 
245
                    {
 
246
                        return result;
 
247
                    }
 
248
                    continue;
 
249
                }
 
250
 
 
251
                // Set objective
 
252
 
 
253
                auto tmp_lp = context.makeBaseLP();
 
254
                if(tmp_lp == nullptr)
 
255
                    return result;
 
256
            
 
257
                // Max the objective
 
258
                glp_set_obj_dir(tmp_lp, GLP_MAX);
 
259
 
 
260
                for(size_t i = 1; i <= nCol; i++) {
 
261
                    glp_set_obj_coef(tmp_lp, i, row[i]);
 
262
                    glp_set_col_kind(tmp_lp, i, GLP_IV);
 
263
                    glp_set_col_bnds(tmp_lp, i, GLP_LO, 0, infty);
 
264
                }
 
265
 
 
266
                auto rs = glp_simplex(tmp_lp, &settings);
 
267
                if (rs == GLP_ETMLIM)
 
268
                {
 
269
                    std::cerr << "glpk: timeout" << std::endl;
 
270
                }
 
271
                else if(rs == 0)
 
272
                {
 
273
                    auto status = glp_get_status(tmp_lp);
 
274
                    if(status == GLP_OPT) {
 
275
                        glp_iocp isettings;
 
276
                        glp_init_iocp(&isettings);
 
277
                        isettings.tm_lim = std::max<int>(((double) timeout * 1000) - (glp_time() - stime), 1);
 
278
                        isettings.msg_lev = 0;
 
279
                        isettings.presolve = GLP_OFF;
 
280
                        auto rs = glp_intopt(tmp_lp, &isettings);
 
281
                        if (rs == GLP_ETMLIM) {
 
282
                            std::cerr << "glpk mip: timeout" << std::endl;
 
283
                        } else if (rs == 0) {
 
284
                            auto status = glp_mip_status(tmp_lp);
 
285
                            if (status == GLP_OPT) {
 
286
                                auto org = p0 + glp_mip_obj_val(tmp_lp);
 
287
                                result[pi].first = std::round(org);
 
288
                                result[pi].second = all_zero;
 
289
                            }
 
290
                            else if (status != GLP_UNBND && status != GLP_FEAS)
 
291
                            {
 
292
                                result[pi].first = p0;
 
293
                                result[pi].second = all_zero;
 
294
                            }
 
295
                        }
 
296
                    }
 
297
                    else if (status == GLP_INFEAS || status == GLP_NOFEAS || status == GLP_UNDEF)
 
298
                    {
 
299
                        result[pi].first = p0;
 
300
                        result[pi].second = all_zero;
 
301
                    }
 
302
                }
 
303
                else if (rs == GLP_ENOPFS || rs == GLP_ENODFS || rs == GLP_ENOFEAS)
 
304
                {
 
305
                    result[pi].first = p0;
 
306
                    result[pi].second = all_zero;
 
307
                }
 
308
                glp_erase_prob(tmp_lp);
 
309
                if(pi == places.size() && result[places.size()].first >= p0)
 
310
                {
 
311
                    return result;
 
312
                }
 
313
                if(pi == places.size() && places.size() == 1)
 
314
                {
 
315
                    result[0] = result[1];
 
316
                    return result;
 
317
                }
 
318
            }
 
319
            return result;
 
320
        }
 
321
 
 
322
 
 
323
        void LinearProgram::make_union(const LinearProgram& other)
 
324
        {
 
325
            if(_result == IMPOSSIBLE || other._result == IMPOSSIBLE)
 
326
            {
 
327
                _result = IMPOSSIBLE;
 
328
                _equations.clear();
 
329
                assert(false);
 
330
                return;
 
331
            }
 
332
 
 
333
            auto it1 = _equations.begin();
 
334
            auto it2 = other._equations.begin();
 
335
 
 
336
            while(it1 != _equations.end() && it2 != other._equations.end())
 
337
            {
 
338
                if(it1->row < it2->row)
 
339
                {
 
340
                    ++it1;
 
341
                }
 
342
                else if(it2->row < it1->row)
 
343
                {
 
344
                    it1 = _equations.insert(it1, *it2);
 
345
                    ++it2;
 
346
                    ++it1;
 
347
                }
 
348
                else
 
349
                {
 
350
                    equation_t& n = *it1;
 
351
                    n.lower = std::max(n.lower, it2->lower);
 
352
                    n.upper = std::min(n.upper, it2->upper);
 
353
                    /*if(n.upper < n.lower)
 
354
                    {
 
355
                        _result = result_t::IMPOSSIBLE;
 
356
                        _equations.clear();
 
357
                        return;
 
358
                    }*/
 
359
                    ++it1;
 
360
                    ++it2;
 
361
                }
 
362
            }
 
363
 
 
364
            if(it2 != other._equations.end())
 
365
                _equations.insert(_equations.end(), it2, other._equations.end());
 
366
        }
 
367
    }
 
368
}
 
369