View of /trunk/gadget/tags.cc

#include "tags.h"
#include "areatime.h"
#include "errorhandler.h"
#include "readfunc.h"
#include "readword.h"
#include "stock.h"
#include "stockprey.h"
#include "gadget.h"
#include "global.h"

Tags::Tags(CommentStream& infile, const char* givenname, const AreaClass* const Area,
  const TimeClass* const TimeInfo, Keeper* const keeper, StockPtrVector stockvec)
  : HasName(givenname) {

  taggingstock = 0;
  numtagtimesteps = 0;
  char text[MaxStrLength];
  strncpy(text, "", MaxStrLength);
  ifstream subfile;
  CommentStream subcomment(subfile);

  keeper->addString("tags");
  keeper->addString(givenname);

  //Currently can only have one stock per tagging experiment
  readWordAndValue(infile, "stock", text);
  stocknames.resize(new char[strlen(text) + 1]);
  strcpy(stocknames[0], text);

  //Currently can only have one area per tagging experiment
  readWordAndVariable(infile, "tagarea", tagarea);
  tagarea = Area->getInnerArea(tagarea);

  infile >> ws;
  char c = infile.peek();
  if ((c == 'e') || (c == 'E'))
    readWordAndVariable(infile, "endyear", endyear);
  else
    endyear = TimeInfo->getLastYear();

  int i, found = 0;
  for (i = 0; i < stockvec.Size(); i++) {
    if (strcasecmp(stockvec[i]->getName(), stocknames[0]) == 0) {
      if (found == 0) {
        taggingstock = stockvec[i];
        LgrpDiv = new LengthGroupDivision(*(taggingstock->getLengthGroupDiv()));
        if (LgrpDiv->Error())
          handle.logMessage(LOGFAIL, "Error in tags - failed to create length group");
        tagStocks.resize(taggingstock);
        taggingstock->setTaggedStock();
      }
      found++;
    }
  }
  if (found == 0)
    handle.logMessage(LOGFAIL, "Error in tags - failed to match stock", stocknames[0]);
  if (found > 1)
    handle.logMessage(LOGFAIL, "Error in tags - repeated stock", stocknames[0]);

  //Now read in the tagloss information
  readWordAndVariable(infile, "tagloss", tagloss);
  tagloss.Inform(keeper);

  //read in the numbers format: tagid - length - number
  readWordAndValue(infile, "numbers", text);
  subfile.open(text, ios::in);
  handle.checkIfFailure(subfile, text);
  handle.Open(text);
  readNumbers(subcomment, givenname, TimeInfo);
  handle.Close();
  subfile.close();
  subfile.clear();
  keeper->clearLast();
  keeper->clearLast();
}

