View of /trunk/gadget/recstatistics.cc

#include "recstatistics.h"
#include "commentstream.h"
#include "readfunc.h"
#include "readword.h"
#include "errorhandler.h"
#include "areatime.h"
#include "fleet.h"
#include "stock.h"
#include "stockprey.h"
#include "mathfunc.h"
#include "readaggregation.h"
#include "gadget.h"
#include "global.h"

RecStatistics::RecStatistics(CommentStream& infile, const AreaClass* const Area,
  const TimeClass* const TimeInfo, double weight, TagPtrVector Tag, const char* name)
  : Likelihood(RECSTATISTICSLIKELIHOOD, weight, name), aggregator(0), alptr(0) {

  char text[MaxStrLength];
  strncpy(text, "", MaxStrLength);
  int i, j, check, numarea = 0;

  char datafilename[MaxStrLength];
  char aggfilename[MaxStrLength];
  strncpy(datafilename, "", MaxStrLength);
  strncpy(aggfilename, "", MaxStrLength);
  ifstream datafile;
  CommentStream subdata(datafile);

  functionname = new char[MaxStrLength];
  strncpy(functionname, "", MaxStrLength);
  readWordAndValue(infile, "datafile", datafilename);
  readWordAndValue(infile, "function", functionname);

  if ((strcasecmp(functionname, "lengthcalcvar") == 0) || (strcasecmp(functionname, "lengthcalcstddev") == 0))
    functionnumber = 1;
  else if ((strcasecmp(functionname, "lengthgivenvar") == 0) || (strcasecmp(functionname, "lengthgivenstddev") == 0))
    functionnumber = 2;
  else if ((strcasecmp(functionname, "lengthnovar") == 0) || (strcasecmp(functionname, "lengthnostddev") == 0))
    functionnumber = 3;
  else
    handle.logFileMessage(LOGFAIL, "\nError in recstatistics - unrecognised function", functionname);

  //read in area aggregation from file
  readWordAndValue(infile, "areaaggfile", aggfilename);
  datafile.open(aggfilename, ios::in);
  handle.checkIfFailure(datafile, aggfilename);
  handle.Open(aggfilename);
  numarea = readAggregation(subdata, areas, areaindex);
  handle.Close();
  datafile.close();
  datafile.clear();

  //Must change from outer areas to inner areas.
  for (i = 0; i < areas.Nrow(); i++)
    for (j = 0; j < areas.Ncol(i); j++)
      areas[i][j] = Area->getInnerArea(areas[i][j]);

  //read in the fleetnames
  i = 0;
  infile >> text >> ws;
  if (strcasecmp(text, "fleetnames") != 0)
    handle.logFileUnexpected(LOGFAIL, "fleetnames", text);
  infile >> text;
  while (!infile.eof() && (strcasecmp(text, "[component]") != 0)) {
    infile >> ws;
    fleetnames.resize(new char[strlen(text) + 1]);
    strcpy(fleetnames[i++], text);
    infile >> text;
  }
  if (fleetnames.Size() == 0)
    handle.logFileMessage(LOGFAIL, "\nError in recstatistics - failed to read fleets");
  handle.logMessage(LOGMESSAGE, "Read fleet data - number of fleets", fleetnames.Size());

  //We have now read in all the data from the main likelihood file
  //But we have to read in the statistics data from datafilename
  datafile.open(datafilename, ios::in);
  handle.checkIfFailure(datafile, datafilename);
  handle.Open(datafilename);
  readStatisticsData(subdata, TimeInfo, numarea, Tag);
  handle.Close();
  datafile.close();
  datafile.clear();

  for (j = 0; j < tagnames.Size(); j++) {
    check = 0;
    for (i = 0; i < Tag.Size(); i++) {
      if (strcasecmp(tagnames[j], Tag[i]->getName()) == 0) {
        check++;
        tagvec.resize(Tag[i]);
      }
    }
    if (check == 0)
      handle.logMessage(LOGFAIL, "Error in recstatistics - failed to match tag", tagnames[j]);
  }
}

