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

Contents of /trunk/pygtk/interface/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 319 - (show annotations) (download) (as text)
Thu Feb 23 12:30:40 2006 UTC (14 years, 9 months ago) by johnpye
File MIME type: text/x-c++src
File size: 18455 byte(s)
Fixed a bug with the diagnose window after a successful solve
Added preference: SolverReporter:close_on_converged
1 #include <iostream>
2 #include <iomanip>
3 #include <stdexcept>
4 #include <sstream>
5 using namespace std;
6
7 extern "C"{
8 #include <utilities/ascConfig.h>
9 #include <utilities/ascSignal.h>
10 #include <utilities/ascMalloc.h>
11 #include <general/dstring.h>
12 #include <general/tm_time.h>
13 #include <compiler/instance_enum.h>
14 #include <compiler/fractions.h>
15 #include <compiler/compiler.h>
16 #include <compiler/dimen.h>
17 #include <compiler/symtab.h>
18 #include <compiler/instance_io.h>
19 #include <compiler/instantiate.h>
20 #include <compiler/bintoken.h>
21 #include <compiler/instance_enum.h>
22 #include <compiler/instquery.h>
23 #include <compiler/check.h>
24 #include <compiler/name.h>
25 #include <compiler/pending.h>
26
27 #include <utilities/readln.h>
28 #include <solver/mtx.h>
29 #include <solver/slv_types.h>
30 #include <solver/var.h>
31 #include <solver/rel.h>
32 #include <solver/discrete.h>
33 #include <solver/conditional.h>
34 #include <solver/logrel.h>
35 #include <solver/bnd.h>
36 #include <solver/calc.h>
37 #include <solver/relman.h>
38 #include <solver/slv_common.h>
39 #include <solver/linsol.h>
40 #include <solver/linsolqr.h>
41 #include <solver/slv_client.h>
42 #include <solver/system.h>
43 #include <solver/slv_interface.h>
44 #include <solver/slvDOF.h>
45 #include <solver/slv3.h>
46 #include <solver/slv_stdcalls.h>
47 }
48
49 #include "simulation.h"
50 #include "solver.h"
51 #include "solverparameters.h"
52 #include "name.h"
53 #include "incidencematrix.h"
54 #include "variable.h"
55 #include "solverstatus.h"
56 #include "solverreporter.h"
57
58 /**
59 Create an instance of a type (call compiler etc)
60
61 @TODO fix mutex on compile command filenames
62 */
63 Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
64 is_built = false;
65 // Create an Instance object for the 'simulation root' (we'll call
66 // it the 'simulation model') and it can be fetched using 'getModel()'
67 // any time later.
68 //simroot = Instanc(GetSimulationRoot(i),name);
69 }
70
71 Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
72 is_built = old.is_built;
73 sys = old.sys;
74 bin_srcname = old.bin_srcname;
75 bin_objname = old.bin_objname;
76 bin_libname = old.bin_libname;
77 bin_cmd = old.bin_cmd;
78 bin_rm = old.bin_rm;
79 }
80
81 Simulation::~Simulation(){
82 //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
83 }
84
85 Instanc &
86 Simulation::getModel(){
87 if(!simroot.getInternalType()){
88 throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
89 }
90 return simroot;
91 }
92
93 void
94 Simulation::checkDoF() const{
95 cerr << "CHECKING DOF..." << endl;
96 int dof, status;
97 if(!sys){
98 throw runtime_error("System not yet built");
99 }
100 slvDOF_status(sys, &status, &dof);
101 switch(status){
102 case 1: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Underspecified; %d degrees of freedom",dof); break;
103 case 2: ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Square"); break;
104 case 3: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally singular"); break;
105 case 4: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Overspecified"); break;
106 case 5:
107 throw runtime_error("Unable to resolve degrees of freedom"); break;
108 default:
109 throw runtime_error("Invalid return status from slvDOF_status");
110 }
111 }
112
113 void
114 Simulation::checkConsistency() const{
115 cerr << "CHECKING CONSISTENCY..." << endl;
116 int *fixedarrayptr;
117
118 int res = consistency_analysis(sys, &fixedarrayptr);
119 struct var_variable **vp = slv_get_master_var_list(sys);
120
121 if(res==1){
122 cerr << "STRUCTURALLY CONSISTENT" << endl;
123 return;
124 }else{
125 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally inconsistent. Free the variables listed on the console\nin order to make system consistent.");
126 cerr << "INCONSISTENT: Free these vars:" << endl;
127 for(int i=0; fixedarrayptr[i]!=-1; ++i){
128 Instanc i1((struct Instance *)var_instance(vp[fixedarrayptr[i]]));
129 cerr << " " << getInstanceName(i1) << endl;
130 }
131 }
132 }
133
134 void
135 Simulation::checkStructuralSingularity() const{
136 cerr << "CHECKING STRUCTURAL SINGULARITY..." << endl;
137
138 int *vil;
139 int *ril;
140 int *fil;
141
142 int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
143 struct var_variable **varlist = slv_get_solvers_var_list(sys);
144 struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
145
146 if(res==0){
147 cerr << "UNABLE TO DETERMINE SINGULARITY LISTS" << endl;
148 return;
149 }else if(res==1){
150 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally singular. Check the listing on the console.");
151 cerr << "STRUCTURALLY SINGULAR: The found singularity involves these relations:" << endl;
152 for(int i=0; ril[i]!=-1; ++i){
153 Instanc i1((struct Instance *)rel_instance(rellist[ril[i]]));
154 cerr << " " << getInstanceName(i1) << endl;
155 }
156
157 cerr << "STRUCTURALLY SINGULAR: ... and these variables:" << endl;
158 for(int i=0; vil[i]!=-1; ++i){
159 Instanc i1((struct Instance *)var_instance(varlist[vil[i]]));
160 cerr << " " << getInstanceName(i1) << endl;
161 }
162
163 cerr << "STRUCTURALLY SINGULAR: ... and may be mitigated by freeing these variables:" << endl;
164 for(int i=0; fil[i]!=-1; ++i){
165 Instanc i1((struct Instance *)var_instance(varlist[fil[i]]));
166 cerr << " " << getInstanceName(i1) << endl;
167 }
168 }else{
169 throw runtime_error("Invalid return from slvDOF_structsing.");
170 }
171 ascfree(vil);
172 ascfree(ril);
173 ascfree(fil);
174 }
175
176 void
177 Simulation::run(const Method &method){
178 cerr << "RUNNING PROCEDURE " << method.getName() << endl;
179 Nam name = Nam(method.getSym());
180 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
181 Proc_enum pe;
182 pe = Initialize(
183 &*(getModel().getInternalType()) ,name.getInternalType(), "__not_named__"
184 ,ASCERR
185 ,0, NULL, NULL
186 );
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
260 ss << " (" << int(pe) << ")";
261 throw runtime_error(ss.str());
262 }
263 }
264
265 const bool
266 Simulation::check(){
267 cerr << "CHECKING SIMULATION" << endl;
268 Instance *i1 = getModel().getInternalType();
269 CheckInstance(stderr, &*i1);
270 cerr << "...DONE CHECKING" << endl;
271 this->checkConsistency();
272 this->checkStructuralSingularity();
273 }
274
275 void
276 Simulation::build(){
277 cerr << "BUILDING SIMULATION..." << endl;
278 Instance *i1 = getModel().getInternalType();
279 sys = system_build(&*i1);
280 if(!sys){
281 throw runtime_error("Unable to build system");
282 }
283 is_built = true;
284 cerr << "...DONE BUILDING" << endl;
285 }
286
287 vector<Variable>
288 Simulation::getFixableVariables(){
289 cerr << "GETTING FIXABLE VARIABLES..." << endl;
290 vector<Variable> vars;
291
292 if(!sys){
293 throw runtime_error("Simulation system not yet built");
294 }
295
296 int32 *vip; /** TODO ensure 32 bit integers are used */
297
298 // Get IDs of elegible variables in array at vip...
299 if(!slvDOF_eligible(sys,&vip)){
300 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
301 }else{
302 //cerr << "FIXABLE VARS FOUND" << endl;
303 struct var_variable **vp = slv_get_solvers_var_list(sys);
304
305 /*struct var_variable *first_var = vp[0];
306 char *first_var_name = var_make_name(sys,first_var);
307 cerr << "FIRST SYS VAR IS NAMED " << var_make_name(s,first_var) << endl;
308 ascfree(first_var_name);*/
309
310 if(vp==NULL){
311 throw runtime_error("Simulation variable list is null");
312 }
313
314 // iterate through this list until we find a -1:
315 int i=0;
316 int var_index = vip[i];
317 while(var_index >= 0){
318 //cerr << "FOUND VARIABLE var_index = " << var_index << endl;
319 struct var_variable *var = vp[var_index];
320 //cerr << "VARIABLE " << var_index << " IS ELIGIBLE" << endl;
321 char *var_name = var_make_name(sys,var);
322 //cerr << "ELIGIBLE VAR: " << var_name << endl;
323 ascfree(var_name);
324 vars.push_back( Variable(this, var) );
325 ++i;
326 var_index = vip[i];
327 }
328 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
329 //cerr << "END ELEGIBLE VARS LIST" << endl;
330 ascfree(vip);
331 //cerr << "FREED VIP LIST" << endl;
332 }
333
334 //cerr << "FINISHED WITH FINDING ELEGIBLE VARIABLES" << endl;
335 return vars;
336 }
337
338
339 void
340 Simulation::solve(Solver solver, SolverReporter &reporter){
341 if(!is_built){
342 throw runtime_error("Simulation::solver: simulation is not yet built, can't start solving.");
343 }
344
345 cerr << "SIMULATION::SOLVE STARTING..." << endl;
346 enum inst_t k = getModel().getKind();
347 if(k!=MODEL_INST)throw runtime_error("Can't solve: not an instance of type MODEL_INST");
348
349 Instance *i1 = getInternalType();
350 int npend = NumberPendingInstances(&*i1);
351 if(npend)throw runtime_error("Can't solve: There are still %d pending instances");
352
353 if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
354
355 cerr << "SIMULATION::SOLVE: SET SOLVER..." << endl;
356 setSolver(solver);
357
358
359 cerr << "PRESOLVING SYSTEM...";
360 slv_presolve(sys);
361 cerr << "DONE" << endl;
362
363 cerr << "SOLVING SYSTEM..." << endl;
364 // Add some stuff here for cleverer iteration....
365 unsigned niter = 1000;
366 double updateinterval = 0.02;
367
368 double starttime = tm_cpu_time();
369 double lastupdate = starttime;
370 SolverStatus status;
371 int solved_vars=0;
372 bool stop=false;
373
374 status.getSimulationStatus(*this);
375 reporter.report(&status);
376
377 for(int iter = 1; iter <= niter && !stop; ++iter){
378
379 if(status.isReadyToSolve()){
380 slv_iterate(sys);
381 }
382
383 status.getSimulationStatus(*this);
384
385 if(reporter.report(&status)){
386 stop = true;
387 }
388 }
389
390 double elapsed = tm_cpu_time() - starttime;
391
392
393 activeblock = status.getCurrentBlockNum();
394
395 reporter.finalise(&status);
396
397 if(status.isOK()){
398 cerr << "... SOLVED, STATUS OK" << endl;
399 }else{
400 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Solver failed");
401 }
402
403 cerr << "SOLVER PERFORMED " << status.getIterationNum() << " ITERATIONS IN " << elapsed << "s" << endl;
404
405 if(status.hasExceededTimeLimit()){
406 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Exceeded interation limit");
407 }
408
409 if(status.isConverged()){
410 ERROR_REPORTER_NOLINE(ASC_USER_SUCCESS,"Solver converged: %d iterations (%.2f s)"
411 ,status.getIterationNum(),elapsed);
412 }else{
413 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Solver not converged in block %d after overall %d iterations (see console)"
414 " (%.2f s).",status.getCurrentBlockNum(),status.getIterationNum(),elapsed);
415 IncidenceMatrix inc = getIncidenceMatrix();
416
417 /*
418 cerr << "VARIABLES IN NON-CONVERGED BLOCK:" << endl;
419 vector<Variable> v = inc.getBlockVars(status.block.current_block);
420 for(vector<Variable>::iterator vi = v.begin(); vi < v.end(); ++vi){
421 cerr << vi->getName() << " = " << vi->getValue() << endl;
422 }
423
424 cerr << "RELATIONS IN NON-CONVERGED BLOCK:" << endl;
425 vector<Relation> r = inc.getBlockRels(status.block.current_block);
426 for(vector<Relation>::iterator ri = r.begin(); ri < r.end(); ++ri){
427 cerr << ri->getName() << endl;
428 }
429 */
430 }
431
432 }
433
434 void
435 Simulation::write(){
436 simroot.write();
437 }
438
439 //------------------------------------------
440 // ASSIGNING SOLVER TO SIMULATION
441
442 void
443 Simulation::setSolver(Solver &solver){
444 cerr << "SETTING SOLVER ON SIMULATION TO " << solver.getName() << endl;
445
446 if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
447 // Update the solver object because sometimes an alternative solver can be returned, apparently.
448
449 int selected = slv_select_solver(sys, solver.getIndex());
450 //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
451
452 if(selected<0){
453 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
454 throw runtime_error("Failed to select solver");
455 }
456
457 if(selected!=solver.getIndex()){
458 solver = Solver(slv_solver_name(selected));
459 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
460 }
461
462 if( slv_eligible_solver(sys) <= 0){
463 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
464 throw runtime_error("Inelegible solver");
465 }
466 }
467
468 const Solver
469 Simulation::getSolver() const{
470 int index = slv_get_selected_solver(sys);
471 //cerr << "Simulation::getSolver: index = " << index << endl;
472 if(index<0)throw runtime_error("No solver selected");
473
474 return Solver(slv_solver_name(index));
475 }
476
477
478 /**
479 Get solver parameters struct wrapped up as a SolverParameters class.
480 */
481 SolverParameters
482 Simulation::getSolverParameters() const{
483 if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
484
485 slv_parameters_t p;
486 slv_get_parameters(sys,&p);
487 return SolverParameters(p);
488 }
489
490 /**
491 Update the solver parameters by passing a new set back
492 */
493 void
494 Simulation::setSolverParameters(SolverParameters &P){
495 if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
496 slv_set_parameters(sys, &(P.getInternalType()));
497 }
498
499 slv_system_structure *
500 Simulation::getSystem(){
501 if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
502 return sys;
503 }
504
505 IncidenceMatrix
506 Simulation::getIncidenceMatrix(){
507 return IncidenceMatrix(*this);
508 }
509
510 const string
511 Simulation::getInstanceName(const Instanc &i) const{
512 char *n;
513 n = WriteInstanceNameString(i.getInternalType(),simroot.getInternalType());
514 string s(n);
515 ascfree(n);
516 return s;
517 }
518
519 const int
520 Simulation::getNumVars(){
521 return slv_get_num_solvers_vars(getSystem());
522 }
523
524 void
525 Simulation::processVarStatus(){
526
527 // this is a cheap function call:
528 const mtx_block_t *bb = slv_get_solvers_blocks(getSystem());
529
530 var_variable **vlist = slv_get_solvers_var_list(getSystem());
531 int nvars = slv_get_num_solvers_vars(getSystem());
532
533 slv_status_t status;
534 slv_get_status(getSystem(), &status);
535
536 if(status.block.number_of == 0){
537 cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
538 return;
539 }
540
541 int activeblock = status.block.current_block;
542 int low = bb->block[activeblock].col.low;
543 int high = bb->block[activeblock].col.high;
544 bool allsolved = status.converged;
545 for(int c=0; c < nvars; ++c){
546 var_variable *v = vlist[c];
547 Instanc i((Instance *)var_instance(v));
548 VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
549 if(i.isFixed()){
550 s = ASCXX_VAR_FIXED;
551 }else if(var_incident(v) && var_active(v)){
552 if(allsolved || c < low){
553 s = ASCXX_VAR_SOLVED;
554 }else if(c <= high){
555 s = ASCXX_VAR_ACTIVE;
556 }else{
557 s = ASCXX_VAR_UNSOLVED;
558 }
559 }
560 i.setVarStatus(s);
561 }
562 }
563
564 const int
565 Simulation::getActiveBlock() const{
566 return activeblock;
567 }

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