void Tags::readNumbers(CommentStream& infile, const char* tagname, const TimeClass* const TimeInfo) {

  int year, step, count, reject;
  int i, lenid, keepdata, timeid;
  int numlen = LgrpDiv->numLengthGroups();
  double tmplength, tmpnumber;
  char tmpname[MaxStrLength];
  strncpy(tmpname, "", MaxStrLength);

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

  year = step = count = reject = 0;
  while (!infile.eof()) {
    keepdata = 1;
    infile >> tmpname >> year >> step >> tmplength >> tmpnumber >> ws;

    //only keep the data if tmpname matches tagname
    if (strcasecmp(tagname, tmpname) != 0)
      keepdata = 0;

    //only keep the data if the length is valid
    lenid = -1;
    for (i = 0; i < numlen; i++)
      if (isEqual(tmplength, LgrpDiv->minLength(i)))
        lenid = i;

    //OK the length doesnt match a minimum length so find the length group it is in
    if ((lenid == -1) && (tmplength > LgrpDiv->minLength()) && (tmplength < LgrpDiv->maxLength())) {
      for (i = 1; i < numlen; i++) {
        if (tmplength < LgrpDiv->minLength(i)) {
          lenid = i - 1;
          break;
        }
      }
      if (lenid == -1)
        lenid = numlen - 1;  //then this must be the last length group
    }

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

    //only keep the data if the number is positive
    if (tmpnumber < 0.0) {
      handle.logMessage(LOGWARN, "Warning in tags - found negative number of tags", tmpnumber);
      keepdata = 0;
    }

    timeid = -1;
    if ((TimeInfo->isWithinPeriod(year, step)) && (keepdata == 1)) {
      for (i = 0; i < Years.Size(); i++)
        if ((Years[i] == year) && (Steps[i] == step))
          timeid = i;

      if (timeid == -1) {
        Years.resize(1, year);
        Steps.resize(1, step);
        timeid = Years.Size() - 1;
        NumberByLength.resize(new DoubleMatrix(1, numlen, 0.0));
      }

    } else
      keepdata = 0;

    if (keepdata == 1) {
      count++;
      (*NumberByLength[timeid])[0][lenid] += tmpnumber;
    } else
      reject++;  //count number of rejected data points read from file
  }

  if (count == 0)
    handle.logMessage(LOGWARN, "Warning in tags - found no data in the data file for", tagname);
  if (reject != 0)
    handle.logMessage(LOGMESSAGE, "Discarded invalid tags data - number of invalid entries", reject);
  handle.logMessage(LOGMESSAGE, "Read tags data file - number of entries", count);

  tagyear = 9999;
  tagstep = 9999;
  for (i = 0; i < Years.Size(); i++)
    if ((Years[i] < tagyear) || (Years[i] == tagyear && Steps[i] < tagstep)) {
      tagyear = Years[i];
      tagstep = Steps[i];
    }

  timeid = -1;
  for (i = 0; i < Years.Size(); i++)
    if ((Years[i] == tagyear) && (Steps[i] == tagstep))
      timeid = i;

  if (timeid == -1)
    handle.logMessage(LOGFAIL, "Error in tags - calculated invalid timestep");

  numtagtimesteps = (TimeInfo->numSteps() * (endyear - tagyear)) - tagstep + 1;
  if (endyear == TimeInfo->getLastYear())
    numtagtimesteps += TimeInfo->getLastStep();
  else
    numtagtimesteps += TimeInfo->numSteps();
}

Tags::~Tags() {
  int i;
  for (i = 0; i < stocknames.Size(); i++)
    delete[] stocknames[i];
  for (i = 0; i < NumberByLength.Size(); i++)
    delete NumberByLength[i];
  for (i = 0; i < CI.Size(); i++)
    delete CI[i];
  while (AgeLengthStock.Nrow() > 0)
    AgeLengthStock.Delete(0);
  while (NumBeforeEating.Nrow() > 0)
    NumBeforeEating.Delete(0);
  delete LgrpDiv;
}

void Tags::Reset() {
  int i;
  while (AgeLengthStock.Nrow() > 0)
    AgeLengthStock.Delete(0);
  while (NumBeforeEating.Nrow() > 0)
    NumBeforeEating.Delete(0);
  while (CI.Size() > 0) {
    delete CI[0];
    CI.Delete(0);
  }
  for (i = 0; i < preyindex.Size(); i++)
    preyindex[i] = -1;
  for (i = 0; i < updated.Size(); i++)
    updated[i] = 0;
}

void Tags::setStock(StockPtrVector& Stocks) {
  int i, j, found;
  StockPtrVector tmpStockVector;
  char* stockname;

  preyindex.resize(1, -1);
  updated.resize(1, 0);
  if (!(taggingstock->isInArea(tagarea)))
    handle.logMessage(LOGFAIL, "Error in tags - stock isnt defined on tagging area");

  if (taggingstock->doesMove()) {
    tmpStockVector = taggingstock->getTransitionStocks();
    for (i = 0; i < tmpStockVector.Size(); i++) {
      transitionStocks.resize(tmpStockVector[i]);
      preyindex.resize(1, -1);
      updated.resize(1, 0);
      tagStocks.resize(tmpStockVector[i]);
      tmpStockVector[i]->setTaggedStock();
    }
  }

  if (taggingstock->doesMature()) {
    tmpStockVector = taggingstock->getMatureStocks();
    for (i = 0; i < tmpStockVector.Size(); i++) {
      matureStocks.resize(tmpStockVector[i]);
      found = 0;
      for (j = 0; j < transitionStocks.Size(); j++)
        if (strcasecmp(transitionStocks[j]->getName(), tmpStockVector[i]->getName()) != 0)
          found++;

      if (found == 0) {
        preyindex.resize(1, -1);
        updated.resize(1, 0);
        tagStocks.resize(tmpStockVector[i]);
        tmpStockVector[i]->setTaggedStock();
      }
    }
  }

  if (taggingstock->doesStray()) {
    tmpStockVector = taggingstock->getStrayStocks();
    for (i = 0; i < tmpStockVector.Size(); i++) {
      strayStocks.resize(tmpStockVector[i]);
      found = 0;
      for (j = 0; j < transitionStocks.Size(); j++)
        if (strcasecmp(transitionStocks[j]->getName(), tmpStockVector[i]->getName()) != 0)
          found++;

      for (j = 0; j < matureStocks.Size(); j++)
        if (strcasecmp(matureStocks[j]->getName(), tmpStockVector[i]->getName()) != 0)
          found++;

      if (found == 0) {
        preyindex.resize(1, -1);
        updated.resize(1, 0);
        tagStocks.resize(tmpStockVector[i]);
        tmpStockVector[i]->setTaggedStock();
      }
    }
  }

  for (i = 1; i < tagStocks.Size(); i++) {
    stockname = new char[strlen(tagStocks[i]->getName()) + 1];
    strcpy(stockname, tagStocks[i]->getName());
    stocknames.resize(stockname);
  }
}

