/[ascend]/branches/relerrorlist/ascxx/simulation.cpp
ViewVC logotype

Contents of /branches/relerrorlist/ascxx/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3197 - (show annotations) (download) (as text)
Sun Apr 23 03:30:25 2017 UTC (15 months, 3 weeks ago) by jpye
File MIME type: text/x-c++src
File size: 29750 byte(s)
test suite almost done, one memory leak still to be fixed.

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

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