void RecStatistics::readStatisticsData(CommentStream& infile,
  const TimeClass* TimeInfo, int numarea, TagPtrVector Tag) {

  double tmpnumber, tmpmean, tmpstddev;
  char tmparea[MaxStrLength], tmptag[MaxStrLength];
  strncpy(tmparea, "", MaxStrLength);
  strncpy(tmptag, "", MaxStrLength);
  int keepdata, needvar, readvar;
  int i, timeid, tagid, areaid, tmpindex;
  int year, step, count, reject;
  char* tagName;

  readvar = 0;
  if (functionnumber == 2)
    readvar = 1;
  needvar = 0;
  if (functionnumber == 1)
    needvar = 1;

  //Check the number of columns in the inputfile
  infile >> ws;
  if ((readvar) && (countColumns(infile) != 7))
    handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 7");
  else if ((!readvar) && (countColumns(infile) != 6))
    handle.logFileMessage(LOGFAIL, "wrong number of columns in inputfile - should be 6");

  year = step = count = reject = 0;
  while (!infile.eof()) {
    keepdata = 1;
    if (readvar)
      infile >> tmptag >> year >> step >> tmparea >> tmpnumber >> tmpmean >> tmpstddev >> ws;
    else
      infile >> tmptag >> year >> step >> tmparea >> tmpnumber >> tmpmean >> ws;

    //crude check to see if something has gone wrong and avoid infinite loops
    if (strlen(tmparea) == 0)
      handle.logFileMessage(LOGFAIL, "failed to read data from file");

    //if tmparea is in areaindex find areaid, else dont keep the data
    areaid = -1;
    for (i = 0; i < areaindex.Size(); i++)
      if (strcasecmp(areaindex[i], tmparea) == 0)
        areaid = i;

    if (areaid == -1)
      keepdata = 0;

    //check if the year and step are in the simulation
    timeid = -1;
    tagid = -1;
    if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
      for (i = 0; i < tagvec.Size(); i++)
        if (strcasecmp(tagvec[i]->getName(), tmptag) == 0)
          tagid = i;

      if (tagid == -1) {
        tmpindex = -1;
        for (i = 0; i < Tag.Size(); i++)
          if (strcasecmp(Tag[i]->getName(), tmptag) == 0)
            tmpindex = i;

        if (tmpindex == -1) {
          keepdata = 0;
        } else {
          tagName = new char[strlen(tmptag) + 1];
          strcpy(tagName, tmptag);
          tagnames.resize(tagName);
          tagid = tagnames.Size() - 1;
          Years.AddRows(1, 1, year);
          Steps.AddRows(1, 1, step);
          timeid = 0;
          numbers.resize(new DoubleMatrix(1, numarea, 0.0));
          obsMean.resize(new DoubleMatrix(1, numarea, 0.0));
          modelMean.resize(new DoubleMatrix(1, numarea, 0.0));
          if (needvar)
            modelStdDev.resize(new DoubleMatrix(1, numarea, 0.0));
          if (readvar)
            obsStdDev.resize(new DoubleMatrix(1, numarea, 0.0));
        }

      } else {
        for (i = 0; i < Years.Ncol(tagid); i++)
          if ((Years[tagid][i] == year) && (Steps[tagid][i] == step))
            timeid = i;

        //if this is a new timestep, resize to store the data
        if (timeid == -1) {
          Years[tagid].resize(1, year);
          Steps[tagid].resize(1, step);
          (*numbers[tagid]).AddRows(1, numarea, 0.0);
          (*obsMean[tagid]).AddRows(1, numarea, 0.0);
          if (readvar)
            (*obsStdDev[tagid]).AddRows(1, numarea, 0.0);
          timeid = Years.Ncol(tagid) - 1;
        }
      }

    } else
      keepdata = 0;

    if (keepdata == 1) {
      //statistics data is required, so store it
      count++;
      (*numbers[tagid])[timeid][areaid] = tmpnumber;
      (*obsMean[tagid])[timeid][areaid] = tmpmean;
      if (readvar)
        (*obsStdDev[tagid])[timeid][areaid] = tmpstddev;
    } else
      reject++;  //count number of rejected data points read from file
  }

  timeindex.resize(tagvec.Size(), -1);
  if (count == 0)
    handle.logMessage(LOGWARN, "Warning in recstatistics - found no data in the data file for", this->getName());
  if (reject != 0)
    handle.logMessage(LOGMESSAGE, "Discarded invalid recstatistics data - number of invalid entries", reject);
  handle.logMessage(LOGMESSAGE, "Read recstatistics data file - number of entries", count);
}

RecStatistics::~RecStatistics() {
  int i;
  for (i = 0; i < fleetnames.Size(); i++)
    delete[] fleetnames[i];
  for (i = 0; i < tagnames.Size(); i++)
    delete[] tagnames[i];
  for (i = 0; i < areaindex.Size(); i++)
    delete[] areaindex[i];
  for (i = 0; i < numbers.Size(); i++) {
    delete numbers[i];
    delete obsMean[i];
    delete modelMean[i];
  }
  for (i = 0; i < obsStdDev.Size(); i++)
    delete obsStdDev[i];
  for (i = 0; i < modelStdDev.Size(); i++)
    delete modelStdDev[i];
  delete[] functionname;
  if (aggregator != 0)  {
    for (i = 0; i < tagvec.Size(); i++)
      delete aggregator[i];
    delete[] aggregator;
    aggregator = 0;
  }
  delete LgrpDiv;
}

