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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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