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

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