Parent Directory | Revision Log

Revision **909** -
(**show annotations**)
(**download**)
(**as text**)

*Thu Oct 26 12:44:41 2006 UTC*
(17 years, 9 months ago)
by *johnpye*

File MIME type: text/x-chdr

File size: 16460 byte(s)

File MIME type: text/x-chdr

File size: 16460 byte(s)

Added finite-difference evaluation of gradients in blackboxes. Some work on J*v evaluation with IDA.

1 | /* ASCEND modelling environment |

2 | Copyright (C) 2006 Carnegie Mellon University |

3 | Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan |

4 | Copyright (C) 1993 Joseph Zaher |

5 | Copyright (C) 1990 Karl Michael Westerberg |

6 | |

7 | This program is free software; you can redistribute it and/or modify |

8 | it under the terms of the GNU General Public License as published by |

9 | the Free Software Foundation; either version 2, or (at your option) |

10 | any later version. |

11 | |

12 | This program is distributed in the hope that it will be useful, |

13 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |

15 | GNU General Public License for more details. |

16 | |

17 | You should have received a copy of the GNU General Public License |

18 | along with this program; if not, write to the Free Software |

19 | Foundation, Inc., 59 Temple Place - Suite 330, |

20 | Boston, MA 02111-1307, USA. |

21 | *//** @file |

22 | Relation manipulator module for the SLV solver. |

23 | |

24 | This module will provide supplemental operations for |

25 | relations such as simplification, evaluation, and |

26 | differentiation. |

27 | |

28 | Dates: 06/90 - original version |

29 | 06/93 - separated out exprman |

30 | 11/93 - added relman_calc_satisfied |

31 | |

32 | Requires: |

33 | #include "utilities/ascConfig.h" |

34 | #include "mtx.h" |

35 | #include "var.h" |

36 | #include "rel.h" |

37 | *//* |

38 | * by Karl Michael Westerberg and Joseph Zaher |

39 | * Created: 2/6/90 |

40 | * Version: $Revision: 1.29 $ |

41 | * Version control file: $RCSfile: relman.h,v $ |

42 | * Date last modified: $Date: 1998/04/23 23:56:24 $ |

43 | * Last modified by: $Author: ballan $ |

44 | */ |

45 | |

46 | #ifndef ASC_RELMAN_H |

47 | #define ASC_RELMAN_H |

48 | |

49 | #define relman_is_linear(a,b) (FALSE) |

50 | /**< |

51 | * Determines whether or not the given relation is linear in |

52 | * all of the variables which pass through the variable filter, treating |

53 | * those variables which fail to pass as constants.<br><br> |

54 | * |

55 | * Example: |

56 | * x1 + x2 >= 3 is a linear relation. |

57 | */ |

58 | |

59 | extern real64 relman_linear_coef(struct rel_relation *rel, |

60 | struct var_variable *var, |

61 | var_filter_t *filter); |

62 | /**< |

63 | * Computes the coefficient of the given variable in a linear |

64 | * relation. If var=NULL, then the "RHS" is returned instead. More |

65 | * nprecisely, "a[var]" is returned where:<br><br> |

66 | * |

67 | * rel : sum[variables v](a[v] * v) COMPARISON a[NULL]<br><br> |

68 | * |

69 | * It is assumed that the relation is linear, otherwise, |

70 | * something will be returned, but it won't be very meaningful. |

71 | */ |

72 | |

73 | extern void relman_simplify(struct rel_relation *rel, int opt_level); |

74 | /**< NOT IMPLEMENTED |

75 | * The left and right hand sides of the given relation are simplified |

76 | * to the extent given by opt_level. The effect of varying values of |

77 | * opt_level is given in the description for exprman_simplify(). |

78 | */ |

79 | |

80 | extern void relman_dissolve_vars(struct rel_relation *rel, var_filter_t *filter); |

81 | /**< NOT IMPLEMENTED |

82 | * Variables which pass through the filter are replaced in the |

83 | * relation by their current values. |

84 | */ |

85 | |

86 | extern void relman_decide_incidence(struct rel_relation *rel); |

87 | /**< |

88 | * Sets the incidence field of each variable which is found to be |

89 | * incident in the relation rel to TRUE. If these variables make |

90 | * up a subset of some larger variable list, it is important to first |

91 | * set the incidence field of all of the variables to FALSE before |

92 | * using this function in order to determine the unattached variables. |

93 | */ |

94 | |

95 | extern void relman_get_incidence(struct rel_relation *rel, |

96 | var_filter_t *filter, |

97 | mtx_matrix_t matrix); |

98 | /**< |

99 | * Upon return, coefficient (rel_n,var_n) (using original row and column |

100 | * numbers) is non-zero if and only if the relation rel with index rel_n |

101 | * depends on a variable with index var_n. |

102 | */ |

