prob.cpp

/*
 * Copyright (C) 2004 Bryant Lee
 *
 * This file is part of FProbPeriod.
 *
 * FProbPeriod is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * FProbPeriod is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with FProbPeriod; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/* This is the main function of the FProbPeriod program
 * FProbPeriod program samples points over N-shift that are shift periodic
 * with period k and finds period and preperiod.  The result gives some
 * indication of what is the average period of all points of shift period k
 * over the N-shift.
 *
 * Written by: Bryant Lee
 * Date: 9/01/05
 */

#define USEARRAYS 0
#define DEBUG 0

//Standard C++
#include <string.h>
#include <iostream>
#include <iomanip>
#include <fstream> //for file ops

//C
#include <stdlib.h> //for atoi
#include <math.h> //for pow, and ceil
#include <time.h> //to get seed for random numbers

//STL
#include <vector>
#include <list>
#include <map>

//local
#include "StringOps.h"
#include "Comp.h"
#include "StorageKey.h"
#include "StorageVal.h"
#include "RecNode.h"
#include "StorageKeyUInt.h"

//mersenne twister, pseudo-random number generator
#include "MTRand.h"

//function definitions
void convertPeriodList(const string &strInput, vector<int> &periodList);
void probExp(unsigned int shiftSize, const vector<int> &periodList,
	     Comp *comp, int numPoints, list<RecNode*> &recList,
	     string &prefix, vector<string> &commandbuffer);
void probExpTwoShift(const vector<int> &periodList,
		     Comp *comp, int numPoints, list<RecNode*> &recList,
		     string &prefix, vector<string> &commandbuffer);
void clearStorageNodes(vector<StorageVal*> &table);
void clearRecList(list<RecNode*> &recList);
void output(unsigned int shiftSize, int numPoints, string &prefix,
	    const vector<int> &periodList, Comp * comp,
	    list<RecNode*> &recList, vector<string> &commandbuffer);
void findDivisors(unsigned int x, vector<int> &divisorsList);
bool leastPeriod(byte *word, unsigned int currPeriod,
		 const vector<int> &divisorsList);
bool leastPeriod(unsigned int word, unsigned int currPeriod,
		 const vector<int> &divisorsList);

//GLOBAL

//this debug variable allows testing the "special case 2-shift" with
//arrays of unsigned ints in which each unsigned int carries UINTSIZE
//bits of info (each bit representing 1 coordinate).  Values up to 32 allowed.
int UINTSIZE = 32;

//boolean that determines whether output is given after each completed experiment (after each K)
//or only at end
bool OUTPUT_AFTER_EACH = false;

