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