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

Contents of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 774 - (show annotations) (download) (as text)
Fri Jul 14 08:01:48 2006 UTC (13 years, 7 months ago) by johnpye
File MIME type: text/x-c++src
File size: 19036 byte(s)
Added the ability to run methods on sub-models within a simulation. Use the right-click context menu.
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
48 #include <utilities/readln.h>
49 #include <solver/mtx.h>
50 #include <solver/slv_types.h>
51 #include <solver/var.h>
52 #include <solver/rel.h>
53 #include <solver/discrete.h>
54 #include <solver/conditional.h>
55 #include <solver/logrel.h>
56 #include <solver/bnd.h>
57 #include <solver/calc.h>
58 #include <solver/relman.h>
59 #include <solver/slv_common.h>
60 #include <solver/linsol.h>
61 #include <solver/linsolqr.h>
62 #include <solver/slv_client.h>
63 #include <solver/system.h>
64 #include <solver/slv_interface.h>
65 #include <solver/slvDOF.h>
66 #include <solver/slv3.h>
67 #include <solver/slv_stdcalls.h>
68 #include <solver/slv_server.h>
69 }
70
71 #include "simulation.h"
72 #include "solver.h"
73 #include "solverparameters.h"
74 #include "name.h"
75 #include "incidencematrix.h"
76 #include "variable.h"
77 #include "solverstatus.h"
78 #include "solverreporter.h"
79
80 /**
81 Create an instance of a type (call compiler etc)
82
83 @TODO fix mutex on compile command filenames
84 */
85 Simulation::Simulation(Instance *i, const SymChar &name) : Instanc(i, name), simroot(GetSimulationRoot(i),SymChar("simroot")){
86 is_built = false;
87 // Create an Instance object for the 'simulation root' (we'll call
88 // it the 'simulation model') and it can be fetched using 'getModel()'
89 // any time later.
90 //simroot = Instanc(GetSimulationRoot(i),name);
91 }
92
93 Simulation::Simulation(const Simulation &old) : Instanc(old), simroot(old.simroot){
94 is_built = old.is_built;
95 sys = old.sys;
96 bin_srcname = old.bin_srcname;
97 bin_objname = old.bin_objname;
98 bin_libname = old.bin_libname;
99 bin_cmd = old.bin_cmd;
100 bin_rm = old.bin_rm;
101 }
102
103 Simulation::~Simulation(){
104 //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
105 }
106
107 Instanc &
108 Simulation::getModel(){
109 if(!simroot.getInternalType()){
110 throw runtime_error("Simulation::getModel: simroot.getInternalType()is NULL");
111 }
112 return simroot;
113 }
114
115 void
116 Simulation::checkDoF() const{
117 cerr << "CHECKING DOF..." << endl;
118 int dof, status;
119 if(!sys){
120 throw runtime_error("System not yet built");
121 }
122 slvDOF_status(sys, &status, &dof);
123 switch(status){
124 case 1: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Underspecified; %d degrees of freedom",dof); break;
125 case 2: ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Square"); break;
126 case 3: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally singular"); break;
127 case 4: ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Overspecified"); break;
128 case 5:
129 throw runtime_error("Unable to resolve degrees of freedom"); break;
130 default:
131 throw runtime_error("Invalid return status from slvDOF_status");
132 }
133 }
134
135 void
136 Simulation::checkConsistency() const{
137 cerr << "CHECKING CONSISTENCY..." << endl;
138 int *fixedarrayptr;
139
140 int res = consistency_analysis(sys, &fixedarrayptr);
141 struct var_variable **vp = slv_get_master_var_list(sys);
142
143 if(res==1){
144 cerr << "STRUCTURALLY CONSISTENT" << endl;
145 return;
146 }else{
147 ERROR_REPORTER_NOLINE(ASC_USER_ERROR,"Structurally inconsistent. Free the variables listed on the console\nin order to make system consistent.");
148 cerr << "INCONSISTENT: Free these vars:" << endl;
149 for(int i=0; fixedarrayptr[i]!=-1; ++i){
150 Instanc i1((struct Instance *)var_instance(vp[fixedarrayptr[i]]));
151 cerr << " " << getInstanceName(i1) << endl;
152 }
153 }
154 }
155
156 /** Returns TRUE if all is OK (not singular) */
157 bool
158 Simulation::checkStructuralSingularity(){
159 cerr << "CHECKING STRUCTURAL SINGULARITY..." << endl;
160
161 int *vil;
162 int *ril;
163 int *fil;
164
165 int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
166 struct var_variable **varlist = slv_get_solvers_var_list(sys);
167 struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
168
169 if(this->sing){
170 delete this->sing;
171 this->sing = NULL;
172 }
173
174 if(res==1){
175 CONSOLE_DEBUG("processing singularity data...");
176 sing = new SingularityInfo;
177
178 // 'singular relations'
179 for(int i=0; ril[i]!=-1; ++i){
180 sing->rels.push_back( Relation( this, rellist[ril[i]] ) );
181 }
182
183 // 'singular variables'
184 for(int i=0; vil[i]!=-1; ++i){
185 sing->vars.push_back( Variable( this, varlist[vil[i]] ) );
186 }
187
188 // 'free these variables'
189 for(int i=0; fil[i]!=-1; ++i){
190 sing->freeablevars.push_back( Variable( this, varlist[fil[i]] ) );
191 }
192
193 // we're done with those lists now
194 ASC_FREE(vil);
195 ASC_FREE(ril);
196 ASC_FREE(fil);
197
198 if(sing->isSingular()){
199 CONSOLE_DEBUG("singularity found");
200 this->sing = sing;
201 return FALSE;
202 }
203 CONSOLE_DEBUG("no singularity");
204 delete sing;
205 return TRUE;
206 }else{
207 if(res==0){
208 throw runtime_error("Unable to determine singularity lists");
209 }else{
210 throw runtime_error("Invalid return from slvDOF_structsing.");
211 }
212 }
213 }
214
215 const SingularityInfo &
216 Simulation::getSingularityInfo() const{
217 if(sing==NULL){
218 throw runtime_error("No singularity info present");
219 }
220 return *sing;
221 }
222
223 void
224 Simulation::run(const Method &method){
225 Instanc &model = getModel();
226 this->run(method,model);
227 }
228
229 void
230 Simulation::run(const Method &method, Instanc &model){
231
232 cerr << "RUNNING PROCEDURE " << method.getName() << endl;
233 Nam name = Nam(method.getSym());
234 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
235 Proc_enum pe;
236 pe = Initialize(
237 &*(model.getInternalType()) ,name.getInternalType(), "__not_named__"
238 ,ASCERR
239 ,0, NULL, NULL
240 );
241
242 if(pe == Proc_all_ok){
243 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Method '%s' was run (check above for errors)\n",method.getName());
244 //cerr << "METHOD " << method.getName() << " COMPLETED OK" << endl;
245 }else{
246 stringstream ss;
247 ss << "Simulation::run: Method '" << method.getName() << "' returned error: ";
248 switch(pe){
249 case Proc_CallOK: ss << "Call OK"; break;
250 case Proc_CallError: ss << "Error occurred in call"; break;
251 case Proc_CallReturn: ss << "Request that caller return (OK)"; break;
252 case Proc_CallBreak: ss << "Break out of enclosing loop"; break;
253 case Proc_CallContinue: ss << "Skip to next iteration"; break;
254
255 case Proc_break: ss << "Break"; break;
256 case Proc_continue: ss << "Continue"; break;
257 case Proc_fallthru: ss << "Fall-through"; break;
258 case Proc_return: ss << "Return"; break;
259 case Proc_stop: ss << "Stop"; break;
260 case Proc_stack_exceeded: ss << "Stack exceeded"; break;
261 case Proc_stack_exceeded_this_frame: ss << "Stack exceeded this frame"; break;
262 case Proc_case_matched: ss << "Case matched"; break;
263 case Proc_case_unmatched: ss << "Case unmatched"; break;
264
265 case Proc_case_undefined_value: ss << "Undefined value in case"; break;
266 case Proc_case_boolean_mismatch: ss << "Boolean mismatch in case"; break;
267 case Proc_case_integer_mismatch: ss << "Integer mismatch in case"; break;
268 case Proc_case_symbol_mismatch: ss << "Symbol mismatch in case"; break;
269 case Proc_case_wrong_index: ss << "Wrong index in case"; break;
270 case Proc_case_wrong_value: ss << "Wrong value in case"; break;
271 case Proc_case_extra_values: ss << "Extra values in case"; break;
272 case Proc_bad_statement: ss << "Bad statement"; break;
273 case Proc_bad_name: ss << "Bad name"; break;
274 case Proc_for_duplicate_index: ss << "Duplicate index"; break;
275 case Proc_for_set_err: ss << "For set error"; break;
276 case Proc_for_not_set: ss << "For not set"; break;
277 case Proc_illegal_name_use: ss << "Illegal name use"; break;
278 case Proc_name_not_found: ss << "Name not found"; break;
279 case Proc_instance_not_found: ss << "Instance not found"; break;
280 case Proc_type_not_found: ss << "Type not found"; break;
281 case Proc_illegal_type_use: ss << "Illegal use"; break;
282 case Proc_proc_not_found: ss << "Method not found"; break;
283 case Proc_if_expr_error_typeconflict: ss << "Type conflict in 'if' expression"; break;
284 case Proc_if_expr_error_nameunfound: ss << "Name not found in 'if' expression"; break;
285 case Proc_if_expr_error_incorrectname: ss << "Incorrect name in 'if' expression"; break;
286 case Proc_if_expr_error_undefinedvalue: ss << "Undefined value in 'if' expression"; break;
287 case Proc_if_expr_error_dimensionconflict: ss << "Dimension conflict in 'if' expression"; break;
288 case Proc_if_expr_error_emptychoice: ss << "Empty choice in 'if' expression"; break;
289 case Proc_if_expr_error_emptyintersection: ss << "Empty intersection in 'if' expression"; break;
290 case Proc_if_expr_error_confused: ss << "Confused in 'if' expression"; break;
291 case Proc_if_real_expr: ss << "Real-valued result in 'if' expression"; break;
292 case Proc_if_integer_expr: ss << "Integeter-valued result in 'if' expression"; break;
293 case Proc_if_symbol_expr: ss << "Symbol-valued result in 'if' expression"; break;
294 case Proc_if_set_expr: ss << "Set-valued result in 'if' expression"; break;
295 case Proc_if_not_logical: ss << "If expression is not logical"; break;
296 case Proc_user_interrupt: ss << "User interrupt"; break;
297 case Proc_infinite_loop: ss << "Infinite loop"; break;
298 case Proc_declarative_constant_assignment: ss << "Declarative constant assignment"; break;
299 case Proc_nonsense_assignment: ss << "Nonsense assginment (bogus)"; break;
300 case Proc_nonconsistent_assignment: ss << "Inconsistent assignment"; break;
301 case Proc_nonatom_assignment: ss << "Non-atom assignment"; break;
302 case Proc_nonboolean_assignment: ss << "Non-boolean assignment"; break;
303 case Proc_noninteger_assignment: ss << "Non-integer assignment"; break;
304 case Proc_nonreal_assignment: ss << "Non-real assignment"; break;
305 case Proc_nonsymbol_assignment: ss << "Non-symbol assignment"; break;
306 case Proc_lhs_error: ss << "Left-hand-side error"; break;
307 case Proc_rhs_error: ss << "Right-hand-side error"; break;
308 case Proc_unknown_error: ss << "Unknown error"; break;
309 default:
310 ss << "Invalid error code";
311 }
312
313 ss << " (" << int(pe) << ")";
314 throw runtime_error(ss.str());
315 }
316 }
317
318 /**
319 @return TRUE if all is OK
320 */
321 const bool
322 Simulation::check(){
323 cerr << "CHECKING SIMULATION" << endl;
324 Instance *i1 = getModel().getInternalType();
325 CheckInstance(stderr, &*i1);
326 cerr << "...DONE CHECKING" << endl;
327 this->checkConsistency();
328
329 return this->checkStructuralSingularity();
330 }
331
332 void
333 Simulation::build(){
334 cerr << "BUILDING SIMULATION..." << endl;
335 Instance *i1 = getModel().getInternalType();
336 sys = system_build(&*i1);
337 if(!sys){
338 throw runtime_error("Unable to build system");
339 }
340 is_built = true;
341 cerr << "...DONE BUILDING" << endl;
342 }
343
344 vector<Variable>
345 Simulation::getFixableVariables(){
346 cerr << "GETTING FIXABLE VARIABLES..." << endl;
347 vector<Variable> vars;
348
349 if(!sys){
350 throw runtime_error("Simulation system not yet built");
351 }
352
353 int32 *vip; /** TODO ensure 32 bit integers are used */
354
355 // Get IDs of elegible variables in array at vip...
356 if(!slvDOF_eligible(sys,&vip)){
357 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
358 }else{
359 struct var_variable **vp = slv_get_solvers_var_list(sys);
360
361 if(vp==NULL){
362 throw runtime_error("Simulation variable list is null");
363 }
364
365 // iterate through this list until we find a -1:
366 int i=0;
367 int var_index = vip[i];
368 while(var_index >= 0){
369 struct var_variable *var = vp[var_index];
370 vars.push_back( Variable(this, var) );
371 ++i;
372 var_index = vip[i];
373 }
374 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
375 ascfree(vip);
376 }
377
378 return vars;
379 }
380
381 vector<Variable>
382 Simulation::getVariablesNearBounds(const double &epsilon){
383 cerr << "GETTING VARIABLES NEAR BOUNDS..." << endl;
384 vector<Variable> vars;
385
386 if(!sys){
387 throw runtime_error("Simulation system not yet built");
388 }
389
390 int *vip;
391 if(slv_near_bounds(sys,epsilon,&vip)){
392 struct var_variable **vp = slv_get_solvers_var_list(sys);
393 struct var_variable *var;
394 cerr << "VARS FOUND NEAR BOUNDS" << endl;
395 int nlow = vip[0];
396 int nhigh = vip[1];
397 int lim1 = 2 + nlow;
398 for(int i=2; i<lim1; ++i){
399 var = vp[vip[i]];
400 char *var_name = var_make_name(sys,var);
401 cerr << "AT LOWER BOUND: " << var_name << endl;
402 ascfree(var_name);
403 vars.push_back(Variable(this,var));
404 };
405 int lim2 = lim1 + nhigh;
406 for(int i=lim1; i<lim2; ++i){
407 var = vp[vip[i]];
408 char *var_name = var_make_name(sys,var);
409 cerr << "AT UPPER BOUND: " << var_name << endl;
410 ascfree(var_name);
411 vars.push_back(Variable(this,var));
412 }
413 }
414 ascfree(vip);
415 return vars;
416 }
417
418 void
419 Simulation::solve(Solver solver, SolverReporter &reporter){
420 if(!is_built){
421 throw runtime_error("Simulation::solver: simulation is not yet built, can't start solving.");
422 }
423
424 cerr << "SIMULATION::SOLVE STARTING..." << endl;
425 enum inst_t k = getModel().getKind();
426 if(k!=MODEL_INST)throw runtime_error("Can't solve: not an instance of type MODEL_INST");
427
428 Instance *i1 = getInternalType();
429 int npend = NumberPendingInstances(&*i1);
430 if(npend)throw runtime_error("Can't solve: There are still %d pending instances");
431
432 if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
433
434 cerr << "SIMULATION::SOLVE: SET SOLVER..." << endl;
435 setSolver(solver);
436
437
438 cerr << "PRESOLVING SYSTEM...";
439 slv_presolve(sys);
440 cerr << "DONE" << endl;
441
442 cerr << "SOLVING SYSTEM..." << endl;
443 // Add some stuff here for cleverer iteration....
444 unsigned niter = 1000;
445 //double updateinterval = 0.02;
446
447 double starttime = tm_cpu_time();
448 //double lastupdate = starttime;
449 SolverStatus status;
450 //int solved_vars=0;
451 bool stop=false;
452
453 status.getSimulationStatus(*this);
454 reporter.report(&status);
455
456 for(unsigned iter = 1; iter <= niter && !stop; ++iter){
457
458 if(status.isReadyToSolve()){
459 slv_iterate(sys);
460 }
461
462 status.getSimulationStatus(*this);
463
464 if(reporter.report(&status)){
465 stop = true;
466 }
467 }
468
469 double elapsed = tm_cpu_time() - starttime;
470
471
472 activeblock = status.getCurrentBlockNum();
473
474 reporter.finalise(&status);
475
476 // Just a little bit of console output:
477
478 if(status.isOK()){
479 cerr << "... SOLVED, STATUS OK" << endl;
480 }else{
481 cerr << "... SOLVER FAILED" << endl;
482 }
483
484 cerr << "SOLVER PERFORMED " << status.getIterationNum() << " ITERATIONS IN " << elapsed << "s" << endl;
485 }
486
487 void
488 Simulation::write(){
489 simroot.write();
490 }
491
492 //------------------------------------------
493 // ASSIGNING SOLVER TO SIMULATION
494
495 void
496 Simulation::setSolver(Solver &solver){
497 cerr << "SETTING SOLVER ON SIMULATION TO " << solver.getName() << endl;
498
499 if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
500 // Update the solver object because sometimes an alternative solver can be returned, apparently.
501
502 int selected = slv_select_solver(sys, solver.getIndex());
503 //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
504
505 if(selected<0){
506 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
507 throw runtime_error("Failed to select solver");
508 }
509
510 if(selected!=solver.getIndex()){
511 solver = Solver(slv_solver_name(selected));
512 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
513 }
514
515 if( slv_eligible_solver(sys) <= 0){
516 ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
517 throw runtime_error("Inelegible solver");
518 }
519 }
520
521 const Solver
522 Simulation::getSolver() const{
523 int index = slv_get_selected_solver(sys);
524 //cerr << "Simulation::getSolver: index = " << index << endl;
525 if(index<0)throw runtime_error("No solver selected");
526
527 return Solver(slv_solver_name(index));
528 }
529
530
531 /**
532 Get solver parameters struct wrapped up as a SolverParameters class.
533 */
534 SolverParameters
535 Simulation::getSolverParameters() const{
536 if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
537
538 slv_parameters_t p;
539 slv_get_parameters(sys,&p);
540 return SolverParameters(p);
541 }
542
543 /**
544 Update the solver parameters by passing a new set back
545 */
546 void
547 Simulation::setSolverParameters(SolverParameters &P){
548 if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
549 slv_set_parameters(sys, &(P.getInternalType()));
550 }
551
552 slv_system_structure *
553 Simulation::getSystem(){
554 if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
555 return sys;
556 }
557
558 IncidenceMatrix
559 Simulation::getIncidenceMatrix(){
560 return IncidenceMatrix(*this);
561 }
562
563 const string
564 Simulation::getInstanceName(const Instanc &i) const{
565 char *n;
566 n = WriteInstanceNameString(i.getInternalType(),simroot.getInternalType());
567 string s(n);
568 ascfree(n);
569 return s;
570 }
571
572 const int
573 Simulation::getNumVars(){
574 return slv_get_num_solvers_vars(getSystem());
575 }
576
577 void
578 Simulation::processVarStatus(){
579
580 // this is a cheap function call:
581 const mtx_block_t *bb = slv_get_solvers_blocks(getSystem());
582
583 var_variable **vlist = slv_get_solvers_var_list(getSystem());
584 int nvars = slv_get_num_solvers_vars(getSystem());
585
586 slv_status_t status;
587 slv_get_status(getSystem(), &status);
588
589 if(status.block.number_of == 0){
590 cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
591 return;
592 }
593
594 int activeblock = status.block.current_block;
595 int low = bb->block[activeblock].col.low;
596 int high = bb->block[activeblock].col.high;
597 bool allsolved = status.converged;
598 for(int c=0; c < nvars; ++c){
599 var_variable *v = vlist[c];
600 Instanc i((Instance *)var_instance(v));
601 VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
602 if(i.isFixed()){
603 s = ASCXX_VAR_FIXED;
604 }else if(var_incident(v) && var_active(v)){
605 if(allsolved || c < low){
606 s = ASCXX_VAR_SOLVED;
607 }else if(c <= high){
608 s = ASCXX_VAR_ACTIVE;
609 }else{
610 s = ASCXX_VAR_UNSOLVED;
611 }
612 }
613 i.setVarStatus(s);
614 }
615 }
616
617 const int
618 Simulation::getActiveBlock() const{
619 return activeblock;
620 }
621
622 bool
623 SingularityInfo::isSingular() const{
624 if(vars.size()||rels.size()){
625 return true;
626 }
627 return false;
628 }

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