//main
int main() {
  unsigned int shiftSize; //# symbols in shift
  string prefix; //prefix for output files
  string strInput, strInput2, strInput3; //temp variables
  vector<int> periodList; //list of periods
  unsigned int i = 0, temp; //counter
  list<RecNode*> recList; //list of records
  list<RecNode*> leastRecList; //list of records for points of least period k
  Comp *comp = NULL; //composition to apply

  int numPoints;
  string outoption;

  //hold the command strings
  vector<string> commandbuffer;
  int commandbufferit;

  //grab all input from cin and put in an array of strings
  while(!cin.eof()) {
    std::getline(cin, strInput);
    commandbuffer.push_back(strInput);
  }

  commandbufferit = 0;

  //get "output after each" or "output at end"
  outoption = commandbuffer[commandbufferit++];
  trim(outoption);
  transform(outoption.begin(), outoption.end(), outoption.begin(), tolower);
  
  if(outoption == "output after each")
    OUTPUT_AFTER_EACH = true;
  else if(outoption == "output at end")
    OUTPUT_AFTER_EACH = false;
  else {
    cout << "Error on line 1: invalid value for output mode\n";
    goto cleanup;
  }

  //get filename prefix
  prefix = commandbuffer[commandbufferit++];
  trim(prefix);
  
  //get N, the number of symbols in the full shift
  strInput = commandbuffer[commandbufferit++];
  trim(strInput);
  shiftSize = atoi(strInput.c_str());

  //get set of K, the periodicities to be tested
  strInput = commandbuffer[commandbufferit++];
  convertPeriodList(strInput, periodList);

  //get m, the number of sample points
  strInput = commandbuffer[commandbufferit++];
  trim(strInput);
  numPoints = atoi(strInput.c_str());

  //Allocate composition object
  comp = new Comp(shiftSize);

  //get "Begin Definitions"
  strInput = commandbuffer[commandbufferit++];

  //make it lowercase
  trim(strInput);
  transform(strInput.begin(), strInput.end(), strInput.begin(), tolower);
  
  if(strInput.compare("begin definitions") != 0) {
    cout << "ERROR: Expecting \"BEGIN DEFINITIONS\"\n";
    goto cleanup;
  }

  //read definitions
  while(1) {
    strInput = commandbuffer[commandbufferit++];

    //strInput2 is lowercased and trimmed
    strInput2.resize(strInput.length());
    transform(strInput.begin(), strInput.end(), strInput2.begin(), tolower);
    trim(strInput2);

    //strInput3 is lowercased and NOT TRIMMED
    strInput3.resize(strInput.length());
    transform(strInput.begin(), strInput.end(), strInput3.begin(), tolower);

    if(strInput2 == "" || (strInput2[0] == '/' && strInput2[1] == '/')) {
      //ignore blank lines and allow comment out
    }
    else if(strInput2 == "begin commands") {
      break; //BREAK
    }
    else {
      i = 0;

      while(i < strInput.length()) {
	if(strInput[i] == '=') {
	  break;
	}

	i++;
      }

      if(strInput2[0] == 'o' && strInput2[1] == 'p' && strInput2[2] == 't') {
	temp = strInput3.find("opt");
	temp += 3;

	comp->addFunc(strInput.substr(i+1, strInput.length() - (i+1)),
		      strInput.substr(temp, i-temp), true);
      }
      else {
	comp->addFunc(strInput.substr(i+1, strInput.length() - (i+1)),
		      strInput.substr(0,i), false);
      }
    }
  }
  
  //read commands
  while(commandbufferit < (int)commandbuffer.size()) { //check for EOF
    strInput = commandbuffer[commandbufferit++];
    trim(strInput);

    //ignore blank lines and allow comment out with "//"
    if(strInput != "" && !(strInput[0] == '/' && strInput[1] == '/')) {

      //set composition
      if(!comp->setComposition(strInput)) {
	cout << "ERROR: Command \"" << strInput
	     << "\" references undefined function\n";
	goto cleanup;
      }

      //do experiment
      if(shiftSize == 2 && sizeof(unsigned int) == 4 && !USEARRAYS)
	probExpTwoShift(periodList, comp, numPoints, recList, prefix,
			commandbuffer);
      else
	probExp(shiftSize, periodList, comp, numPoints, recList, prefix,
		commandbuffer);

      if(!OUTPUT_AFTER_EACH)
	output(shiftSize, numPoints, prefix, periodList, comp, recList,
	       commandbuffer);
    	
      //clear records (delete the nodes with ptrs in lists)
      clearRecList(recList);
      clearRecList(leastRecList);
      recList.clear();
      leastRecList.clear();
    }
  }
  
  //cleanup and exit
  cleanup:
  //delete the nodes with ptrs in lists
  clearRecList(recList);
  clearRecList(leastRecList);
  
  //delete composition (and the functions in the Comp object's storage)
  delete comp;

  return 0;
}


/*
 * Full output
 */

