/[ascend]/trunk/base/generic/solver/rel.c
ViewVC logotype

Contents of /trunk/base/generic/solver/rel.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 771 - (show annotations) (download) (as text)
Fri Jul 14 04:31:54 2006 UTC (18 years ago) by johnpye
File MIME type: text/x-csrc
File size: 35854 byte(s)
Removed some debug noise.
Minor changes to dsg models.
1 /* ASCEND modelling environment
2 Copyright (C) 1990 Karl Michael Westerberg
3 Copyright (C) 1993 Joseph Zaher
4 Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
5 Copyright (C) 2006 Carnegie Mellon University
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 *//*
22 by Karl Michael Westerberg
23 Created: 2/6/90
24 Last in CVS $Revision: 1.32 $ $Date: 1998/01/29 00:42:28 $ $Author: ballan $
25 */
26
27 #include <math.h>
28 #include <stdarg.h>
29 #include <utilities/ascConfig.h>
30 #include <utilities/ascPanic.h>
31 #include <utilities/ascMalloc.h>
32 #include <utilities/mem.h>
33 #include <general/list.h>
34 #include <general/dstring.h>
35 #include <compiler/compiler.h>
36 #include <compiler/symtab.h>
37 #include <compiler/fractions.h>
38 #include <compiler/instance_enum.h>
39 #include <compiler/instance_io.h>
40 #include <compiler/symtab.h>
41 #include <compiler/extfunc.h>
42 #include <compiler/extcall.h>
43 #include <compiler/functype.h>
44 #include <compiler/safe.h>
45 #include <compiler/dimen.h>
46 #include <compiler/expr_types.h>
47 #include <compiler/find.h>
48 #include <compiler/atomvalue.h>
49 #include <compiler/instquery.h>
50 #include <compiler/mathinst.h>
51 #include <compiler/parentchild.h>
52 #include <compiler/instance_io.h>
53 #include <compiler/relation_type.h>
54 #include <compiler/relation.h>
55 #include <compiler/relation_util.h>
56 #include <compiler/packages.h>
57 #define _SLV_SERVER_C_SEEN_ /* for the extrel stuff in header */
58 #include "mtx.h"
59 #include "slv_types.h"
60 #include "var.h"
61 #include "rel.h"
62 #include "discrete.h"
63 #include "conditional.h"
64 #include "logrel.h"
65 #include "bnd.h"
66 #include "slv_server.h"
67
68 /*------------------------------------------------------------------------------
69 forward declarations, constants, typedefs
70 */
71
72 #ifdef DEBUG
73 # define REL_DEBUG CONSOLE_DEBUG
74 #else
75 # define REL_DEBUG(MSG,...) ((void)0)
76 #endif
77
78 #define IPTR(i) ((struct Instance *)(i))
79 #define REIMPLEMENT 0 /* if set to 1, compiles code tagged with it. */
80
81 /* define symchar names needed */
82 static symchar *g_strings[1];
83 #define INCLUDED_R g_strings[0]
84
85 static const struct rel_relation g_rel_defaults = {
86 NULL, /* instance */
87 NULL, /* extnode */
88 NULL, /* incidence */
89 e_undefined, /* e_token */
90 0, /* n_incidences */
91 -1, /* mindex */
92 -1, /* sindex */
93 -1, /* model index */
94 (REL_INCLUDED) /* flags */
95 };
96 /*
97 Don't forget to update the initialization when the structure is modified.
98 */
99
100
101 /* forward decls */
102
103 static struct rel_relation *rel_create_extnode(struct rel_relation * rel
104 , struct ExtCallNode *ext
105 );
106
107 /*------------------------------------------------------------------------------
108 GENERAL 'RELATION' ROUTINES
109 */
110
111 static struct rel_relation *rel_copy(const struct rel_relation *rel){
112 struct rel_relation *newrel;
113 newrel = ASC_NEW(struct rel_relation);
114 REL_DEBUG("Copying REL_RELATION from %p to %p",rel,newrel);
115 *newrel = *rel;
116 return(newrel);
117 }
118
119 struct rel_relation *rel_create(SlvBackendToken instance
120 , struct rel_relation *newrel
121 ){
122 CONST struct relation *instance_relation;
123 struct ExtCallNode *ext;
124 enum Expr_enum ctype;
125
126 REL_DEBUG("instance = %p",IPTR(instance));
127 REL_DEBUG("REL_RELATION newrel = %p",newrel);
128
129 if(newrel==NULL){
130 /* if newrel was not provided, create new copy of a 'default relation' */
131 newrel = rel_copy(&g_rel_defaults);
132 REL_DEBUG("CREATED NEW REL_RELATION at %p",newrel);
133 }else{
134 /* else copy the default relation into the memory space we were given */
135 *newrel = g_rel_defaults;
136 REL_DEBUG("CLEARED REL_RELATION at %p, SETTING DEFAULTS", newrel);
137 }
138 assert(newrel!=NULL);
139
140 /* the rel_relation points to the instance */
141 newrel->instance = instance;
142
143 /* get the 'struct relation' object for this relation */
144 instance_relation = GetInstanceRelation(IPTR(instance),&ctype);
145
146 REL_DEBUG("The 'relation' struct is at %p",instance_relation);
147
148 REL_DEBUG("Instance %p --> RELATION = %p",IPTR(instance),instance_relation);
149 switch (ctype) {
150 case e_token:
151 newrel->type = e_rel_token;
152 break;
153 case e_opcode:
154 Asc_Panic(2,__FUNCTION__, "switching on e_opcode");
155 break;
156 case e_glassbox:
157 newrel->type = e_rel_glassbox;
158 break;
159 case e_blackbox:
160 REL_DEBUG("Blackbox...");
161 newrel->type = e_rel_blackbox;
162 ext = BlackBoxExtCall(instance_relation);
163
164 REL_DEBUG("Subject instance at %p",ExternalCallVarInstance(ext));
165 REL_DEBUG("Subject instance type '%s'",instance_typename(ExternalCallVarInstance(ext)));
166
167 if(ext){
168 REL_DEBUG("REL_EXTNODE FOUND, ATTACHING REL_RELATION TO EXT at %p",ext);
169 newrel = rel_create_extnode(newrel,ext);
170 }else{
171 REL_DEBUG("SET NODEINFO TO NULL IN NEWREL AT %p",newrel);
172 newrel->nodeinfo = NULL;
173 }
174
175 REL_DEBUG("Subject instance is at %p",ExternalCallVarInstance(ext));
176 break;
177 default:
178 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Unknown relation type in rel_create");
179 break;
180 }
181
182 return(newrel);
183 }
184
185 SlvBackendToken rel_instance(struct rel_relation *rel){
186 if (rel==NULL) return NULL;
187 return (SlvBackendToken) rel->instance;
188 }
189
190 /**
191 This is an evil routine to determine a variable pointer inside a
192 relation, given the variable's Instance (backend) pointer.
193
194 Good design will eliminate the need for this function.
195
196 Aborts with Asc_Panic in case of var not found.
197
198 @param rel relation in which we're looking
199 @param inst token that we need to match
200 @return var_variable corresponding to the inst parameter
201 */
202 static struct var_variable *rel_instance_to_var(struct rel_relation *rel,
203 SlvBackendToken inst
204 ){
205 int j, nincid;
206 struct var_variable *var;
207 struct var_variable **incid;
208
209 incid = rel_incidence_list_to_modify(rel);
210 nincid = rel_n_incidences(rel);
211
212 REL_DEBUG("Looking for var in list of %d incident on rel %p",nincid,rel);
213
214 var = NULL;
215 for(j=0;j<nincid;++j){
216 if(( var_instance(incid[j]) )==inst){
217 var = incid[j];
218 break;
219 }
220 }
221 if(var==NULL){
222 Asc_Panic(2,__FUNCTION__,"Var not found");
223 }
224
225 return var;
226 }
227
228 void rel_write_name(slv_system_t sys,struct rel_relation *rel,FILE *fp){
229 if(rel == NULL || fp==NULL) return;
230 if(sys!=NULL) {
231 WriteInstanceName(fp,rel_instance(rel),slv_instance(sys));
232 }else{
233 WriteInstanceName(fp,rel_instance(rel),NULL);
234 }
235 }
236
237 void rel_destroy(struct rel_relation *rel){
238 struct Instance *inst;
239 switch (rel->type) {
240 case e_rel_token:
241 break;
242 default:
243 break;
244 }
245 if (rel->nodeinfo) {
246 rel->nodeinfo->cache = NULL;
247 ascfree(rel->nodeinfo);
248 rel->nodeinfo = NULL;
249 }
250 ascfree((POINTER)rel->incidence);
251 inst = IPTR(rel->instance);
252 if (inst) {
253 if (GetInterfacePtr(inst)==rel) { /* we own the pointer -- reset it. */
254 SetInterfacePtr(inst,NULL);
255 }
256 }
257 ascfree((POINTER)rel);
258 }
259
260 uint32 rel_flags( struct rel_relation *rel){
261 return rel->flags;
262 }
263
264 void rel_set_flags(struct rel_relation *rel, uint32 flags){
265 rel->flags = flags;
266 }
267
268 uint32 rel_flagbit(struct rel_relation *rel, uint32 one){
269 if (rel==NULL || rel->instance == NULL) {
270 FPRINTF(stderr,"ERROR: rel_flagbit called with bad var.\n");
271 return 0;
272 }
273 return (rel->flags & one);
274 }
275
276 void rel_set_flagbit(struct rel_relation *rel, uint32 field,uint32 one){
277 if (one) {
278 rel->flags |= field;
279 } else {
280 rel->flags &= ~field;
281 }
282 }
283
284 /* this needs to be reimplemented properly. */
285 boolean rel_less(struct rel_relation *rel){
286 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
287 case e_notequal:
288 case e_less:
289 case e_lesseq:
290 return TRUE;
291 default:
292 return FALSE;
293 }
294 }
295
296 /**
297 this needs to be reimplemented properly.
298 */
299 boolean rel_equal( struct rel_relation *rel){
300 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
301 case e_equal:
302 case e_lesseq:
303 case e_greatereq:
304 return TRUE;
305 default:
306 return FALSE;
307 }
308 }
309
310 boolean rel_greater(struct rel_relation *rel){
311 switch( RelationRelop(GetInstanceRelationOnly(IPTR(rel->instance))) ) {
312 case e_notequal:
313 case e_greater:
314 case e_greatereq:
315 return TRUE;
316 default:
317 return FALSE;
318 }
319 }
320
321 static enum rel_enum compenum2relenum(enum Expr_enum t){
322 switch (t) {
323 case e_equal: return e_rel_equal;
324 case e_less: return e_rel_less;
325 case e_greater: return e_rel_greater;
326 case e_lesseq: return e_rel_lesseq;
327 case e_greatereq:return e_rel_greatereq;
328 default:
329 FPRINTF(ASCERR,"ERROR (rel.c): compenum2relenum miscalled.\n");
330 return e_rel_not_equal;
331 }
332 }
333 enum rel_enum rel_relop( struct rel_relation *rel){
334 return
335 compenum2relenum(RelationRelop(
336 GetInstanceRelationOnly(IPTR(rel->instance))));
337 }
338
339 char *rel_make_name(slv_system_t sys,struct rel_relation *rel){
340 return WriteInstanceNameString(IPTR(rel->instance),IPTR(slv_instance(sys)));
341 }
342
343 int32 rel_mindex( struct rel_relation *rel){
344 return( rel->mindex );
345 }
346
347 void rel_set_mindex( struct rel_relation *rel, int32 index){
348 rel->mindex = index;
349 }
350
351 int32 rel_sindex( const struct rel_relation *rel){
352 return( rel->sindex );
353 }
354
355 void rel_set_sindex( struct rel_relation *rel, int32 index){
356 rel->sindex = index;
357 }
358
359 int32 rel_model(const struct rel_relation *rel){
360 return((const int32) rel->model );
361 }
362
363 void rel_set_model( struct rel_relation *rel, int32 index){
364 rel->model = index;
365 }
366
367 real64 rel_residual(struct rel_relation *rel){
368 return( RelationResidual(GetInstanceRelationOnly(IPTR(rel->instance))));
369 }
370
371 void rel_set_residual( struct rel_relation *rel, real64 residual){
372 struct relation *reln;
373 reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
374 SetRelationResidual(reln,residual);
375 }
376
377 real64 rel_multiplier(struct rel_relation *rel){
378 return( RelationMultiplier(GetInstanceRelationOnly(IPTR(rel->instance))));
379 }
380
381 void rel_set_multiplier(struct rel_relation *rel, real64 multiplier){
382 struct relation *reln;
383 reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
384 SetRelationMultiplier(reln,multiplier);
385 }
386
387 real64 rel_nominal( struct rel_relation *rel){
388 return( RelationNominal(GetInstanceRelationOnly(IPTR(rel->instance))));
389 }
390
391 void rel_set_nominal( struct rel_relation *rel, real64 nominal){
392 struct relation *reln;
393 reln = (struct relation *)GetInstanceRelationOnly(IPTR(rel->instance));
394 SetRelationNominal(reln,nominal);
395 }
396
397 /**
398 too bad there's no entry point that rel must call before being used
399 generally, like the FindType checking stuff in var.c
400 */
401 static void check_included_flag(void){
402 if (INCLUDED_R == NULL || AscFindSymbol(INCLUDED_R) == NULL) {
403 INCLUDED_R = AddSymbolL("included",8);
404 }
405 }
406 uint32 rel_included( struct rel_relation *rel){
407 struct Instance *c;
408 check_included_flag();
409 c = ChildByChar(IPTR(rel->instance),INCLUDED_R);
410 if( c == NULL ) {
411 FPRINTF(stderr,"ERROR: (rel) rel_included\n");
412 FPRINTF(stderr," No 'included' field in relation.\n");
413 WriteInstance(stderr,IPTR(rel->instance));
414 return FALSE;
415 }
416 rel_set_flagbit(rel,REL_INCLUDED,GetBooleanAtomValue(c));
417 return( GetBooleanAtomValue(c) );
418 }
419
420 void rel_set_included( struct rel_relation *rel, uint32 included){
421 struct Instance *c;
422 check_included_flag();
423 c = ChildByChar(IPTR(rel->instance),INCLUDED_R);
424 if( c == NULL ) {
425 FPRINTF(stderr,"ERROR: (rel) rel_set_included\n");
426 FPRINTF(stderr," No 'included' field in relation.\n");
427 WriteInstance(stderr,IPTR(rel->instance));
428 return;
429 }
430 SetBooleanAtomValue(c,included,(unsigned)0);
431 rel_set_flagbit(rel,REL_INCLUDED,included);
432 }
433
434 int32 rel_apply_filter( const struct rel_relation *rel, rel_filter_t *filter){
435 if (rel==NULL || filter==NULL) {
436 FPRINTF(stderr,"rel_apply_filter miscalled with NULL\n");
437 return FALSE;
438 }
439 /* AND to mask off irrelevant bits in flags and match value, then compare */
440 return ( (filter->matchbits & rel->flags) ==
441 (filter->matchbits & filter->matchvalue)
442 );
443 }
444
445 /**
446 Implementation function for rel_n_incidences(). Do not call
447 this function directly - use rel_n_incidences() instead.
448 */
449 int32 rel_n_incidencesF(struct rel_relation *rel){
450 if (rel!=NULL) return rel->n_incidences;
451 FPRINTF(stderr,"rel_n_incidences miscalled with NULL\n");
452 return 0;
453 }
454
455 void rel_set_incidencesF(struct rel_relation *rel,int32 n,
456 struct var_variable **i
457 ){
458 if(rel!=NULL && n >=0) {
459 if (n && i==NULL) {
460 FPRINTF(stderr,"rel_set_incidence miscalled with NULL ilist\n");
461 }
462 rel->n_incidences = n;
463 rel->incidence = i;
464 return;
465 }
466 FPRINTF(stderr,"rel_set_incidence miscalled with NULL or n < 0\n");
467 }
468
469 const struct var_variable **rel_incidence_list( struct rel_relation *rel){
470 if (rel==NULL) return NULL;
471 return( (const struct var_variable **)rel->incidence );
472 }
473
474 struct var_variable **rel_incidence_list_to_modify( struct rel_relation *rel){
475 if (rel==NULL) return NULL;
476 return( (struct var_variable **)rel->incidence );
477 }
478
479 #if KILL
480 expr_t rel_lhs( struct rel_relation *rel){
481 if (rel==NULL) return NULL;
482 return( rel->lhs );
483 }
484
485 expr_t rel_rhs( struct rel_relation *rel){
486 if (rel==NULL) return NULL;
487 return( rel->rhs );
488 }
489 #endif /* KILL */
490
491 /*==============================================================================
492 EXTERNAL RELATIONS CACHE
493
494 External Relations Cache for solvers.
495 by Kirk Andre Abbott
496 Created: 08/10/94
497 */
498
499 double g_external_tolerance = 1.0e-12;
500
501 /* - - - - - - - -
502 REL->EXTREL ACCESS FUNCTIONS
503 */
504
505 static struct rel_relation *rel_create_extnode(struct rel_relation * rel
506 , struct ExtCallNode *ext
507 ){
508 struct rel_extnode *nodeinfo;
509 /* struct Instance *inst; */
510
511 /* REL_DEBUG("Creating rel_extnode"); */
512 /* REL_DEBUG("REL = %p",rel); */
513 nodeinfo = ASC_NEW(struct rel_extnode);
514 nodeinfo->whichvar = (int)ExternalCallVarIndex(ext);
515 asc_assert(nodeinfo->whichvar >= 1);
516 nodeinfo->cache = NULL;
517 rel->nodeinfo = nodeinfo;
518
519 /* inst = ExternalCallVarInstance(ext); */
520 /* REL_DEBUG("rel_extnode whichvar IS INSTANCE AT %p",inst); */
521 /* REL_DEBUG("INSTANCE type is %s",instance_typename(inst)); */
522
523 /* REL_DEBUG("REL NODEINFO = %p",rel->nodeinfo); */
524 return rel;
525 }
526
527 struct rel_extnode *rel_extnodeinfo( struct rel_relation *rel){
528 /* REL_DEBUG("REL NODEINFO = %p",rel->nodeinfo); */
529 return(rel->nodeinfo);
530 }
531
532 unsigned long rel_extwhichvar( struct rel_relation *rel){
533 /* REL_DEBUG("REL NODEINFO = %p",rel->nodeinfo); */
534 if (rel->nodeinfo) {
535 return(rel->nodeinfo->whichvar);
536 } else {
537 return 0; /* never a valid index */
538 }
539 }
540
541 struct ExtRelCache *rel_extcache( struct rel_relation *rel){
542 /* REL_DEBUG("REL NODEINFO = %p",rel->nodeinfo); */
543 if(rel->nodeinfo!=NULL){
544 return(rel->nodeinfo->cache);
545 }else{
546 return NULL;
547 }
548 }
549
550 /**
551 This function is naughty!
552 */
553 struct Instance *rel_extsubject(struct rel_relation *rel){
554 unsigned long subject = rel_extwhichvar(rel);
555 struct ExtRelCache *cache = rel_extcache(rel);
556 return GetSubjectInstance(cache->arglist,subject);
557 }
558
559 void rel_set_extnodeinfo( struct rel_relation *rel
560 , struct rel_extnode *nodeinfo
561 ){
562 rel->nodeinfo = nodeinfo;
563 /* REL_DEBUG("REL NODEINFO = %p",rel->nodeinfo); */
564 }
565
566 void rel_set_extwhichvar(struct rel_relation *rel, int32 whichvar){
567 rel->nodeinfo->whichvar = whichvar;
568 }
569
570 void rel_set_extcache( struct rel_relation *rel,struct ExtRelCache * cache){
571 rel->nodeinfo->cache = cache;
572 }
573
574 /*------------------------------------------------------------------------------
575 EXTERNAL RELATION CACHE (EXTRELCACHE)
576 */
577
578 struct ExtRelCache *CreateExtRelCache(struct ExtCallNode *ext){
579 struct ExtRelCache *cache;
580 unsigned long n_input_args, n_output_args;
581 int32 ninputs, noutputs;
582
583 assert(ext!=NULL);
584
585 cache = ASC_NEW(struct ExtRelCache);
586 cache->user_data = NULL; /* it's vital that this is initialized to NULL ! */
587
588 /* Copy various pointers from the ExtCallNode to our cache object */
589 cache->nodestamp = ExternalCallNodeStamp(ext);
590 cache->efunc = ExternalCallExtFunc(ext);
591 cache->data = ExternalCallDataInstance(ext);
592 cache->arglist = ExternalCallArgList(ext);
593 REL_DEBUG("ASSIGNED efunc %p to ExtRelCache %p",cache->efunc,cache);
594
595 /* Fetch the size of the input/output argument lists */
596 n_input_args = NumberInputArgs(cache->efunc);
597 n_output_args = NumberOutputArgs(cache->efunc);
598
599 /* Linearise the argument lists (for fast access, presumably) */
600 cache->inputlist = LinearizeArgList(cache->arglist,1,n_input_args);
601 /* Note: we own the cache of the LinearizeArgList call. */
602
603 ninputs = (int32)gl_length(cache->inputlist);
604 noutputs = (int32)CountNumberOfArgs(cache->arglist,n_input_args+1,
605 n_input_args+n_output_args);
606 cache->ninputs = ninputs;
607 cache->noutputs = noutputs;
608
609 /*
610 Create the 'invars' and 'outvars' lists so that we can insert stuff
611 into the solverside matrix correctly
612 */
613
614 cache->invars = NULL;
615 cache->outvars = ASC_NEW_ARRAY_CLEAR(struct var_variable *,noutputs);
616
617 REL_DEBUG("ALLOCATED FOR %d OUTPUTS",noutputs);
618 cache->inputs = ASC_NEW_ARRAY_CLEAR(double,ninputs);
619 cache->outputs = ASC_NEW_ARRAY_CLEAR(double,noutputs);
620 cache->jacobian = ASC_NEW_ARRAY_CLEAR(double,ninputs*noutputs);
621
622 /* Setup default flags for controlling calculations. */
623 cache->evaluation_required = 1;
624 cache->first_func_eval = 1;
625 cache->first_deriv_eval = 1;
626
627 REL_DEBUG("NEW CACHE = %p",cache);
628 return cache;
629 }
630
631 struct ExtRelCache *CreateCacheFromInstance(SlvBackendToken relinst){
632 struct ExtCallNode *ext;
633 struct ExtRelCache *cache;
634 CONST struct relation *reln;
635 enum Expr_enum type;
636
637 REL_DEBUG("CREATE CACHE FROM INSTANCE");
638
639 assert(relinst != NULL && InstanceKind(IPTR(relinst))==REL_INST);
640 reln = GetInstanceRelation(IPTR(relinst),&type);
641 if (type!=e_blackbox) {
642 FPRINTF(stderr,"Invalid relation type in (CreateCacheFromInstance)\n");
643 return NULL;
644 }
645 ext = BlackBoxExtCall(reln);
646 cache = CreateExtRelCache(ext);
647 return cache;
648 }
649
650 void extrel_store_output_var(struct rel_relation *rel){
651 struct ExtRelCache *cache;
652 int whichvar;
653 struct var_variable *var;
654 struct Instance *inst;
655
656 cache = rel_extcache(rel);
657 inst = rel_extsubject(rel);
658 var = rel_instance_to_var(rel,inst);
659 whichvar = rel_extwhichvar(rel);
660
661 REL_DEBUG("outvar[%d] at %p",whichvar-cache->ninputs-1,var);
662 asc_assert(cache->outvars!=NULL);
663 cache->outvars[whichvar - cache->ninputs - 1] = var;
664 }
665
666 static struct var_variable *extrel_outvar(
667 struct ExtRelCache *cache
668 ,int whichvar
669 ){
670 asc_assert(cache->outvars!=NULL);
671 return cache->outvars[whichvar-cache->ninputs-1];
672 }
673
674 void ExtRel_DestroyCache(struct ExtRelCache *cache)
675 {
676 if (cache) {
677 cache->nodestamp=-1;
678 cache->efunc = NULL;
679 cache->data = NULL;
680 gl_destroy(cache->inputlist); /* we own the list */
681 ascfree(cache->inputs);
682 ascfree(cache->outputs);
683 ascfree(cache->jacobian);
684 cache->inputlist = NULL;
685 cache->inputs = NULL;
686 cache->outputs = NULL;
687 cache->jacobian = NULL;
688 ascfree(cache);
689 }
690 }
691
692 /**
693 This function creates the lists of input and output var_variables so that
694 we can insert variables where required in the overall matrix (mtx).
695
696 This might form the core of a new implementation, if we decided to keep
697 the external relation caching stuff on the solver side.
698
699 At the moment the implementation is pretty naive, because it's the only
700 part in the ExtRel stuff where we have anything to do with var_variable
701 pointers. Ideally, all references to Instances in the ExtRelCache will
702 be switched over to var_variable, so we don't have the border trafficking
703 in struct Instance pointers, and so it all goes via var_variable...?
704 */
705 void extrel_store_input_vars(struct rel_relation *rel){
706 struct ExtRelCache *cache;
707 int i, n;
708 struct var_variable *var;
709
710 cache = rel_extcache(rel);
711 n = rel_n_incidences(rel);
712
713 /* new stuff: create the 'invars' and 'outvars' lists */
714 cache->invars = ASC_NEW_ARRAY_CLEAR(struct var_variable *,cache->ninputs);
715
716 for(i=0;i<cache->ninputs;++i){
717 var = rel_instance_to_var(rel, gl_fetch(cache->inputlist,i+1));
718 cache->invars[i] = var;
719 REL_DEBUG("invar[%d] at %p",i,var);
720 }
721 }
722
723
724 /*- - - - - - -
725 'INIT' FUNCTION CALLS
726 */
727
728 int32 ExtRel_PreSolve(struct ExtRelCache *cache, int32 setup){
729 struct ExternalFunc *efunc;
730 struct Slv_Interp slv_interp;
731
732 ExtBBoxInitFunc *init_func;
733 int32 nok = 0;
734
735 if(cache==NULL){
736 ERROR_REPORTER_HERE(ASC_PROG_ERR,"cache is NULL");
737 return 1;
738 }
739
740 /* prepare parameters to pass to the init function */
741 efunc = cache->efunc;
742 REL_DEBUG("CACHE = %p",cache);
743 init_func = GetInitFunc(efunc);
744 Init_Slv_Interp(&slv_interp);
745 slv_interp.nodestamp = cache->nodestamp;
746 slv_interp.user_data = cache->user_data;
747
748 if(setup){
749 slv_interp.task = bb_first_call;
750 }else{
751 slv_interp.task = bb_last_call;
752 }
753
754 /* call the init function */
755 if(init_func!=NULL){
756 nok = (*init_func)(&slv_interp,cache->data,cache->arglist);
757
758 if (nok) {
759 ERROR_REPORTER_HERE(ASC_PROG_ERR,"Error running init function (%d)",nok);
760 return 1;
761 }
762 REL_DEBUG("Ran init function");
763 }
764
765 /* Save the user's data and update our status. */
766 cache->user_data = slv_interp.user_data;
767 cache->evaluation_required = (unsigned)1; /* force at least one calculation */
768 cache->first_func_eval = (unsigned)0;
769 return 0;
770 }
771
772 /* - - - - - - -
773 EVALUATION FUNCTIONS
774 */
775
776 /**
777 Comparison for two values being sent as inputs to an external relation.
778 If they are sufficently similar, we will avoid re-evaluating the external
779 relation.
780 */
781 static int ArgsDifferent(double new, double old){
782 if(fabs(new - old) >= g_external_tolerance){
783 return 1;
784 }else{
785 return 0;
786 }
787 }
788
789 real64 ExtRel_Evaluate_Residual(struct rel_relation *rel){
790 double value;
791 REL_DEBUG("EVALUATING RELATION %p",rel);
792 value = ExtRel_Evaluate_RHS(rel) - ExtRel_Evaluate_LHS(rel);
793 REL_DEBUG("RESIDUAL = %f",value);
794 return value;
795 }
796
797 real64 ExtRel_Evaluate_RHS(struct rel_relation *rel){
798 struct Slv_Interp slv_interp;
799 struct ExtRelCache *cache;
800 struct ExternalFunc *efunc;
801 struct Instance *arg;
802 struct gl_list_t *inputlist;
803 double value;
804 long unsigned c;
805 int32 ninputs;
806 int32 nok;
807 unsigned long whichvar;
808 int32 newcalc_reqd=0;
809
810 /* REL_DEBUG("REL_RELATION = %p",rel); */
811
812 ExtBBoxFunc *eval_func;
813
814 assert(rel_extnodeinfo(rel));
815
816 /*
817 struct rel_extnode *nodeinfo = rel_extnodeinfo(rel);
818 REL_DEBUG("REL NODEINFO = %p",nodeinfo);
819 */
820
821 cache = rel_extcache(rel);
822 efunc = cache->efunc;
823 /*
824 REL_DEBUG("CACHE = %p",cache);
825 REL_DEBUG("efunc = %p",efunc);
826 REL_DEBUG("efunc->etype = %u",efunc->etype);
827 REL_DEBUG("efunc->value = %p",efunc->u.black.value);
828 */
829
830 eval_func = GetValueFunc(efunc);
831 inputlist = cache->inputlist;
832 ninputs = cache->ninputs;
833 whichvar = rel_extwhichvar(rel);
834
835 /*
836 The determination of whether a new calculation is required needs
837 some more thought. One thing we should insist upon is that all
838 the relations for an external relation are forced into the same
839 block.
840 */
841 for (c=0;c<ninputs;c++) {
842 arg = (struct Instance *)gl_fetch(inputlist,c+1);
843 value = RealAtomValue(arg);
844 REL_DEBUG("FOR INPUT %lu, value=%f (arg at %p), CACHED=%f"
845 ,c+1, value, arg,cache->inputs[c]
846 );
847 if(ArgsDifferent(value, cache->inputs[c])){
848 REL_DEBUG("ARGS ARE DIFFERENT, new calc will be required");
849 newcalc_reqd = 1;
850 cache->inputs[c] = value;
851 }
852 }
853 value = 0;
854 /*
855 Do the calculations, if necessary. Note that we have to *ensure*
856 that we send the user the information that he provided to us.
857 We have to update our user_data after each call of the user's function
858 as he might move information around (not smart but possible), on us.
859 If a function call is made, mark a new calculation as having been,
860 done, otherwise dont.
861 */
862 if (newcalc_reqd) {
863 REL_DEBUG("NEW CALCULATION REQUIRED");
864 Init_Slv_Interp(&slv_interp);
865 slv_interp.nodestamp = cache->nodestamp;
866 slv_interp.user_data = cache->user_data;
867
868 slv_interp.task = bb_func_eval;
869
870 for (c=0;c<ninputs;c++){
871 REL_DEBUG("input %lu: value = %f",c+1, cache->inputs[c]);
872 }
873
874 nok = (*eval_func)(&slv_interp, ninputs, cache->noutputs,
875 cache->inputs, cache->outputs, cache->jacobian);
876 if(nok){
877 REL_DEBUG("EXTERNAL CALCULATION ERROR (%d)",nok);
878 /* return, don't change the output values, don't update flags */
879 return 0;
880 }else{
881 REL_DEBUG("EVAL FUNC OK");
882 }
883
884 for (c=0;c<cache->noutputs;c++){
885 REL_DEBUG("output %lu: value = %f",c+1, cache->outputs[c]);
886 }
887
888 value = cache->outputs[whichvar - ninputs - 1];
889 /* REL_DEBUG("CALCULATED VALUE IS %f",value); */
890 cache->evaluation_required = (unsigned)1; /* newcalc done */
891 cache->user_data = slv_interp.user_data; /* update user_data */
892 }
893 else{
894 value = cache->outputs[whichvar - ninputs - 1];
895 cache->evaluation_required = 0; /* a result was simply returned */
896 }
897
898 REL_DEBUG("RHS VALUE = %f",value);
899 /*
900 FIXME need to get the value of the output and return the difference from
901 the computed value and the current value as the residual
902 */
903 return value;
904 }
905
906 real64 ExtRel_Evaluate_LHS(struct rel_relation *rel){
907 struct Instance *inst;
908 double value;
909
910 REL_DEBUG("...");
911
912 assert(rel_extnodeinfo(rel));
913
914 inst = rel_extsubject(rel);
915
916 /* REL_DEBUG("VAR IS INSTANCE AT %p",inst); */
917
918 /* REL_DEBUG("INSTANCE TYPE = %s",instance_typename(inst)); */
919
920 value = RealAtomValue(inst);
921 REL_DEBUG("LHS VALUE = %f",value);
922 return value;
923 }
924
925 /*- - - - -
926 GRADIENT EVALUATION
927
928 The following code implements gradient evaluation routines for
929 external relations. The routines here will prepare the arguements
930 and call a user supplied derivative routine, if same is non-NULL.
931
932 If it's NULL, the user supplied function evaluation routine will be
933 used to compute the gradients via finite differencing.
934
935 The current solver will not necessarily call for the derivative
936 all at once. This makes it necessary to do the gradient
937 computations (user supplied or finite difference), and to cache
938 away the results, as for the residuals. Based on calculation flags, the
939 appropriate *row* of this cached jacobian will be extracted and mapped to
940 the main solve matrix.
941
942 The cached jacobian is a contiguous vector ninputs*noutputs long
943 and is loaded row wise. Indexing starts from 0. Each row corresponds
944 to the partial derivatives of the output variable (associated with
945 that row, wrt to all its incident input variables.
946
947 Careful attention needs to be paid to the way this jacobian is
948 loaded/unloaded, because of the multiple indexing schemes in use.
949 i.e, arglist's and inputlists index 1..nvars, whereas all numeric
950 vectors number from 0.
951 */
952
953 struct deriv_data {
954 var_filter_t *filter;
955 mtx_matrix_t mtx;
956 mtx_coord_t nz;
957 };
958
959
960 /**
961 This function attempts to map the information from the contiguous
962 jacobian back into the main matrix, based on the current row <=>
963 whichvar. The jacobian assumed to have been calculated.
964 Since we are operating a relation at a time, we have to find
965 out where to index into our jacobian. This index is computed as
966 follows:
967
968 @param inputlist Instance objects corresponding to the blackbox inputs, in order
969 @param whichvar Output Instance index (for use with GetSubjectInstance(cache->arglist,whichvar))
970 @param ninputs The length of the inputlist
971 @param jacobian Dense-matrix Jacobian data returned from the blackbox func.
972 @param deriv_data Data cache (See rel.c)
973
974 index = (whichvar - ninputs - 1) * ninputs
975
976 Example: a problem with 4 inputs, 3 outputs and whichvar = 6.
977 with the counting for vars 1..nvars, but the jacobian indexing
978 starting from 0 (c-wise).
979
980 v-------- first output variable
981 I I I I O O O
982 1 2 3 4 5 6 7
983 ^--------- whichvar
984
985 ------------------ grads for whichvar = 6
986 | | | |
987 row 1 1 1 1 2 2 2 2 3 3 3 3
988 col 1 2 3 4 1 2 3 4 1 2 3 4
989
990 index = 0 1 2 3 4 5 6 7 8 9 10 11
991 jacobian =2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
992
993 Hence jacobian index = (6 - 4 - 1) * 4 = 4
994
995 @NOTE This only corresponds to jacobian elements from the RHS of
996 blackbox equations. There remain minus-ones to put down the diagonal for
997 each output variable, since
998
999 RESID[j] = Y(x[i],..x[n]) - y[j]
1000
1001 (noting that it's RHS - LHS, see ExtRel_Evaluate_Residual), so
1002
1003 dRESID[j]/dy[i] = -1
1004 dRESID[j]/dx[i] = dy[j]/dx[i]
1005
1006 The latter is what's being filled in here.
1007 */
1008 static void ExtRel_MapDataToMtx(struct ExtRelCache *cache,
1009 unsigned long whichvar,
1010 struct deriv_data *d
1011 ){
1012 struct var_variable *var = NULL;
1013 double value, *ptr;
1014 //boolean used;
1015 unsigned long c;
1016 int32 index;
1017 unsigned long ninputs;
1018
1019 ninputs = cache->ninputs;
1020
1021 REL_DEBUG("whichvar = %lu, ninputs = %lu",whichvar, ninputs);
1022 index = ((int)whichvar - ninputs - 1) * ninputs;
1023 REL_DEBUG("JACOBIAN INDEX = %d",index);
1024 ptr = &(cache->jacobian[index]);
1025
1026 asc_assert(ninputs >= 0);
1027
1028 /*
1029 REL_DEBUG("Filter matchbits %x, matchvalue %x"
1030 ,d->filter->matchbits
1031 ,d->filter->matchvalue
1032 );
1033 */
1034
1035 /* for input variables, the residual is taken from the matrix */
1036 for (c=0;c<(unsigned long)ninputs;c++) {
1037 var = cache->invars[c];
1038 REL_DEBUG("invar[%lu] at %p",c+1,var);
1039 /*
1040 // this is perhaps conditional modelling stuff, broken for the moment
1041 // for this first crack, all input variables are active and in the block
1042 used = var_apply_filter(var,d->filter);
1043 if (used) {
1044 */
1045 d->nz.col = mtx_org_to_col(d->mtx,var_sindex(var));
1046 REL_DEBUG("column = %d",d->nz.col);
1047 value = ptr[c] + mtx_value(d->mtx,&(d->nz));
1048 REL_DEBUG("input %lu is used, value = %f",c,value);
1049 mtx_set_value(d->mtx,&(d->nz), value);
1050 /*
1051 // disused, continued
1052 }else{
1053 var_filter_t f2 = {VAR_ACTIVE,VAR_ACTIVE};
1054 if(!var_apply_filter(var,&f2)){
1055 REL_DEBUG("var is not active");
1056 }else{
1057 REL_DEBUG("var not used...???");
1058 }
1059 }
1060 */
1061 }
1062 }
1063
1064
1065 static double CalculateInterval(double varvalue){
1066 UNUSED_PARAMETER(varvalue);
1067
1068 return (1.0e-05);
1069 }
1070
1071 /**
1072 Evaluate jacobian elements for an external relation
1073 using finite difference (peturbation of each input variable).
1074
1075 ExtRel Finite Differencing.
1076
1077 This routine actually does the finite differencing.
1078 The jacobian is a single contiguous vector. We load information
1079 in it *row* wise. If we have noutputs x ninputs = 3 x 4, variables,
1080 then jacobian entry 4,
1081 would correspond to jacobian[1][0], i.e., = 0.5 for this eg.
1082
1083 2.0 9.0 4.0 6.0
1084 0.5 1.3 0.0 9.7
1085 80 7.0 1.0 2.5
1086
1087 2.0 9.0 4.0 6.0 0.5 1.3 0.0 9.7 80 7.0 1.0 2.5
1088 [0][0] [1][0] [2][1]
1089
1090 When we are finite differencing variable c, we will be loading
1091 jacobian positions c, c+ninputs, c+2*ninputs ....
1092 */
1093 static int32 ExtRel_FDiff(struct Slv_Interp *slv_interp,
1094 ExtBBoxFunc *eval_func,
1095 int32 ninputs, int32 noutputs,
1096 double *inputs, double *outputs,
1097 double *jacobian
1098 ){
1099 int32 c1,c2, nok = 0;
1100 double *tmp_vector;
1101 double *ptr;
1102 double old_x,interval,value;
1103
1104 REL_DEBUG("NUMERICAL DERIVATIVE...");
1105
1106 tmp_vector = ASC_NEW_ARRAY_CLEAR(double,noutputs);
1107 for (c1=0;c1<ninputs;c1++){
1108 /* perturb x */
1109 old_x = inputs[c1];
1110 interval = CalculateInterval(old_x);
1111 inputs[c1] = old_x + interval;
1112 REL_DEBUG("PETURBATION WITH input[%d]=%f",c1+1,inputs[c1]);
1113
1114 /* call routine */
1115 nok = (*eval_func)(slv_interp, ninputs, noutputs, inputs, tmp_vector, jacobian);
1116 if(nok){
1117 REL_DEBUG("External evaluation error (%d)",nok);
1118 break;
1119 }
1120
1121 /* load jacobian */
1122 ptr = &jacobian[c1];
1123 for (c2=0;c2<noutputs;c2++) {
1124 value = (tmp_vector[c2] - outputs[c2])/interval;
1125 REL_DEBUG("output[%d]: value = %f, gradient = %f",c2+1,tmp_vector[c2],value);
1126 *ptr = value;
1127 ptr += ninputs;
1128 }
1129 inputs[c1] = old_x;
1130 }
1131 ASC_FREE(tmp_vector);
1132 if(nok){
1133 REL_DEBUG("External evaluation error");
1134 }
1135 return nok;
1136 }
1137
1138 static int32 ExtRel_CalcDeriv(struct rel_relation *rel, struct deriv_data *d){
1139 int32 nok = 0;
1140 struct Slv_Interp slv_interp;
1141 struct ExtRelCache *cache;
1142 struct ExternalFunc *efunc;
1143 unsigned long whichvar;
1144 ExtBBoxFunc *eval_func;
1145 ExtBBoxFunc *deriv_func;
1146
1147 REL_DEBUG("...");
1148
1149 assert(rel_extnodeinfo(rel));
1150 cache = rel_extcache(rel);
1151 whichvar = rel_extwhichvar(rel);
1152 efunc = cache->efunc;
1153
1154 /*
1155 * Check and deal with the special case of the first
1156 * computation.
1157 */
1158 if(cache->first_deriv_eval) {
1159 REL_DEBUG("FIRST DERIV EVAL");
1160 cache->evaluation_required = (unsigned)1;
1161 cache->first_deriv_eval = (unsigned)0;
1162 }
1163
1164 /*
1165 * If a function evaluation was not recently done, then we
1166 * can return the results from the cached jacobian.
1167 */
1168 if(!cache->evaluation_required){
1169 REL_DEBUG("NO NEW CALC DONE, RETURN CACHED JACOBIAN");
1170 ExtRel_MapDataToMtx(cache, whichvar, d);
1171 return 0;
1172 }
1173
1174 /*
1175 * If we are here, we have to do the derivative calculation.
1176 * The only thing to determine is whether we do analytical
1177 * derivatives (deriv_func != NULL) or finite differencing.
1178 * In any case init the interpreter.
1179 */
1180 Init_Slv_Interp(&slv_interp);
1181 slv_interp.task = bb_deriv_eval;
1182 slv_interp.user_data = cache->user_data;
1183
1184 deriv_func = GetDerivFunc(efunc);
1185 if(deriv_func){
1186 REL_DEBUG("USING EXTERNAL DERIVATIVE FUNCTION");
1187 nok = (*deriv_func)(&slv_interp, cache->ninputs, cache->noutputs,
1188 cache->inputs, cache->outputs, cache->jacobian);
1189 if(nok){
1190 cache->evaluation_required = 1;
1191 return nok;
1192 }
1193 }else{
1194 REL_DEBUG("USING NUMERICAL DERIVATIVE");
1195 eval_func = GetValueFunc(efunc);
1196 nok = ExtRel_FDiff(&slv_interp, eval_func,
1197 cache->ninputs, cache->noutputs,
1198 cache->inputs, cache->outputs, cache->jacobian);
1199 if(nok){
1200 cache->evaluation_required = 1;
1201 return nok;
1202 }
1203 }
1204
1205 /*
1206 * Cleanup. Ensure that we update the users data, and load
1207 * the main matrix with the derivative information.
1208 */
1209 cache->user_data = slv_interp.user_data; /* save user info */
1210 ExtRel_MapDataToMtx(cache, whichvar, d);
1211 return 0;
1212 }
1213
1214 /**
1215 ExtRel Deriv routines.
1216
1217 This is the entry point for most cases. ExtRel_CalcDeriv depends
1218 on ExtRel_Evaluate being called immediately before it.
1219 */
1220
1221 double ExtRel_Diffs_RHS(struct rel_relation *rel,
1222 var_filter_t *filter,
1223 int32 row,
1224 mtx_matrix_t mtx
1225 ){
1226 int32 nok;
1227 double rhs;
1228 struct deriv_data data;
1229
1230 data.filter = filter;
1231 data.mtx = mtx;
1232 data.nz.row = row;
1233
1234 rhs = ExtRel_Evaluate_RHS(rel);
1235 nok = ExtRel_CalcDeriv(rel,&data);
1236 if (nok)
1237 return 0.0;
1238 else
1239 return rhs;
1240 }
1241
1242
1243 double ExtRel_Diffs_LHS(struct rel_relation *rel, var_filter_t *filter,
1244 int32 row, mtx_matrix_t mtx
1245 ){
1246 double lhs;
1247 struct ExtRelCache *cache;
1248 int whichvar;
1249 struct var_variable *var;
1250 mtx_coord_t nz;
1251
1252 lhs = ExtRel_Evaluate_LHS(rel);
1253
1254 cache = rel_extcache(rel);
1255 whichvar = rel_extwhichvar(rel);
1256 var = extrel_outvar(cache,whichvar);
1257
1258 /* enter '-1' into the matrix where the output var belongs */
1259 nz.row = row;
1260 nz.col = mtx_org_to_col(mtx,var_sindex(var));
1261 mtx_set_value(mtx,&(nz), -1.0 );
1262
1263 /*
1264 REL_DEBUG("OUTPUTING MATRIX");
1265 mtx_region_t r;
1266 mtx_region(&r, 0, 2, 0, 4);
1267 mtx_write_region_human(ASCERR,mtx,&r);
1268 REL_DEBUG("Mapping LHS jacobian entry -1.0 to var at %p",var);
1269 */
1270
1271 return lhs;
1272 }
1273
1274 double extrel_resid_and_jacobian(struct rel_relation *rel
1275 , var_filter_t *filter, int32 row, mtx_matrix_t mtx
1276 ){
1277 return ExtRel_Diffs_RHS(rel,filter,row,mtx)
1278 - ExtRel_Diffs_LHS(rel,filter,row,mtx);
1279 }

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