103 | |

104 | ASC_DLLSPEC(real64 ) relman_eval(struct rel_relation *rel, int32 *calc_okp, int safe); |

105 | /**< |

106 | The residual of the relation is calculated and returned. In addition |

107 | to returning the residual, the residual field of the relation is |

108 | updated. The residual field of the relation is not updated when an error |

109 | occurs. |

110 | |

111 | @param safe If set nonzero, "safe" functions are used for evaluation (means |

112 | that overflows, divide by zero, etc are avoided) |

113 | @param calc_ok (returned) status of the calculation. 0=error, else ok. |

114 | @return residual (= LHS - RHS, regardless of comparison) |

115 | |

116 | @NOTE |

117 | This function should be surrounded by Asc_SignalHandlerPush/Pop both |

118 | with arguments (SIGFPE,SIG_IGN). If it is being called in a loop, |

119 | the push/pop should be _outside_ the loop. |

120 | */ |

121 | |

122 | extern int32 relman_obj_direction(struct rel_relation *rel); |

123 | /**< |

124 | * Returns: |

125 | * - direction = -1 if objective is minimization |

126 | * - direction = 1 if objective is maximization |

127 | * - direction = 0 otherwise. (ie. if not an objective) |

128 | */ |

129 | |

130 | extern real64 relman_scale(struct rel_relation *rel); |

131 | /**< |

132 | * Calculates relation nominal scaling factor for |

133 | * current values stored in the relations variables. |

134 | * Fills the relations nominal field and also returns |

135 | * the relations nominal. |

136 | */ |

137 | |

138 | #define relman_diff(a,b,c,d) (abort(),1) |

139 | /**< |

140 | Calculates the derivative of the relation residual with respect to |

141 | the specified variable and stuffs it in pd. if problem with |

142 | calculation, returns 1, else 0. |

143 | If the value of safe is nonzero, "safe" functions will be used to |

144 | calculate the residual. |

145 | |

146 | @TODO relman_diff() needs to be reimplemented - needs compiler-side work. |

147 | */ |

148 | |

149 | extern int relman_diff2(struct rel_relation *rel, |

150 | var_filter_t *filter, |

151 | real64 *derivatives, |

152 | int32 *variables, |

153 | int32 *count, |

154 | int32 safe); |

155 | /**< |

156 | Calculates the row of the jacobian matrix (the transpose gradient of |

157 | the relation residual, $ \grad^{T}(f) $) corresponding to the relation |

158 | rel. The filter determines which variables actually contribute to the |

159 | jacobian. |

160 | |

161 | derivatives[i] will contain the derivative of the relation with |

162 | respect to the variable whose solver index is stored in |

163 | variables[i]. |

164 | |

165 | @param rel Relation being differentiated |

166 | @param filter Filter for variables for which derivs are desired |

167 | @param safe If nonzero, "safe" functions are used to for the calculations |

168 | @param derivatives output vector (allocated by the calling function) |

169 | @param variables output vector (allocated by the calling function) |

170 | @param count output value, will be set to the number of elements assigned upon exit. |

171 | |

172 | @return 0 on success, non-zero if an error is encountered in the calculation |

173 | */ |

174 | |

175 | extern int relman_diff_grad(struct rel_relation *rel, |

176 | var_filter_t *filter, |

177 | real64 *derivatives, |

178 | int32 *variables_master, |

179 | int32 *variables_solver, |

180 | int32 *count, |

181 | real64 *resid, |

182 | int32 safe); |

183 | /**< |

184 | * Calculates the row of the jacobian matrix (the transpose gradient of |

185 | * the relation residual grad^T(f) ) corresponding to the relation |

186 | * rel. The filter determines which variables actually contribute to the |

187 | * jacobian. The residual of the relation is also computed and returned. |

188 | * If an error is encountered in the calculation, the status returned is |

189 | * 1. Status = 0 is OK. |

190 | * If the value of safe is nonzero, "safe" functions are used to for |

191 | * the calculations.<br><br> |

192 | * The calling function should allocate the output vectors 'derivatives', |

193 | * 'variables_master' and 'variables_solver'. 'count' will be set to |

194 | * the number of elements assigned upon exit. |

195 | * derivative(i) will contain the derivative of the relation with |

196 | * respect to the variable whose master index is stored in |

197 | * variables_master(i). The solver index of each variable is stored in |

198 | * variables_solver(i). |

199 | * |

200 | * There are two differences wrt to relman_diff2: |

201 | * - the master index (solver independent) is obtained |

202 | * - the residual is evaluated |

203 | */ |

204 | |

