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

Annotation of /trunk/pygtk/simulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

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