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

Annotation of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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