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

Annotation of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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