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

Contents of /trunk/ascxx/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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