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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 918 - (show annotations) (download) (as text)
Sat Nov 4 03:58:06 2006 UTC (15 years, 7 months ago) by johnpye
File MIME type: text/x-c++src
File size: 22714 byte(s)
Little bit of added debug output
1 /* ASCEND modelling environment
2 Copyright (C) 2006 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 */
19 #include <iostream>
20 #include <iomanip>
21 #include <stdexcept>
22 #include <sstream>
23 using namespace std;
24
25 #include "config.h"
26
27 extern "C"{
28 #include <utilities/ascConfig.h>
29 #include <utilities/error.h>
30 #include <utilities/ascSignal.h>
31 #include <utilities/ascMalloc.h>
32 #include <general/dstring.h>
33 #include <general/tm_time.h>
34 #include <compiler/instance_enum.h>
35 #include <compiler/fractions.h>
36 #include <compiler/compiler.h>
37 #include <compiler/dimen.h>
38 #include <compiler/symtab.h>
39 #include <compiler/instance_io.h>
40 #include <compiler/instantiate.h>
41 #include <compiler/bintoken.h>
42 #include <compiler/instance_enum.h>
43 #include <compiler/instquery.h>
44 #include <compiler/check.h>
45 #include <compiler/name.h>
46 #include <compiler/pending.h>
47 #include <compiler/importhandler.h>
48
49 #include <utilities/readln.h>
50 #include <solver/mtx.h>
51 #include <solver/slv_types.h>
52 #include <solver/var.h>
53 #include <solver/rel.h>
54 #include <solver/discrete.h>
55 #include <solver/conditional.h>
56 #include <solver/logrel.h>
57 #include <solver/bnd.h>
58 #include <solver/calc.h>
59 #include <solver/relman.h>
60 #include <solver/slv_common.h>
61 #include <solver/linsol.h>
62 #include <solver/linsolqr.h>
63 #include <solver/slv_client.h>
64 #include <solver/system.h>
65 #include <solver/slv_interface.h>
66 #include <solver/slvDOF.h>
67 #include <solver/slv3.h>
68 #include <solver/slv_stdcalls.h>
69 #include <solver/slv_server.h>
70 }
71
72 #include "simulation.h"
73 #include "solver.h"
74 #include "solverparameters.h"
75 #include "name.h"
76 #include "incidencematrix.h"
77 #include "variable.h"
78 #include "solverstatus.h"
79 #include "solverreporter.h"
80
81 /**
82 Create an instance of a type (call compiler etc)
83
84 @TODO fix mutex on compile command filenames
85 */
86 Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
87 is_built = false;
88 // Create an Instance object for the 'simulation root' (we'll call
89 // it the 'simulation model') and it can be fetched using 'getModel()'
90 // any time later.
91 //simroot = Instanc(GetSimulationRoot(i),name);
92 }
93
94 Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
95 is_built = old.is_built;
96 sys = old.sys;
97 bin_srcname = old.bin_srcname;
98 bin_objname = old.bin_objname;
99 bin_libname = old.bin_libname;
100 bin_cmd = old.bin_cmd;
101 bin_rm = old.bin_rm;
102 sing = NULL;
103 }
104
105 Simulation::~Simulation(){
106 //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
107 }
108
109 Instanc &
110 Simulation::getModel(){
111 if(!simroot.getInternalType()){
112 throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
113 }
114 return simroot;
115 }
116
117
118 slv_system_structure *
119 Simulation::getSystem(){
120 if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
121 return sys;
122 }
123
124
125 const string
126 Simulation::getInstanceName(const Instanc &i) const{
127 char *n;
128 n = WriteInstanceNameString(i.getInternalType(),simroot.getInternalType());
129 string s(n);
130 ascfree(n);
131 return s;
132 }
133
134 const int
135 Simulation::getNumVars(){
136 return slv_get_num_solvers_vars(getSystem());
137 }
138
139
140 void
141 Simulation::write(){
142 simroot.write();
143 }
144
145 //------------------------------------------------------------------------------
146 // RUNNING MODEL 'METHODS'
147
148 void
149 Simulation::run(const Method &method){
150 Instanc &model = getModel();
151 this->run(method,model);
152 }
153
154 void
155 Simulation::runDefaultMethod(){
156 Method m = getType().getMethod(SymChar("on_load"));
157 run(m);
158 }
159
160 void
161 Simulation::run(const Method &method, Instanc &model){
162
163 // set the 'sim' pointer to our local variable...
164 CONSOLE_DEBUG("Setting shared pointer 'sim'");
165 importhandler_setsharedpointer("sim",this);
166
167 if(not is_built){
168 CONSOLE_DEBUG("WARNING, SIMULATION NOT YET BUILT");
169 }
170
171 CONSOLE_DEBUG("Running method %s...", method.getName());
172
173 Nam name = Nam(method.getSym());
174 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
175
176
177 Proc_enum pe;
178 pe = Initialize(
179 &*(model.getInternalType()) ,name.getInternalType(), "__not_named__"
180 ,ASCERR
181 ,0, NULL, NULL
182 );
183
184 // clear out the 'sim' pointer (soon it will be invalid)
185 importhandler_setsharedpointer("sim",NULL);
186 CONSOLE_DEBUG("Cleared shared pointer 'sim'");
187
188 if(pe == Proc_all_ok){
189 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Method '%s' was run (check above for errors)\n",method.getName());
190 //cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
191 }else{
192 stringstream ss;
193 ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
194 switch(pe){
195 case Proc_CallOK: ss << "Call OK"; break;
196 case Proc_CallError: ss << "Error occurred in call"; break;
197 case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
198 case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
199 case Proc_CallContinue: ss << "Skip to next iteration"; break;
200
201 case Proc_break: ss << "Break"; break;
202 case Proc_continue: ss << "Continue"; break;
203 case Proc_fallthru: ss << "Fall-through"; break;
204 case Proc_return: ss << "Return"; break;
205 case Proc_stop: ss << "Stop"; break;
206 case Proc_stack_exceeded: ss << "Stack exceeded"; break;
207 case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
208 case Proc_case_matched: ss << "Case matched"; break;
209 case Proc_case_unmatched: ss << "Case unmatched"; break;
210
211 case Proc_case_undefined_value: ss << "Undefined value in case"; break;
212 case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
213 case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
214 case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
215 case Proc_case_wrong_index: ss << "Wrong index in case"; break;
216 case Proc_case_wrong_value: ss << "Wrong value in case"; break;
217 case Proc_case_extra_values: ss << "Extra values in case"; break;
218 case Proc_bad_statement: ss << "Bad statement"; break;
219 case Proc_bad_name: ss << "Bad name"; break;
220 case Proc_for_duplicate_index: ss << "Duplicate index"; break;
221 case Proc_for_set_err: ss << "For set error"; break;
222 case Proc_for_not_set: ss << "For not set"; break;
223 case Proc_illegal_name_use: ss << "Illegal name use"; break;
224 case Proc_name_not_found: ss << "Name not found"; break;
225 case Proc_instance_not_found: ss << "Instance not found"; break;
226 case Proc_type_not_found: ss << "Type not found"; break;
227 case Proc_illegal_type_use: ss << "Illegal use"; break;
228 case Proc_proc_not_found: ss << "Method not found"; break;
229 case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
230 case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
231 case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
232 case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
233 case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
234 case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
235 case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
236 case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
237 case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
238 case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
239 case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
240 case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
241 case Proc_if_not_logical: ss << "If expression is not logical"; break;
242 case Proc_user_interrupt: ss << "User interrupt"; break;
243 case Proc_infinite_loop: ss << "Infinite loop"; break;
244 case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
245 case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
246 case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
247 case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
248 case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
249 case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
250 case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
251 case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
252 case Proc_lhs_error: ss << "Left-hand-side error"; break;
253 case Proc_rhs_error: ss << "Right-hand-side error"; break;
254 case Proc_unknown_error: ss << "Unknown error"; break;
255 default:
256 ss << "Invalid error code";
257 }
258
259 ss << " (" << int(pe) << ")";
260 throw runtime_error(ss.str());
261 }
262 }
263
264 //-----------------------------------------------------------------------------
265 // CHECKING METHODS
266
267 /**
268 Check that all the analysis went OK: solver lists are all there, etc...?
269
270 Can't return anything here because of limitations in the C API
271
272 @TODO there's something wrong with this at the moment: even after 'FIX'
273 methods are run, check shows them as not fixed, up until the point that 'SOLVE'
274 successfully completes. Something's not being synchronised properly...
275 */
276 void
277 Simulation::checkInstance(){
278 //cerr << "CHECKING SIMULATION INSTANCE" << endl;
279 if(!is_built){
280 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Simulation has not been built");
281 return;
282 }
283 Instance *i1 = getModel().getInternalType();
284 CheckInstance(stderr, &*i1);
285 //cerr << "DONE CHECKING INSTANCE" << endl;
286 }
287
288 /**
289 @return 1 = underspecified, 2 = square, 3 = structurally singular, 4 = overspecified
290 */
291 enum StructuralStatus
292 Simulation::checkDoF() const{
293 int dof, status;
294
295 if(!is_built){
296 throw runtime_error("System not yet built");
297 }
298 CONSOLE_DEBUG("Calling slvDOF_status...");
299 slvDOF_status(sys, &status, &dof);
300 switch(status){
301 case ASCXX_DOF_UNDERSPECIFIED:
302 case ASCXX_DOF_SQUARE:
303 case ASCXX_DOF_OVERSPECIFIED:
304 case ASCXX_DOF_STRUCT_SINGULAR:
305 return (enum StructuralStatus)status;
306 case 5:
307 throw runtime_error("Unable to resolve degrees of freedom"); break;
308 default:
309 throw runtime_error("Invalid return status from slvDOF_status");
310 }
311 }
312
313 /**
314 Check consistency
315
316 @TODO what is the difference between this and checkStructuralSingularity?
317
318 @return list of freeable variables. List will be empty if sys is consistent.
319 */
320 vector<Variable>
321 Simulation::getFreeableVariables(){
322 vector<Variable> v;
323
324 //cerr << "CHECKING CONSISTENCY..." << endl;
325 int *fixedarrayptr=NULL;
326
327 if(!sys){
328 throw runtime_error("System not yet built");
329 }
330
331 int res = consistency_analysis(sys, &fixedarrayptr);
332
333 if(res==1){
334 cerr << "STRUCTURALLY CONSISTENT" << endl;
335 }else{
336 if(fixedarrayptr ==NULL){
337 ERROR_REPORTER_HERE(ASC_USER_ERROR,"STRUCTURALLY INCONSISTENT");
338 throw runtime_error("Invalid consistency analysis result returned!");
339 }
340
341 struct var_variable **vp = slv_get_master_var_list(sys);
342 for(int i=0; fixedarrayptr[i]!=-1; ++i){
343 v.push_back( Variable(this, vp[fixedarrayptr[i]]) );
344 }
345 }
346 return v;
347 }
348
349 /** Returns TRUE if all is OK (not singular) */
350 bool
351 Simulation::checkStructuralSingularity(){
352 int *vil;
353 int *ril;
354 int *fil;
355
356 cerr << "RETRIEVING slfDOF_structsing INFO" << endl;
357
358 int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
359 struct var_variable **varlist = slv_get_solvers_var_list(sys);
360 struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
361
362 if(this->sing){
363 cerr << "DELETING OLD SINGULATING INFO" << endl;
364 delete this->sing;
365 this->sing = NULL;
366 }
367
368 if(res==1){
369 CONSOLE_DEBUG("processing singularity data...");
370 sing = new SingularityInfo();
371
372 // pull in the lists of vars and rels, and the freeable vars:
373 for(int i=0; ril[i]!=-1; ++i){
374 sing->rels.push_back( Relation(this, rellist[ril[i]]) );
375 }
376
377 for(int i=0; vil[i]!=-1; ++i){
378 sing->vars.push_back( Variable(this, varlist[vil[i]]) );
379 }
380
381 for(int i=0; fil[i]!=-1; ++i){
382 sing->freeablevars.push_back( Variable(this, varlist[fil[i]]) );
383 }
384
385 // we're done with those lists now
386 ASC_FREE(vil);
387 ASC_FREE(ril);
388 ASC_FREE(fil);
389
390 if(sing->isSingular()){
391 CONSOLE_DEBUG("singularity found");
392 this->sing = sing;
393 return FALSE;
394 }
395 CONSOLE_DEBUG("no singularity");
396 delete sing;
397 return TRUE;
398 }else{
399 if(res==0){
400 throw runtime_error("Unable to determine singularity lists");
401 }else{
402 throw runtime_error("Invalid return from slvDOF_structsing.");
403 }
404 }
405 }
406
407 /**
408 If the checkStructuralSingularity analysis has been done,
409 this funciton will let you access the SingularityInfo data that was
410 stored.
411 */
412 const SingularityInfo &
413 Simulation::getSingularityInfo() const{
414 if(sing==NULL){
415 throw runtime_error("No singularity info present");
416 }
417 return *sing;
418 }
419
420 //------------------------------------------
421 // ASSIGNING SOLVER TO SIMULATION
422
423 void
424 Simulation::setSolver(Solver &solver){
425 /* CONSOLE_DEBUG("Setting solver on sim %p, root inst %p",this,this->simroot.getInternalType()); */
426 if(!is_built)throw runtime_error("Can't setSolver: simulation has not been built yet");
427 if(!sys)throw runtime_error("Can't setSolver: 'sys' is not assigned.");
428
429 // Update the solver object because sometimes an alternative solver can be returned, apparently.
430
431 CONSOLE_DEBUG("Calling slv_select_solver...");
432 int selected = slv_select_solver(sys, solver.getIndex());
433 //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
434
435 if(selected<0){
436 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
437 throw runtime_error("Failed to select solver");
438 }
439
440 if(selected!=solver.getIndex()){
441 solver = Solver(slv_solver_name(selected));
442 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
443 }
444
445 if( slv_eligible_solver(sys) <= 0){
446 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
447 throw runtime_error("Inelegible solver");
448 }
449 }
450
451 const Solver
452 Simulation::getSolver() const{
453 int index = slv_get_selected_solver(sys);
454 //cerr << "Simulation::getSolver: index = " << index << endl;
455 if(index<0)throw runtime_error("No solver selected");
456
457 return Solver(slv_solver_name(index));
458 }
459
460 //------------------------------------------------------------------------------
461 // BUILD THE SYSTEM (SEND IT TO THE SOLVER)
462
463 /**
464 Build the system (send it to the solver)
465 */
466 void
467 Simulation::build(){
468 //cerr << "BUILDING SIMULATION..." << endl;
469 if(is_built){
470 CONSOLE_DEBUG("Note: rebuilding system (was already built)");
471 }
472
473 CONSOLE_DEBUG("Calling system_build...");
474 sys = system_build(simroot.getInternalType());
475 if(!sys){
476 throw runtime_error("Unable to build system");
477 }
478 is_built = true;
479 //cerr << "...DONE BUILDING" << endl;
480 }
481
482
483 //------------------------------------------------------------------------------
484 // SOLVER CONFIGURATION PARAMETERS
485
486 /**
487 Get solver parameters struct wrapped up as a SolverParameters class.
488 */
489 SolverParameters
490 Simulation::getSolverParameters() const{
491 if(!is_built)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
492 if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has no 'sys' assigned.");
493
494 slv_parameters_t p;
495 slv_get_parameters(sys,&p);
496 return SolverParameters(p);
497 }
498
499 /**
500 Update the solver parameters by passing a new set back
501 */
502 void
503 Simulation::setSolverParameters(SolverParameters &P){
504 if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
505 CONSOLE_DEBUG("Calling slv_set_parameters");
506 slv_set_parameters(sys, &(P.getInternalType()));
507 }
508
509 //------------------------------------------------------------------------------
510 // PRE-SOLVE DIAGNOSTICS
511
512 /**
513 Get a list of variables to fix to make an underspecified system
514 become square. Also seems to return stuff when you have a structurally
515 singuler system.
516 */
517 vector<Variable>
518 Simulation::getFixableVariables(){
519 //cerr << "GETTING FIXABLE VARIABLES..." << endl;
520 vector<Variable> vars;
521
522 if(!sys){
523 throw runtime_error("Simulation system not yet built");
524 }
525
526 int32 *vip; /** TODO ensure 32 bit integers are used */
527
528 // Get IDs of elegible variables in array at vip...
529 CONSOLE_DEBUG("Calling slvDOF_eligible");
530 if(!slvDOF_eligible(sys,&vip)){
531 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
532 }else{
533 struct var_variable **vp = slv_get_solvers_var_list(sys);
534
535 if(vp==NULL){
536 throw runtime_error("Simulation variable list is null");
537 }
538
539 // iterate through this list until we find a -1:
540 int i=0;
541 int var_index = vip[i];
542 while(var_index >= 0){
543 struct var_variable *var = vp[var_index];
544 vars.push_back( Variable(this, var) );
545 ++i;
546 var_index = vip[i];
547 }
548 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
549 ascfree(vip);
550 }
551
552 return vars;
553 }
554
555 /**
556 Get the list of variables near their bounds. Helps to indentify why
557 you might be having non-convergence problems.
558 */
559 vector<Variable>
560 Simulation::getVariablesNearBounds(const double &epsilon){
561 //cerr << "GETTING VARIABLES NEAR BOUNDS..." << endl;
562 vector<Variable> vars;
563
564 if(!sys){
565 throw runtime_error("Simulation system not yet built");
566 }
567
568 int *vip;
569 CONSOLE_DEBUG("Calling slv_near_bounds...");
570 if(slv_near_bounds(sys,epsilon,&vip)){
571 struct var_variable **vp = slv_get_solvers_var_list(sys);
572 struct var_variable *var;
573 cerr << "VARS FOUND NEAR BOUNDS" << endl;
574 int nlow = vip[0];
575 int nhigh = vip[1];
576 int lim1 = 2 + nlow;
577 for(int i=2; i<lim1; ++i){
578 var = vp[vip[i]];
579 char *var_name = var_make_name(sys,var);
580 cerr << "AT LOWER BOUND: " << var_name << endl;
581 ascfree(var_name);
582 vars.push_back(Variable(this,var));
583 };
584 int lim2 = lim1 + nhigh;
585 for(int i=lim1; i<lim2; ++i){
586 var = vp[vip[i]];
587 char *var_name = var_make_name(sys,var);
588 cerr << "AT UPPER BOUND: " << var_name << endl;
589 ascfree(var_name);
590 vars.push_back(Variable(this,var));
591 }
592 }
593 ascfree(vip);
594 return vars;
595 }
596
597
598 bool
599 SingularityInfo::isSingular() const{
600 if(vars.size()||rels.size()){
601 return true;
602 }
603 return false;
604 }
605
606 //------------------------------------------------------------------------------
607 // SOLVING
608
609 /**
610 Solve the system through to convergence. This function is hardwired with
611 a maximum of 1000 iterations, but will interrupt itself when the 'stop'
612 condition comes back from the SolverReporter.
613 */
614 void
615 Simulation::solve(Solver solver, SolverReporter &reporter){
616 if(!is_built){
617 throw runtime_error("Simulation::solver: simulation is not yet built, can't start solving.");
618 }
619
620 //cerr << "SIMULATION::SOLVE STARTING..." << endl;
621 enum inst_t k = getModel().getKind();
622 if(k!=MODEL_INST)throw runtime_error("Can't solve: not an instance of type MODEL_INST");
623
624 Instance *i1 = getInternalType();
625 int npend = NumberPendingInstances(&*i1);
626 if(npend)throw runtime_error("Can't solve: There are still %d pending instances");
627
628 if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
629
630 //cerr << "SIMULATION::SOLVE: SET SOLVER..." << endl;
631 setSolver(solver);
632
633
634 //cerr << "PRESOLVING SYSTEM...";
635 CONSOLE_DEBUG("Calling slv_presolve...");
636 slv_presolve(sys);
637 //cerr << "DONE" << endl;
638
639 //cerr << "SOLVING SYSTEM..." << endl;
640 // Add some stuff here for cleverer iteration....
641 unsigned niter = 1000;
642 //double updateinterval = 0.02;
643
644 double starttime = tm_cpu_time();
645 //double lastupdate = starttime;
646 SolverStatus status;
647 //int solved_vars=0;
648 bool stop=false;
649
650 status.getSimulationStatus(*this);
651 reporter.report(&status);
652
653 for(unsigned iter = 1; iter <= niter && !stop; ++iter){
654
655 if(status.isReadyToSolve()){
656 CONSOLE_DEBUG("Calling slv_iterate...");
657 slv_iterate(sys);
658 }
659
660 status.getSimulationStatus(*this);
661
662 if(reporter.report(&status)){
663 stop = true;
664 }
665 }
666
667 double elapsed = tm_cpu_time() - starttime;
668
669
670 activeblock = status.getCurrentBlockNum();
671
672 reporter.finalise(&status);
673
674 // Just a little bit of console output:
675
676 if(status.isOK()){
677 //cerr << "... SOLVED, STATUS OK" << endl;
678 }else{
679 cerr << "... SOLVER FAILED" << endl;
680 }
681
682 //cerr << "SOLVER PERFORMED " << status.getIterationNum() << " ITERATIONS IN " << elapsed << "s" << endl;
683
684 // communicate solver variable status back to the instance tree
685 processVarStatus();
686 }
687
688 //------------------------------------------------------------------------------
689 // POST-SOLVE DIAGNOSTICS
690
691 const int
692 Simulation::getActiveBlock() const{
693 return activeblock;
694 }
695
696 /**
697 Return an IncidenceMatrix built from the current state of the solver system.
698
699 This will actually return something meaningful even before solve.
700 */
701 IncidenceMatrix
702 Simulation::getIncidenceMatrix(){
703 return IncidenceMatrix(*this);
704 }
705
706 /**
707 This function looks at all the variables in the solve's list and updates
708 the variable status for the corresponding instances.
709
710 It does this by using the 'interface pointer' in the Instance, see
711 the C-API function GetInterfacePtr.
712
713 This is used to display visually which variables have been solved, which
714 ones have not yet been attempted, and which ones were active when the solver
715 failed (ASCXX_VAR_ACTIVE).
716 */
717 void
718 Simulation::processVarStatus(){
719
720 // this is a cheap function call:
721 const mtx_block_t *bb = slv_get_solvers_blocks(getSystem());
722
723 var_variable **vlist = slv_get_solvers_var_list(getSystem());
724 int nvars = slv_get_num_solvers_vars(getSystem());
725
726 slv_status_t status;
727 slv_get_status(getSystem(), &status);
728
729 if(status.block.number_of == 0){
730 cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
731 return;
732 }
733
734 int activeblock = status.block.current_block;
735 int low = bb->block[activeblock].col.low;
736 int high = bb->block[activeblock].col.high;
737 bool allsolved = status.converged;
738 for(int c=0; c < nvars; ++c){
739 var_variable *v = vlist[c];
740 Instanc i((Instance *)var_instance(v));
741 VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
742 if(i.isFixed()){
743 s = ASCXX_VAR_FIXED;
744 }else if(var_incident(v) && var_active(v)){
745 if(allsolved || c < low){
746 s = ASCXX_VAR_SOLVED;
747 }else if(c <= high){
748 s = ASCXX_VAR_ACTIVE;
749 }else{
750 s = ASCXX_VAR_UNSOLVED;
751 }
752 }
753 i.setVarStatus(s);
754 }
755 }
756

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