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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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