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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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