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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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