void output(unsigned int shiftSize, int numPoints, string &prefix,
const vector<int> &periodList, Comp * comp, list<RecNode*> &recList,
vector<string> &commandbuffer) {
  list<RecNode*>::const_iterator it;

  unsigned int commandbufferit = 0;

  if(DEBUG)
    cout << "Output\n";

  //output
  ofstream fout;
  string fname = prefix + comp->returnName() + ".txt";
  fout.open(fname.c_str(), ios::trunc | ios::out);
  
  if(!fout.is_open()) {
    cout << "ERROR: could not create file " << fname << "\n";
    return;
  }
    
  //HEADER
  fout << "FProbPeriod 3.67\n";

  /*
  fout << "N = " << shiftSize << "\n";
  fout << "K = ";
  for(i = 0; i < periodList.size(); i++) {
    fout << periodList[i];
    if(i != periodList.size() - 1)
      fout << ",";
  }
  fout << "\n";
  fout << "m = " << numPoints << "\n";
  */

  fout << "\n";

  //copy of the input
  fout << "INPUT\n";
  fout << "--------------------\n";
  for(commandbufferit = 0; commandbufferit < commandbuffer.size(); commandbufferit++) {
    string test;
    test = commandbuffer[commandbufferit];
    trim(test);
    if(test != "") {
      fout << commandbuffer[commandbufferit] << "\n";
    }
  }
  fout << "--------------------\n";
  fout << "\n";

  //file for preperiods
  ofstream fout_preperiod;
  string fname_preperiod = prefix + comp->returnName() +
    "_preperiods.txt";
  fout_preperiod.open(fname_preperiod.c_str(), ios::trunc | ios::out);
  
  if(!fout_preperiod.is_open()) {
    cout << "ERROR: could not create file " << fname_preperiod << "\n";
    return;
  }
 
  //HEADER
  fout_preperiod << "FProbPeriod 3.67\n";
  /*
  fout_preperiod << "N = " << shiftSize << "\n";
  fout_preperiod << "K = ";
  for(i = 0; i < periodList.size(); i++) {
    fout_preperiod << periodList[i];
    if(i != periodList.size() - 1)
      fout_preperiod << ",";
  }
  fout_preperiod << "\n";
  fout_preperiod << "m = " << numPoints << "\n";
  */

  fout_preperiod << "\n";

  //copy of the input
  fout_preperiod << "INPUT\n";
  fout_preperiod << "--------------------\n";
  for(commandbufferit = 0; commandbufferit < commandbuffer.size(); commandbufferit++) {
    string test;
    test = commandbuffer[commandbufferit];
    trim(test);
    if(test != "") {
      fout_preperiod << commandbuffer[commandbufferit] << "\n";
    }
  }
  fout_preperiod << "--------------------\n";
  fout_preperiod << "\n";

  /* Output file */

  //fill on right side, show in fixed notation
  fout << setiosflags(ios::left | ios::fixed);

  fout << "COMMAND:\n";
  fout << comp->returnName() << "\n";
  fout << "WITH DEFINITIONS:\n";
  comp->printDefinitions(fout);
  fout << "\n";

  /* Preperiod file */

  //fill on right side, show in fixed notation
  fout_preperiod << setiosflags(ios::left | ios::fixed);

  fout_preperiod << "Preperiod output\n\n";

  fout_preperiod << "COMMAND:\n";
  fout_preperiod << comp->returnName() << "\n";
  fout_preperiod << "WITH DEFINITIONS:\n";
  comp->printDefinitions(fout_preperiod);
  fout_preperiod << "\n";

  /* output file */
  fout << "PERIOD K (NOT LEAST)\n\n";

  /* preperiod file */
  fout_preperiod << "PERIOD K (NOT LEAST)\n\n";

  /* output file */
  //Summary info

  fout << "K     (AvgPd)^(1/K)   (MaxPd)^(1/K)   (AvgPrepd)^(1/K)   (MaxPrepd)^(1/K)\n";

  for(it = recList.begin(); it != recList.end(); it++) {
    fout << setw(3) << (*it)->currPeriod << "   ";

    if((*it)->pending) {
      fout << "PENDING" << "\n";
    }
    else {
      fout << setw(13) << setprecision(4) <<
	pow((*it)->avgPeriod(),((double)1/((double)(*it)->currPeriod))) << "   ";
      fout << setw(13) << setprecision(4) <<
	pow((*it)->maxPeriod(),((double)1/(double)(*it)->currPeriod)) << "   ";
      fout << setw(16) << setprecision(4) <<
	pow((*it)->avgPreperiod(),((double)1/((double)(*it)->currPeriod))) << "   ";
      fout << setw(16) << setprecision(4) <<
	pow((*it)->maxPreperiod(), ((double)1/((double)(*it)->currPeriod))) <<"\n";
    }
  }
  fout <<"\n";

  fout << "K     AvgPd           MaxPd           AvgPrepd           MaxPrepd\n";

  for(it = recList.begin(); it != recList.end(); it++) {
    fout << setw(3) << (*it)->currPeriod << "   ";
    
    if((*it)->pending) {
      fout << "PENDING" << "\n";
    }
    else {
      fout << setw(13) << setprecision(2) << (*it)->avgPeriod() << "   ";
      fout << setw(13) << setprecision(2) << (*it)->maxPeriod() << "   ";
      fout << setw(16) << setprecision(2) << (*it)->avgPreperiod() << "   ";
      fout << setw(16) << setprecision(2) << (*it)->maxPreperiod() <<"\n";
    }
  }

  //orbit and period multiplicity
  fout << "\n";
  fout << "K     Period        Mult(#points)\n";
  for(it = recList.begin(); it != recList.end(); it++) {
    fout << setw(3) << (*it)->currPeriod << "   ";

    if((*it)->pending) {
      fout << "PENDING" << "\n";
    }
    else {
      map<unsigned long long, unsigned long long>::const_iterator pitpoint;
      bool pfirst = true;
      
      //note that the period and orbit maps have the same periods in the same orders
      // so we can iterate through both at once
      for(pitpoint = (*it)->periodMap.begin();
	  pitpoint != (*it)->periodMap.end(); pitpoint++) {
	if(!pfirst) {
	  fout << "      ";
	}
	else {
	pfirst = false;
	}
	
	fout << setw(11) << pitpoint->first << "   ";
	fout << setw(13) << pitpoint->second << "\n"; //points
      }
    }
  }

  /* Preperiod file   */
  
  //Preperiod multiplicity
  fout_preperiod << "K     Preperiod     Mult\n";
  for(it = recList.begin(); it != recList.end(); it++) {
    fout_preperiod << setw(3) << (*it)->currPeriod << "   ";
    
    if((*it)->pending) {
      fout_preperiod << "PENDING" << "\n";
    }
    else {
      map<unsigned long long, unsigned long long>::const_iterator pit;
      bool pfirst = true;
      
      for(pit = (*it)->preperiodMap.begin(); pit != (*it)->preperiodMap.end();
	  pit++) {
	if(!pfirst) {
	  fout_preperiod << "      ";
	}
	else {
	pfirst = false;
	}
	
	fout_preperiod << setw(11) << pit->first << "   ";
	fout_preperiod << setw(13) << pit->second << "\n";
      }
    }
  }
 
  //DONE
  fout_preperiod.close();
  fout.close();
}

