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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1133 - (show annotations) (download) (as text)
Sun Jan 14 11:51:48 2007 UTC (13 years, 4 months ago) by johnpye
File MIME type: text/x-c++src
File size: 25453 byte(s)
Removed the 'REX' and 'IEX' array aliases (added the solver parameter comments directly in the slv_param_* calls).
Renamed Simulation::[gs]etSolverParameters to Simulation::[gs]etParameters.
Added [gs]etParameter (singular) methods to SWIG wrapper.
Fixed missing NULL in IDA.
Fixed bug with BDF/AM wrong way round in LSODE.
Removed some debug output from slv.c
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/ascConfig.h>
29 #include <utilities/error.h>
30 #include <utilities/ascSignal.h>
31 #include <utilities/ascMalloc.h>
32 #include <general/dstring.h>
33 #include <general/tm_time.h>
34 #include <compiler/instance_enum.h>
35 #include <compiler/fractions.h>
36 #include <compiler/compiler.h>
37 #include <compiler/dimen.h>
38 #include <compiler/symtab.h>
39 #include <compiler/instance_io.h>
40 #include <compiler/instantiate.h>
41 #include <compiler/bintoken.h>
42 #include <compiler/instance_enum.h>
43 #include <compiler/instquery.h>
44 #include <compiler/check.h>
45 #include <compiler/name.h>
46 #include <compiler/pending.h>
47 #include <compiler/importhandler.h>
48
49 #include <utilities/readln.h>
50 #include <solver/mtx.h>
51 #include <solver/slv_types.h>
52 #include <solver/var.h>
53 #include <solver/rel.h>
54 #include <solver/discrete.h>
55 #include <solver/conditional.h>
56 #include <solver/logrel.h>
57 #include <solver/bnd.h>
58 #include <solver/calc.h>
59 #include <solver/relman.h>
60 #include <solver/slv_common.h>
61 #include <solver/linsol.h>
62 #include <solver/linsolqr.h>
63 #include <solver/slv_client.h>
64 #include <solver/system.h>
65 #include <solver/slv_interface.h>
66 #include <solver/slvDOF.h>
67 #include <solver/slv3.h>
68 #include <solver/slv_stdcalls.h>
69 #include <solver/slv_server.h>
70 }
71
72 #include "simulation.h"
73 #include "solver.h"
74 #include "solverparameters.h"
75 #include "name.h"
76 #include "incidencematrix.h"
77 #include "variable.h"
78 #include "solverstatus.h"
79 #include "solverreporter.h"
80 #include "matrix.h"
81
82 Simulation::Simulation(){
83 throw runtime_error("Can't create Simulation without arguments");
84 }
85
86 /**
87 Create an instance of a type (call compiler etc)
88
89 @TODO fix mutex on compile command filenames
90 */
91 Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
92 CONSOLE_DEBUG("Created simulation");
93 sys = NULL;
94 //is_built = false;
95 // Create an Instance object for the 'simulation root' (we'll call
96 // it the 'simulation model') and it can be fetched using 'getModel()'
97 // any time later.
98 //simroot = Instanc(GetSimulationRoot(i),name);
99 }
100
101 Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
102 //is_built = old.is_built;
103 CONSOLE_DEBUG("Copying Simulation...");
104 sys = old.sys;
105 sing = NULL;
106 }
107
108 Simulation::~Simulation(){
109 CONSOLE_DEBUG("Destroying Simulation...");
110 /*
111 // FIXME removing this here, because Python overzealously seems to delete simulations
112
113 CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
114 system_free_reused_mem();
115 if(sys){
116 CONSOLE_DEBUG("Destroying simulation system...");
117 system_destroy(sys);
118 }
119 */
120 sys = NULL;
121 }
122
123 Instanc &
124 Simulation::getModel(){
125 if(!simroot.getInternalType()){
126 throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
127 }
128 if(InstanceKind(simroot.getInternalType())!=MODEL_INST){
129 throw runtime_error("Simulation::getModel: simroot is not a MODEL instance");
130 }
131 return simroot;
132 }
133
134
135 slv_system_structure *
136 Simulation::getSystem(){
137 if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
138 return sys;
139 }
140
141
142 const string
143 Simulation::getInstanceName(const Instanc &i) const{
144 char *n;
145 n = WriteInstanceNameString(i.getInternalType(),simroot.getInternalType());
146 string s(n);
147 ascfree(n);
148 return s;
149 }
150
151 const int
152 Simulation::getNumVars(){
153 return slv_get_num_solvers_vars(getSystem());
154 }
155
156
157 void
158 Simulation::write(){
159 simroot.write();
160 }
161
162 //------------------------------------------------------------------------------
163 // RUNNING MODEL 'METHODS'
164
165 void
166 Simulation::run(const Method &method){
167 Instanc &model = getModel();
168 this->run(method,model);
169 }
170
171 void
172 Simulation::runDefaultMethod(){
173 const Type &type = getType();
174 Method m;
175 try{
176 m = type.getMethod(SymChar("on_load"));
177 }catch(runtime_error &e){
178 ERROR_REPORTER_NOLINE(ASC_USER_WARNING,"There is no 'on_load' method defined for type '%s'",type.getName().toString());
179 return;
180 }
181 run(m);
182 }
183
184 void
185 Simulation::run(const Method &method, Instanc &model){
186
187 // set the 'sim' pointer to our local variable...
188 CONSOLE_DEBUG("Setting shared pointer 'sim' = %p",this);
189 importhandler_setsharedpointer("sim",this);
190
191 /*if(not is_built){
192 CONSOLE_DEBUG("WARNING, SIMULATION NOT YET BUILT");
193 }*/
194
195 CONSOLE_DEBUG("Running method %s...", method.getName());
196
197 Nam name = Nam(method.getSym());
198 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
199
200 error_reporter_tree_start();
201
202 CONSOLE_DEBUG("sys = %p",sys);
203 CONSOLE_DEBUG("simroot = %p",simroot.getInternalType());
204
205 Proc_enum pe;
206 pe = Initialize(
207 &*(model.getInternalType()) ,name.getInternalType(), "__not_named__"
208 ,ASCERR
209 ,0, NULL, NULL
210 );
211
212 int haserror=0;
213 if(error_reporter_tree_has_error()){
214 haserror=1;
215 }
216 error_reporter_tree_end();
217
218 // clear out the 'sim' pointer (soon it will be invalid)
219 importhandler_setsharedpointer("sim",NULL);
220 CONSOLE_DEBUG("Cleared shared pointer 'sim'");
221
222 if(pe == Proc_all_ok){
223 if(haserror){
224 ERROR_REPORTER_NOLINE(ASC_PROG_ERR,"Method '%s' had error(s).",method.getName());
225 stringstream ss;
226 ss << "Method '"<<method.getName()<<"' returned 'all_ok' status but output error(s)";
227 throw runtime_error(ss.str());
228 }else{
229 ERROR_REPORTER_NOLINE(ASC_USER_SUCCESS,"Method '%s' returned 'all_ok' and output no errors.\n",method.getName());
230 }
231 //cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
232 }else{
233 stringstream ss;
234 ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
235 switch(pe){
236 case Proc_CallOK: ss << "Call OK"; break;
237 case Proc_CallError: ss << "Error occurred in call"; break;
238 case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
239 case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
240 case Proc_CallContinue: ss << "Skip to next iteration"; break;
241
242 case Proc_break: ss << "Break"; break;
243 case Proc_continue: ss << "Continue"; break;
244 case Proc_fallthru: ss << "Fall-through"; break;
245 case Proc_return: ss << "Return"; break;
246 case Proc_stop: ss << "Stop"; break;
247 case Proc_stack_exceeded: ss << "Stack exceeded"; break;
248 case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
249 case Proc_case_matched: ss << "Case matched"; break;
250 case Proc_case_unmatched: ss << "Case unmatched"; break;
251
252 case Proc_case_undefined_value: ss << "Undefined value in case"; break;
253 case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
254 case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
255 case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
256 case Proc_case_wrong_index: ss << "Wrong index in case"; break;
257 case Proc_case_wrong_value: ss << "Wrong value in case"; break;
258 case Proc_case_extra_values: ss << "Extra values in case"; break;
259 case Proc_bad_statement: ss << "Bad statement"; break;
260 case Proc_bad_name: ss << "Bad name"; break;
261 case Proc_for_duplicate_index: ss << "Duplicate index"; break;
262 case Proc_for_set_err: ss << "For set error"; break;
263 case Proc_for_not_set: ss << "For not set"; break;
264 case Proc_illegal_name_use: ss << "Illegal name use"; break;
265 case Proc_name_not_found: ss << "Name not found"; break;
266 case Proc_instance_not_found: ss << "Instance not found"; break;
267 case Proc_type_not_found: ss << "Type not found"; break;
268 case Proc_illegal_type_use: ss << "Illegal use"; break;
269 case Proc_proc_not_found: ss << "Method not found"; break;
270 case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
271 case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
272 case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
273 case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
274 case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
275 case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
276 case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
277 case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
278 case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
279 case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
280 case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
281 case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
282 case Proc_if_not_logical: ss << "If expression is not logical"; break;
283 case Proc_user_interrupt: ss << "User interrupt"; break;
284 case Proc_infinite_loop: ss << "Infinite loop"; break;
285 case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
286 case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
287 case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
288 case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
289 case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
290 case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
291 case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
292 case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
293 case Proc_lhs_error: ss << "Left-hand-side error"; break;
294 case Proc_rhs_error: ss << "Right-hand-side error"; break;
295 case Proc_unknown_error: ss << "Unknown error"; break;
296 default:
297 ss << "Invalid error code";
298 }
299
300 ss << " (" << int(pe) << ")";
301 throw runtime_error(ss.str());
302 }
303 }
304
305 //-----------------------------------------------------------------------------
306 // CHECKING METHODS
307
308 /**
309 Check that all the analysis went OK: solver lists are all there, etc...?
310
311 Can't return anything here because of limitations in the C API
312
313 @TODO there's something wrong with this at the moment: even after 'FIX'
314 methods are run, check shows them as not fixed, up until the point that 'SOLVE'
315 successfully completes. Something's not being synchronised properly...
316 */
317 void
318 Simulation::checkInstance(){
319 //cerr << "CHECKING SIMULATION INSTANCE" << endl;
320 /*if(!is_built){
321 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Simulation has not been built");
322 return;
323 }*/
324 Instance *i1 = getModel().getInternalType();
325 CheckInstance(stderr, &*i1);
326 //cerr << "DONE CHECKING INSTANCE" << endl;
327 }
328
329 /**
330 @return 1 = underspecified, 2 = square, 3 = structurally singular, 4 = overspecified
331 */
332 enum StructuralStatus
333 Simulation::checkDoF() const{
334 int dof, status;
335
336 if(!sys){
337 throw runtime_error("System is not built");
338 }
339
340 /*if(!is_built){
341 throw runtime_error("System not yet built");
342 }*/
343 CONSOLE_DEBUG("Calling slvDOF_status...");
344 slvDOF_status(sys, &status, &dof);
345 switch(status){
346 case ASCXX_DOF_UNDERSPECIFIED:
347 case ASCXX_DOF_SQUARE:
348 case ASCXX_DOF_OVERSPECIFIED:
349 case ASCXX_DOF_STRUCT_SINGULAR:
350 return (enum StructuralStatus)status;
351 case 5:
352 throw runtime_error("Unable to resolve degrees of freedom"); break;
353 default:
354 throw runtime_error("Invalid return status from slvDOF_status");
355 }
356 }
357
358 /**
359 Check consistency
360
361 @TODO what is the difference between this and checkStructuralSingularity?
362
363 @return list of freeable variables. List will be empty if sys is consistent.
364 */
365 vector<Variable>
366 Simulation::getFreeableVariables(){
367 vector<Variable> v;
368
369 //cerr << "CHECKING CONSISTENCY..." << endl;
370 int *fixedarrayptr=NULL;
371
372 if(!sys){
373 throw runtime_error("System not yet built");
374 }
375
376 int res = consistency_analysis(sys, &fixedarrayptr);
377
378 if(res==1){
379 cerr << "STRUCTURALLY CONSISTENT" << endl;
380 }else{
381 if(fixedarrayptr ==NULL){
382 ERROR_REPORTER_HERE(ASC_USER_ERROR,"STRUCTURALLY INCONSISTENT");
383 throw runtime_error("Invalid consistency analysis result returned!");
384 }
385
386 struct var_variable **vp = slv_get_master_var_list(sys);
387 for(int i=0; fixedarrayptr[i]!=-1; ++i){
388 v.push_back( Variable(this, vp[fixedarrayptr[i]]) );
389 }
390 }
391 return v;
392 }
393
394 /** Returns TRUE if all is OK (not singular) */
395 bool
396 Simulation::checkStructuralSingularity(){
397 int *vil;
398 int *ril;
399 int *fil;
400
401 if(this->sing){
402 cerr << "DELETING OLD SINGULATING INFO" << endl;
403 delete this->sing;
404 this->sing = NULL;
405 }
406
407 cerr << "RETRIEVING slfDOF_structsing INFO" << endl;
408
409 int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
410
411
412 if(res==1){
413 throw runtime_error("Unable to determine singularity lists");
414 }
415
416 if(res!=0){
417 throw runtime_error("Invalid return from slvDOF_structsing.");
418 }
419
420
421 CONSOLE_DEBUG("processing singularity data...");
422 sing = new SingularityInfo();
423
424 struct var_variable **varlist = slv_get_solvers_var_list(sys);
425 struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
426
427 // pull in the lists of vars and rels, and the freeable vars:
428 for(int i=0; ril[i]!=-1; ++i){
429 sing->rels.push_back( Relation(this, rellist[ril[i]]) );
430 }
431
432 for(int i=0; vil[i]!=-1; ++i){
433 sing->vars.push_back( Variable(this, varlist[vil[i]]) );
434 }
435
436 for(int i=0; fil[i]!=-1; ++i){
437 sing->freeablevars.push_back( Variable(this, varlist[fil[i]]) );
438 }
439
440 // we're done with those lists now
441 ASC_FREE(vil);
442 ASC_FREE(ril);
443 ASC_FREE(fil);
444
445 if(sing->isSingular()){
446 CONSOLE_DEBUG("singularity found");
447 this->sing = sing;
448 return FALSE;
449 }
450 CONSOLE_DEBUG("no singularity");
451 delete sing;
452 return TRUE;
453 }
454
455 /**
456 If the checkStructuralSingularity analysis has been done,
457 this funciton will let you access the SingularityInfo data that was
458 stored.
459 */
460 const SingularityInfo &
461 Simulation::getSingularityInfo() const{
462 if(sing==NULL){
463 throw runtime_error("No singularity info present");
464 }
465 return *sing;
466 }
467
468 //------------------------------------------
469 // ASSIGNING SOLVER TO SIMULATION
470
471 void
472 Simulation::setSolver(Solver &solver){
473 /* CONSOLE_DEBUG("Setting solver on sim %p, root inst %p",this,this->simroot.getInternalType()); */
474
475 try{
476 // build the system (if not built already)
477 build();
478 }catch(runtime_error &e){
479 stringstream ss;
480 ss << "Couldn't prepare system for solving:";
481 ss << e.what();
482 throw runtime_error(ss.str());
483 }
484
485 CONSOLE_DEBUG("Selecting solver '%s'",solver.getName().c_str());
486 int selected = slv_select_solver(sys, solver.getIndex());
487 //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
488
489 if(selected<0){
490 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
491 throw runtime_error("Failed to select solver");
492 }
493
494 if(selected!=solver.getIndex()){
495 solver = Solver(slv_solver_name(selected));
496 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
497 }
498
499 if( slv_eligible_solver(sys) <= 0){
500 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
501 throw runtime_error("Inelegible solver");
502 }
503 }
504
505 const Solver
506 Simulation::getSolver() const{
507 int index = slv_get_selected_solver(sys);
508 //cerr << "Simulation::getSolver: index = " << index << endl;
509 if(index<0)throw runtime_error("No solver selected");
510
511 return Solver(slv_solver_name(index));
512 }
513
514 //------------------------------------------------------------------------------
515 // BUILD THE SYSTEM
516
517 /**
518 Build the system (send it to the solver)
519 */
520 void
521 Simulation::build(){
522 if(sys){
523 CONSOLE_DEBUG("System is already built (%p)",sys);
524 return;
525 }
526
527 if(simroot.getKind() != MODEL_INST){
528 throw runtime_error("Simulation does not contain a MODEL_INST");
529 }
530
531 if(NumberPendingInstances(simroot.getInternalType())){
532 throw runtime_error("System has pending instances; can't yet send to solver.");
533 }
534
535 CONSOLE_DEBUG("============== REALLY building system...");
536 sys = system_build(simroot.getInternalType());
537 if(!sys){
538 throw runtime_error("Unable to build system");
539 }
540
541 CONSOLE_DEBUG("System built OK");
542 }
543
544
545 //------------------------------------------------------------------------------
546 // SOLVER CONFIGURATION PARAMETERS
547
548 /**
549 Get solver parameters struct wrapped up as a SolverParameters class.
550 */
551 SolverParameters
552 Simulation::getParameters() const{
553 //if(!is_built)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
554 if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has no 'sys' assigned.");
555
556 slv_parameters_t p;
557 slv_get_parameters(sys,&p);
558 return SolverParameters(p);
559 }
560
561 /**
562 Update the solver parameters by passing a new set back
563 */
564 void
565 Simulation::setParameters(SolverParameters &P){
566 if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
567 CONSOLE_DEBUG("Calling slv_set_parameters");
568 slv_set_parameters(sys, &(P.getInternalType()));
569 }
570
571 //------------------------------------------------------------------------------
572 // PRE-SOLVE DIAGNOSTICS
573
574 /**
575 Get a list of variables to fix to make an underspecified system
576 become square. Also seems to return stuff when you have a structurally
577 singuler system.
578 */
579 vector<Variable>
580 Simulation::getFixableVariables(){
581 //cerr << "GETTING FIXABLE VARIABLES..." << endl;
582 vector<Variable> vars;
583
584 if(!sys){
585 throw runtime_error("Simulation system not yet built");
586 }
587
588 int32 *vip; /** TODO ensure 32 bit integers are used */
589
590 // Get IDs of elegible variables in array at vip...
591 CONSOLE_DEBUG("Calling slvDOF_eligible");
592 if(!slvDOF_eligible(sys,&vip)){
593 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
594 }else{
595 struct var_variable **vp = slv_get_solvers_var_list(sys);
596
597 if(vp==NULL){
598 throw runtime_error("Simulation variable list is null");
599 }
600
601 // iterate through this list until we find a -1:
602 int i=0;
603 int var_index = vip[i];
604 while(var_index >= 0){
605 struct var_variable *var = vp[var_index];
606 vars.push_back( Variable(this, var) );
607 ++i;
608 var_index = vip[i];
609 }
610 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
611 ascfree(vip);
612 }
613
614 return vars;
615 }
616
617 /**
618 Return a list of ALL the fixed variables in the solver's variable list
619 */
620 vector<Variable>
621 Simulation::getFixedVariables(){
622 if(!sys)throw runtime_error("Simulation system not build yet");
623 vector<Variable> vars;
624 var_variable **vlist = slv_get_solvers_var_list(sys);
625 unsigned long nvars = slv_get_num_solvers_vars(sys);
626 for(unsigned long i=0;i<nvars;++i){
627 if(!var_fixed(vlist[i]))continue;
628 vars.push_back(Variable(this,vlist[i]));
629 }
630 return vars;
631 }
632
633 /**
634 For solvers that store a big matrix for the system, return a pointer to that
635 matrix (struct mtx_header*) as a C++-wrapped object of class Matrix.
636 */
637 Matrix
638 Simulation::getMatrix(){
639 if(!sys)throw runtime_error("Simulation system not built yet");
640 mtx_matrix_t M = slv_get_sys_mtx(sys);
641 if(M==NULL)throw runtime_error("Simulation system does not possess a matrix");
642 return Matrix(M);
643 }
644
645 /**
646 Get the list of variables near their bounds. Helps to indentify why
647 you might be having non-convergence problems.
648 */
649 vector<Variable>
650 Simulation::getVariablesNearBounds(const double &epsilon){
651 //cerr << "GETTING VARIABLES NEAR BOUNDS..." << endl;
652 vector<Variable> vars;
653
654 if(!sys){
655 throw runtime_error("Simulation system not yet built");
656 }
657
658 int *vip;
659 CONSOLE_DEBUG("Calling slv_near_bounds...");
660 if(slv_near_bounds(sys,epsilon,&vip)){
661 struct var_variable **vp = slv_get_solvers_var_list(sys);
662 struct var_variable *var;
663 cerr << "VARS FOUND NEAR BOUNDS" << endl;
664 int nlow = vip[0];
665 int nhigh = vip[1];
666 int lim1 = 2 + nlow;
667 for(int i=2; i<lim1; ++i){
668 var = vp[vip[i]];
669 char *var_name = var_make_name(sys,var);
670 cerr << "AT LOWER BOUND: " << var_name << endl;
671 ascfree(var_name);
672 vars.push_back(Variable(this,var));
673 };
674 int lim2 = lim1 + nhigh;
675 for(int i=lim1; i<lim2; ++i){
676 var = vp[vip[i]];
677 char *var_name = var_make_name(sys,var);
678 cerr << "AT UPPER BOUND: " << var_name << endl;
679 ascfree(var_name);
680 vars.push_back(Variable(this,var));
681 }
682 }
683 ascfree(vip);
684 return vars;
685 }
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 slv_get_status(sys, &status);
826
827 if(status.block.number_of == 0){
828 cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
829 return;
830 }
831
832 int activeblock = status.block.current_block;
833 int low = bb->block[activeblock].col.low;
834 int high = bb->block[activeblock].col.high;
835 bool allsolved = status.converged;
836 for(int c=0; c < nvars; ++c){
837 var_variable *v = vlist[c];
838 Instanc i((Instance *)var_instance(v));
839 VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
840 if(i.isFixed()){
841 s = ASCXX_VAR_FIXED;
842 }else if(var_incident(v) && var_active(v)){
843 if(allsolved || c < low){
844 s = ASCXX_VAR_SOLVED;
845 }else if(c <= high){
846 s = ASCXX_VAR_ACTIVE;
847 }else{
848 s = ASCXX_VAR_UNSOLVED;
849 }
850 }
851 i.setVarStatus(s);
852 }
853
854 CONSOLE_DEBUG(" ...done var status");
855
856 }
857

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