//Must have set stocks according to stocknames using setStock before calling Update()
//Now we need to distribute the tagged fish to the same age/length groups as the tagged stock.
void Tags::Update(int timeid) {
  int i, j;
  PopInfoVector NumberInArea;
  NumberInArea.resizeBlank(LgrpDiv->numLengthGroups());

  const AgeBandMatrix* stockPopInArea;
  const LengthGroupDivision* tmpLgrpDiv;

  stockPopInArea = &(taggingstock->getCurrentALK(tagarea));
  stockPopInArea->sumColumns(NumberInArea);

  //Now we have total number of stock per length in tagarea, N(., l) (NumberInArea) and
  //number of stock per age/length, N(a, l) (stockPopInArea) so we must initialise
  //AgeLengthStock so that it can hold all information of number of tagged stock
  //per area/age/length after endtime. We must make AgeBandMatrixPtrVector same size as
  //the one in stock even though have only one area entry at the beginning
  IntVector stockareas = taggingstock->getAreas();
  int numareas = stockareas.Size();
  int tagareaindex = -1;
  j = 0;
  while (j <= numareas && tagareaindex == -1) {
    if (tagarea == stockareas[j])
      tagareaindex = j;
    j++;
  }
  if (tagareaindex == -1)
    handle.logMessage(LOGFAIL, "Error in tags - invalid area for tagged stock");

  int maxage = stockPopInArea->maxAge();
  int minage = stockPopInArea->minAge();
  int numberofagegroups = maxage - minage + 1;
  int upperlgrp, minl, maxl, age, length, stockid;
  double numfishinarea, numstockinarea;
  IntVector lgrpsize(numberofagegroups, 0);
  IntVector lgrpmin(numberofagegroups, 0);

  for (i = 0; i < numberofagegroups; i++) {
    lgrpmin[i]= stockPopInArea->minLength(i + minage);
    upperlgrp = stockPopInArea->maxLength(i + minage);
    lgrpsize[i] = upperlgrp - lgrpmin[i];
  }

  AgeLengthStock.resize(new AgeBandMatrixPtrVector(numareas, minage, lgrpmin, lgrpsize));
  for (age = minage; age <= maxage; age++) {
    minl = stockPopInArea->minLength(age);
    maxl = stockPopInArea->maxLength(age);
    for (length = minl; length < maxl; length++) {
      numfishinarea = NumberInArea[length].N;
      numstockinarea = (*stockPopInArea)[age][length].N;
      if (numfishinarea > verysmall && numstockinarea > verysmall)
        (*AgeLengthStock[0])[tagareaindex][age][length].N = (*NumberByLength[timeid])[0][length - minl] * numstockinarea / numfishinarea;
      else
        (*AgeLengthStock[0])[tagareaindex][age][length].N = 0.0;
    }
  }
  taggingstock->addTags(AgeLengthStock[0], this, exp(-tagloss));
  updated[0] = 1;

  if (taggingstock->isEaten()) {
    tmpLgrpDiv = taggingstock->getPrey()->getLengthGroupDiv();
    lgrpmin.Reset();
    lgrpsize.Reset();
    lgrpmin.resize(numberofagegroups, 0);
    lgrpsize.resize(numberofagegroups, tmpLgrpDiv->numLengthGroups());
    NumBeforeEating.resize(new AgeBandMatrixPtrVector(numareas, minage, lgrpmin, lgrpsize));
    CI.resize(new ConversionIndex(LgrpDiv, tmpLgrpDiv));
    if (CI[CI.Size() - 1]->Error())
      handle.logMessage(LOGFAIL, "Error in tags - error when checking length structure");

    stockid = stockIndex(taggingstock->getName());
    if (stockid < 0 || stockid >= preyindex.Size())
      handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

    preyindex[stockid] = NumBeforeEating.Nrow() - 1;
  }

  for (i = 1; i < tagStocks.Size(); i++) {
    stockPopInArea = &tagStocks[i]->getCurrentALK(tagarea);
    stockareas = tagStocks[i]->getAreas();
    numareas = stockareas.Size();
    maxage = stockPopInArea->maxAge();
    minage = stockPopInArea->minAge();
    numberofagegroups = maxage - minage + 1;
    lgrpmin.Reset();
    lgrpsize.Reset();
    lgrpmin.resize(numberofagegroups, 0);
    lgrpsize.resize(numberofagegroups, 0);
    for (j = 0; j < numberofagegroups; j++) {
      lgrpmin[j] = stockPopInArea->minLength(j + minage);
      upperlgrp = stockPopInArea->maxLength(j + minage);
      lgrpsize[j] = upperlgrp - lgrpmin[j];
    }

    AgeLengthStock.resize(new AgeBandMatrixPtrVector(numareas, minage, lgrpmin, lgrpsize));
    if (tagStocks[i]->isEaten()) {
      tmpLgrpDiv = tagStocks[i]->getPrey()->getLengthGroupDiv();
      lgrpmin.Reset();
      lgrpsize.Reset();
      lgrpmin.resize(numberofagegroups, 0);
      lgrpsize.resize(numberofagegroups, tmpLgrpDiv->numLengthGroups());
      NumBeforeEating.resize(new AgeBandMatrixPtrVector(numareas, minage, lgrpmin, lgrpsize));
      CI.resize(new ConversionIndex(LgrpDiv, tmpLgrpDiv));
      if (CI[CI.Size() - 1]->Error())
        handle.logMessage(LOGFAIL, "Error in tags - error when checking length structure");

      stockid = stockIndex(tagStocks[i]->getName());
      if (stockid < 0 || stockid >= preyindex.Size())
        handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

      preyindex[stockid] = NumBeforeEating.Nrow() - 1;
    }
  }
}

