blockmath.c

Go to the documentation of this file.

/***************************************************
Apache License Version 2.0

Copyright 2007-2011 EcoEquity and Stockholm Environment Institute

Licensed under the Apache License, Version 2.0 (the "License"); you may not use this
file except in compliance with the License. You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under
the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
KIND, either express or implied. See the License for the specific language governing
permissions and limitations under the License.
****************************************************/

#include <math.h>
#include <stdlib.h>
#include "datastruct.h"
#include "blockmath.h"
#ifdef DMALLOC
#include <dmalloc.h>
#endif


struct datablock transform(struct datablock dbin, double (*f)(double), struct appdata ad) {
   struct datablock dbout;
   int row, year;

   dbout.y_start = dbin.y_start;
   dbout.y_end = dbin.y_end;
   // Constrain by time range of overall analysis
   if (ad.minyear > dbout.y_start) dbout.y_start = ad.minyear;
   if (ad.maxyear < dbout.y_end) dbout.y_end = ad.maxyear;
   dbout.data = NULL;

   allocdata(&dbout, ad);

   for (row = 0; row < ad.numrecs; row++) {
      for (year = dbout.y_start; year <= dbout.y_end; year++) {
         setval(dbout, (*f)(getval(dbin, row, year, ad)), row, year, ad);
      }
   }

   return dbout;
}

struct datablock binop(struct datablock a, struct datablock b, char op, struct appdata ad) {
   struct datablock result;
   int row, year;

   // Constrain by time range of arguments
   result.y_start = a.y_start > b.y_start ? a.y_start : b.y_start;
   result.y_end = a.y_end < b.y_end ? a.y_end : b.y_end;
   // Constrain by time range of overall analysis
   if (ad.minyear > result.y_start) result.y_start = ad.minyear;
   if (ad.maxyear < result.y_end) result.y_end = ad.maxyear;
   allocdata(&result, ad);

   for (row = 0; row < ad.numrecs; row++) {
      for (year = result.y_start; year <= result.y_end; year++) {
         switch (op) {
            case '+': setval(result, getval(a, row, year, ad) + getval(b, row, year, ad), row, year, ad); break;
            case '-': setval(result, getval(a, row, year, ad) - getval(b, row, year, ad), row, year, ad); break;
            case '*': setval(result, getval(a, row, year, ad) * getval(b, row, year, ad), row, year, ad); break;
            case '/': setval(result, getval(a, row, year, ad) / getval(b, row, year, ad), row, year, ad); break;
            default: setval(result, NAN, row, year, ad);
         }
      }
   }
   return result;
}

struct datablock scalmult(double scalar, struct datablock dbin, struct appdata ad) {
   struct datablock result;
   int row, year;

   result.y_start = dbin.y_start;
   result.y_end = dbin.y_end;
   // Constrain by time range of overall analysis
   if (ad.minyear > result.y_start) result.y_start = ad.minyear;
   if (ad.maxyear < result.y_end) result.y_end = ad.maxyear;
   allocdata(&result, ad);

   for (row = 0; row < ad.numrecs; row++) {
      for (year = result.y_start; year <= result.y_end; year++) {
         setval(result, scalar * getval(dbin, row, year, ad), row, year, ad);
      }
   }
   return result;
}

double pnorm(double x) {
   return 0.5 * erfc(-M_SQRT1_2 * x);
}

double quick_qnorm(double p) {
   double a0 = 2.30753, a1 = 0.27061, b1 = 0.99229, b2 = 0.04481;
   double num, denom;
   double sign, t;

   sign = -1;
   if (p > 0.5) {
      p = 1 - p;
      sign = 1;
   }

   t = sqrt(-2.0 * log(p));

   num = a0 + a1 * t;
   denom = 1.0 + b1 * t + b2 * t * t;

   return sign * (t - num/denom);
}

double qnorm(double p) {
   int niter;
   double pstar, pderiv, x;

   if (p < EPSILON) {
      return quick_qnorm(EPSILON);
   }
   if (p > 1 - EPSILON) {
      return quick_qnorm(1 - EPSILON);
   }

   // Get initial guess (to within 3e-3)
   x = quick_qnorm(p);

   pstar = pnorm(x);
   for (niter = 0; (fabs(pstar - p) > EPSILON) && (niter < MAXITER); niter++) {
      pderiv = dnorm(x);
      x = x + (p - pstar)/pderiv;
      pstar = pnorm(x);
   }

   if  (niter == MAXITER) {
      return NAN;
   } else {
      return x;
   }
}

double dnorm(double x) {
    return INV_SQRT_2PI * exp(-0.5 * x * x);
}

double ztransf(double y, double ybar, double sigma) {
    double z;
    
    if (ybar == 0 || sigma == 0) {
        z = NAN;
    } else if (y == 0) {
        z = -INFINITY;
    } else if (isinf(y) == 1) {
        z = INFINITY;
    } else {
        z = (1/sigma) * log(y/ybar) + 0.5 * sigma;
    }
    
    return z;
}

double gini2sdlog(double gini) {
   if (gini <= 0 || gini >= 100) {
      return NAN;
   }
   return M_SQRT2 * qnorm(0.5 * (1.0 + gini/100.0));
}

struct timeseries filter(struct datablock data, struct fixedfactor filter, struct appdata ad) {
   struct timeseries result;
   int year, i;
   double accum;

   result.y_start = data.y_start;
   result.y_end = data.y_end;
   allocts(&result, ad);

   for (year = result.y_start; year <= result.y_end; year++) {
      accum = 0.0;
      for (i = 0; i < ad.numrecs; i++) {
         accum += filter.data[i] * getval(data, i, year, ad);
      }
      settsval(result, accum, year, ad);
   }

   return result;
}

double filteryr(int year, struct datablock data, struct fixedfactor filter, struct appdata ad) {
   int i;
   double accum;

   accum = 0.0;
   for (i = 0; i < ad.numrecs; i++) {
      accum += filter.data[i] * getval(data, i, year, ad);
   }

   return accum;
}