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

Contents of /trunk/ascxx/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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