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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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