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

Annotation of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 912 - (hide annotations) (download) (as text)
Fri Oct 27 07:18:21 2006 UTC (13 years, 3 months ago) by johnpye
File MIME type: text/x-c++src
File size: 22326 byte(s)
Removed BBOXWHINE (replaced with some one-time-only warnings for the moment)
Added ExtMethodDestroyFn to allow 'user_data' associated with external methods to be destroyed.
Implemented the destroy fn through to 'extpy' module.
Added 'name' as an extra parameter in the user_data for extpy, to help with debug msgs.
Moved 'solvernotes' to a file of its own (was part of listnotes.py)
Added 'repaint' to GTK 'tools' menu (for debugging)
Added 'python.h' to top of library, type files (pygtk) to stop silly warnings.
Working on some diagnosing of problems as noted in Simulation::checkInstance.
Removed some old comments from namio.h and others.
Renamed 'blsys' to 'sys' in integrator.c.
Some work on fixing up the J*v function for IDA (not yet complete).
Added new 'destroyfn' parameter (as NULL) to all calls to 'CreateUserFunctionMethod'.
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     #include <utilities/ascConfig.h>
29 johnpye 506 #include <utilities/error.h>
30 johnpye 132 #include <utilities/ascSignal.h>
31     #include <utilities/ascMalloc.h>
32     #include <general/dstring.h>
33 johnpye 196 #include <general/tm_time.h>
34 johnpye 132 #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 johnpye 900 #include <compiler/importhandler.h>
48 johnpye 132
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 johnpye 328 #include <solver/slv_server.h>
70 johnpye 132 }
71    
72     #include "simulation.h"
73     #include "solver.h"
74 johnpye 208 #include "solverparameters.h"
75 johnpye 132 #include "name.h"
76 johnpye 233 #include "incidencematrix.h"
77 johnpye 237 #include "variable.h"
78 johnpye 307 #include "solverstatus.h"
79 johnpye 310 #include "solverreporter.h"
80 johnpye 132
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 johnpye 153 is_built = false;
88 johnpye 132 // Create an Instance object for the 'simulation root' (we'll call
89 johnpye 190 // it the 'simulation model') and it can be fetched using 'getModel()'
90 johnpye 132 // any time later.
91     //simroot = Instanc(GetSimulationRoot(i),name);
92     }
93    
94 johnpye 164 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 johnpye 190 bin_rm = old.bin_rm;
102 johnpye 775 sing = NULL;
103 johnpye 164 }
104    
105     Simulation::~Simulation(){
106     //CONSOLE_DEBUG("Deleting simulation %s", getName().toString());
107     }
108    
109 johnpye 132 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 johnpye 775
118     slv_system_structure *
119     Simulation::getSystem(){
120     if(!sys)throw runtime_error("Can't getSystem: simulation not yet built");
121     return sys;
122 johnpye 132 }
123    
124 johnpye 295
125 johnpye 775 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 johnpye 295
134 johnpye 775 const int
135     Simulation::getNumVars(){
136     return slv_get_num_solvers_vars(getSystem());
137 johnpye 295 }
138    
139    
140 johnpye 775 void
141     Simulation::write(){
142     simroot.write();
143 johnpye 295 }
144    
145 johnpye 775 //------------------------------------------------------------------------------
146     // RUNNING MODEL 'METHODS'
147 johnpye 772
148 johnpye 295 void
149 johnpye 132 Simulation::run(const Method &method){
150 johnpye 774 Instanc &model = getModel();
151     this->run(method,model);
152     }
153    
154     void
155 johnpye 900 Simulation::runDefaultMethod(){
156     Method m = getType().getMethod(SymChar("on_load"));
157     run(m);
158     }
159    
160     void
161 johnpye 774 Simulation::run(const Method &method, Instanc &model){
162    
163 johnpye 900 // set the 'sim' pointer to our local variable...
164     CONSOLE_DEBUG("Setting shared pointer 'sim'");
165     importhandler_setsharedpointer("sim",this);
166    
167 johnpye 912 if(not is_built){
168     CONSOLE_DEBUG("WARNING, SIMULATION NOT YET BUILT");
169     }
170    
171     CONSOLE_DEBUG("Running method %s...", method.getName());
172    
173 johnpye 132 Nam name = Nam(method.getSym());
174 johnpye 154 //cerr << "CREATED NAME '" << name.getName() << "'" << endl;
175 johnpye 900
176    
177 johnpye 132 Proc_enum pe;
178     pe = Initialize(
179 johnpye 774 &*(model.getInternalType()) ,name.getInternalType(), "__not_named__"
180 johnpye 132 ,ASCERR
181     ,0, NULL, NULL
182     );
183    
184 johnpye 900 // clear out the 'sim' pointer (soon it will be invalid)
185     importhandler_setsharedpointer("sim",NULL);
186     CONSOLE_DEBUG("Cleared shared pointer 'sim'");
187    
188 johnpye 132 if(pe == Proc_all_ok){
189 johnpye 190 ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Method '%s' was run (check above for errors)\n",method.getName());
190 johnpye 132 //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 johnpye 775 //-----------------------------------------------------------------------------
265     // CHECKING METHODS
266    
267 johnpye 772 /**
268 johnpye 775 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 johnpye 912
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 johnpye 772 */
276 johnpye 775 void
277     Simulation::checkInstance(){
278 johnpye 897 //cerr << "CHECKING SIMULATION INSTANCE" << endl;
279 johnpye 912 if(!is_built){
280     ERROR_REPORTER_HERE(ASC_PROG_ERR,"Simulation has not been built");
281     return;
282     }
283 johnpye 132 Instance *i1 = getModel().getInternalType();
284 johnpye 164 CheckInstance(stderr, &*i1);
285 johnpye 897 //cerr << "DONE CHECKING INSTANCE" << endl;
286 johnpye 775 }
287 johnpye 772
288 johnpye 775 /**
289     @return 1 = underspecified, 2 = square, 3 = structurally singular, 4 = overspecified
290     */
291     enum StructuralStatus
292     Simulation::checkDoF() const{
293     cerr << "CHECKING DOF..." << endl;
294     int dof, status;
295     if(!sys){
296     throw runtime_error("System not yet built");
297     }
298     slvDOF_status(sys, &status, &dof);
299     switch(status){
300     case ASCXX_DOF_UNDERSPECIFIED:
301     case ASCXX_DOF_SQUARE:
302     case ASCXX_DOF_OVERSPECIFIED:
303     case ASCXX_DOF_STRUCT_SINGULAR:
304     return (enum StructuralStatus)status;
305     case 5:
306     throw runtime_error("Unable to resolve degrees of freedom"); break;
307     default:
308     throw runtime_error("Invalid return status from slvDOF_status");
309     }
310 johnpye 190 }
311 johnpye 132
312 johnpye 775 /**
313     Check consistency
314    
315     @TODO what is the difference between this and checkStructuralSingularity?
316    
317     @return list of freeable variables. List will be empty if sys is consistent.
318     */
319     vector<Variable>
320     Simulation::getFreeableVariables(){
321     vector<Variable> v;
322    
323 johnpye 897 //cerr << "CHECKING CONSISTENCY..." << endl;
324 johnpye 841 int *fixedarrayptr=NULL;
325 johnpye 775
326 johnpye 841 if(!sys){
327     throw runtime_error("System not yet built");
328     }
329    
330 johnpye 775 int res = consistency_analysis(sys, &fixedarrayptr);
331    
332     if(res==1){
333     cerr << "STRUCTURALLY CONSISTENT" << endl;
334     }else{
335 johnpye 841 if(fixedarrayptr ==NULL){
336     ERROR_REPORTER_HERE(ASC_USER_ERROR,"STRUCTURALLY INCONSISTENT");
337     throw runtime_error("Invalid constistency analysis result returned!");
338     }
339    
340     struct var_variable **vp = slv_get_master_var_list(sys);
341 johnpye 775 for(int i=0; fixedarrayptr[i]!=-1; ++i){
342     v.push_back( Variable(this, vp[fixedarrayptr[i]]) );
343     }
344     }
345     return v;
346     }
347    
348     /** Returns TRUE if all is OK (not singular) */
349     bool
350     Simulation::checkStructuralSingularity(){
351     int *vil;
352     int *ril;
353     int *fil;
354    
355     cerr << "RETREIVING slfDOF_structsing INFO" << endl;
356    
357     int res = slvDOF_structsing(sys, mtx_FIRST, &vil, &ril, &fil);
358     struct var_variable **varlist = slv_get_solvers_var_list(sys);
359     struct rel_relation **rellist = slv_get_solvers_rel_list(sys);
360    
361     if(this->sing){
362     cerr << "DELETING OLD SINGULATING INFO" << endl;
363     delete this->sing;
364     this->sing = NULL;
365     }
366    
367     if(res==1){
368     CONSOLE_DEBUG("processing singularity data...");
369     sing = new SingularityInfo();
370    
371     // pull in the lists of vars and rels, and the freeable vars:
372     for(int i=0; ril[i]!=-1; ++i){
373     sing->rels.push_back( Relation(this, rellist[ril[i]]) );
374     }
375    
376     for(int i=0; vil[i]!=-1; ++i){
377     sing->vars.push_back( Variable(this, varlist[vil[i]]) );
378     }
379    
380     for(int i=0; fil[i]!=-1; ++i){
381     sing->freeablevars.push_back( Variable(this, varlist[fil[i]]) );
382     }
383    
384     // we're done with those lists now
385     ASC_FREE(vil);
386     ASC_FREE(ril);
387     ASC_FREE(fil);
388    
389     if(sing->isSingular()){
390     CONSOLE_DEBUG("singularity found");
391     this->sing = sing;
392     return FALSE;
393     }
394     CONSOLE_DEBUG("no singularity");
395     delete sing;
396     return TRUE;
397     }else{
398     if(res==0){
399     throw runtime_error("Unable to determine singularity lists");
400     }else{
401     throw runtime_error("Invalid return from slvDOF_structsing.");
402     }
403     }
404     }
405    
406     /**
407     If the checkStructuralSingularity analysis has been done,
408     this funciton will let you access the SingularityInfo data that was
409     stored.
410     */
411     const SingularityInfo &
412     Simulation::getSingularityInfo() const{
413     if(sing==NULL){
414     throw runtime_error("No singularity info present");
415     }
416     return *sing;
417     }
418    
419     //------------------------------------------
420     // ASSIGNING SOLVER TO SIMULATION
421    
422 johnpye 132 void
423 johnpye 775 Simulation::setSolver(Solver &solver){
424 johnpye 910 /* CONSOLE_DEBUG("Setting solver on sim %p, root inst %p",this,this->simroot.getInternalType()); */
425 johnpye 900 if(!is_built)throw runtime_error("Can't setSolver: simulation has not been built yet");
426     if(!sys)throw runtime_error("Can't setSolver: 'sys' is not assigned.");
427 johnpye 775
428     // Update the solver object because sometimes an alternative solver can be returned, apparently.
429    
430     int selected = slv_select_solver(sys, solver.getIndex());
431     //cerr << "Simulation::setSolver: slv_select_solver returned " << selected << endl;
432    
433     if(selected<0){
434     ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Failed to select solver");
435     throw runtime_error("Failed to select solver");
436     }
437    
438     if(selected!=solver.getIndex()){
439     solver = Solver(slv_solver_name(selected));
440     ERROR_REPORTER_NOLINE(ASC_PROG_NOTE,"Substitute solver '%s' (index %d) selected.\n", solver.getName().c_str(), selected);
441     }
442    
443     if( slv_eligible_solver(sys) <= 0){
444     ERROR_REPORTER_NOLINE(ASC_PROG_ERROR,"Inelegible solver '%s'", solver.getName().c_str() );
445     throw runtime_error("Inelegible solver");
446     }
447     }
448    
449     const Solver
450     Simulation::getSolver() const{
451     int index = slv_get_selected_solver(sys);
452     //cerr << "Simulation::getSolver: index = " << index << endl;
453     if(index<0)throw runtime_error("No solver selected");
454    
455     return Solver(slv_solver_name(index));
456     }
457    
458     //------------------------------------------------------------------------------
459     // BUILD THE SYSTEM (SEND IT TO THE SOLVER)
460    
461     /**
462     Build the system (send it to the solver)
463     */
464     void
465 johnpye 132 Simulation::build(){
466 johnpye 897 //cerr << "BUILDING SIMULATION..." << endl;
467 johnpye 900 if(is_built){
468     throw runtime_error("Simulation is already built");
469     }
470    
471 johnpye 164 Instance *i1 = getModel().getInternalType();
472     sys = system_build(&*i1);
473     if(!sys){
474 johnpye 132 throw runtime_error("Unable to build system");
475     }
476 johnpye 153 is_built = true;
477 johnpye 897 //cerr << "...DONE BUILDING" << endl;
478 johnpye 132 }
479    
480 johnpye 775
481     //------------------------------------------------------------------------------
482     // SOLVER CONFIGURATION PARAMETERS
483    
484     /**
485     Get solver parameters struct wrapped up as a SolverParameters class.
486     */
487     SolverParameters
488     Simulation::getSolverParameters() const{
489 johnpye 900 if(!is_built)throw runtime_error("Can't getSolverParameters: Simulation system has not been built yet.");
490     if(!sys)throw runtime_error("Can't getSolverParameters: Simulation system has no 'sys' assigned.");
491 johnpye 775
492     slv_parameters_t p;
493     slv_get_parameters(sys,&p);
494     return SolverParameters(p);
495     }
496    
497     /**
498     Update the solver parameters by passing a new set back
499     */
500     void
501     Simulation::setSolverParameters(SolverParameters &P){
502     if(!sys)throw runtime_error("Can't set solver parameters: simulation has not been built yet.");
503     slv_set_parameters(sys, &(P.getInternalType()));
504     }
505    
506     //------------------------------------------------------------------------------
507     // PRE-SOLVE DIAGNOSTICS
508    
509     /**
510     Get a list of variables to fix to make an underspecified system
511     become square. Also seems to return stuff when you have a structurally
512     singuler system.
513     */
514 johnpye 132 vector<Variable>
515     Simulation::getFixableVariables(){
516 johnpye 897 //cerr << "GETTING FIXABLE VARIABLES..." << endl;
517 johnpye 132 vector<Variable> vars;
518    
519 johnpye 164 if(!sys){
520 johnpye 132 throw runtime_error("Simulation system not yet built");
521     }
522    
523 johnpye 316 int32 *vip; /** TODO ensure 32 bit integers are used */
524 johnpye 132
525     // Get IDs of elegible variables in array at vip...
526 johnpye 316 if(!slvDOF_eligible(sys,&vip)){
527 johnpye 190 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"No fixable variables found.");
528 johnpye 132 }else{
529     struct var_variable **vp = slv_get_solvers_var_list(sys);
530    
531     if(vp==NULL){
532     throw runtime_error("Simulation variable list is null");
533     }
534 johnpye 190
535 johnpye 132 // iterate through this list until we find a -1:
536     int i=0;
537 johnpye 316 int var_index = vip[i];
538 johnpye 132 while(var_index >= 0){
539     struct var_variable *var = vp[var_index];
540 johnpye 237 vars.push_back( Variable(this, var) );
541 johnpye 132 ++i;
542 johnpye 316 var_index = vip[i];
543 johnpye 132 }
544 johnpye 190 ERROR_REPORTER_NOLINE(ASC_USER_NOTE,"Found %d fixable variables.",i);
545 johnpye 316 ascfree(vip);
546 johnpye 132 }
547    
548     return vars;
549     }
550    
551 johnpye 775 /**
552     Get the list of variables near their bounds. Helps to indentify why
553     you might be having non-convergence problems.
554     */
555 johnpye 328 vector<Variable>
556     Simulation::getVariablesNearBounds(const double &epsilon){
557 johnpye 897 //cerr << "GETTING VARIABLES NEAR BOUNDS..." << endl;
558 johnpye 328 vector<Variable> vars;
559    
560     if(!sys){
561     throw runtime_error("Simulation system not yet built");
562     }
563    
564     int *vip;
565     if(slv_near_bounds(sys,epsilon,&vip)){
566     struct var_variable **vp = slv_get_solvers_var_list(sys);
567     struct var_variable *var;
568     cerr << "VARS FOUND NEAR BOUNDS" << endl;
569     int nlow = vip[0];
570     int nhigh = vip[1];
571     int lim1 = 2 + nlow;
572     for(int i=2; i<lim1; ++i){
573     var = vp[vip[i]];
574     char *var_name = var_make_name(sys,var);
575     cerr << "AT LOWER BOUND: " << var_name << endl;
576     ascfree(var_name);
577     vars.push_back(Variable(this,var));
578     };
579     int lim2 = lim1 + nhigh;
580     for(int i=lim1; i<lim2; ++i){
581     var = vp[vip[i]];
582     char *var_name = var_make_name(sys,var);
583     cerr << "AT UPPER BOUND: " << var_name << endl;
584     ascfree(var_name);
585     vars.push_back(Variable(this,var));
586     }
587     }
588     ascfree(vip);
589     return vars;
590     }
591 johnpye 502
592 johnpye 775
593     bool
594     SingularityInfo::isSingular() const{
595     if(vars.size()||rels.size()){
596     return true;
597     }
598     return false;
599     }
600    
601     //------------------------------------------------------------------------------
602     // SOLVING
603    
604     /**
605     Solve the system through to convergence. This function is hardwired with
606     a maximum of 1000 iterations, but will interrupt itself when the 'stop'
607     condition comes back from the SolverReporter.
608     */
609 johnpye 132 void
610 johnpye 310 Simulation::solve(Solver solver, SolverReporter &reporter){
611 johnpye 153 if(!is_built){
612     throw runtime_error("Simulation::solver: simulation is not yet built, can't start solving.");
613     }
614    
615 johnpye 890 //cerr << "SIMULATION::SOLVE STARTING..." << endl;
616 johnpye 132 enum inst_t k = getModel().getKind();
617     if(k!=MODEL_INST)throw runtime_error("Can't solve: not an instance of type MODEL_INST");
618 johnpye 190
619 johnpye 164 Instance *i1 = getInternalType();
620     int npend = NumberPendingInstances(&*i1);
621 johnpye 190 if(npend)throw runtime_error("Can't solve: There are still %d pending instances");
622 johnpye 132
623     if(!sys)throw runtime_error("Can't solve: Simulation system has not been built yet.");
624 johnpye 190
625 johnpye 890 //cerr << "SIMULATION::SOLVE: SET SOLVER..." << endl;
626 johnpye 132 setSolver(solver);
627 johnpye 190
628    
629 johnpye 890 //cerr << "PRESOLVING SYSTEM...";
630 johnpye 132 slv_presolve(sys);
631 johnpye 890 //cerr << "DONE" << endl;
632 johnpye 190
633 johnpye 890 //cerr << "SOLVING SYSTEM..." << endl;
634 johnpye 196 // Add some stuff here for cleverer iteration....
635 johnpye 207 unsigned niter = 1000;
636 johnpye 714 //double updateinterval = 0.02;
637 johnpye 132
638 johnpye 196 double starttime = tm_cpu_time();
639 johnpye 714 //double lastupdate = starttime;
640 johnpye 307 SolverStatus status;
641 johnpye 714 //int solved_vars=0;
642 johnpye 196 bool stop=false;
643    
644 johnpye 307 status.getSimulationStatus(*this);
645 johnpye 314 reporter.report(&status);
646 johnpye 307
647 johnpye 746 for(unsigned iter = 1; iter <= niter && !stop; ++iter){
648 johnpye 285
649 johnpye 307 if(status.isReadyToSolve()){
650 johnpye 196 slv_iterate(sys);
651     }
652 johnpye 307
653     status.getSimulationStatus(*this);
654 johnpye 502
655 johnpye 317 if(reporter.report(&status)){
656     stop = true;
657     }
658 johnpye 196 }
659 johnpye 317
660 johnpye 230 double elapsed = tm_cpu_time() - starttime;
661 johnpye 196
662 johnpye 319
663 johnpye 307 activeblock = status.getCurrentBlockNum();
664 johnpye 285
665 johnpye 319 reporter.finalise(&status);
666    
667 johnpye 322 // Just a little bit of console output:
668    
669 johnpye 307 if(status.isOK()){
670 johnpye 890 //cerr << "... SOLVED, STATUS OK" << endl;
671 johnpye 132 }else{
672 johnpye 321 cerr << "... SOLVER FAILED" << endl;
673 johnpye 132 }
674    
675 johnpye 890 //cerr << "SOLVER PERFORMED " << status.getIterationNum() << " ITERATIONS IN " << elapsed << "s" << endl;
676 johnpye 132 }
677    
678 johnpye 775 //------------------------------------------------------------------------------
679     // POST-SOLVE DIAGNOSTICS
680 johnpye 132
681 johnpye 775 const int
682     Simulation::getActiveBlock() const{
683     return activeblock;
684 johnpye 132 }
685    
686 johnpye 208 /**
687 johnpye 775 Return an IncidenceMatrix built from the current state of the solver system.
688 johnpye 132
689 johnpye 775 This will actually return something meaningful even before solve.
690 johnpye 225 */
691 johnpye 233 IncidenceMatrix
692     Simulation::getIncidenceMatrix(){
693     return IncidenceMatrix(*this);
694     }
695 johnpye 252
696 johnpye 775 /**
697     This function looks at all the variables in the solve's list and updates
698     the variable status for the corresponding instances.
699 johnpye 317
700 johnpye 775 It does this by using the 'interface pointer' in the Instance, see
701     the C-API function GetInterfacePtr.
702 johnpye 317
703 johnpye 775 This is used to display visually which variables have been solved, which
704     ones have not yet been attempted, and which ones were active when the solver
705     failed (ASCXX_VAR_ACTIVE).
706     */
707 johnpye 255 void
708     Simulation::processVarStatus(){
709    
710     // this is a cheap function call:
711     const mtx_block_t *bb = slv_get_solvers_blocks(getSystem());
712 johnpye 297
713 johnpye 255 var_variable **vlist = slv_get_solvers_var_list(getSystem());
714     int nvars = slv_get_num_solvers_vars(getSystem());
715    
716     slv_status_t status;
717     slv_get_status(getSystem(), &status);
718    
719 johnpye 297 if(status.block.number_of == 0){
720     cerr << "Variable statuses can't be set: block structure not yet determined." << endl;
721     return;
722     }
723    
724 johnpye 255 int activeblock = status.block.current_block;
725     int low = bb->block[activeblock].col.low;
726     int high = bb->block[activeblock].col.high;
727 johnpye 258 bool allsolved = status.converged;
728 johnpye 255 for(int c=0; c < nvars; ++c){
729     var_variable *v = vlist[c];
730     Instanc i((Instance *)var_instance(v));
731     VarStatus s = ASCXX_VAR_STATUS_UNKNOWN;
732 johnpye 258 if(i.isFixed()){
733     s = ASCXX_VAR_FIXED;
734     }else if(var_incident(v) && var_active(v)){
735     if(allsolved || c < low){
736 johnpye 255 s = ASCXX_VAR_SOLVED;
737     }else if(c <= high){
738     s = ASCXX_VAR_ACTIVE;
739     }else{
740     s = ASCXX_VAR_UNSOLVED;
741     }
742     }
743     i.setVarStatus(s);
744     }
745     }
746 johnpye 285

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