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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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