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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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