void Tags::updateTags(int year, int step) {
  int i, timeid;

  timeid = -1;
  for (i = 0; i < Years.Size(); i++)
    if (Years[i] == year && Steps[i] == step)
      timeid = i;

  if (timeid != -1) {
    if (tagyear == year && tagstep == step)
      this->Update(timeid);
    else
      this->addToTagStock(timeid);
  }
}

void Tags::deleteStockTags() {
  int i;
  for (i = 0; i < tagStocks.Size(); i++) {
    if (updated[i] == 1) {
      tagStocks[i]->deleteTags(this->getName());
      updated[i] = 2;
    }
  }
}

void Tags::updateMatureStock(const TimeClass* const TimeInfo) {
  int i, id;

  if (endyear <= TimeInfo->getYear())
    handle.logMessage(LOGWARN, "Warning in tags - tagging experiment has finished");
  else {
    for (i = 0; i < matureStocks.Size(); i++) {
      id = stockIndex(matureStocks[i]->getName());
      if (id < 0 || id >= AgeLengthStock.Nrow())
        handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

      if (updated[id] == 0) {
        matureStocks[i]->addTags(AgeLengthStock[id], this, exp(-tagloss));
        updated[id] = 1;
      }
    }
  }
}

void Tags::updateTransitionStock(const TimeClass* const TimeInfo) {
  int i, id;

  if (endyear <= TimeInfo->getYear())
    handle.logMessage(LOGWARN, "Warning in tags - tagging experiment has finished");
  else {
    for (i = 0; i < transitionStocks.Size(); i++) {
      id = stockIndex(transitionStocks[i]->getName());
      if (id < 0 || id >= AgeLengthStock.Nrow())
        handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

      if (updated[id] == 0) {
        transitionStocks[i]->addTags(AgeLengthStock[id], this, exp(-tagloss));
        updated[id] = 1;
      }
    }
  }
}

