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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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