void probExp(unsigned int shiftSize, const vector<int> &periodList,
	     Comp *comp, int numPoints, list<RecNode*> &recList,
	     string &prefix, vector<string> &commandbuffer) {
  unsigned int i = 0; //k = 0; //counters
  unsigned long long j;
  unsigned int currPeriod;
  int point_it;

  map<StorageKey,unsigned long long> tree; //records points that we have seen

  list<RecNode*>::const_iterator recIterator;

  //initialize random number generator
  MTRand random;

  if(DEBUG)
    cout << "Arrays\n";

  if(DEBUG)
     random.seed(2005); //hand seed

  //make a list of recNodes, one recNode per K value
  for(i =0; i < periodList.size(); i++) {
    recList.push_back(new RecNode(periodList[i], numPoints));
  }

  //initialize recIterator
  recIterator = recList.begin();

  //do for every K value
  for(i = 0; i < periodList.size(); i++) {
    currPeriod = periodList[i];

    //record keeping
    RecNode *record = *recIterator;
    
    //start the sampled test
    for(point_it = 0; point_it < numPoints; point_it++) {
      byte cWord[currPeriod]; //curr word
      byte output[currPeriod]; //helper word
      unsigned long long period = 0, preperiod = 0;
      map<StorageKey,unsigned long long>::iterator rt;

      tree.clear(); //clear the tree

      //Set current word to random point
      for(j = 0; j < currPeriod; j++) {
	cWord[j] = random.randInt(shiftSize - 1);
      }

      //insert the first word into tree
      StorageKey sk1 = StorageKey(cWord, currPeriod);
      tree.insert(pair<StorageKey,unsigned long long>(sk1,0));

      //DEBUG
      if(DEBUG >= 5) {
	cout << "\n\n-- Start orbit --\n";
	printArray(cWord, currPeriod);
      }

      //Perform experiment
      j = 1;
      while(true) {
	comp->image(cWord, currPeriod, output);
	copyWord(cWord, output, currPeriod);
	
	//DEBUG
	if(DEBUG >= 5)
	  printArray(cWord, currPeriod);
	
	StorageKey sk2 = StorageKey(cWord, currPeriod);
	
	rt = tree.find(sk2);
	if(rt != tree.end()) {
	  break;
	}
	else {
	  tree.insert(pair<StorageKey,unsigned long long>(sk2,j));
	}
	
	j++;
      }

      period = j - rt->second;
      preperiod = rt->second;
      
      record->recordPeriod(period, preperiod);
    }

    record->pending = false;

    //advance recIterator to next in list
    recIterator++;

    if(OUTPUT_AFTER_EACH) {
      output(shiftSize, numPoints, prefix, periodList, comp, recList,
	     commandbuffer);
    }
  }
}

