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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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