/[ascend]/trunk/pygtk/simulation.cpp
ViewVC logotype

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 953 - (show annotations) (download) (as text)
Thu Dec 7 14:47:15 2006 UTC (17 years, 10 months ago) by johnpye
File MIME type: text/x-c++src
File size: 23030 byte(s)
Added test for C99 FPE handling
Fixing mess-up of ChildByChar in arrayinst.h header.
Added 'safeeval' config option to IDA.
Changed 'SigHandler' to 'SigHandlerFn *' in line with other function pointer datatypes being used in ASCEND.
Moved processVarStatus *after* 'Failed integrator' exception (ongoing issue).
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 */
19 #include <iostream>
20 #include <iomanip>
21 #include <stdexcept>
22 #include <sstream>
23 using namespace std;
24
25 #include "config.h"
26
27 extern "C"{
28 #include <utilities/ascConfig.h>
29 #include <utilities/error.h>
30 #include <utilities/ascSignal.h>
31 #include <utilities/ascMalloc.h>
32 #include <general/dstring.h>
33 #include <general/tm_time.h>
34 #include <compiler/instance_enum.h>
35 #include <compiler/fractions.h>
36 #include <compiler/compiler.h>
37 #include <compiler/dimen.h>
38 #include <compiler/symtab.h>
39 #include <compiler/instance_io.h>
40 #include <compiler/instantiate.h>
41 #include <compiler/bintoken.h>
42 #include <compiler/instance_enum.h>
43 #include <compiler/instquery.h>
44 #include <compiler/check.h>
45 #include <compiler/name.h>
46 #include <compiler/pending.h>
47 #include <compiler/importhandler.h>
48
49 #include <utilities/readln.h>
50 #include <solver/mtx.h>
51 #include <solver/slv_types.h>
52 #include <solver/var.h>
53 #include <solver/rel.h>
54 #include <solver/discrete.h>
55 #include <solver/conditional.h>
56 #include <solver/logrel.h>
57 #include <solver/bnd.h>
58 #include <solver/calc.h>
59 #include <solver/relman.h>
60 #include <solver/slv_common.h>
61 #include <solver/linsol.h>
62 #include <solver/linsolqr.h>
63 #include <solver/slv_client.h>
64 #include <solver/system.h>
65 #include <solver/slv_interface.h>
66 #include <solver/slvDOF.h>
67 #include <solver/slv3.h>
68 #include <solver/slv_stdcalls.h>
69 #include <solver/slv_server.h>
70 }
71
72 #include "simulation.h"
73 #include "solver.h"
74 #include "solverparameters.h"
75 #include "name.h"
76 #include "incidencematrix.h"
77 #include "variable.h"
78 #include "solverstatus.h"
79 #include "solverreporter.h"
80
81 /**
82 Create an instance of a type (call compiler etc)
83
84 @TODO fix mutex on compile command filenames
85 */
86 Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
87 sys = NULL;
88 //is_built = false;
89 // Create an Instance object for the 'simulation root' (we'll call
90 // it the 'simulation model') and it can be fetched using 'getModel()'
91 // any time later.
92 //simroot = Instanc(GetSimulationRoot(i),name);
93 }
94
95 Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
96 //is_built = old.is_built;
97 sys = old.sys;
98 bin_srcname = old.bin_srcname;
99 bin_objname = old.bin_objname;
100 bin_libname = old.bin_libname;
101 bin_cmd = old.bin_cmd;
102 bin_rm = old.bin_rm;
103 sing = NULL;
104 }
105
106 Simulation::~Simulation(){
107 //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
108 }
109
110 Instanc &
111 Simulation::getModel(){
112 if(!simroot.getInternalType()){
113 throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
114 }
115 if(InstanceKind(simroot.getInternalType())!=MODEL_INST){
116 throw runtime_error("Simulation::getModel: simroot is not a MODEL instance");
117 }
118 return simroot;
119 }
120
121
122 slv_system_structure *
123 Simulation::getSystem(){
124 if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
125 return sys;
126 }
127
128
129 const string
130 Simulation::getInstanceName(const Instanc &i) const{
131 char *n;
132 n = WriteInstanceNameString(i.getInternalType(),simroot.getInternalType());
133 string s(n);
134 ascfree(n);
135 return s;
136 }
137
138 const int
139 Simulation::getNumVars(){
140 return slv_get_num_solvers_vars(getSystem());
141 }
142
143
144 void
145 Simulation::write(){
146 simroot.write();
147 }
148
149 //------------------------------------------------------------------------------
150 // RUNNING MODEL 'METHODS'
151
152 void
153 Simulation::run(const Method &method){
154 Instanc &model = getModel();
155 this->run(method,model);
156 }
157
158 void
159 Simulation::runDefaultMethod(){
160 Method m = getType().getMethod(SymChar("on_load"));
161 run(m);
162 }
163
164 void
165 Simulation::run(const Method &method, Instanc &model){
166
167 // set the 'sim' pointer to our local variable...
168 CONSOLE_DEBUG("Setting shared pointer 'sim'");
169 importhandler_setsharedpointer("sim",this);
170
171 /*if(not is_built){
172 CONSOLE_DEBUG("WARNING, SIMULATION NOT YET BUILT");
173 }*/
174
175 CONSOLE_DEBUG("Running method %s...", method.getName());
176
177 Nam name = Nam(method.getSym());
178 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
179
180 error_reporter_tree_start();
181
182 Proc_enum pe;
183 pe = Initialize(
184 &*(model.getInternalType()) ,name.getInternalType(), "__not_named__"
185 ,ASCERR
186 ,0, NULL, NULL
187 );
188
189 int haserror=0;
190 if(error_reporter_tree_has_error()){
191 haserror=1;
192 }
193 error_reporter_tree_end();
194
195 // clear out the 'sim' pointer (soon it will be invalid)
196 importhandler_setsharedpointer("sim",NULL);
197 CONSOLE_DEBUG("Cleared shared pointer 'sim'");
198
199 if(pe == Proc_all_ok){
200 if(haserror){
201 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Method '%s' had error(s).",method.getName());
202 stringstream ss;
203 ss << "Method '"<<method.getName()<<"' returned 'all_ok' status but output error(s)";
204 throw runtime_error(ss.str());
205 }else{
206 ERROR_REPORTER_NOLINE(ASC_USER_SUCCESS,"Method '%s' returned 'all_ok' and output no errors.\n",method.getName());
207 }
208 //cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
209 }else{
210 stringstream ss;
211 ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
212 switch(pe){
213 case Proc_CallOK: ss << "Call OK"; break;
214 case Proc_CallError: ss << "Error occurred in call"; break;
215 case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
216 case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
217 case Proc_CallContinue: ss << "Skip to next iteration"; break;
218
219 case Proc_break: ss << "Break"; break;
220 case Proc_continue: ss << "Continue"; break;
221 case Proc_fallthru: ss << "Fall-through"; break;
222 case Proc_return: ss << "Return"; break;
223 case Proc_stop: ss << "Stop"; break;
224 case Proc_stack_exceeded: ss << "Stack exceeded"; break;
225 case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
226 case Proc_case_matched: ss << "Case matched"; break;
227 case Proc_case_unmatched: ss << "Case unmatched"; break;
228
229 case Proc_case_undefined_value: ss << "Undefined value in case"; break;
230 case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
231 case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
232 case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
233 case Proc_case_wrong_index: ss << "Wrong index in case"; break;
234 case Proc_case_wrong_value: ss << "Wrong value in case"; break;
235 case Proc_case_extra_values: ss << "Extra values in case"; break;
236 case Proc_bad_statement: ss << "Bad statement"; break;
237 case Proc_bad_name: ss << "Bad name"; break;
238 case Proc_for_duplicate_index: ss << "Duplicate index"; break;
239 case Proc_for_set_err: ss << "For set error"; break;
240 case Proc_for_not_set: ss << "For not set"; break;
241 case Proc_illegal_name_use: ss << "Illegal name use"; break;
242 case Proc_name_not_found: ss << "Name not found"; break;
243 case Proc_instance_not_found: ss << "Instance not found"; break;
244 case Proc_type_not_found: ss << "Type not found"; break;
245 case Proc_illegal_type_use: ss << "Illegal use"; break;
246 case Proc_proc_not_found: ss << "Method not found"; break;
247 case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
248 case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
249 case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
250 case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
251 case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
252 case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
253 case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
254 case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
255 case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
256 case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
257 case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
258 case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
259 case Proc_if_not_logical: ss << "If expression is not logical"; break;
260 case Proc_user_interrupt: ss << "User interrupt"; break;
261 case Proc_infinite_loop: ss << "Infinite loop"; break;
262 case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
263 case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
264 case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
265 case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
266 case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
267 case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
268 case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
269 case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
270 case Proc_lhs_error: ss << "Left-hand-side error"; break;
271 case Proc_rhs_error: ss << "Right-hand-side error"; break;
272 case Proc_unknown_error: ss << "Unknown error"; break;
273 default:
274 ss << "Invalid error code";
275 }
276
277 ss << " (" << int(pe) << ")";
278 throw runtime_error(ss.str());
279 }
280 }
281
282 //-----------------------------------------------------------------------------
283 // CHECKING METHODS
284
285 /**
286 Check that all the analysis went OK: solver lists are all there, etc...?
287
288 Can't return anything here because of limitations in the C API
289
290 @TODO there's something wrong with this at the moment: even after 'FIX'
291 methods are run, check shows them as not fixed, up until the point that 'SOLVE'
292 successfully completes. Something's not being synchronised properly...
293 */
294 void
295 Simulation::checkInstance(){
296 //cerr << "CHECKING SIMULATION INSTANCE" << endl;
297 /*if(!is_built){
298 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Simulation has not been built");
299 return;
300 }*/
301 Instance *i1 = getModel().getInternalType();
302 CheckInstance(stderr, &*i1);
303 //cerr << "DONE CHECKING INSTANCE" << endl;
304 }
305
306 /**
307 @return 1 = underspecified, 2 = square, 3 = structurally singular, 4 = overspecified
308 */
309 enum StructuralStatus
310 Simulation::checkDoF() const{
311 int dof, status;
312
313 if(!sys){
314 throw runtime_error("System is not built");
315 }
316
317 /*if(!is_built){
318 throw runtime_error("System not yet built");
319 }*/
320 CONSOLE_DEBUG("Calling slvDOF_status...");
321 slvDOF_status(sys, &status, &dof);
322 switch(status){
323 case ASCXX_DOF_UNDERSPECIFIED:
324 case ASCXX_DOF_SQUARE:
325 case ASCXX_DOF_OVERSPECIFIED:
326 case ASCXX_DOF_STRUCT_SINGULAR:
327 return (enum StructuralStatus)status;
328 case 5:
329 throw runtime_error("Unable to resolve degrees of freedom"); break;
330 default:
331 throw runtime_error("Invalid return status from slvDOF_status");
332 }
333 }
334
335 /**
336 Check consistency
337
338 @TODO what is the difference between this and checkStructuralSingularity?
339
340 @return list of freeable variables. List will be empty if sys is consistent.
341 */
342 vector<Variable>
343 Simulation::getFreeableVariables(){
344 vector<Variable> v;
345
346 //cerr << "CHECKING CONSISTENCY..." << endl;
347 int *fixedarrayptr=NULL;
348
349 if(!sys){
350 throw runtime_error("System not yet built");
351 }
352
353 int res = consistency_analysis(sys, &fixedarrayptr);
354
355 if(res==1){
356 cerr << "STRUCTURALLY CONSISTENT" << endl;
357 }else{
358 if(fixedarrayptr ==NULL){
359 ERROR_REPORTER_HERE(ASC_USER_ERROR,"STRUCTURALLY INCONSISTENT");
360 throw runtime_error("Invalid consistency analysis result returned!");
361 }
362
363 struct var_variable **vp = slv_get_master_var_list(sys);
364 for(int i=0; fixedarrayptr[i]!=-1; ++i){
365 v.push_back( Variable(this, vp[fixedarrayptr[i]]) );
366 }
367 }
368 return v;
369 }
370
371 /** Returns TRUE if all is OK (not singular) */
372 bool
373 Simulation::checkStructuralSingularity(){
374 int *vil;
375 int *ril;
376 int *fil;
377
378 cerr << "RETRIEVING slfDOF_structsing INFO" << endl;
379
380 int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
381 struct var_variable **varlist = slv_get_solvers_var_list(sys);
382 struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
383
384 if(this->sing){
385 cerr << "DELETING OLD SINGULATING INFO" << endl;
386 delete this->sing;
387 this->sing = NULL;
388 }
389
390 if(res==1){
391 CONSOLE_DEBUG("processing singularity data...");
392 sing = new SingularityInfo();
393
394 // pull in the lists of vars and rels, and the freeable vars:
395 for(int i=0; ril[i]!=-1; ++i){
396 sing->rels.push_back( Relation(this, rellist[ril[i]]) );
397 }
398
399 for(int i=0; vil[i]!=-1; ++i){
400 sing->vars.push_back( Variable(this, varlist[vil[i]]) );
401 }
402
403 for(int i=0; fil[i]!=-1; ++i){
404 sing->freeablevars.push_back( Variable(this, varlist[fil[i]]) );
405 }
406
407 // we're done with those lists now
408 ASC_FREE(vil);
409 ASC_FREE(ril);
410 ASC_FREE(fil);
411
412 if(sing->isSingular()){
413 CONSOLE_DEBUG("singularity found");
414 this->sing = sing;
415 return FALSE;
416 }
417 CONSOLE_DEBUG("no singularity");
418 delete sing;
419 return TRUE;
420 }else{
421 if(res==0){
422 throw runtime_error("Unable to determine singularity lists");
423 }else{
424 throw runtime_error("Invalid return from slvDOF_structsing.");
425 }
426 }
427 }
428
429 /**
430 If the checkStructuralSingularity analysis has been done,
431 this funciton will let you access the SingularityInfo data that was
432 stored.
433 */
434 const SingularityInfo &
435 Simulation::getSingularityInfo() const{
436 if(sing==NULL){
437 throw runtime_error("No singularity info present");
438 }
439 return *sing;
440 }
441
442 //------------------------------------------
443 // ASSIGNING SOLVER TO SIMULATION
444
445 void
446 Simulation::setSolver(Solver &solver){
447 /* CONSOLE_DEBUG("Setting solver on sim %p, root inst %p",this,this->simroot.getInternalType()); */
448
449 try{
450 build();
451 }catch(runtime_error &e){
452 stringstream ss;
453 ss << "Couldn't prepare system for solving:";
454 ss << e.what();
455 throw runtime_error(ss.str());
456 }
457
458 CONSOLE_DEBUG("Calling slv_select_solver...");
459 int selected = slv_select_solver(sys, solver.getIndex());
460 //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
461
462 if(selected<0){
463 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
464 throw runtime_error("Failed to select solver");
465 }
466
467 if(selected!=solver.getIndex()){
468 solver = Solver(slv_solver_name(selected));
469 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
470 }
471
472 if( slv_eligible_solver(sys) <= 0){
473 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
474 throw runtime_error("Inelegible solver");
475 }
476 }
477
478 const Solver
479 Simulation::getSolver() const{
480 int index = slv_get_selected_solver(sys);
481 //cerr << "Simulation::getSolver: index = " << index << endl;
482 if(index<0)throw runtime_error("No solver selected");
483
484 return Solver(slv_solver_name(index));
485 }
486
487 //------------------------------------------------------------------------------
488 // BUILD THE SYSTEM
489
490 /**
491 Build the system (send it to the solver)
492 */
493 void
494 Simulation::build(){
495 if(sys){
496 CONSOLE_DEBUG("System is already built");
497 return;
498 }
499
500 if(simroot.getKind() != MODEL_INST){
501 throw runtime_error("Simulation does not contain a MODEL_INST");
502 }
503
504 if(NumberPendingInstances(simroot.getInternalType())){
505 throw runtime_error("System has pending instances; can't yet send to solver.");
506 }
507
508 sys = system_build(simroot.getInternalType());
509 if(!sys){
510 throw runtime_error("Unable to build system");
511 }
512
513 CONSOLE_DEBUG("System built OK");
514 }
515
516
517 //------------------------------------------------------------------------------
518 // SOLVER CONFIGURATION PARAMETERS
519
520 /**
521 Get solver parameters struct wrapped up as a SolverParameters class.
522 */
523 SolverParameters
524 Simulation::getSolverParameters() const{
525 //if(!is_built)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
526 if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has no 'sys' assigned.");
527
528 slv_parameters_t p;
529 slv_get_parameters(sys,&p);
530 return SolverParameters(p);
531 }
532
533 /**
534 Update the solver parameters by passing a new set back
535 */
536 void
537 Simulation::setSolverParameters(SolverParameters &P){
538 if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
539 CONSOLE_DEBUG("Calling slv_set_parameters");
540 slv_set_parameters(sys, &(P.getInternalType()));
541 }
542
543 //------------------------------------------------------------------------------
544 // PRE-SOLVE DIAGNOSTICS
545
546 /**
547 Get a list of variables to fix to make an underspecified system
548 become square. Also seems to return stuff when you have a structurally
549 singuler system.
550 */
551 vector<Variable>
552 Simulation::getFixableVariables(){
553 //cerr << "GETTING FIXABLE VARIABLES..." << endl;
554 vector<Variable> vars;
555
556 if(!sys){
557 throw runtime_error("Simulation system not yet built");
558 }
559
560 int32 *vip; /** TODO ensure 32 bit integers are used */
561
562 // Get IDs of elegible variables in array at vip...
563 CONSOLE_DEBUG("Calling slvDOF_eligible");
564 if(!slvDOF_eligible(sys,&vip)){
565 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
566 }else{
567 struct var_variable **vp = slv_get_solvers_var_list(sys);
568
569 if(vp==NULL){
570 throw runtime_error("Simulation variable list is null");
571 }
572
573 // iterate through this list until we find a -1:
574 int i=0;
575 int var_index = vip[i];
576 while(var_index >= 0){
577 struct var_variable *var = vp[var_index];
578 vars.push_back( Variable(this, var) );
579 ++i;
580 var_index = vip[i];
581 }
582 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
583 ascfree(vip);
584 }
585
586 return vars;
587 }
588
589 /**
590 Get the list of variables near their bounds. Helps to indentify why
591 you might be having non-convergence problems.
592 */
593 vector<Variable>
594 Simulation::getVariablesNearBounds(const double &epsilon){
595 //cerr << "GETTING VARIABLES NEAR BOUNDS..." << endl;
596 vector<Variable> vars;
597
598 if(!sys){
599 throw runtime_error("Simulation system not yet built");
600 }
601
602 int *vip;
603 CONSOLE_DEBUG("Calling slv_near_bounds...");
604 if(slv_near_bounds(sys,epsilon,&vip)){
605 struct var_variable **vp = slv_get_solvers_var_list(sys);
606 struct var_variable *var;
607 cerr << "VARS FOUND NEAR BOUNDS" << endl;
608 int nlow = vip[0];
609 int nhigh = vip[1];
610 int lim1 = 2 + nlow;
611 for(int i=2; i<lim1; ++i){
612 var = vp[vip[i]];
613 char *var_name = var_make_name(sys,var);
614 cerr << "AT LOWER BOUND: " << var_name << endl;
615 ascfree(var_name);
616 vars.push_back(Variable(this,var));
617 };
618 int lim2 = lim1 + nhigh;
619 for(int i=lim1; i<lim2; ++i){
620 var = vp[vip[i]];
621 char *var_name = var_make_name(sys,var);
622 cerr << "AT UPPER BOUND: " << var_name << endl;
623 ascfree(var_name);
624 vars.push_back(Variable(this,var));
625 }
626 }
627 ascfree(vip);
628 return vars;
629 }
630
631
632 bool
633 SingularityInfo::isSingular() const{
634 if(vars.size()||rels.size()){
635 return true;
636 }
637 return false;
638 }
639
640 //------------------------------------------------------------------------------
641 // SOLVING
642
643 /**
644 Solve the system through to convergence. This function is hardwired with
645 a maximum of 1000 iterations, but will interrupt itself when the 'stop'
646 condition comes back from the SolverReporter.
647 */
648 void
649 Simulation::solve(Solver solver, SolverReporter &reporter){
650
651 if(!sys){
652 try{
653 build();
654 }catch(runtime_error &e){
655 stringstream ss;
656 ss << "Unable to build system: ";
657 ss << e.what();
658 throw runtime_error(ss.str());
659 }
660 }
661
662 setSolver(solver);
663
664 //cerr << "PRESOLVING SYSTEM...";
665 CONSOLE_DEBUG("Calling slv_presolve...");
666 slv_presolve(sys);
667 //cerr << "DONE" << endl;
668
669 //cerr << "SOLVING SYSTEM..." << endl;
670 // Add some stuff here for cleverer iteration....
671 unsigned niter = 1000;
672 //double updateinterval = 0.02;
673
674 double starttime = tm_cpu_time();
675 //double lastupdate = starttime;
676 SolverStatus status;
677 //int solved_vars=0;
678 bool stop=false;
679
680 status.getSimulationStatus(*this);
681 reporter.report(&status);
682
683 for(unsigned iter = 1; iter <= niter && !stop; ++iter){
684
685 if(status.isReadyToSolve()){
686 /* CONSOLE_DEBUG("Calling slv_iterate..."); */
687 slv_iterate(sys);
688 }
689
690 status.getSimulationStatus(*this);
691
692 if(reporter.report(&status)){
693 stop = true;
694 }
695 }
696
697 double elapsed = tm_cpu_time() - starttime;
698 CONSOLE_DEBUG("Elapsed time: %0.3f", elapsed);
699
700 activeblock = status.getCurrentBlockNum();
701
702 reporter.finalise(&status);
703
704 // Just a little bit of console output:
705
706 if(status.isOK()){
707 //cerr << "... SOLVED, STATUS OK" << endl;
708 }else{
709 cerr << "... SOLVER FAILED" << endl;
710 }
711
712 //cerr << "SOLVER PERFORMED " << status.getIterationNum() << " ITERATIONS IN " << elapsed << "s" << endl;
713
714 // communicate solver variable status back to the instance tree
715 processVarStatus();
716 }
717
718 //------------------------------------------------------------------------------
719 // POST-SOLVE DIAGNOSTICS
720
721 const int
722 Simulation::getActiveBlock() const{
723 return activeblock;
724 }
725
726 /**
727 Return an IncidenceMatrix built from the current state of the solver system.
728
729 This will actually return something meaningful even before solve.
730 */
731 IncidenceMatrix
732 Simulation::getIncidenceMatrix(){
733 return IncidenceMatrix(*this);
734 }
735
736 /**
737 This function looks at all the variables in the solve's list and updates
738 the variable status for the corresponding instances.
739
740 It does this by using the 'interface pointer' in the Instance, see
741 the C-API function GetInterfacePtr.
742
743 This is used to display visually which variables have been solved, which
744 ones have not yet been attempted, and which ones were active when the solver
745 failed (ASCXX_VAR_ACTIVE).
746 */
747 void
748 Simulation::processVarStatus(){
749 if(!sys)throw runtime_error("No system built");
750
751 // this is a cheap function call:
752 const mtx_block_t *bb = slv_get_solvers_blocks(getSystem());
753
754 var_variable **vlist = slv_get_solvers_var_list(getSystem());
755 int nvars = slv_get_num_solvers_vars(getSystem());
756
757 slv_status_t status;
758 slv_get_status(getSystem(), &status);
759
760 if(status.block.number_of == 0){
761 cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
762 return;
763 }
764
765 int activeblock = status.block.current_block;
766 int low = bb->block[activeblock].col.low;
767 int high = bb->block[activeblock].col.high;
768 bool allsolved = status.converged;
769 for(int c=0; c < nvars; ++c){
770 var_variable *v = vlist[c];
771 Instanc i((Instance *)var_instance(v));
772 VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
773 if(i.isFixed()){
774 s = ASCXX_VAR_FIXED;
775 }else if(var_incident(v) && var_active(v)){
776 if(allsolved || c < low){
777 s = ASCXX_VAR_SOLVED;
778 }else if(c <= high){
779 s = ASCXX_VAR_ACTIVE;
780 }else{
781 s = ASCXX_VAR_UNSOLVED;
782 }
783 }
784 i.setVarStatus(s);
785 }
786 }
787

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22