void probExpTwoShift(const vector<int> &periodList,
		     Comp *comp, int numPoints, list<RecNode*> &recList,
		     string &prefix, vector<string> &commandbuffer) {
  unsigned int i = 0;
  unsigned long long j = 0; //k = 0; //counters
  unsigned int currPeriod;
  int point_it;

  map<StorageKeyUInt, unsigned long long> tree;
  
  list<RecNode*>::const_iterator recIterator;

  MTRand random;

  if(DEBUG)
     random.seed(2005); //hand seed

  if(DEBUG)
    cout << "Special case 2-shift\n";

  for(i = 0; i < periodList.size(); i++) {
    recList.push_back(new RecNode(periodList[i], numPoints));
  }

  recIterator = recList.begin();

  for(i = 0; i < periodList.size(); i++) {
    unsigned int blocks;
    currPeriod = periodList[i];

    blocks = (int) ceil( ((double)currPeriod) / ((double)UINTSIZE) );

    if(DEBUG) {
      cout << "K = " << currPeriod << " blocks = " << blocks << "\n";
    }

    RecNode *record = *recIterator;

    for(point_it = 0; point_it < numPoints; point_it++) {

      unsigned int cWord[blocks]; //curr word
      unsigned int output[blocks]; //helper word
      unsigned long long period = 0, preperiod = 0;
      map<StorageKeyUInt,unsigned long long>::iterator rt;

      //clear the tree
      tree.clear();

      //set cWord to a random word
      //NOTE: it is very important to clear the array first or
      //else you may start out with junk
      for(j = 0; j < blocks; j++) {
	cWord[j] = 0;
      }

      int blocki = 0, offseti = 0;
      for(j = 0; j < currPeriod; j++) {
	cWord[blocki] |= (random.randInt(1) << (offseti));

	offseti++;
	if(offseti >= UINTSIZE) {
	  blocki++;
	  offseti -= UINTSIZE;
	}
      }

      tree.insert(pair<StorageKeyUInt, unsigned long long>(StorageKeyUInt(cWord, blocks), 0));

      //DEBUG
      if(DEBUG >= 5) {
	cout << "\n-- Start Orbit --\n";
	printUIntArray(cWord, currPeriod);
	cout << "\n";
      }

      //Perform experiment
      j = 1;
      while(true) {
	//NOTE: no need to clear output[].  Comp::image will automatically
	//initialize output[] to 0's and then fill it in.

	comp->image(cWord, currPeriod, output, blocks);

	copyUIntArray(cWord, output, blocks);

	//DEBUG
	if(DEBUG >= 5) {
	  printUIntArray(cWord, currPeriod);
	  cout << "\n";
	}
	
	StorageKeyUInt sk(cWord,blocks);

	rt = tree.find(sk);
	if(rt != tree.end()) {
	  break;
	}
	else {
	  tree.insert(pair<StorageKeyUInt, unsigned long long>(sk, j));
	}
	
	j++;
      }

      period = j - rt->second;
      preperiod = rt->second;
      
      record->recordPeriod(period, preperiod);
    }

    record->pending = false;

    //advance recIterator to next in list
    recIterator++;

    if(OUTPUT_AFTER_EACH) {
      output(2, numPoints, prefix, periodList, comp, recList, commandbuffer);
    }
  }
}