void RecStatistics::Reset(const Keeper* const keeper) {
  int i;
  Likelihood::Reset(keeper);
  for (i = 0; i < timeindex.Size(); i++)
    timeindex[i] = -1;
  for (i = 0; i < modelMean.Size(); i++)
    (*modelMean[i]).setToZero();
  if (handle.getLogLevel() >= LOGMESSAGE)
    handle.logMessage(LOGMESSAGE, "Reset recstatistics component", this->getName());
}

void RecStatistics::Print(ofstream& outfile) const {
  int i;

  outfile << "\nRecapture Statistics " << this->getName() << " - likelihood value " << likelihood
    << "\n\tFunction " << functionname;
  outfile << "\n\tFleet names:";
  for (i = 0; i < fleetnames.Size(); i++)
    outfile << sep << fleetnames[i];
  outfile << endl;
  for (i = 0; i < tagnames.Size(); i++) {
    outfile << "\tTagging experiment:\t" << tagnames[i] << endl;
    aggregator[i]->Print(outfile);
    outfile << endl;
  }
  outfile.flush();
}

void RecStatistics::printSummary(ofstream& outfile) {
  //JMB to start with just print a summary of all likelihood scores
  //This needs to be changed to print something for each timestep and area
  //but we need to sort out how to get the information in the same format as other components
  if (!(isZero(likelihood))) {
    outfile << "all   all        all" << sep << setw(largewidth) << this->getName() << sep
      << setprecision(smallprecision) << setw(smallwidth) << weight << sep
      << setprecision(largeprecision) << setw(largewidth) << likelihood << endl;
    outfile.flush();
  }
}

void RecStatistics::setFleetsAndStocks(FleetPtrVector& Fleets, StockPtrVector& Stocks) {
  int i, j, t, found, minage, maxage;
  FleetPtrVector fleets;
  StockPtrVector stocks;
  CharPtrVector stocknames;

  for (i = 0; i < fleetnames.Size(); i++) {
    found = 0;
    for (j = 0; j < Fleets.Size(); j++)
      if (strcasecmp(fleetnames[i], Fleets[j]->getName()) == 0) {
        found++;
        fleets.resize(Fleets[j]);
      }

    if (found == 0)
      handle.logMessage(LOGFAIL, "Error in recstatistics - unrecognised fleet", fleetnames[i]);
  }

  for (i = 0; i < fleets.Size(); i++)
    for (j = 0; j < fleets.Size(); j++)
      if ((strcasecmp(fleets[i]->getName(), fleets[j]->getName()) == 0) && (i != j))
        handle.logMessage(LOGFAIL, "Error in recstatistics - repeated fleet", fleets[i]->getName());

  aggregator = new RecAggregator*[tagvec.Size()];
  for (t = 0; t < tagvec.Size(); t++) {
    stocks.Reset();
    stocknames = tagvec[t]->getStockNames();
    for (i = 0; i < stocknames.Size(); i++) {
      found = 0;
      for (j = 0; j < Stocks.Size(); j++)
        if (Stocks[j]->isEaten())
          if (strcasecmp(stocknames[i], Stocks[j]->getName()) == 0) {
            found++;
            stocks.resize(Stocks[j]);
          }

      if (found == 0)
        handle.logMessage(LOGFAIL, "Error in recstatistics - unrecognised stock", stocknames[i]);
    }

    //Check the stock has been tagged
    for (j = 0; j < stocks.Size(); j++)
      if (!stocks[j]->isTagged())
        handle.logMessage(LOGFAIL, "Error in recstatistics - stocks hasnt been tagged", stocks[j]->getName());

    LgrpDiv = new LengthGroupDivision(*(stocks[0]->getPrey()->getLengthGroupDiv()));
    for (i = 1; i < stocks.Size(); i++)
      if (!LgrpDiv->Combine(stocks[i]->getPrey()->getLengthGroupDiv()))
        handle.logMessage(LOGFAIL, "Error in recstatistics - length groups not compatible");

    if (LgrpDiv->Error())
      handle.logMessage(LOGFAIL, "Error in recstatistics - failed to create length group");

    minage = 999;
    maxage = 0;
    for (i = 0; i < stocks.Size(); i++) {
      minage = min(stocks[i]->minAge(), minage);
      maxage = max(stocks[i]->maxAge(), maxage);
    }

    IntMatrix ages(1, 0, 0);
    for (i = 0; i <= maxage - minage; i++)
      ages[0].resize(1, minage + i);
    aggregator[t] = new RecAggregator(fleets, stocks, LgrpDiv, areas, ages, tagvec[t]);
  }
}