205 | ASC_DLLSPEC(int ) relman_diffs(struct rel_relation *rel, |

206 | var_filter_t *filter, mtx_matrix_t mtx, |

207 | real64 *resid, int safe); |

208 | /**< |

209 | Calculates the row of the jacobian matrix (the transpose gradient of |

210 | the relation residual grad^T(f) ) corresponding to the relation |

211 | rel. The filter determines which variables actually contribute to the |

212 | jacobian. The residual of the relation is also computed and returned. |

213 | If an error is encountered in the calculation, the status returned is |

214 | 1 and the residual is set to some number we managed to calculate, |

215 | while the gradient is discarded. status = 0 is OK. |

216 | |

217 | @param rel relation for which jacobian entries are required |

218 | @param filter filter for which variables should actually contribute to the jacobian |

219 | @param mtx matrix into which the row (corresponding to rel) is written |

220 | @param safe if non-zero, "safe" functions are used to for the calucaltions. |

221 | |

222 | It doesn't matter how you have permuted the columns and rows: |

223 | for the vars which pass the filter you send we |

224 | fill the org row determined by rel_sindex and the org cols |

225 | determined by var_sindex. |

226 | |

227 | @return 0 on success, 1 on calculation error (residual will be returned, grad discarded) |

228 | |

229 | @NOTE The row of the mtx corresponding to rel should be cleared |

230 | before calling this function, since this FILLS with the gradient.<br><br> |

231 | |

232 | @NOTE *changed* -- This operator used to just ADD on top of any incidence |

233 | already in the row. This is not TRUE now. |

234 | |

235 | @TODO This operator really needs to be redesigned so it can deal with |

236 | harwellian matrices, glassbox rels and blackbox. |

237 | */ |

238 | |

239 | extern int32 relman_diff_harwell(struct rel_relation **rlist, |

240 | var_filter_t *vfilter, rel_filter_t *rfilter, |

241 | int32 rlen, int32 bias, int32 mors, |

242 | real64 *avec, int32 *ivec, int32 *jvec); |

243 | /**< |

244 | * This fills an "a-i-j" sparse matrix in the avec/ivec/jvec given. |

245 | * @param rlist struct rel_relation **, list of relations rlen long. |

246 | * @param vfilter var_filter_t *, stuffs gradient for matching variables only. |

247 | * @param rlen int32, length of list of relations. |

248 | * @param bias int32, 0 = row grouped together, 1 = column grouped together. |

249 | * There is a substantial penalty for bias = 1. we MODEL by row. |

250 | * @param mors int32, 0 = master var index of columns, master rel index of rows |

251 | * 1 = solver var index of columns, master rel index of rows |

252 | * 2 = master var index of columns, solver rel index of rows |

253 | * 3 = solver var index of columns, solver rel index of rows |

254 | * Size of avec,ivec,jvec given is assumed big enough. |

255 | * big_enough = relman_jacobian_count(rlist,rlen,vfilter,rfilter,&dummy); |

256 | * If ivec or jvec given is NULL, then neither is stuffed, though avec is. |

257 | @return 0 on success, <0 on floating point errors, 1 on unrecoverable error. |

258 | |

259 | err = 1 --> unrecoverable error/bad input. caller should probably punt. |

260 | err < 0 --> -(number of floating point errors in evaluation). |

261 | The matrix will contain an approximation only. |

262 | |

263 | @todo relman_diff_harwell() bias == 1 is not yet implemented. |

264 | */ |

265 | |

266 | extern int32 relman_jacobian_count(struct rel_relation **rlist, |

267 | int32 rlen, |

268 | var_filter_t *vfilter, |

269 | rel_filter_t *rfilter, |

270 | int32 *rhomax); |

271 | /**< |

272 | * Return the number of nonzero gradient entries in the equations |

273 | * given. Only equations passing rfilter and entries passing vfilter |

274 | * are counted. rlen is the length of the relation list. |

275 | * *rhomax is the largest row count on return. |

276 | */ |

277 | |

278 | extern boolean relman_calc_satisfied_scaled(struct rel_relation *rel, |

279 | real64 tolerance); |

280 | /**< |

281 | * This definition of satisfaction includes the notion |

282 | * of scaling by the relation nominal before comparison. |

283 | */ |

284 | extern boolean relman_calc_satisfied(struct rel_relation *rel, |

285 | real64 tolerance); |

286 | /**< |

287 | * Returns TRUE or FALSE depending on whether the relation whose residual |

288 | * has been previously calculated is satisfied based on the value stored |

289 | * in the residual field. The satisfied field of the relation is also |

290 | * updated. A tolerance specification allows equalities to be declared |

291 | * satisfied as long as their residuals are close to zero. |

292 | * |

293 | * <!-- relman_calc_satisfied_scaled: --> |

294 | * <!-- This definition of satisfaction includes the notion --> |

295 | * <!-- of scaling by the relation nominal before comparison. --> |

296 | */ |

