GridLAB-D: core/aggregate.c Source File

00001 
00030 #include <ctype.h>
00031 #include <math.h>
00032 #include "platform.h"
00033 #include "aggregate.h"
00034 #include "output.h"
00035 #include "find.h"
00036 
00040 AGGREGATION *aggregate_mkgroup(char *aggregator, 
00041                                char *group_expression) 
00042 {
00043     AGGREGATOR op = AGGR_NOP;
00044     AGGREGATION *result=NULL;
00045     char aggrop[9], aggrval[33], aggrpart[9]="";
00046     unsigned char flags=0x00;
00047     if (sscanf(aggregator,"%8[A-Za-z0-9_](%32[A-Za-z0-9_].%8[A-Za-z0-9_])",aggrop,aggrval,aggrpart)!=3 &&
00048         sscanf(aggregator,"%8[A-Za-z0-9_](%32[A-Za-z0-9_])",aggrop,aggrval)!=2 &&
00049         (flags|=AF_ABS,
00050         sscanf(aggregator,"%8[A-Za-z0-9_]|%32[A-Za-z0-9_].%8[A-Za-z0-9_]|",aggrop,aggrval,aggrpart)!=3) &&
00051         sscanf(aggregator,"%8[A-Za-z0-9_]|%32[A-Za-z0-9_]|",aggrop,aggrval)!=2 
00052         )
00053     {
00054         output_error("aggregate group '%s' is not valid", aggregator);
00055         errno = EINVAL;
00056         return NULL;
00057     }
00058 
00059     if (stricmp(aggrop,"min")==0) op=AGGR_MIN;
00060     else if (stricmp(aggrop,"max")==0) op=AGGR_MAX;
00061     else if (stricmp(aggrop,"avg")==0) op=AGGR_AVG;
00062     else if (stricmp(aggrop,"std")==0) op=AGGR_STD;
00063     else if (stricmp(aggrop,"sum")==0) op=AGGR_SUM;
00064     else if (stricmp(aggrop,"prod")==0) op=AGGR_SUM;
00065     else if (stricmp(aggrop,"mbe")==0) op=AGGR_MBE;
00066     else if (stricmp(aggrop,"mean")==0) op=AGGR_MEAN;
00067     else if (stricmp(aggrop,"var")==0) op=AGGR_VAR;
00068     else if (stricmp(aggrop,"kur")==0) op=AGGR_KUR;
00069     else if (stricmp(aggrop,"count")==0) op=AGGR_COUNT;
00070     else if (stricmp(aggrop,"gamma")==0) op=AGGR_GAMMA;
00071     else
00072     {
00073         output_error("aggregate group '%s' does not use a known aggregator", aggregator);
00074         errno = EINVAL;
00075         return NULL;
00076     }
00077     if (op!=AGGR_NOP)
00078     {
00079         PROPERTY *pinfo=NULL;
00080         AGGRPART part=AP_NONE;
00081         FINDLIST *list=NULL;
00082 
00083         /* compile and check search program */
00084         FINDPGM *pgm = find_mkpgm(group_expression);
00085         if (pgm==NULL)
00086         {
00087             output_error("aggregate group expression '%s' is not valid", group_expression);
00088             errno = EINVAL;
00089             return NULL;
00090         }
00091         else
00092         {
00093             PGMCONSTFLAGS flags = find_pgmconstants(pgm); 
00094             
00095             /* the search must be over the same class so that the property offset is known in advance */
00096             if ((flags&CF_CLASS)!=CF_CLASS)
00097             {
00098                 output_error("aggregate group expression '%s' does not result in a set with a fixed class", group_expression);
00099                 errno = EINVAL;
00100                 free(pgm);
00101                 return NULL;
00102             }
00103             else
00104             {
00105                 OBJECT *obj;
00106                 list = find_runpgm(NULL,pgm);
00107                 if (list==NULL)
00108                 {
00109                     output_error("aggregate group expression '%s' does not result is a usable object list", group_expression);
00110                     free(pgm);
00111                     errno=EINVAL;
00112                     return NULL;
00113                 }
00114                 obj = find_first(list);
00115                 if (obj==NULL)
00116                 {
00117                     output_error("aggregate group expression '%s' results is an empty object list", group_expression);
00118                     free(pgm);
00119                     free(list);
00120                     errno=EINVAL;
00121                     return NULL;
00122                 }
00123                 pinfo = class_find_property(obj->oclass,aggrval);
00124                 if (pinfo==NULL)
00125                 {
00126                     output_error("aggregate group property '%s' is not found in the objects satisfying search criteria '%s'", aggrval, group_expression);
00127                     errno = EINVAL;
00128                     free(pgm);
00129                     free(list);
00130                     return NULL;
00131                 }
00132                 else if (pinfo->ptype==PT_double)
00133                 {
00134                     if (strcmp(aggrpart,"")!=0)
00135                     {   /* doubles cannot have parts */
00136                         output_error("aggregate group property '%s' cannot have part '%s'", aggrval, aggrpart);
00137                         errno = EINVAL;
00138                         free(pgm);
00139                         free(list);
00140                         return NULL;
00141                     }
00142                     part = AP_NONE;
00143                 }
00144                 else if (pinfo->ptype==PT_complex)
00145                 {   /* complex must have parts */
00146                     if (strcmp(aggrpart,"real")==0)
00147                         part = AP_REAL;
00148                     else if (strcmp(aggrpart,"imag")==0)
00149                         part = AP_IMAG;
00150                     else if (strcmp(aggrpart,"mag")==0)
00151                         part = AP_MAG;
00152                     else if (strcmp(aggrpart,"ang")==0)
00153                         part = AP_ANG;
00154                     else if (strcmp(aggrpart,"arg")==0)
00155                         part = AP_ARG;
00156                     else
00157                     {
00158                         output_error("aggregate group property '%s' cannot have part '%s'", aggrval, aggrpart);
00159                         errno = EINVAL;
00160                         free(pgm);
00161                         free(list);
00162                         return NULL;
00163                     }
00164                 }
00165                 else
00166                 {
00167                     output_error("aggregate group property '%s' cannot be aggregated", aggrval);
00168                     errno = EINVAL;
00169                     free(pgm);
00170                     free(list);
00171                     return NULL;
00172                 }
00173             }
00174         }
00175 
00176         /* build aggregation unit */
00177         result = (AGGREGATION*)malloc(sizeof(AGGREGATION));
00178         if (result!=NULL)
00179         {
00180             result->op = op;
00181             result->group = pgm;
00182             result->pinfo = pinfo;
00183             result->part = part;
00184             result->last = list;
00185             result->next = NULL;
00186             result->flags = flags;
00187         }
00188         else
00189         {
00190             errno=ENOMEM;
00191             free(pgm);
00192             free(list);
00193             return NULL;
00194         }
00195     }
00196 
00197     return result;
00198 }
00199 
00200 double mag(complex *x)
00201 {
00202     return sqrt(x->r*x->r + x->i*x->i);
00203 }
00204 
00205 double arg(complex *x)
00206 {
00207     return (x->r==0) ? (x->i>0 ? PI/2 : (x->i==0 ? 0 : -PI/2)) : ((x->i>0) ? (x->r>0 ? atan(x->i/x->r) : PI-atan(x->i/x->r)) : (x->r>0 ? -atan(x->i/x->r) : PI+atan(x->i/x->r)));
00208 }
00209 
00212 double aggregate_value(AGGREGATION *aggr) 
00213 {
00214     OBJECT *obj;
00215     double numerator=0, denominator=0, secondary=0;
00216 
00217     /* non-constant groups need search program rerun */
00218     if ((aggr->group->constflags&CF_CONSTANT)!=CF_CONSTANT)
00219         aggr->last = find_runpgm(NULL,aggr->group); 
00220     for (obj=find_first(aggr->last); obj!=NULL; obj=find_next(aggr->last,obj))
00221     {
00222         double value=0;
00223         double *pdouble = NULL;
00224         complex *pcomplex = NULL;
00225         switch (aggr->pinfo->ptype) {
00226         case PT_complex:
00227             pcomplex = object_get_complex(obj,aggr->pinfo);
00228             if (pcomplex!=NULL)
00229             {
00230                 switch (aggr->part) {
00231                 case AP_REAL: value=pcomplex->r; break;
00232                 case AP_IMAG: value=pcomplex->i; break;
00233                 case AP_MAG: value=mag(pcomplex); break;
00234                 case AP_ARG: value=arg(pcomplex); break;
00235                 case AP_ANG: value=arg(pcomplex)*180/PI;  break;
00236                 default: pcomplex = NULL; break; /* invalidate the result */
00237                 }
00238             }
00239             break;
00240         case PT_double:
00241             pdouble = object_get_double(obj,aggr->pinfo);
00242             if (pdouble!=NULL)
00243                 value = *pdouble;
00244             break;
00245         default:
00246             break;
00247         }
00248         if (pdouble!=NULL || pcomplex!=NULL) /* valid value */
00249         {
00250             if ((aggr->flags&AF_ABS)==AF_ABS) value=fabs(value);
00251             switch (aggr->op) {
00252             case AGGR_MIN:
00253                 if (value<numerator || denominator==0) numerator=value;
00254                 denominator = 1;
00255                 break;
00256             case AGGR_MAX:
00257                 if (value>numerator || denominator==0) numerator=value;
00258                 denominator = 1;
00259                 break;
00260             case AGGR_COUNT:
00261                 numerator++;
00262                 denominator=1;
00263                 break;
00264             case AGGR_MBE:
00265             case AGGR_AVG:
00266             case AGGR_MEAN:
00267                 numerator+=value;
00268                 denominator++;
00269                 break;
00270             case AGGR_SUM:
00271                 numerator+=value;
00272                 denominator = 1;
00273                 break;
00274             case AGGR_PROD:
00275                 numerator*=value;
00276                 denominator = 1;
00277                 break;
00278             case AGGR_GAMMA:
00279                 denominator+=log(value);
00280                 if (numerator==0 || secondary>value)
00281                     secondary = value;
00282                 numerator++;
00283                 break;
00284             case AGGR_STD:
00285             case AGGR_VAR:
00286                 denominator++;
00287                 numerator += value;
00288                 secondary += value*value;
00289                 break;
00290             case AGGR_KUR:
00291                 /* not yet supported (see below for todos) */
00292             default:
00293                 break;
00294             }
00295         }
00296     }
00297     switch (aggr->op) {
00298         double v = 0.0, t = 0.0, m = 0.0;
00299     case AGGR_GAMMA:
00300         return 1 + numerator/(denominator-numerator*log(secondary));
00301     case AGGR_STD:
00302         return sqrt((secondary + numerator*numerator/denominator)/(denominator-1));
00303     case AGGR_VAR:
00304         return (secondary + numerator*numerator/denominator) / (denominator-1);
00305     case AGGR_MBE:
00306         v = 0.0;
00307         m = numerator/denominator;
00308         for (obj=find_first(aggr->last); obj!=NULL; obj=find_next(aggr->last,obj)){
00309             t = *(object_get_double(obj,aggr->pinfo)) - m;
00310             v += ((t > 0) ? t : -t);
00311         }
00312         return v/denominator;
00313     case AGGR_KUR:
00315         throw_exception("kurtosis aggregation is not implemented");
00316     default:
00317         return numerator/denominator;
00318     }
00319 }