void convertPeriodList(const string &strInput,
		       vector<int> &periodList) {
  vector<string> strList; //all terms
  vector<string> subList; //vector representation of an "[a,b]" term
  unsigned int strListLength = 0, i = 0;

  split(strInput,';', strList);
  strListLength = strList.size();

  for(i = 0; i < strListLength; i++) {
    string *temp = &(strList[i]);

    trim(*temp);

    if((*temp)[0] == '[') {
      int j, low, high;

      temp->erase(0,1); //remove '['
      temp->erase(temp->length() - 1, 1); //remove ']'
      
      split(*temp, ',', subList);
      
      trim(subList[0]);
      trim(subList[1]);
      
      low = atoi(subList[0].c_str());
      high = atoi(subList[1].c_str());
      
      for(j = low; j <= high; j++) {
	periodList.push_back(j);
      }

      subList.clear();
    }
    else {
      periodList.push_back(atoi(temp->c_str()));
    }
  }
}

/*
 * Delete the nodes pointed to by ptrs in recList
 */
void clearRecList(list<RecNode*> &recList) {
  list<RecNode*>::iterator it;
  
  for(it = recList.begin(); it != recList.end(); it++) {
    delete (*it);
  }
}

/*
 * Delete the nodes pointed to by the values in the table
 */
void clearStorageNodes(vector<StorageVal*> &table) {
  vector<StorageVal*>::iterator it;

  for(it = table.begin(); it != table.end(); it++) {
    delete (*it);
  }
}

/*
 * Find divisors of X, excluding X itself
 */
void findDivisors(unsigned int x, vector<int> &divisorsList) {
  unsigned int i = 2;
  unsigned int halfX = x/2;

  divisorsList.clear();

  if(x > 1) {
    divisorsList.push_back(1);
  }

  for(i = 2; i <= halfX; i++) {
    if(x % i == 0) {
      divisorsList.push_back(i);
    }
  }
}

/*
 * returns true if this word has least period p
 */
bool leastPeriod(byte *word, unsigned int currPeriod,
		 const vector<int> &divisorsList) {
  bool result = false;
  unsigned int i, j, k;

  if(currPeriod == 1)
    return true;
  
  //for every divisor d
  for(i = 0; i < divisorsList.size(); i++) {
    unsigned int d = divisorsList[i];

    result = false;

    //check in chunks of length d
    for(j = 0; j < currPeriod; j+= d) {
      //check that this chunk is same as other chunks
      for(k = 0; k < d; k++) {
	if(word[k + j] != word[k]) {
	  result = true;
	  break;
	}
      }

      if(result == true)
	break;
    }

    if(result == false)
      break;
  }
  
  return result;
}

/*
 * returns true if this word has least period p
 */
bool leastPeriod(unsigned int word, unsigned int currPeriod,
		 const vector<int> &divisorsList) {
  bool result = false;
  unsigned int i, j, k;
  
  if(currPeriod == 1)
    return true;

  //for every divisor d
  for(i = 0; i < divisorsList.size(); i++) {
    unsigned int d = divisorsList[i];

    result = false;

    //check in chunks of length d
    for(j = 0; j < currPeriod; j+= d) {
      //check that this chunk is same as other chunks
      for(k = 0; k < d; k++) {
	if(((word >> (k + j)) & 1) != ((word >> k) & 1)) {
	  result = true;
	  break;
	}
      }

      if(result == true)
	break;
    }

    if(result == false)
      break;
  }
  
  return result;
}