void RecStatistics::addLikelihood(const TimeClass* const TimeInfo) {
  int t, i, check;
  double l = 0.0;

  check = 0;
  for (t = 0; t < tagvec.Size(); t++) {
    timeindex[t] = -1;
    if (tagvec[t]->isWithinPeriod(TimeInfo->getYear(), TimeInfo->getStep())) {
      for (i = 0; i < Years.Ncol(t); i++) {
        if (Years[t][i] == TimeInfo->getYear() && Steps[t][i] == TimeInfo->getStep()) {
          if (check == 0)
            if (handle.getLogLevel() >= LOGMESSAGE)
              handle.logMessage(LOGMESSAGE, "Calculating likelihood score for recstatistics component", this->getName());
          timeindex[t] = i;
          aggregator[t]->Sum();
          check++;
        }
      }
    }
  }

  if (check > 0) {
    l = calcLikSumSquares();
    likelihood += l;
    if (handle.getLogLevel() >= LOGMESSAGE)
      handle.logMessage(LOGMESSAGE, "The likelihood score for this component on this timestep is", l);
  }
}

double RecStatistics::calcLikSumSquares() {
  double lik = 0.0;
  int t, area;
  double simvar, simdiff;

  for (t = 0; t < tagvec.Size(); t++) {
    if (timeindex[t] > -1) {
      alptr = &aggregator[t]->getSum();
      for (area = 0; area < alptr->Size(); area++) {
        simvar = 0.0;
        ps.calcStatistics((*alptr)[area][0], aggregator[t]->getLengthGroupDiv(), 0);
        (*modelMean[t])[timeindex[t]][area] = ps.meanLength();

        switch (functionnumber) {
          case 1:
            (*modelStdDev[t])[timeindex[t]][area] = ps.sdevLength();
            simvar = (*modelStdDev[t])[timeindex[t]][area] * (*modelStdDev[t])[timeindex[t]][area];
            break;
          case 2:
            simvar = (*obsStdDev[t])[timeindex[t]][area] * (*obsStdDev[t])[timeindex[t]][area];
            break;
          case 3:
            simvar = 1.0;
            break;
          default:
            handle.logMessage(LOGWARN, "Warning in recstatistics - unrecognised function", functionname);
            break;
        }

        if (!(isZero(simvar))) {
          simdiff = (*modelMean[t])[timeindex[t]][area] - (*obsMean[t])[timeindex[t]][area];
          lik += simdiff * simdiff * (*numbers[t])[timeindex[t]][area] / simvar;
        }
      }
    }
  }
  return lik;
}

void RecStatistics::printLikelihood(ofstream& outfile, const TimeClass* const TimeInfo) {
  int year = TimeInfo->getYear();
  int step = TimeInfo->getStep();
  int t, ti, timeid, area;

  for (t = 0; t < tagvec.Size(); t++) {
    if (tagvec[t]->isWithinPeriod(year, step)) {
      timeid = -1;
      for (ti = 0; ti < Years.Ncol(t); ti++)
        if (Years[t][ti] == year && Steps[t][ti] == step)
          timeid = ti;

      if (timeid > -1) {
        for (area = 0; area < areaindex.Size(); area++) {
          outfile << setw(printwidth) << tagnames[t] << sep << setw(lowwidth)
            << year << sep << setw(lowwidth) << step << sep << setw(printwidth)
            << areaindex[area] << sep << setprecision(printprecision)
            << setw(printwidth) << (*numbers[t])[timeid][area] << sep
            << setprecision(largeprecision) << setw(largewidth);

          //JMB crude filter to remove the 'silly' values from the output
          if ((*modelMean[t])[timeid][area] < rathersmall)
            outfile << 0;
          else
            outfile << (*modelMean[t])[timeid][area];

          switch (functionnumber) {
            case 1:
              outfile << sep << setprecision(printprecision) << setw(printwidth)
                << (*modelStdDev[t])[timeid][area] << endl;
              break;
            case 2:
              outfile << sep << setprecision(printprecision) << setw(printwidth)
                << (*obsStdDev[t])[timeid][area] << endl;
              break;
            case 3:
              outfile << endl;
              break;
            default:
              handle.logMessage(LOGWARN, "Warning in recstatistics - unrecognised function", functionname);
              break;
          }
        }
      }
    }
  }
}