void Tags::updateStrayStock(const TimeClass* const TimeInfo) {
  int i, id;

  if (endyear <= TimeInfo->getYear())
    handle.logMessage(LOGWARN, "Warning in tags - tagging experiment has finished");
  else {
    for (i = 0; i < strayStocks.Size(); i++) {
      id = stockIndex(strayStocks[i]->getName());
      if (id < 0 || id >= AgeLengthStock.Nrow())
        handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

      if (updated[id] == 0) {
        strayStocks[i]->addTags(AgeLengthStock[id], this, exp(-tagloss));
        updated[id] = 1;
      }
    }
  }
}

int Tags::stockIndex(const char* stockname) {
  int i;
  for (i = 0; i < tagStocks.Size(); i++)
    if (strcasecmp(stockname, tagStocks[i]->getName()) == 0)
      return i;

  return -1;
}

int Tags::areaIndex(const char* stockname, int area) {
  int i, j;
  for (i = 0; i < tagStocks.Size(); i++) {
    if (strcasecmp(stockname, tagStocks[i]->getName()) == 0) {
      IntVector stockareas = tagStocks[i]->getAreas();
      for (j = 0; j < stockareas.Size(); j++)
        if (stockareas[j] == area)
          return j;
      return -1;
    }
  }
  return -1;
}

void Tags::storeConsumptionALK(int area, const char* stockname) {
  int stockid, preyid, areaid;
  stockid = stockIndex(stockname);
  if (stockid < 0)
    handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

  preyid = preyindex[stockid];
  if (preyid > NumBeforeEating.Nrow() || preyid < 0)
    handle.logMessage(LOGFAIL, "Error in tags - invalid prey identifier");

  areaid = areaIndex(stockname, area);
  if (areaid < 0)
    handle.logMessage(LOGFAIL, "Error in tags - invalid area identifier");

  (*NumBeforeEating[preyid])[areaid].setToZero();
  (*NumBeforeEating[preyid])[areaid].Add((*AgeLengthStock[stockid])[areaid], *CI[preyid]);
}

const AgeBandMatrix& Tags::getConsumptionALK(int area, const char* stockname) {
  int stockid, preyid, areaid;
  stockid = stockIndex(stockname);
  if (stockid < 0)
    handle.logMessage(LOGFAIL, "Error in tags - invalid stock identifier");

  preyid = preyindex[stockid];
  if (preyid > NumBeforeEating.Nrow() || preyid < 0)
    handle.logMessage(LOGFAIL, "Error in tags - invalid prey identifier");

  areaid = areaIndex(stockname, area);
  if (areaid < 0)
    handle.logMessage(LOGFAIL, "Error in tags - invalid area identifier");

  return (*NumBeforeEating[preyid])[areaid];
}

int Tags::isWithinPeriod(int year, int step) {
  if ((year > tagyear || (year == tagyear && step >= tagstep)) && (year <= endyear))
    return 1;
  return 0;
}

void Tags::addToTagStock(int timeid) {
  int i, tagareaindex;
  PopInfoVector NumberInArea;
  NumberInArea.resizeBlank(LgrpDiv->numLengthGroups());
  const AgeBandMatrix* stockPopInArea;
  stockPopInArea = &(taggingstock->getCurrentALK(tagarea));
  stockPopInArea->sumColumns(NumberInArea);

  i = 0;
  tagareaindex = -1;
  IntVector stockareas = taggingstock->getAreas();
  while (i <= stockareas.Size() && tagareaindex == -1) {
    if (tagarea == stockareas[i])
      tagareaindex = i;
    i++;
  }
  if (tagareaindex == -1)
    handle.logMessage(LOGFAIL, "Error in tags - invalid area for tagged stock");

  int maxage = stockPopInArea->maxAge();
  int minage = stockPopInArea->minAge();
  int minl, maxl, age, length;
  double numfishinarea, numstockinarea;

  for (age = minage; age <= maxage; age++) {
    minl = stockPopInArea->minLength(age);
    maxl = stockPopInArea->maxLength(age);
    for (length = minl; length < maxl; length++) {
      numfishinarea = NumberInArea[length].N;
      numstockinarea = (*stockPopInArea)[age][length].N;
      if (numfishinarea > verysmall && numstockinarea > verysmall)
        (*AgeLengthStock[0])[tagareaindex][age][length].N += (*NumberByLength[timeid])[0][length - minl] * numstockinarea / numfishinarea;
    }
  }
}

void Tags::Print(ofstream& outfile) const {
//  outfile << "\nTag\nName" << sep << this->getName() << endl;
}