297 | |

298 | #define relman_directly_solve(r,s,a,n) \ |

299 | relman_directly_solve_new(r,s,a,n,1.0e-8) |

300 | /**< @see relman_directly_solve_new(). */ |

301 | extern real64 *relman_directly_solve_new(struct rel_relation *rel, |

302 | struct var_variable *solvefor, |

303 | int *able, |

304 | int *nsolns, |

305 | real64 tol); |

306 | /**< |

307 | * Attempts to solve the given equation for the given variable. If this |

308 | * function is able to determine the solution set, then *able is set to |

309 | * 1 and a newly allocated solution list is returned: *nsolns will be |

310 | * set to the length of this array. Otherwise *able is 0 and NULL |

311 | * is returned. NULL *may* also be returned if the solution set is empty. |

312 | * A return of able == 1, solution_list != NULL, and nsolns == 0 is |

313 | * possible for certain classes of floating point exceptions. |

314 | * It is assumed that the relation is a condition of equality.<br><br> |

315 | * |

316 | * relman_directly_solve_new() handles passing in a tolerance for glassbox |

317 | * relations so a rootfinder can do the work rather than leaving it to |

318 | * someone else. The rootfinder is based on Brent's algorithm. Old clients of |

319 | * relman_directly solve are grandfathered at a default tolerance of 1e-8. |

320 | * Once all these clients are gone, go back to the old name space for |

321 | * the new function.<br><br> |

322 | * |

323 | * Do NOT free or keep a persistent pointer to the solution_list this |

324 | * function returns. |

325 | */ |

326 | |

327 | #define relman_make_string_infix(sys,rel) \ |

328 | relman_make_vstring_infix((sys),(rel),TRUE) |

329 | /**< @see relman_make_xstring_infix() */ |

330 | #define relman_make_string_postfix(sys,rel) \ |

331 | relman_make_vstring_postfix((sys),(rel),TRUE) |

332 | /**< @see relman_make_xstring_infix() */ |

333 | #define relman_make_xstring_infix(sys,rel) \ |

334 | relman_make_vstring_infix((sys),(rel),FALSE) |

335 | /**< |

336 | @return string |

337 | @param sys |

338 | @param rel struct rel_relation *rel |

339 | |

340 | Creates a sufficiently large string (you must free it when you're |

341 | done with it) into which it converts a relation. The string will be |

342 | terminated with a '\0' character. |

343 | |

344 | The xstring versions of this call make strings where all the variables |

345 | are written as x<varindex> (e.g. x23) rather than as object names. |

346 | The MASTER index (var_mindex) is used, not the solver's sindex. |

347 | Currently the compiler does not support xstring postfix format, |

348 | but could easily do so if needed. |

349 | |

350 | The name of a var is context dependent, so you have to provide the |

351 | slv_system_t from which you got the relation. |

352 | */ |

353 | #if 0 /* needs compiler side work */ |

354 | #define relman_make_xstring_postfix(sys,rel) \ |

355 | relman_make_vstring_postfix((sys),(rel),FALSE) |

356 | #else |

357 | #define relman_make_xstring_postfix(sys,rel) \ |

358 | dummy_rel_string((sys),(rel),0) |

359 | #endif |

360 | /**< |

361 | * Not suppported. |

362 | * @TODO Consider adding support for xstring postfix format. |

363 | */ |

364 | |

365 | ASC_DLLSPEC(char *) relman_make_vstring_infix(slv_system_t sys, |

366 | struct rel_relation *rel, |

367 | int style); |

368 | /**< |

369 | * Implementation function for relman_make_string_infix() and |

370 | * relman_make_xstring_infix(). Do not call this function |

371 | * directly - use one of the macros instead. |

372 | */ |

373 | extern char *relman_make_vstring_postfix(slv_system_t sys, |

374 | struct rel_relation *rel, |

375 | int style); |

376 | /**< |

377 | * Implementation function for relman_make_string_postfix(). and |

378 | * Do not call this function directly - use the macro instead. |

379 | */ |

380 | |

381 | extern char *dummyrelstring(slv_system_t sys, |

382 | struct rel_relation *rel, |

383 | int style); |

384 | /**< Temporary no-op function to placehold unimplemented io functions. */ |

385 | |

386 | extern void relman_free_reused_mem(void); |

387 | /**< |

388 | * Call when desired to free memory cached internally. |

389 | */ |

390 | |

391 | #endif /* ASC_RELMAN_H */ |

392 |

john.pye@anu.edu.au | ViewVC Help |

Powered by ViewVC 1.1.22 |