Home Research Papers Posters Manuals Talks Download People Help Internal

ck-cp/controlPoints.C

Go to the documentation of this file.

#include <charm++.h>

// This file is compiled twice to make a version that is capable of not needing the tracing to be turned on. 

#include "controlPoints.h"
#include "trace-controlPoints.h"
#include "LBDatabase.h"
#include "controlPoints.h"
#include "charm++.h"
#include "trace-projections.h"
#include <pathHistory.h>
#include "cp_effects.h"
#include <iostream>
#include <math.h>
#include <climits>

#if CMK_WITH_CONTROLPOINT

#define roundDouble(x)        ((long)(x+0.5))
#define myAbs(x)   (((x)>=0.0)?(x):(-1.0*(x)))
#define isInRange(v,a,b) ( ((v)<=(a)&&(v)>=(b)) || ((v)<=(b)&&(v)>=(a)) )

inline double closestInRange(double v, double a, double b){
  return (v<a) ? a : ((v>b)?b:v);
}


//  A framework for tuning "control points" exposed by an application. Tuning decisions are based upon observed performance measurements.
 

using namespace std;

#define DEFAULT_CONTROL_POINT_SAMPLE_PERIOD  10000000


//#undef DEBUGPRINT
//#define DEBUGPRINT 4

static void periodicProcessControlPoints(void* ptr, double currWallTime);


// A pointer to this PE's controlpoint manager Proxy
/* readonly */ CProxy_controlPointManager controlPointManagerProxy;
/* readonly */ int random_seed;
/* readonly */ long controlPointSamplePeriod;
/* readonly */ int whichTuningScheme;
/* readonly */ bool writeDataFileAtShutdown;
/* readonly */ bool shouldFilterOutputData;
/* readonly */ bool loadDataFileAtStartup;
/* readonly */ bool shouldGatherMemoryUsage;
/* readonly */ bool shouldGatherUtilization;
/* readonly */ bool shouldGatherAll;
/* readonly */ char CPDataFilename[512];

extern bool enableCPTracing;

std::map<std::string, int> defaultControlPointValues;



typedef enum tuningSchemeEnum {RandomSelection, SimulatedAnnealing, ExhaustiveSearch, CriticalPathAutoPrioritization, UseBestKnownTiming, UseSteering, MemoryAware, Simplex, DivideAndConquer, AlwaysDefaults, LDBPeriod, LDBPeriodLinear, LDBPeriodQuadratic, LDBPeriodOptimal}  tuningScheme;



void printTuningScheme(){
  switch(whichTuningScheme){
  case RandomSelection:
    CkPrintf("Tuning Scheme: RandomSelection\n");
    break;
  case AlwaysDefaults:
    CkPrintf("Tuning Scheme: AlwaysDefaults\n");
    break;
  case SimulatedAnnealing:
    CkPrintf("Tuning Scheme: SimulatedAnnealing\n");
    break;
  case ExhaustiveSearch:
    CkPrintf("Tuning Scheme: ExhaustiveSearch\n");
    break;
  case CriticalPathAutoPrioritization:
    CkPrintf("Tuning Scheme: CriticalPathAutoPrioritization\n");
    break;
  case UseBestKnownTiming:
    CkPrintf("Tuning Scheme: UseBestKnownTiming\n");
    break;
  case UseSteering:
    CkPrintf("Tuning Scheme: UseSteering\n");
    break;
  case MemoryAware:
    CkPrintf("Tuning Scheme: MemoryAware\n");
    break;
  case Simplex:
    CkPrintf("Tuning Scheme: Simplex Algorithm\n");
    break;
  case DivideAndConquer:
    CkPrintf("Tuning Scheme: Divide & Conquer Algorithm\n");
    break;
  case LDBPeriod:
    CkPrintf("Tuning Scheme: Load Balancing Period Steering (Constant Prediction)\n");
    break;
  case LDBPeriodLinear:
    CkPrintf("Tuning Scheme: Load Balancing Period Steering (Linear Prediction)\n");
    break;
  case LDBPeriodQuadratic:
    CkPrintf("Tuning Scheme: Load Balancing Period Steering (Quadratic Prediction)\n");
    break;
  case LDBPeriodOptimal:
    CkPrintf("Tuning Scheme: Load Balancing Period Steering (Optimal Prediction)\n");
    break;
  default:
    CkPrintf("Unknown tuning scheme\n");
    break;
  }
  fflush(stdout);
}



CkReduction::reducerType idleTimeReductionType;
CkReductionMsg *idleTimeReduction(int nMsg,CkReductionMsg **msgs){
  double ret[3];
  if(nMsg > 0){
    CkAssert(msgs[0]->getSize()==3*sizeof(double));
    double *m=(double *)msgs[0]->getData();
    ret[0]=m[0];
    ret[1]=m[1];
    ret[2]=m[2];
  }
  for (int i=1;i<nMsg;i++) {
    CkAssert(msgs[i]->getSize()==3*sizeof(double));
    double *m=(double *)msgs[i]->getData();
    ret[0]=min(ret[0],m[0]);
    ret[1]+=m[1];
    ret[2]=max(ret[2],m[2]);
  }  
  return CkReductionMsg::buildNew(3*sizeof(double),ret);   
}



CkReduction::reducerType allMeasuresReductionType;
#define ALL_REDUCTION_SIZE 12
CkReductionMsg *allMeasuresReduction(int nMsg,CkReductionMsg **msgs){
  double ret[ALL_REDUCTION_SIZE];
  if(nMsg > 0){
    CkAssert(msgs[0]->getSize()==ALL_REDUCTION_SIZE*sizeof(double));
    double *m=(double *)msgs[0]->getData();
    memcpy(ret, m, ALL_REDUCTION_SIZE*sizeof(double) );
  }
  for (int i=1;i<nMsg;i++) {
    CkAssert(msgs[i]->getSize()==ALL_REDUCTION_SIZE*sizeof(double));
    double *m=(double *)msgs[i]->getData();
    // idle time (min/sum/max)
    ret[0]=min(ret[0],m[0]);
    ret[1]+=m[1];
    ret[2]=max(ret[2],m[2]);
    // overhead time (min/sum/max)
    ret[3]=min(ret[3],m[3]);
    ret[4]+=m[4];
    ret[5]=max(ret[5],m[5]);
    // mem usage (max)
    ret[6]=max(ret[6],m[6]);
    // bytes per invocation for two types of entry methods
    ret[7]+=m[7];
    ret[8]+=m[8];
    ret[9]+=m[9];
    ret[10]+=m[10];
    // Grain size (avg)
    ret[11]+=m[11];
  }  
  return CkReductionMsg::buildNew(ALL_REDUCTION_SIZE*sizeof(double),ret);   
}


/*initproc*/ void registerCPReductions(void) {
  idleTimeReductionType=CkReduction::addReducer(idleTimeReduction, false, "idleTimeReduction");
  allMeasuresReductionType=CkReduction::addReducer(allMeasuresReduction, false, "allMeasuresReduction");
}






unsigned int randInt(unsigned int num, const char* name, int seed=0){
  CkAssert(num > 0);

  unsigned long hash = 0;
  unsigned int c;
  unsigned char * str = (unsigned char*)name;

  while ((c = *str++)){
    unsigned int c2 = (c+64)%128;
    unsigned int c3 = (c2*5953)%127;
    hash = c3 + (hash << 6) + (hash << 16) - hash;
  }

  unsigned long shuffled1 = (hash*2083)%7907;
  unsigned long shuffled2 = (seed*4297)%2017;

  unsigned long shuffled3 = (random_seed*4799)%7919;

  unsigned int namehash = shuffled3 ^ shuffled1 ^ shuffled2;

  unsigned int result = ((namehash * 6029) % 1117) % num;

  CkAssert(result >=0 && result < num);
  return result;
}



controlPointManager::controlPointManager() {
  generatedPlanForStep = -1;

    exitWhenReady = false;
    alreadyRequestedMemoryUsage = false;   
    alreadyRequestedIdleTime = false;
    alreadyRequestedAll = false;
    
    instrumentedPhase * newPhase = new instrumentedPhase();
    allData.phases.push_back(newPhase);   
    
    frameworkShouldAdvancePhase = false;
    haveControlPointChangeCallback = false;
//    CkPrintf("[%d] controlPointManager() Constructor Initializing control points, and loading data file\n", CkMyPe());
    
    phase_id = 0;

    if(loadDataFileAtStartup){    
      loadDataFile();
    }

    
    if(CkMyPe() == 0){
      CcdCallFnAfterOnPE((CcdVoidFn)periodicProcessControlPoints, (void*)NULL, controlPointSamplePeriod, CkMyPe());
    }

    traceRegisterUserEvent("No currently executing message", 5000);
    traceRegisterUserEvent("Zero time along critical path", 5010);
    traceRegisterUserEvent("Positive total time along critical path", 5020);
    traceRegisterUserEvent("env->setPathHistory()", 6000);
    traceRegisterUserEvent("Critical Path", 5900);
    traceRegisterUserEvent("Table Entry", 5901);


#if USER_EVENT_FOR_REGISTERTERMINALPATH
    traceRegisterUserEvent("registerTerminalPath", 100); 
#endif

  }
  

 controlPointManager::~controlPointManager(){
//    CkPrintf("[%d] controlPointManager() Destructor\n", CkMyPe());
  }

  void controlPointManager::pup(PUP::er &p)
  {
      // FIXME: does not work when control point is actually used,
      // just minimal pup so that it allows exit function to work (exitIfReady).
    p|generatedPlanForStep;
    p|exitWhenReady;
    p|alreadyRequestedMemoryUsage;
    p|alreadyRequestedIdleTime;
    p|alreadyRequestedAll;
    p|frameworkShouldAdvancePhase;
    p|haveControlPointChangeCallback;
    p|phase_id;
  }
  

  void controlPointManager::loadDataFile(){
    ifstream infile(CPDataFilename);
    vector<std::string> names;
    std::string line;
  
    while(getline(infile,line)){
      if(line[0] != '#')
    break;
    }
  
    int numTimings = 0;
    std::istringstream n(line);
    n >> numTimings;
  
    while(getline(infile,line)){ 
      if(line[0] != '#') 
    break; 
    }

    int numControlPointNames = 0;
    std::istringstream n2(line);
    n2 >> numControlPointNames;
  
    for(int i=0; i<numControlPointNames; i++){
      getline(infile,line);
      names.push_back(line);
    }

    for(int i=0;i<numTimings;i++){
      getline(infile,line);
      while(line[0] == '#')
    getline(infile,line); 

      instrumentedPhase * ips = new instrumentedPhase();

      std::istringstream iss(line);

      // Read memory usage for phase
      iss >> ips->memoryUsageMB;
      //     CkPrintf("Memory usage loaded from file: %d\n", ips.memoryUsageMB);

      // Read idle time
      iss >> ips->idleTime.min;
      iss >> ips->idleTime.avg;
      iss >> ips->idleTime.max;

      // Read overhead time
      iss >> ips->overheadTime.min;
      iss >> ips->overheadTime.avg;
      iss >> ips->overheadTime.max;

      // read bytePerInvoke
      iss >> ips->bytesPerInvoke;

      // read grain size
      iss >> ips->grainSize;

      // Read control point values
      for(int cp=0;cp<numControlPointNames;cp++){
    int cpvalue;
    iss >> cpvalue;
    ips->controlPoints.insert(make_pair(names[cp],cpvalue));
      }

      // ignore median time
      double mt;
      iss >> mt;

      // Read all times

      double time;

      while(iss >> time){
    ips->times.push_back(time);
#if DEBUGPRINT > 5
    CkPrintf("read time %lf from file\n", time);
#endif
      }

      allData.phases.push_back(ips);

    }

    infile.close();
  }



  void controlPointManager::writeDataFile(){
    CkPrintf("============= writeDataFile() to %s  ============\n", CPDataFilename);
    ofstream outfile(CPDataFilename);
    allData.cleanupNames();

//    string s = allData.toString();
//    CkPrintf("At end: \n %s\n", s.c_str());

    if(shouldFilterOutputData){
      allData.verify();
      allData.filterOutIncompletePhases();
    }

//    string s2 = allData.toString();
//    CkPrintf("After filtering: \n %s\n", s2.c_str());
    if(allData.toString().length() > 10){
      outfile << allData.toString();
    } else {
      outfile << " No data available to save to disk " << endl;
    }
    outfile.close();
  }

  void controlPointManager::setCPCallback(CkCallback cb, bool _frameworkShouldAdvancePhase){
    frameworkShouldAdvancePhase = _frameworkShouldAdvancePhase;
    controlPointChangeCallback = cb;
    haveControlPointChangeCallback = true;
  }


void controlPointManager::setFrameworkAdvancePhase(bool _frameworkShouldAdvancePhase){
  frameworkShouldAdvancePhase = _frameworkShouldAdvancePhase;
}

  void controlPointManager::processControlPoints(){

#if DEBUGPRINT
    CkPrintf("[%d] processControlPoints() haveControlPointChangeCallback=%d frameworkShouldAdvancePhase=%d\n", CkMyPe(), (int)haveControlPointChangeCallback, (int)frameworkShouldAdvancePhase);
#endif

    //==========================================================================================
    // Print the data for each phase

    const int s = allData.phases.size();

#if DEBUGPRINT
    CkPrintf("\n\nExamining critical paths and priorities and idle times (num phases=%d)\n", s );
    for(int p=0;p<s;++p){
      const instrumentedPhase &phase = allData.phases[p];
      const idleTimeContainer &idle = phase.idleTime;
      //      vector<MergeablePathHistory> const &criticalPaths = phase.criticalPaths;
      vector<double> const &times = phase.times;

      CkPrintf("Phase %d:\n", p);
      idle.print();
     //   CkPrintf("critical paths: (* affected by control point)\n");
    //   for(int i=0;i<criticalPaths.size(); i++){
    // If affected by a control point
    //  criticalPaths[i].print();

      //    criticalPaths[i].print();
      //    CkPrintf("\n");
    

      //   }
      CkPrintf("Timings:\n");
      for(int i=0;i<times.size(); i++){
    CkPrintf("%lf ", times[i]);
      }
      CkPrintf("\n");

    }
    
    CkPrintf("\n\n");


#endif



    //==========================================================================================
    // If this is a phase during which we try to adapt control point values based on critical path
#if 0
    if( s%5 == 4) {

      // Find the most recent phase with valid critical path data and idle time measurements      
      int whichPhase=-1;
      for(int p=0; p<s; ++p){
    const instrumentedPhase &phase = allData.phases[p];
    const idleTimeContainer &idle = phase.idleTime;
    if(idle.isValid() && phase.criticalPaths.size()>0 ){
      whichPhase = p;
    }
      }
      
      
      CkPrintf("Examining phase %d which has a valid idle time and critical paths\n", whichPhase);
      const instrumentedPhase &phase = allData.phases[whichPhase];
      const idleTimeContainer &idle = phase.idleTime;
      
      if(idle.min > 0.1){
    CkPrintf("Min PE idle is HIGH. %.2lf%% > 10%%\n", idle.min*100.0);
    
    // Determine the median critical path for this phase
    int medianCriticalPathIdx = phase.medianCriticalPathIdx();
    const PathHistory &path = phase.criticalPaths[medianCriticalPathIdx];
    CkAssert(phase.criticalPaths.size() > 0);
        CkAssert(phase.criticalPaths.size() > medianCriticalPathIdx); 
    
    CkPrintf("Median Critical Path has time %lf\n", path.getTotalTime());
    
    if(phase.times[medianCriticalPathIdx] > 1.2 * path.getTotalTime()){
      CkPrintf("The application step(%lf) is taking significantly longer than the critical path(%lf). BAD\n",phase.times[medianCriticalPathIdx], path.getTotalTime() );


      CkPrintf("Finding control points related to the critical path\n");
      int cpcount = 0;
      std::set<std::string> controlPointsAffectingCriticalPath;

      
      for(int e=0;e<path.getNumUsed();e++){
        if(path.getUsedCount(e)>0){
          int ep = path.getUsedEp(e);

          std::map<std::string, std::set<int> >::iterator iter;
          for(iter=affectsPrioritiesEP.begin(); iter!= affectsPrioritiesEP.end(); ++iter){
        if(iter->second.count(ep)>0){
          controlPointsAffectingCriticalPath.insert(iter->first);
          CkPrintf("Control Point \"%s\" affects the critical path\n", iter->first.c_str());
          cpcount++;
        }
          }
          
        }
      }
      

      if(cpcount==0){
        CkPrintf("No control points are known to affect the critical path\n");
      } else {
        double beginTime = CmiWallTimer();

        CkPrintf("Attempting to modify control point values for %d control points that affect the critical path\n", controlPointsAffectingCriticalPath.size());
        
        newControlPoints = phase.controlPoints;
        
        if(frameworkShouldAdvancePhase){
          gotoNextPhase();  
        }
        
        if(haveControlPointChangeCallback){ 
#if DEBUGPRINT
          CkPrintf("Calling control point change callback\n");
#endif
          // Create a high priority message and send it to the callback
          controlPointMsg *msg = new (8*sizeof(int)) controlPointMsg; 
          *((int*)CkPriorityPtr(msg)) = -INT_MAX;
          CkSetQueueing(msg, CK_QUEUEING_IFIFO);
          controlPointChangeCallback.send(msg);
          
        }
        
        
        // adjust the control points that can affect the critical path

        char textDescription[4096*2];
        textDescription[0] = '\0';

        std::map<std::string,int>::iterator newCP;
        for(newCP = newControlPoints.begin(); newCP != newControlPoints.end(); ++ newCP){
          if( controlPointsAffectingCriticalPath.count(newCP->first) > 0 ){
        // decrease the value (increase priority) if within range
        int lowerbound = controlPointSpace[newCP->first].first;
        if(newCP->second > lowerbound){
          newControlPointsAvailable = true;
          newCP->second --;
        }
          }
        }
       
        // Create a string for a projections user event
        if(1){
          std::map<std::string,int>::iterator newCP;
          for(newCP = newControlPoints.begin(); newCP != newControlPoints.end(); ++ newCP){
        sprintf(textDescription+strlen(textDescription), "<br>\"%s\"=%d", newCP->first.c_str(), newCP->second);
          }
        }

        char *userEventString = new char[4096*2];
        sprintf(userEventString, "Adjusting control points affecting critical path: %s ***", textDescription);
        
        static int userEventCounter = 20000;
        CkPrintf("USER EVENT: %s\n", userEventString);
        
        traceRegisterUserEvent(userEventString, userEventCounter); 
        traceUserBracketEvent(userEventCounter, beginTime, CmiWallTimer());
        userEventCounter++;
        
      }
      
    } else {
      CkPrintf("The application step(%lf) is dominated mostly by the critical path(%lf). GOOD\n",phase.times[medianCriticalPathIdx], path.getTotalTime() );
    }
    
    
      }
      
    } else {
      
      
    }
   
#endif



    if(frameworkShouldAdvancePhase){
      gotoNextPhase();  
    }
    
    if(haveControlPointChangeCallback){ 
      // Create a high priority message and send it to the callback
      controlPointMsg *msg = new (8*sizeof(int)) controlPointMsg; 
      *((int*)CkPriorityPtr(msg)) = -INT_MAX;
      CkSetQueueing(msg, CK_QUEUEING_IFIFO);
      controlPointChangeCallback.send(msg);
    }
    
  }
  
  bool controlPointManager::controlPointAffectsThisEP(int ep){
    std::map<std::string, std::set<int> >::iterator iter;
    for(iter=affectsPrioritiesEP.begin(); iter!= affectsPrioritiesEP.end(); ++iter){
      if(iter->second.count(ep)>0){
    return true;
      }
    }
    return false;    
  }
  
  bool controlPointManager::controlPointAffectsThisArray(int array){
    std::map<std::string, std::set<int> >::iterator iter;
    for(iter=affectsPrioritiesArray.begin(); iter!= affectsPrioritiesArray.end(); ++iter){
      if(iter->second.count(array)>0){
    return true;
      }
    }
    return false;   
  }
  

  instrumentedPhase * controlPointManager::currentPhaseData(){
    int s = allData.phases.size();
    CkAssert(s>=1);
    return allData.phases[s-1];
  }
 

  instrumentedPhase * controlPointManager::previousPhaseData(){
    int s = allData.phases.size();
    if(s >= 2 && phase_id > 0) {
      return allData.phases[s-2];
    } else {
      return NULL;
    }
  }
 
  instrumentedPhase * controlPointManager::twoAgoPhaseData(){
    int s = allData.phases.size();
    if(s >= 3 && phase_id > 0) {
      return allData.phases[s-3];
    } else {
      return NULL;
    }
  }
  

  void controlPointManager::gotoNextPhase(){
    CkPrintf("gotoNextPhase shouldGatherAll=%d enableCPTracing=%d\n", (int)shouldGatherAll, (int)enableCPTracing);
    fflush(stdout);
      
    if(enableCPTracing){
      if(shouldGatherAll && CkMyPe() == 0 && !alreadyRequestedAll){
    alreadyRequestedAll = true;
    CkCallback *cb = new CkCallback(CkIndex_controlPointManager::gatherAll(NULL), 0, thisProxy);
    CkPrintf("Requesting all measurements\n");
    thisProxy.requestAll(*cb);
    delete cb;
    
      } else {
    
    if(shouldGatherMemoryUsage && CkMyPe() == 0 && !alreadyRequestedMemoryUsage){
      alreadyRequestedMemoryUsage = true;
      CkCallback *cb = new CkCallback(CkIndex_controlPointManager::gatherMemoryUsage(NULL), 0, thisProxy);
      thisProxy.requestMemoryUsage(*cb);
      delete cb;
    }
    
    if(shouldGatherUtilization && CkMyPe() == 0 && !alreadyRequestedIdleTime){
      alreadyRequestedIdleTime = true;
      CkCallback *cb = new CkCallback(CkIndex_controlPointManager::gatherIdleTime(NULL), 0, thisProxy);
      thisProxy.requestIdleTime(*cb);
      delete cb;
    }
      }
    }



#if CMK_LBDB_ON && 0
    LBDatabase * myLBdatabase = LBDatabaseObj();
    LBDB * myLBDB = myLBdatabase->getLBDB();       // LBDB is Defined in LBDBManager.h
    const CkVec<LBObj*> objs = myLBDB->getObjs();
    
    LBRealType maxObjWallTime = -1.0;
    
    for(int i=0;i<objs.length();i++){
      LBObj* o = objs[i];
      const LDObjData d = o->ObjData();
      LBRealType cpuTime = d.cpuTime;
      LBRealType wallTime = d.wallTime;
      // can also get object handles from the LDObjData struct
      CkPrintf("[%d] LBDB Object[%d]: cpuTime=%f wallTime=%f\n", CkMyPe(), i, cpuTime, wallTime);
      if(wallTime > maxObjWallTime){

      }
      
    }

    myLBDB->ClearLoads(); // BUG: Probably very dangerous if we are actually using load balancing
    
#endif    


    
    // increment phase id
    phase_id++;
    

    // Create new entry for the phase we are starting now
    instrumentedPhase * newPhase = new instrumentedPhase();
    allData.phases.push_back(newPhase);
    
    CkPrintf("Now in phase %d allData.phases.size()=%d\n", phase_id, allData.phases.size());

  }

  void controlPointManager::setTiming(double time){
    currentPhaseData()->times.push_back(time);

#if USE_CRITICAL_PATH_HEADER_ARRAY
       
    // First we should register this currently executing message as a path, because it is likely an important one to consider.
    //    registerTerminalEntryMethod();
    
    // save the critical path for this phase
    //   thisPhaseData.criticalPaths.push_back(maxTerminalPathHistory);
    //    maxTerminalPathHistory.reset();
    
    
    // Reset the counts for the currently executing message
    //resetThisEntryPath();
        
#endif
    
  }
  
  void controlPointManager::requestIdleTime(CkCallback cb){
    CkAssert(enableCPTracing);
   
    double i = localControlPointTracingInstance()->idleRatio();
    double idle[3];
    idle[0] = i;
    idle[1] = i;
    idle[2] = i;
    
    //    CkPrintf("[%d] idleRatio=%f\n", CkMyPe(), i);
    
    localControlPointTracingInstance()->resetTimings();
    
    contribute(3*sizeof(double),idle,idleTimeReductionType, cb);
  }
  
  void controlPointManager::gatherIdleTime(CkReductionMsg *msg){
    CkAssert(enableCPTracing);

    int size=msg->getSize() / sizeof(double);
    CkAssert(size==3);
    double *r=(double *) msg->getData();
        
    instrumentedPhase* prevPhase = previousPhaseData();
    if(prevPhase != NULL){
      prevPhase->idleTime.min = r[0];
      prevPhase->idleTime.avg = r[1]/CkNumPes();
      prevPhase->idleTime.max = r[2];
      prevPhase->idleTime.print();
      CkPrintf("Stored idle time min=%lf avg=%lf max=%lf in prevPhase=%p\n", prevPhase->idleTime.min, prevPhase->idleTime.avg, prevPhase->idleTime.max, prevPhase);
    } else {
      CkPrintf("There is no previous phase to store the idle time measurements\n");
    }
    
    alreadyRequestedIdleTime = false;
    checkForShutdown();
    delete msg;
  }






  void controlPointManager::requestAll(CkCallback cb){
    CkAssert(enableCPTracing);

    TraceControlPoints *t = localControlPointTracingInstance();

    double data[ALL_REDUCTION_SIZE];

    double *idle = data;
    double *over = data+3;
    double *mem = data+6;
    double *msgBytes = data+7;  
    double *grainsize = data+11;  

    const double i = t->idleRatio();
    idle[0] = i;
    idle[1] = i;
    idle[2] = i;

    const double o = t->overheadRatio();
    over[0] = o;
    over[1] = o;
    over[2] = o;

    const double m = t->memoryUsageMB();
    mem[0] = m;

    msgBytes[0] = t->b2;
    msgBytes[1] = t->b3;
    msgBytes[2] = t->b2mlen;
    msgBytes[3] = t->b3mlen;

    grainsize[0] = t->grainSize();
    
    localControlPointTracingInstance()->resetAll();

    contribute(ALL_REDUCTION_SIZE*sizeof(double),data,allMeasuresReductionType, cb);
  }
  
  void controlPointManager::gatherAll(CkReductionMsg *msg){
    CkAssert(enableCPTracing);

    CkAssert(msg->getSize()==ALL_REDUCTION_SIZE*sizeof(double));
    int size=msg->getSize() / sizeof(double);
    double *data=(double *) msg->getData();
        
    double *idle = data;
    double *over = data+3;
    double *mem = data+6;
    double *msgBytes = data+7;
    double *granularity = data+11;


    //    std::string b = allData.toString();

    instrumentedPhase* prevPhase = previousPhaseData();
    if(prevPhase != NULL){
      prevPhase->idleTime.min = idle[0];
      prevPhase->idleTime.avg = idle[1]/CkNumPes();
      prevPhase->idleTime.max = idle[2];
      
      prevPhase->overheadTime.min = over[0];
      prevPhase->overheadTime.avg = over[1]/CkNumPes();
      prevPhase->overheadTime.max = over[2];
      
      prevPhase->memoryUsageMB = mem[0];
      
      CkPrintf("Stored idle time min=%lf avg=%lf max=%lf  mem=%lf in prevPhase=%p\n", (double)prevPhase->idleTime.min, prevPhase->idleTime.avg, prevPhase->idleTime.max, (double)prevPhase->memoryUsageMB, prevPhase);

      double bytesPerInvoke2 = msgBytes[2] / msgBytes[0];
      double bytesPerInvoke3 = msgBytes[3] / msgBytes[1];

      prevPhase->grainSize = (granularity[0] / (double)CkNumPes());

      CkPrintf("Bytes Per Invokation 2 = %f\n", bytesPerInvoke2);
      CkPrintf("Bytes Per Invokation 3 = %f\n", bytesPerInvoke3);

      CkPrintf("Bytes Per us of work 2 = %f\n", bytesPerInvoke2/prevPhase->grainSize);
      CkPrintf("Bytes Per us of work 3 = %f\n", bytesPerInvoke3/prevPhase->grainSize);

      if(bytesPerInvoke2 > bytesPerInvoke3)
    prevPhase->bytesPerInvoke = bytesPerInvoke2;
      else
    prevPhase->bytesPerInvoke = bytesPerInvoke3;

      //      prevPhase->print();
      // CkPrintf("prevPhase=%p number of timings=%d\n", prevPhase, prevPhase->times.size() );

      //      std::string a = allData.toString();

      //CkPrintf("Before:\n%s\nAfter:\n%s\n\n", b.c_str(), a.c_str());
      
    } else {
      CkPrintf("There is no previous phase to store measurements\n");
    }
    
    alreadyRequestedAll = false;
    checkForShutdown();
    delete msg;
  }




  void controlPointManager::checkForShutdown(){
    if( exitWhenReady && !alreadyRequestedAll && !alreadyRequestedMemoryUsage && !alreadyRequestedIdleTime && CkMyPe()==0){
      doExitNow();
    }
  }


  void controlPointManager::exitIfReady(){
     if( !alreadyRequestedMemoryUsage && !alreadyRequestedAll && !alreadyRequestedIdleTime && CkMyPe()==0){
       //  CkPrintf("controlPointManager::exitIfReady exiting immediately\n");
       doExitNow();
     } else {
       // CkPrintf("controlPointManager::exitIfReady Delaying exiting\n");
       exitWhenReady = true;
     }
  }



  void controlPointManager::doExitNow(){
          _TRACE_BEGIN_EXECUTE_DETAILED(-1, -1, _threadEP,CkMyPe(), 0, NULL, this);
      writeOutputToDisk();
    // CkPrintf("[%d] Control point manager calling CkContinueExit()\n", CkMyPe());
    CkContinueExit();
  }

  void controlPointManager::writeOutputToDisk(){
      if(writeDataFileAtShutdown){
          controlPointManagerProxy.ckLocalBranch()->writeDataFile();
      }
  }


  void controlPointManager::requestMemoryUsage(CkCallback cb){
    int m = CmiMaxMemoryUsage() / 1024 / 1024;
    CmiResetMaxMemory();
    //    CkPrintf("PE %d Memory Usage is %d MB\n",CkMyPe(), m);
    contribute(sizeof(int),&m,CkReduction::max_int, cb);
  }

  void controlPointManager::gatherMemoryUsage(CkReductionMsg *msg){
    int size=msg->getSize() / sizeof(int);
    CkAssert(size==1);
    int *m=(int *) msg->getData();

    CkPrintf("[%d] Max Memory Usage for all processors is %d MB\n", CkMyPe(), m[0]);
    
    instrumentedPhase* prevPhase = previousPhaseData();
    if(prevPhase != NULL){
      prevPhase->memoryUsageMB = m[0];
    } else {
      CkPrintf("No place to store memory usage");
    }

    alreadyRequestedMemoryUsage = false;
    checkForShutdown();
    delete msg;
  }


//   /// Inform the control point framework that a named control point affects the priorities of some array  
//   void controlPointManager::associatePriorityArray(const char *name, int groupIdx){
//     CkPrintf("Associating control point \"%s\" affects priority of array id=%d\n", name, groupIdx );
    
//     if(affectsPrioritiesArray.count(std::string(name)) > 0 ) {
//       affectsPrioritiesArray[std::string(name)].insert(groupIdx);
//     } else {
//       std::set<int> s;
//       s.insert(groupIdx);
//       affectsPrioritiesArray[std::string(name)] = s;
//     }
    
// #if DEBUGPRINT   
//     std::map<std::string, std::set<int> >::iterator f;
//     for(f=affectsPrioritiesArray.begin(); f!=affectsPrioritiesArray.end();++f){
//       std::string name = f->first;
//       std::set<int> &vals = f->second;
//       cout << "Control point " << name << " affects arrays: ";
//       std::set<int>::iterator i;
//       for(i=vals.begin(); i!=vals.end();++i){
//  cout << *i << " ";
//       }
//       cout << endl;
//     }
// #endif
    
//   }
  
//   /// Inform the control point framework that a named control point affects the priority of some entry method
//   void controlPointManager::associatePriorityEntry(const char *name, int idx){
//     CkPrintf("Associating control point \"%s\" with EP id=%d\n", name, idx);

//       if(affectsPrioritiesEP.count(std::string(name)) > 0 ) {
//       affectsPrioritiesEP[std::string(name)].insert(idx);
//     } else {
//       std::set<int> s;
//       s.insert(idx);
//       affectsPrioritiesEP[std::string(name)] = s;
//     }
    
// #if DEBUGPRINT
//     std::map<std::string, std::set<int> >::iterator f;
//     for(f=affectsPrioritiesEP.begin(); f!=affectsPrioritiesEP.end();++f){
//       std::string name = f->first;
//       std::set<int> &vals = f->second;
//       cout << "Control point " << name << " affects EP: ";
//       std::set<int>::iterator i;
//       for(i=vals.begin(); i!=vals.end();++i){
//  cout << *i << " ";
//       }
//       cout << endl;
//     }
// #endif


//   }
  


void gotoNextPhase(){
  controlPointManagerProxy.ckLocalBranch()->gotoNextPhase();
}

FLINKAGE void FTN_NAME(GOTONEXTPHASE,gotonextphase)()
{
  gotoNextPhase();
}


class controlPointMain : public CBase_controlPointMain {
public:
  controlPointMain(CkArgMsg* args){
#if OLDRANDSEED
    struct timeval tp;
    gettimeofday(& tp, NULL);
    random_seed = (int)tp.tv_usec ^ (int)tp.tv_sec;
#else
    double time = CmiWallTimer();
    int sec = (int)time;
    int usec = (int)((time - (double)sec)*1000000.0);
    random_seed =  sec ^ usec;
#endif
    
    
    double period, periodms;
    bool haveSamplePeriod = CmiGetArgDoubleDesc(args->argv,"+CPSamplePeriod", &period,"The time between Control Point Framework samples (in seconds)");
    bool haveSamplePeriodMs = CmiGetArgDoubleDesc(args->argv,"+CPSamplePeriodMs", &periodms,"The time between Control Point Framework samples (in milliseconds)");
    if(haveSamplePeriod){
      CkPrintf("controlPointSamplePeriod = %lf sec\n", period);
      controlPointSamplePeriod =  (int)(period * 1000); 
    } else if(haveSamplePeriodMs){
      CkPrintf("controlPointSamplePeriodMs = %lf ms\n", periodms);
      controlPointSamplePeriod = periodms; 
    } else {
      controlPointSamplePeriod =  DEFAULT_CONTROL_POINT_SAMPLE_PERIOD;
    }

  
    
    whichTuningScheme = RandomSelection;


    if( CmiGetArgFlagDesc(args->argv,"+CPSchemeRandom","Randomly Select Control Point Values") ){
      whichTuningScheme = RandomSelection;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPExhaustiveSearch","Exhaustive Search of Control Point Values") ){
      whichTuningScheme = ExhaustiveSearch;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPAlwaysUseDefaults","Always Use The Provided Default Control Point Values") ){
      whichTuningScheme = AlwaysDefaults;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPSimulAnneal","Simulated Annealing Search of Control Point Values") ){
      whichTuningScheme = SimulatedAnnealing;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPCriticalPathPrio","Use Critical Path to adapt Control Point Values") ){
      whichTuningScheme = CriticalPathAutoPrioritization;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPBestKnown","Use BestKnown Timing for Control Point Values") ){
      whichTuningScheme = UseBestKnownTiming;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPSteering","Use Steering to adjust Control Point Values") ){
      whichTuningScheme = UseSteering;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPMemoryAware", "Adjust control points to approach available memory") ){
      whichTuningScheme = MemoryAware;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPSimplex", "Nelder-Mead Simplex Algorithm") ){
      whichTuningScheme = Simplex;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPDivideConquer", "A divide and conquer program specific steering scheme") ){
      whichTuningScheme = DivideAndConquer;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPLDBPeriod", "Adjust the load balancing period (Constant Predictor)") ){
      whichTuningScheme = LDBPeriod;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPLDBPeriodLinear", "Adjust the load balancing period (Linear Predictor)") ){
      whichTuningScheme = LDBPeriodLinear;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPLDBPeriodQuadratic", "Adjust the load balancing period (Quadratic Predictor)") ){
      whichTuningScheme = LDBPeriodQuadratic;
    } else if ( CmiGetArgFlagDesc(args->argv,"+CPLDBPeriodOptimal", "Adjust the load balancing period (Optimal Predictor)") ){
      whichTuningScheme = LDBPeriodOptimal;
    }

    char *defValStr = NULL;
    if( CmiGetArgStringDesc(args->argv, "+CPDefaultValues", &defValStr, "Specify the default control point values used for the first couple phases") ){
      CkPrintf("You specified default value string: %s\n", defValStr);
      
      // Parse the string, looking for commas
     

      char *tok = strtok(defValStr, ",");
      while (tok) {
    // split on the equal sign
    int len = strlen(tok);
    char *eqsign = strchr(tok, '=');
    if(eqsign != NULL && eqsign != tok){
      *eqsign = '\0';
      char *cpname = tok;
      std::string cpName(tok);
      char *cpDefVal = eqsign+1;      
      int v=-1;
      if(sscanf(cpDefVal, "%d", &v) == 1){
        CkPrintf("Command Line Argument Specifies that Control Point \"%s\" defaults to %d\n", cpname, v);
        CkAssert(CkMyPe() == 0); // can only access defaultControlPointValues on PE 0
        defaultControlPointValues[cpName] = v;
      }
    }
    tok = strtok(NULL, ",");
      }

    }

    shouldGatherAll = false;
    shouldGatherMemoryUsage = false;
    shouldGatherUtilization = false;
    
    if ( CmiGetArgFlagDesc(args->argv,"+CPGatherAll","Gather all types of measurements for each phase") ){
      shouldGatherAll = true;
    } else {
      if ( CmiGetArgFlagDesc(args->argv,"+CPGatherMemoryUsage","Gather memory usage after each phase") ){
    shouldGatherMemoryUsage = true;
      }
      if ( CmiGetArgFlagDesc(args->argv,"+CPGatherUtilization","Gather utilization & Idle time after each phase") ){
    shouldGatherUtilization = true;
      }
    }
    
    writeDataFileAtShutdown = false;   
    if( CmiGetArgFlagDesc(args->argv,"+CPSaveData","Save Control Point timings & configurations at completion") ){
      writeDataFileAtShutdown = true;
    }

    shouldFilterOutputData = true;
    if( CmiGetArgFlagDesc(args->argv,"+CPNoFilterData","Don't filter phases from output data") ){
      shouldFilterOutputData = false;
    }


   loadDataFileAtStartup = false;   
    if( CmiGetArgFlagDesc(args->argv,"+CPLoadData","Load Control Point timings & configurations at startup") ){
      loadDataFileAtStartup = true;
    }

    char *cpdatafile;
    if( CmiGetArgStringDesc(args->argv, "+CPDataFilename", &cpdatafile, "Specify control point data file to save/load") ){
      sprintf(CPDataFilename, "%s", cpdatafile);
    } else {
      sprintf(CPDataFilename, "controlPointData.txt");
    }


    controlPointManagerProxy = CProxy_controlPointManager::ckNew();

    delete args;
  }
  ~controlPointMain(){}
};

void registerCPChangeCallback(CkCallback cb, bool frameworkShouldAdvancePhase){
  CkAssert(CkMyPe() == 0);
  CkPrintf("Application has registered a control point change callback\n");
  controlPointManagerProxy.ckLocalBranch()->setCPCallback(cb, frameworkShouldAdvancePhase);
}

void setFrameworkAdvancePhase(bool frameworkShouldAdvancePhase){
  if(CkMyPe() == 0){
    CkPrintf("Application has specified that framework should %sadvance phase\n", frameworkShouldAdvancePhase?"":"not ");
    controlPointManagerProxy.ckLocalBranch()->setFrameworkAdvancePhase(frameworkShouldAdvancePhase);
  }
}

void registerControlPointTiming(double time){
  CkAssert(CkMyPe() == 0);
#if DEBUGPRINT>0
  CkPrintf("Program registering its own timing with registerControlPointTiming(time=%lf)\n", time);
#endif
  controlPointManagerProxy.ckLocalBranch()->setTiming(time);
}

void controlPointTimingStamp() {
  CkAssert(CkMyPe() == 0);
#if DEBUGPRINT>0
  CkPrintf("Program registering its own timing with controlPointTimingStamp()\n", time);
#endif
  
  static double prev_time = 0.0;
  double t = CmiWallTimer();
  double duration = t - prev_time;
  prev_time = t;
    
  controlPointManagerProxy.ckLocalBranch()->setTiming(duration);
}

FLINKAGE void FTN_NAME(CONTROLPOINTTIMINGSTAMP,controlpointtimingstamp)()
{
  controlPointTimingStamp();
}


FLINKAGE void FTN_NAME(SETFRAMEWORKADVANCEPHASEF,setframeworkadvancephasef)(CMK_TYPEDEF_INT4 *value)
{
  setFrameworkAdvancePhase(*value);
}




extern "C" void controlPointShutdown(){
  if(CkMyPe() == 0){

    if (!controlPointManagerProxy.ckGetGroupID().isZero()) {
      // wait for gathering of idle time & memory usage to complete
      controlPointManagerProxy.ckLocalBranch()->exitIfReady();
    } else {
      CkContinueExit();
    }
  }
}

void controlPointInitNode(){
  registerExitFn(controlPointShutdown);
}

static void periodicProcessControlPoints(void* ptr, double currWallTime){
#ifdef DEBUGPRINT
  CkPrintf("[%d] periodicProcessControlPoints()\n", CkMyPe());
#endif
  controlPointManagerProxy.ckLocalBranch()->processControlPoints();
  CcdCallFnAfterOnPE((CcdVoidFn)periodicProcessControlPoints, (void*)NULL, controlPointSamplePeriod, CkMyPe());
}





void controlPointManager::generatePlan() {
  const double startGenerateTime = CmiWallTimer();
  const int phase_id = this->phase_id;
  const int effective_phase = allData.phases.size();

  // Only generate a plan once per phase
  if(generatedPlanForStep == phase_id) 
    return;
  generatedPlanForStep = phase_id;
 
  CkPrintf("Generating Plan for phase %d\n", phase_id); 
  printTuningScheme();

  // By default lets put the previous phase data into newControlPoints
  instrumentedPhase *prevPhase = previousPhaseData();
  for(std::map<std::string, int >::const_iterator cpsIter=prevPhase->controlPoints.begin(); cpsIter != prevPhase->controlPoints.end(); ++cpsIter){
      newControlPoints[cpsIter->first] = cpsIter->second;
  }


  if( whichTuningScheme == RandomSelection ){
    std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
    for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
      const std::string &name = cpsIter->first;
      const std::pair<int,int> &bounds = cpsIter->second;
      const int lb = bounds.first;
      const int ub = bounds.second;
      newControlPoints[name] = lb + randInt(ub-lb+1, name.c_str(), phase_id);
    }
  } else if( whichTuningScheme == CriticalPathAutoPrioritization) {
    // -----------------------------------------------------------
    //  USE CRITICAL PATH TO ADJUST PRIORITIES
    //   
    // This scheme will return the median value for the range 
    // early on, and then will transition over to the new control points
    // determined by the critical path adapting code

    // This won't work until the periodic function is fixed up a bit

    std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
    for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
      const std::string &name = cpsIter->first;
      const std::pair<int,int> &bounds = cpsIter->second;
      const int lb = bounds.first;
      const int ub = bounds.second;
      newControlPoints[name] =  (lb+ub)/2;
    }

  } else if ( whichTuningScheme == MemoryAware ) {

    // -----------------------------------------------------------
    //  STEERING BASED ON MEMORY USAGE

    instrumentedPhase *twoAgoPhase = twoAgoPhaseData();
    instrumentedPhase *prevPhase = previousPhaseData();
 
    if(phase_id%2 == 0){
      CkPrintf("Steering (memory based) based on 2 phases ago:\n");
      twoAgoPhase->print();
      CkPrintf("\n");
      fflush(stdout);
      
      // See if memory usage is low:
      double memUsage = twoAgoPhase->memoryUsageMB;
      CkPrintf("Steering (memory based) encountered memory usage of (%f MB)\n", memUsage);
      fflush(stdout);
      if(memUsage < 1100.0 && memUsage > 0.0){ // Kraken has about 16GB and 12 cores per node
    CkPrintf("Steering (memory based) encountered low memory usage (%f) < 1200 \n", memUsage);
    CkPrintf("Steering (memory based) controlPointSpace.size()=\n", controlPointSpace.size());
    
    // Initialize plan to be the values from two phases ago (later we'll adjust this)
    newControlPoints = twoAgoPhase->controlPoints;


    CkPrintf("Steering (memory based) initialized plan\n");
    fflush(stdout);

    // look for a possible control point knob to turn
    std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > > &possibleCPsToTune = CkpvAccess(cp_effects)["MemoryConsumption"];
    
    // FIXME: assume for now that we just have one control point with the effect, and one direction to turn it
    bool found = false;
    std::string cpName;
    std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > *info;
    std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > >::iterator iter;
    for(iter = possibleCPsToTune.begin(); iter != possibleCPsToTune.end(); iter++){
      cpName = iter->first;
      info = &iter->second;
      found = true;
      break;
    }

    // Adapt the control point value
    if(found){
      CkPrintf("Steering found knob to turn that should increase memory consumption\n");
      fflush(stdout);
      const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
      const int maxValue = controlPointSpace[cpName].second;
      
      if(twoAgoValue+1 <= maxValue){
        newControlPoints[cpName] = twoAgoValue+1; // increase from two phases back
      }
    }
    
      }
    }

    CkPrintf("Steering (memory based) done for this phase\n");
    fflush(stdout);

  } else if ( whichTuningScheme == UseBestKnownTiming ) {

    // -----------------------------------------------------------
    //  USE BEST KNOWN TIME  
    static bool firstTime = true;
    if(firstTime){
      firstTime = false;
      instrumentedPhase *best = allData.findBest(); 
      CkPrintf("Best known phase is:\n");
      best->print();
      CkPrintf("\n");
      
      std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
      for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter) {
    const std::string &name = cpsIter->first;
    newControlPoints[name] =  best->controlPoints[name];
      }
    }
  } else if( whichTuningScheme == LDBPeriod) {
    // Assume this is used in this manner:
    //  1) go to next phase
    //  2) request control point
    //  3) load balancing
    //  4) computation
    
    
    instrumentedPhase *twoAgoPhase = twoAgoPhaseData();
    instrumentedPhase *prevPhase = previousPhaseData();
    
    
    const std::vector<double> &times = twoAgoPhase->times;
    const int oldNumTimings = times.size();


    const std::vector<double> &timesNew = prevPhase->times;
    const int newNumTimings = timesNew.size();


    if(oldNumTimings > 4 && newNumTimings > 4){
      
      // Build model of execution time based on two phases ago
      // Compute the average times for each 1/3 of the steps, except for the 2 first steps where load balancing occurs
      
      double oldSum = 0;
      
      for(int i=2; i<oldNumTimings; i++){
    oldSum += times[i];
      }
      
      const double oldAvg = oldSum / (oldNumTimings-2);
      
      
      
      
      // Computed as an integral from 0.5 to the number of bins of the same size as two ago phase + 0.5
      const double expectedTotalTime = oldAvg * newNumTimings;
      
      
      // Measure actual time
      double newSum = 0.0;
      for(int i=2; i<newNumTimings; ++i){
    newSum += timesNew[i];
      }
      
      const double newAvg = newSum / (newNumTimings-2);
      const double newTotalTimeExcludingLBSteps = newAvg * ((double)newNumTimings); // excluding the load balancing abnormal steps
      
      const double benefit = expectedTotalTime - newTotalTimeExcludingLBSteps;
      
      // Determine load balance cost
      const double lbcost = timesNew[0] + timesNew[1] - 2.0*newAvg;
      
      const double benefitAfterLB = benefit - lbcost;
    
    
      // Determine whether LB cost outweights the estimated benefit
      CkPrintf("Constant Model: lbcost = %f, expected = %f, actual = %f\n", lbcost, expectedTotalTime, newTotalTimeExcludingLBSteps);
    
    
      std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
      for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
    const std::string &name = cpsIter->first;
    const std::pair<int,int> &bounds = cpsIter->second;
    const int lb = bounds.first;
    const int ub = bounds.second;
      
    if(benefitAfterLB > 0){
      CkPrintf("Constant Model: Beneficial LB\n");
      int newval = newControlPoints[name] / 2;
      if(newval > lb)
        newControlPoints[name] = newval;
      else 
        newControlPoints[name] = lb;
    } else {
      CkPrintf("Constant Model: Detrimental LB\n");
      int newval = newControlPoints[name] * 2;
      if(newval < ub)
        newControlPoints[name] = newval;
      else
        newControlPoints[name] = ub;
    }
      }
    }
    
    
  }  else if( whichTuningScheme == LDBPeriodLinear) {
    // Assume this is used in this manner:
    //  1) go to next phase
    //  2) request control point
    //  3) load balancing
    //  4) computation


    instrumentedPhase *twoAgoPhase = twoAgoPhaseData();
    instrumentedPhase *prevPhase = previousPhaseData();
    
    const std::vector<double> &times = twoAgoPhase->times;
    const int oldNumTimings = times.size();

    const std::vector<double> &timesNew = prevPhase->times;
    const int newNumTimings = timesNew.size();
    

    if(oldNumTimings > 4 && newNumTimings > 4){

      // Build model of execution time based on two phases ago
      // Compute the average times for each 1/3 of the steps, except for the 2 first steps where load balancing occurs
      const int b1 = 2 + (oldNumTimings-2)/2;
      double s1 = 0;
      double s2 = 0;
    
      const double ldbStepsTime = times[0] + times[1];
    
      for(int i=2; i<b1; i++){
    s1 += times[i];
      }
      for(int i=b1; i<oldNumTimings; i++){
    s2 += times[i];
      }
      
      
      // Compute the estimated time for the last phase's data
    
      const double a1 = s1 / (double)(b1-2);
      const double a2 = s2 / (double)(oldNumTimings-b1);
      const double aold = (a1+a2)/2.0;

      const double expectedTotalTime = newNumTimings*(aold+(oldNumTimings+newNumTimings)*(a2-a1)/oldNumTimings);
        
    
      // Measure actual time
      double sum = 0.0;
      for(int i=2; i<newNumTimings; ++i){
    sum += timesNew[i];
      }

      const double avg = sum / ((double)(newNumTimings-2));
      const double actualTotalTime = avg * ((double)newNumTimings); // excluding the load balancing abnormal steps

      const double benefit = expectedTotalTime - actualTotalTime;

      // Determine load balance cost
      const double lbcost = timesNew[0] + timesNew[1] - 2.0*avg;

      const double benefitAfterLB = benefit - lbcost;

    
      // Determine whether LB cost outweights the estimated benefit
      CkPrintf("Linear Model: lbcost = %f, expected = %f, actual = %f\n", lbcost, expectedTotalTime, actualTotalTime);
    
    
    
      std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
      for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
    const std::string &name = cpsIter->first;
    const std::pair<int,int> &bounds = cpsIter->second;
    const int lb = bounds.first;
    const int ub = bounds.second;
      
    if(benefitAfterLB > 0){
      CkPrintf("Linear Model: Beneficial LB\n");
      int newval = newControlPoints[name] / 2;
      if(newval > lb)
        newControlPoints[name] = newval;
      else 
        newControlPoints[name] = lb;
    } else {
      CkPrintf("Linear Model: Detrimental LB\n");
      int newval = newControlPoints[name] * 2;
      if(newval < ub)
        newControlPoints[name] = newval;
      else 
        newControlPoints[name] = ub;
    }
      }
    }

  }

  else if( whichTuningScheme == LDBPeriodQuadratic) {
    // Assume this is used in this manner:
    //  1) go to next phase
    //  2) request control point
    //  3) load balancing
    //  4) computation


    instrumentedPhase *twoAgoPhase = twoAgoPhaseData();
    instrumentedPhase *prevPhase = previousPhaseData();
        
    const std::vector<double> &times = twoAgoPhase->times;
    const int oldNumTimings = times.size();

    const std::vector<double> &timesNew = prevPhase->times;
    const int newNumTimings = timesNew.size();

    
    if(oldNumTimings > 4 && newNumTimings > 4){


      // Build model of execution time based on two phases ago
      // Compute the average times for each 1/3 of the steps, except for the 2 first steps where load balancing occurs
      const int b1 = 2 + (oldNumTimings-2)/3;
      const int b2 = 2 + (2*(oldNumTimings-2))/3;
      double s1 = 0;
      double s2 = 0;
      double s3 = 0;

      const double ldbStepsTime = times[0] + times[1];
    
      for(int i=2; i<b1; i++){
    s1 += times[i];
      }
      for(int i=b1; i<b2; i++){
    s2 += times[i];
      }
      for(int i=b2; i<oldNumTimings; i++){
    s3 += times[i];
      }

    
      // Compute the estimated time for the last phase's data
    
      const double a1 = s1 / (double)(b1-2);
      const double a2 = s2 / (double)(b2-b1);
      const double a3 = s3 / (double)(oldNumTimings-b2);
    
      const double a = (a1-2.0*a2+a3)/2.0;
      const double b = (a1-4.0*a2+3.0*a3)/2.0;
      const double c = a3;
    
      // Computed as an integral from 0.5 to the number of bins of the same size as two ago phase + 0.5
      const double x1 = (double)newNumTimings / (double)oldNumTimings * 3.0 + 0.5;  // should be 3.5 if ratio is same
      const double x2 = 0.5;
      const double expectedTotalTime = a*x1*x1*x1/3.0 + b*x1*x1/2.0 + c*x1 - (a*x2*x2*x2/3.0 + b*x2*x2/2.0 + c*x2);
   
    
      // Measure actual time
      double sum = 0.0;
      for(int i=2; i<newNumTimings; ++i){
    sum += timesNew[i];
      }

      const double avg = sum / ((double)(newNumTimings-2));
      const double actualTotalTime = avg * ((double)newNumTimings); // excluding the load balancing abnormal steps

      const double benefit = expectedTotalTime - actualTotalTime;

      // Determine load balance cost
      const double lbcost = timesNew[0] + timesNew[1] - 2.0*avg;

      const double benefitAfterLB = benefit - lbcost;

    
      // Determine whether LB cost outweights the estimated benefit
      CkPrintf("Quadratic Model: lbcost = %f, expected = %f, actual = %f, x1=%f, a1=%f, a2=%f, a3=%f, b1=%d, b2=%d, a=%f, b=%f, c=%f\n", lbcost, expectedTotalTime, actualTotalTime, x1, a1, a2, a3, b1, b2, a, b, c);
    
    
    
      std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
      for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
    const std::string &name = cpsIter->first;
    const std::pair<int,int> &bounds = cpsIter->second;
    const int lb = bounds.first;
    const int ub = bounds.second;
      
    if(benefitAfterLB > 0){
      CkPrintf("QuadraticModel: Beneficial LB\n");
      int newval = newControlPoints[name] / 2;
      if(newval > lb)
        newControlPoints[name] = newval;
      else 
        newControlPoints[name] = lb;
    } else {
      CkPrintf("QuadraticModel: Detrimental LB\n");
      int newval = newControlPoints[name] * 2;
      if(newval < ub)
        newControlPoints[name] = newval;
      else 
        newControlPoints[name] = ub;
    }
      
      }
    }
    

  }  else if( whichTuningScheme == LDBPeriodOptimal) {
    // Assume this is used in this manner:
    //  1) go to next phase
    //  2) request control point
    //  3) load balancing
    //  4) computation



    instrumentedPhase *prevPhase = previousPhaseData();
    
    const std::vector<double> &times = prevPhase->times;
    const int numTimings = times.size();
    
    if( numTimings > 4){

      const int b1 = 2 + (numTimings-2)/2;
      double s1 = 0;
      double s2 = 0;
    
    
      for(int i=2; i<b1; i++){
    s1 += times[i];
      }
      for(int i=b1; i<numTimings; i++){
    s2 += times[i];
      }
      
    
      const double a1 = s1 / (double)(b1-2);
      const double a2 = s2 / (double)(numTimings-b1);
      const double avg = (a1+a1) / 2.0;

      const double m = (a2-a1)/((double)(numTimings-2)/2.0); // An approximation of the slope of the execution times    

      const double ldbStepsTime = times[0] + times[1];
      const double lbcost = ldbStepsTime - 2.0*avg; // An approximation of the 
      

      int newval = roundDouble(sqrt(2.0*lbcost/m));
      
      // We don't really know what to do if m<=0, so we'll just double the period
      if(m<=0)
    newval = 2*numTimings;     
      
      CkPrintf("Optimal Model (double when negative): lbcost = %f, m = %f, new ldbperiod should be %d\n", lbcost, m, newval);    
    
    
      std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
      for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
    // TODO: lookup only control points that are relevant instead of all of them
    const std::string &name = cpsIter->first;
    const std::pair<int,int> &bounds = cpsIter->second;
    const int lb = bounds.first;
    const int ub = bounds.second;
    
    if(newval < lb){
      newControlPoints[name] = lb;
    } else if(newval > ub){
      newControlPoints[name] = ub;
    } else {
      newControlPoints[name] = newval;
    } 
    
      }
      
      
    }
    
 

  } else if ( whichTuningScheme == UseSteering ) {
      // -----------------------------------------------------------
      //  STEERING BASED ON KNOWLEDGE

      // after 3 phases (and only on even steps), do steering performance. Otherwise, just use previous phase's configuration
      // plans are only generated after 3 phases

      instrumentedPhase *twoAgoPhase = twoAgoPhaseData();
      instrumentedPhase *prevPhase = previousPhaseData();

      if(phase_id%4 == 0){
          CkPrintf("Steering based on 2 phases ago:\n");
          twoAgoPhase->print();
          CkPrintf("\n");
          fflush(stdout);

          std::vector<std::map<std::string,int> > possibleNextStepPlans;


          // ========================================= Concurrency =============================================
          // See if idle time is high:
          {
              double idleTime = twoAgoPhase->idleTime.avg;
              CkPrintf("Steering encountered idle time (%f)\n", idleTime);
              fflush(stdout);
              if(idleTime > 0.10){
                  CkPrintf("Steering encountered high idle time(%f) > 10%%\n", idleTime);
                  CkPrintf("Steering controlPointSpace.size()=\n", controlPointSpace.size());

                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > > &possibleCPsToTune = CkpvAccess(cp_effects)["Concurrency"];

                  bool found = false;
                  std::string cpName;
                  std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > *info;
                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > >::iterator iter;
                  for(iter = possibleCPsToTune.begin(); iter != possibleCPsToTune.end(); iter++){
                      cpName = iter->first;
                      info = &iter->second;

                      // Initialize a new plan based on two phases ago
                      std::map<std::string,int> aNewPlan = twoAgoPhase->controlPoints;

                      CkPrintf("Steering found knob to turn\n");
                      fflush(stdout);

                      if(info->first == ControlPoint::EFF_INC){
                          const int maxValue = controlPointSpace[cpName].second;
                          const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                          if(twoAgoValue+1 <= maxValue){
                              aNewPlan[cpName] = twoAgoValue+1; // increase from two phases back
                          }
                      } else {
                          const int minValue = controlPointSpace[cpName].second;
                          const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                          if(twoAgoValue-1 >= minValue){
                              aNewPlan[cpName] = twoAgoValue-1; // decrease from two phases back
                          }
                      }

                      possibleNextStepPlans.push_back(aNewPlan);

                  }
              }
          }

          // ========================================= Grain Size =============================================
          // If the grain size is too small, there may be tons of messages and overhead time associated with scheduling
          {
              double overheadTime = twoAgoPhase->overheadTime.avg;
              CkPrintf("Steering encountered overhead time (%f)\n", overheadTime);
              fflush(stdout);
              if(overheadTime > 0.10){
                  CkPrintf("Steering encountered high overhead time(%f) > 10%%\n", overheadTime);
                  CkPrintf("Steering controlPointSpace.size()=\n", controlPointSpace.size());

                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > > &possibleCPsToTune = CkpvAccess(cp_effects)["GrainSize"];

                  bool found = false;
                  std::string cpName;
                  std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > *info;
                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > >::iterator iter;
                  for(iter = possibleCPsToTune.begin(); iter != possibleCPsToTune.end(); iter++){
                      cpName = iter->first;
                      info = &iter->second;

                      // Initialize a new plan based on two phases ago
                      std::map<std::string,int> aNewPlan = twoAgoPhase->controlPoints;



                      CkPrintf("Steering found knob to turn\n");
                      fflush(stdout);

                      if(info->first == ControlPoint::EFF_INC){
                          const int maxValue = controlPointSpace[cpName].second;
                          const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                          if(twoAgoValue+1 <= maxValue){
                              aNewPlan[cpName] = twoAgoValue+1; // increase from two phases back
                          }
                      } else {
                          const int minValue = controlPointSpace[cpName].second;
                          const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                          if(twoAgoValue-1 >= minValue){
                              aNewPlan[cpName] = twoAgoValue-1; // decrease from two phases back
                          }
                      }

                      possibleNextStepPlans.push_back(aNewPlan);

                  }

              }
          }
          // ========================================= GPU Offload =============================================
          // If the grain size is too small, there may be tons of messages and overhead time associated with scheduling
          {
              double idleTime = twoAgoPhase->idleTime.avg;
              CkPrintf("Steering encountered idle time (%f)\n", idleTime);
              fflush(stdout);
              if(idleTime > 0.10){
                  CkPrintf("Steering encountered high idle time(%f) > 10%%\n", idleTime);
                  CkPrintf("Steering controlPointSpace.size()=\n", controlPointSpace.size());

                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > > &possibleCPsToTune = CkpvAccess(cp_effects)["GPUOffloadedWork"];

                  bool found = false;
                  std::string cpName;
                  std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > *info;
                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > >::iterator iter;
                  for(iter = possibleCPsToTune.begin(); iter != possibleCPsToTune.end(); iter++){
                      cpName = iter->first;
                      info = &iter->second;

                      // Initialize a new plan based on two phases ago
                      std::map<std::string,int> aNewPlan = twoAgoPhase->controlPoints;


                      CkPrintf("Steering found knob to turn\n");
                      fflush(stdout);

                      if(info->first == ControlPoint::EFF_DEC){
                          const int maxValue = controlPointSpace[cpName].second;
                          const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                          if(twoAgoValue+1 <= maxValue){
                              aNewPlan[cpName] = twoAgoValue+1; // increase from two phases back
                          }
                      } else {
                          const int minValue = controlPointSpace[cpName].second;
                          const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                          if(twoAgoValue-1 >= minValue){
                              aNewPlan[cpName] = twoAgoValue-1; // decrease from two phases back
                          }
                      }

                      possibleNextStepPlans.push_back(aNewPlan);

                  }

              }
          }

          // ========================================= Done =============================================


          if(possibleNextStepPlans.size() > 0){
              newControlPoints = possibleNextStepPlans[0];
          }


          CkPrintf("Steering done for this phase\n");
          fflush(stdout);

      }  else {
          // This is not a phase to do steering, so stick with previously used values (one phase ago)
          CkPrintf("not a phase to do steering, so stick with previously planned values (one phase ago)\n");
          fflush(stdout);
      }



  } else if( whichTuningScheme == SimulatedAnnealing ) {

    // -----------------------------------------------------------
    //  SIMULATED ANNEALING
    //  Simulated Annealing style hill climbing method
    //
    //  Find the best search space configuration, and try something
    //  nearby it, with a radius decreasing as phases increase

    CkPrintf("Finding best phase\n");
    instrumentedPhase *bestPhase = allData.findBest();  
    CkPrintf("best found:\n"); 
    bestPhase->print(); 
    CkPrintf("\n"); 
    

    const int convergeByPhase = 100;
    // Determine from 0.0 to 1.0 how far along we are
    const double progress = (double) min(effective_phase, convergeByPhase) / (double)convergeByPhase;

    std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
    for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
      const std::string &name = cpsIter->first;
      const std::pair<int,int> &bounds = cpsIter->second;
      const int minValue = bounds.first;
      const int maxValue = bounds.second;
      
      const int before = bestPhase->controlPoints[name];   
  
      const int range = (int)((maxValue-minValue+1)*(1.0-progress));

      int high = min(before+range, maxValue);
      int low = max(before-range, minValue);
      
      newControlPoints[name] = low;
      if(high-low > 0){
    newControlPoints[name] += randInt(high-low, name.c_str(), phase_id); 
      } 
      
    }

  } else if ( whichTuningScheme == DivideAndConquer ) {

      // -----------------------------------------------------------
      //  STEERING FOR Divide & Conquer Programs
      //  This scheme uses no timing information. It just tries to converge to the point where idle time = overhead time.
      //  For a Fibonacci example, this appears to be a good heurstic for finding the best performing program.
      //  The scheme can be applied within a single program tree computation, if the tree is being traversed depth first.

      // after 3 phases (and only on even steps), do steering performance. Otherwise, just use previous phase's configuration
      // plans are only generated after 3 phases

      instrumentedPhase *twoAgoPhase = twoAgoPhaseData();
      instrumentedPhase *prevPhase = previousPhaseData();

      if(phase_id%4 == 0){
          CkPrintf("Divide & Conquer Steering based on 2 phases ago:\n");
          twoAgoPhase->print();
          CkPrintf("\n");
          fflush(stdout);

          std::vector<std::map<std::string,int> > possibleNextStepPlans;


          // ========================================= Concurrency =============================================
          // See if idle time is high:
          {
              double idleTime = twoAgoPhase->idleTime.avg;
              double overheadTime = twoAgoPhase->overheadTime.avg;


              CkPrintf("Divide & Conquer Steering encountered overhead time (%f) idle time (%f)\n",overheadTime, idleTime);
              fflush(stdout);
              if(idleTime+overheadTime > 0.10){
                  CkPrintf("Steering encountered high idle+overheadTime time(%f) > 10%%\n", idleTime+overheadTime);
                  CkPrintf("Steering controlPointSpace.size()=\n", controlPointSpace.size());

                  int direction = -1;
                  if (idleTime>overheadTime){
                      // We need to decrease the grain size, or increase the available concurrency
                      direction = 1;
                  }

                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > > &possibleCPsToTune = CkpvAccess(cp_effects)["Concurrency"];

                  bool found = false;
                  std::string cpName;
                  std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > *info;
                  std::map<std::string, std::pair<int, std::vector<ControlPoint::ControlPointAssociation> > >::iterator iter;
                  for(iter = possibleCPsToTune.begin(); iter != possibleCPsToTune.end(); iter++){
                      cpName = iter->first;
                      info = &iter->second;
                      
                      // Initialize a new plan based on two phases ago
                      std::map<std::string,int> aNewPlan = twoAgoPhase->controlPoints;

                      CkPrintf("Divide & Conquer Steering found knob to turn\n");
                      fflush(stdout);

                      int adjustByAmount = (int)(myAbs(idleTime-overheadTime)*5.0);
                      
                      if(info->first == ControlPoint::EFF_INC){
                        const int minValue = controlPointSpace[cpName].first;
                        const int maxValue = controlPointSpace[cpName].second;
                        const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                        const int newVal = closestInRange(twoAgoValue+adjustByAmount*direction, minValue, maxValue);                      
                        CkAssert(newVal <= maxValue && newVal >= minValue);
                        aNewPlan[cpName] = newVal;
                      } else {
                        const int minValue = controlPointSpace[cpName].first;
                        const int maxValue = controlPointSpace[cpName].second;
                        const int twoAgoValue =  twoAgoPhase->controlPoints[cpName];
                        const int newVal = closestInRange(twoAgoValue-adjustByAmount*direction, minValue, maxValue);
                        CkAssert(newVal <= maxValue && newVal >= minValue);
                        aNewPlan[cpName] = newVal;
                      }
                      
                      possibleNextStepPlans.push_back(aNewPlan);
                  }
              }
          }

          if(possibleNextStepPlans.size() > 0){
            CkPrintf("Divide & Conquer Steering found %d possible next phases, using first one\n", possibleNextStepPlans.size());
            newControlPoints = possibleNextStepPlans[0];
          } else {
            CkPrintf("Divide & Conquer Steering found no possible next phases\n");
          }
      }

  } else if( whichTuningScheme == Simplex ) {

      // -----------------------------------------------------------
      //  Nelder Mead Simplex Algorithm
      //
      //  A scheme that takes a simplex (n+1 points) and moves it
      //  toward the minimum, eventually converging there.

      s.adapt(controlPointSpace, newControlPoints, phase_id, allData);

  } else if( whichTuningScheme == ExhaustiveSearch ) {
    
    // -----------------------------------------------------------
    // EXHAUSTIVE SEARCH
   
    int numDimensions = controlPointSpace.size();
    CkAssert(numDimensions > 0);
  
    vector<int> lowerBounds(numDimensions);
    vector<int> upperBounds(numDimensions); 
  
    int d=0;
    std::map<std::string, pair<int,int> >::iterator iter;
    for(iter = controlPointSpace.begin(); iter != controlPointSpace.end(); iter++){
      //    CkPrintf("Examining dimension %d\n", d);
      lowerBounds[d] = iter->second.first;
      upperBounds[d] = iter->second.second;
      d++;
    }
   
    // get names for each dimension (control point)
    vector<std::string> names(numDimensions);
    d=0;
    for(std::map<std::string, pair<int,int> >::iterator niter=controlPointSpace.begin(); niter!=controlPointSpace.end(); niter++){
      names[d] = niter->first;
      d++;
    }
  
  
    // Create the first possible configuration
    vector<int> config = lowerBounds;
    config.push_back(0);
  
    // Increment until finding an unused configuration
    allData.cleanupNames(); // put -1 values in for any control points missing
    std::vector<instrumentedPhase*> &phases = allData.phases;     

    while(true){
    
      std::vector<instrumentedPhase*>::iterator piter; 
      bool testedConfiguration = false; 
      for(piter = phases.begin(); piter != phases.end(); piter++){ 
      
    // Test if the configuration matches this phase
    bool match = true;
    for(int j=0;j<numDimensions;j++){
      match &= (*piter)->controlPoints[names[j]] == config[j];
    }
      
    if(match){
      testedConfiguration = true; 
      break;
    } 
      } 
    
      if(testedConfiguration == false){ 
    break;
      } 
    
      // go to next configuration
      config[0] ++;
      // Propagate "carrys"
      for(int i=0;i<numDimensions;i++){
    if(config[i] > upperBounds[i]){
      config[i+1]++;
      config[i] = lowerBounds[i];
    }
      }
      // check if we have exhausted all possible values
      if(config[numDimensions] > 0){
    break;
      }
       
    }
  
    if(config[numDimensions] > 0){
      for(int i=0;i<numDimensions;i++){ 
    config[i]= (lowerBounds[i]+upperBounds[i]) / 2;
      }
    }

    // results are now in config[i]
    
    for(int i=0; i<numDimensions; i++){
      newControlPoints[names[i]] = config[i];
      CkPrintf("Exhaustive search chose:   %s   -> %d\n", names[i].c_str(), config[i]);
    }


  } else {
    CkAbort("Some Control Point tuning strategy must be enabled.\n");
  }


  const double endGenerateTime = CmiWallTimer();
  
  CkPrintf("Time to generate next control point configuration(s): %f sec\n", (endGenerateTime - startGenerateTime) );
  
}





int controlPoint(const char *name, int lb, int ub){
  instrumentedPhase *thisPhaseData = controlPointManagerProxy.ckLocalBranch()->currentPhaseData();
  const int phase_id = controlPointManagerProxy.ckLocalBranch()->phase_id;
  std::map<std::string, pair<int,int> > &controlPointSpace = controlPointManagerProxy.ckLocalBranch()->controlPointSpace;
  int result;

  // if we already have control point values for phase, return them
  if( thisPhaseData->controlPoints.count(std::string(name))>0 && thisPhaseData->controlPoints[std::string(name)]>=0 ){
    CkPrintf("Already have control point values for phase. %s -> %d\n", name, (int)(thisPhaseData->controlPoints[std::string(name)]) );
    return thisPhaseData->controlPoints[std::string(name)];
  }
  

  if( phase_id < 4 || whichTuningScheme == AlwaysDefaults){
    // For the first few phases, just use the lower bound, or the default if one was provided 
    // This ensures that the ranges for all the control points are known before we do anything fancy
    result = lb;

    if(defaultControlPointValues.count(std::string(name)) > 0){
      int v = defaultControlPointValues[std::string(name)];
      CkPrintf("Startup phase using default value of %d for  \"%s\"\n", v, name);   
      result = v;
    }

  } else if(controlPointSpace.count(std::string(name)) == 0){
    // if this is the first time we've seen the CP, then return the lower bound
    result = lb;
    
  }  else {
    // Generate a plan one time for each phase
    controlPointManagerProxy.ckLocalBranch()->generatePlan();
    
    // Use the planned value:
    result = controlPointManagerProxy.ckLocalBranch()->newControlPoints[std::string(name)];
    
  }

  if(!isInRange(result,ub,lb)){
    std::cerr << "control point = " << result << " is out of range: " << lb << " " << ub << std::endl;
    fflush(stdout);
    fflush(stderr);
  }
  CkAssert(isInRange(result,ub,lb));
  thisPhaseData->controlPoints[std::string(name)] = result; // was insert() 

  controlPointSpace.insert(std::make_pair(std::string(name),std::make_pair(lb,ub))); 

  CkPrintf("Control Point \"%s\" for phase %d is: %d\n", name, phase_id, result);
  //  thisPhaseData->print();
  
  return result;
}


FLINKAGE int FTN_NAME(CONTROLPOINT, controlpoint)(CMK_TYPEDEF_INT4 *lb, CMK_TYPEDEF_INT4 *ub){
  CkAssert(CkMyPe() == 0);
  return controlPoint("FortranCP", *lb, *ub);
}




// void controlPointPriorityArray(const char *name, CProxy_ArrayBase &arraybase){
//   CkGroupID aid = arraybase.ckGetArrayID();
//   int groupIdx = aid.idx;
//   controlPointManagerProxy.ckLocalBranch()->associatePriorityArray(name, groupIdx);
//   //  CkPrintf("Associating control point \"%s\" with array id=%d\n", name, groupIdx );
// }


// /// Inform the control point framework that a named control point affects the priorities of some entry method  
// void controlPointPriorityEntry(const char *name, int idx){
//   controlPointManagerProxy.ckLocalBranch()->associatePriorityEntry(name, idx);
//   //  CkPrintf("Associating control point \"%s\" with EP id=%d\n", name, idx);
// }





void simplexScheme::adapt(std::map<std::string, std::pair<int,int> > & controlPointSpace, std::map<std::string,int> &newControlPoints, const int phase_id, instrumentedData &allData){

    if(useBestKnown){
        CkPrintf("Simplex Tuning: Simplex algorithm is done, using best known phase:\n");
        return;
    }


    if(firstSimplexPhase< 0){
        firstSimplexPhase = allData.phases.size()-1;
        CkPrintf("First simplex phase is %d\n", firstSimplexPhase);
    }

    int n = controlPointSpace.size();

    CkAssert(n>=2);


    if(simplexState == beginning){
        // First we evaluate n+1 random points, then we go to a different state
        if(allData.phases.size() < firstSimplexPhase + n+2  ){
            CkPrintf("Simplex Tuning: chose random configuration\n");

            // Choose random values from the middle third of the range in each dimension
            std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
            for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
                const std::string &name = cpsIter->first;
                const std::pair<int,int> &bounds = cpsIter->second;
                const int lb = bounds.first;
                const int ub = bounds.second;
                newControlPoints[name] = lb + randInt(ub-lb+1, name.c_str(), phase_id);
            }
        } else {
            // Set initial simplex:
            for(int i=0; i<n+1; i++){
                simplexIndices.insert(firstSimplexPhase+i);
            }
            CkAssert(simplexIndices.size() == n+1);

            // Transition to reflecting state
            doReflection(controlPointSpace, newControlPoints, phase_id, allData);

        }

    } else if (simplexState == reflecting){
        const double recentPhaseTime = allData.phases[allData.phases.size()-2]->medianTime();
        const double previousWorstPhaseTime = allData.phases[worstPhase]->medianTime();

        // Find the highest time from other points in the simplex
        double highestTimeForOthersInSimplex = 0.0;
        for(std::set<int>::iterator iter = simplexIndices.begin(); iter != simplexIndices.end(); ++iter){
            double t = allData.phases[*iter]->medianTime();
            if(*iter != worstPhase && t > highestTimeForOthersInSimplex) {
                highestTimeForOthersInSimplex = t;
            }
        }

        CkPrintf("After reflecting, the median time for the phase is %f, previous worst phase %d time was %f\n", recentPhaseTime, worstPhase, previousWorstPhaseTime);

        if(recentPhaseTime < highestTimeForOthersInSimplex){
            // if y* < yl,  transition to "expanding"
            doExpansion(controlPointSpace, newControlPoints, phase_id, allData);

        } else if (recentPhaseTime <= highestTimeForOthersInSimplex){
            // else if y* <= yi replace ph with p* and transition to "evaluatingOne"
            CkAssert(simplexIndices.size() == n+1);
            simplexIndices.erase(worstPhase);
            CkAssert(simplexIndices.size() == n);
            simplexIndices.insert(pPhase);
            CkAssert(simplexIndices.size() == n+1);
            // Transition to reflecting state
            doReflection(controlPointSpace, newControlPoints, phase_id, allData);

        } else {
            // if y* > yh
            if(recentPhaseTime <= worstTime){
                // replace Ph with P*
                CkAssert(simplexIndices.size() == n+1);
                simplexIndices.erase(worstPhase);
                CkAssert(simplexIndices.size() == n);
                simplexIndices.insert(pPhase);
                CkAssert(simplexIndices.size() == n+1);
                // Because we later will possibly replace Ph with P**, and we just replaced it with P*, we need to update our Ph reference
                worstPhase = pPhase;
                // Just as a sanity check, make sure we don't use the non-updated values.
                worst.clear();
            }

            // Now, form P** and do contracting phase
            doContraction(controlPointSpace, newControlPoints, phase_id, allData);

        }

    } else if (simplexState == doneExpanding){
        const double recentPhaseTime = allData.phases[allData.phases.size()-2]->medianTime();
        const double previousWorstPhaseTime = allData.phases[worstPhase]->medianTime();
        // A new configuration has been evaluated

        // Check to see if y** < y1
        if(recentPhaseTime < bestTime){
            // replace Ph by P**
            CkAssert(simplexIndices.size() == n+1);
            simplexIndices.erase(worstPhase);
            CkAssert(simplexIndices.size() == n);
            simplexIndices.insert(p2Phase);
            CkAssert(simplexIndices.size() == n+1);
        } else {
            //  replace Ph by P*
            CkAssert(simplexIndices.size() == n+1);
            simplexIndices.erase(worstPhase);
            CkAssert(simplexIndices.size() == n);
            simplexIndices.insert(pPhase);
            CkAssert(simplexIndices.size() == n+1);
        }

        // Transition to reflecting state
        doReflection(controlPointSpace, newControlPoints, phase_id, allData);

    }  else if (simplexState == contracting){
        const double recentPhaseTime = allData.phases[allData.phases.size()-2]->medianTime();
        const double previousWorstPhaseTime = allData.phases[worstPhase]->medianTime();
        // A new configuration has been evaluated

        // Check to see if y** > yh
        if(recentPhaseTime <= worstTime){
            // replace Ph by P**
            CkPrintf("Replacing phase %d with %d\n", worstPhase, p2Phase);
            CkAssert(simplexIndices.size() == n+1);
            simplexIndices.erase(worstPhase);
            CkAssert(simplexIndices.size() == n);
            simplexIndices.insert(p2Phase);
            CkAssert(simplexIndices.size() == n+1);
            // Transition to reflecting state
            doReflection(controlPointSpace, newControlPoints, phase_id, allData);

        } else {
            //  conceptually we will replace all Pi by (Pi+Pl)/2, but there is nothing to store this until after we have tried all of them
            simplexState = stillContracting;

            // A set of phases for which (P_i+P_l)/2 ought to be evaluated
            stillMustContractList = simplexIndices;

            CkPrintf("Simplex Tuning: Switched to state: stillContracting\n");
        }

    } else if (simplexState == stillContracting){
        CkPrintf("Simplex Tuning: stillContracting found %d configurations left to try\n", stillMustContractList.size());

        if(stillMustContractList.size()>0){
            int c = *stillMustContractList.begin();
            stillMustContractList.erase(c);
            CkPrintf("Simplex Tuning: stillContracting evaluating configuration derived from phase %d\n", c);

            std::vector<double> cPhaseConfig = pointCoords(allData, c);

            // Evaluate point P by storing new configuration in newControlPoints, and by transitioning to "reflecting" state
            int v = 0;
            for(std::map<std::string, std::pair<int,int> >::iterator cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
                const std::string &name = cpsIter->first;
                const std::pair<int,int> &bounds = cpsIter->second;
                const int lb = bounds.first;
                const int ub = bounds.second;

                double val = (cPhaseConfig[v] + best[v])/2.0;

                newControlPoints[name] = keepInRange((int)val,lb,ub);
                CkPrintf("Simplex Tuning: v=%d Reflected worst %d %s -> %f (ought to be %f )\n", (int)v, (int)worstPhase, (char*)name.c_str(), (double)newControlPoints[name], (double)P[v]);
                v++;
            }
        } else {
            // We have tried all configurations. We should update the simplex to refer to all the newly tried configurations, and start over
            CkAssert(stillMustContractList.size() == 0);
            simplexIndices.clear();
            CkAssert(simplexIndices.size()==0);
            for(int i=0; i<n+1; i++){
                simplexIndices.insert(allData.phases.size()-2-i);
            }
            CkAssert(simplexIndices.size()==n+1);

            // Transition to reflecting state
            doReflection(controlPointSpace, newControlPoints, phase_id, allData);

        }


    } else if (simplexState == expanding){
        const double recentPhaseTime = allData.phases[allData.phases.size()-2]->medianTime();
        const double previousWorstPhaseTime = allData.phases[worstPhase]->medianTime();
        // A new configuration has been evaluated

        // determine if y** < yl
        if(recentPhaseTime < bestTime){
            // replace Ph by P**
            CkAssert(simplexIndices.size() == n+1);
            simplexIndices.erase(worstPhase);
            CkAssert(simplexIndices.size() == n);
            simplexIndices.insert(p2Phase);
            CkAssert(simplexIndices.size() == n+1);
            // Transition to reflecting state
            doReflection(controlPointSpace, newControlPoints, phase_id, allData);

        } else {
            // else, replace ph with p*
            CkAssert(simplexIndices.size() == n+1);
            simplexIndices.erase(worstPhase);
            CkAssert(simplexIndices.size() == n);
            simplexIndices.insert(pPhase);
            CkAssert(simplexIndices.size() == n+1);
            // Transition to reflecting state
            doReflection(controlPointSpace, newControlPoints, phase_id, allData);
        }


    } else {
        CkAbort("Unknown simplexState");
    }

}



void simplexScheme::doReflection(std::map<std::string, std::pair<int,int> > & controlPointSpace, std::map<std::string,int> &newControlPoints, const int phase_id, instrumentedData &allData){

    int n = controlPointSpace.size();

    printSimplex(allData);

    computeCentroidBestWorst(controlPointSpace, newControlPoints, phase_id, allData);


    // Quit if the diameter of our simplex is small
    double maxr = 0.0;
    for(int i=0; i<n+1; i++){
        //      Compute r^2 of this simplex point from the centroid
        double r2 = 0.0;
        std::vector<double> p = pointCoords(allData, i);
        for(int d=0; d<p.size(); d++){
            double r1 = (p[d] * centroid[d]);
            r2 += r1*r1;
        }
        if(r2 > maxr)
            maxr = r2;
    }

#if 0
    // At some point we quit this tuning
    if(maxr < 10){
        useBestKnown = true;
        instrumentedPhase *best = allData.findBest();
        CkPrintf("Simplex Tuning: Simplex diameter is small, so switching over to best known phase:\n");

        std::map<std::string, std::pair<int,int> >::const_iterator cpsIter;
        for(cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter) {
            const std::string &name = cpsIter->first;
            newControlPoints[name] =  best->controlPoints[name];
        }
    }
#endif

    // Compute new point P* =(1+alpha)*centroid - alpha(worstPoint)

    pPhase = allData.phases.size()-1;
    P.resize(n);
    for(int i=0; i<n; i++){
        P[i] = (1.0+alpha) * centroid[i] - alpha * worst[i] ;
    }

    for(int i=0; i<P.size(); i++){
        CkPrintf("Simplex Tuning: P dimension %d is %f\n", i, P[i]);
    }


    // Evaluate point P by storing new configuration in newControlPoints, and by transitioning to "reflecting" state
    int v = 0;
    for(std::map<std::string, std::pair<int,int> >::iterator cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
        const std::string &name = cpsIter->first;
        const std::pair<int,int> &bounds = cpsIter->second;
        const int lb = bounds.first;
        const int ub = bounds.second;
        newControlPoints[name] = keepInRange((int)P[v],lb,ub);
        CkPrintf("Simplex Tuning: v=%d Reflected worst %d %s -> %f (ought to be %f )\n", (int)v, (int)worstPhase, (char*)name.c_str(), (double)newControlPoints[name], (double)P[v]);
        v++;
    }


    // Transition to "reflecting" state
    simplexState = reflecting;
    CkPrintf("Simplex Tuning: Switched to state: reflecting\n");

}




void simplexScheme::doExpansion(std::map<std::string, std::pair<int,int> > & controlPointSpace, std::map<std::string,int> &newControlPoints, const int phase_id, instrumentedData &allData){
    int n = controlPointSpace.size();
    printSimplex(allData);

    // Note that the original Nelder Mead paper has an error when it displays the equation for P** in figure 1.
    // I believe the equation for P** in the text on page 308 is correct.

    // Compute new point P2 = (1+gamma)*P - gamma(centroid)


    p2Phase = allData.phases.size()-1;
    P2.resize(n);
    for(int i=0; i<n; i++){
        P2[i] = ( (1.0+gamma) * P[i] - gamma * centroid[i] );
    }

    for(int i=0; i<P2.size(); i++){
        CkPrintf("P2 aka P** dimension %d is %f\n", i, P2[i]);
    }


    // Evaluate point P** by storing new configuration in newControlPoints, and by transitioning to "reflecting" state
    int v = 0;
    for(std::map<std::string, std::pair<int,int> >::iterator cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
        const std::string &name = cpsIter->first;
        const std::pair<int,int> &bounds = cpsIter->second;
        const int lb = bounds.first;
        const int ub = bounds.second;
        newControlPoints[name] = keepInRange((int)P2[v],lb,ub);
        CkPrintf("Simplex Tuning: v=%d worstPhase=%d Expanding %s -> %f (ought to be %f )\n", (int)v, (int)worstPhase, (char*)name.c_str(), (double)newControlPoints[name], (double)P[v]);
        v++;
    }


    // Transition to "doneExpanding" state
    simplexState = doneExpanding;
    CkPrintf("Simplex Tuning: Switched to state: doneExpanding\n");

}




void simplexScheme::doContraction(std::map<std::string, std::pair<int,int> > & controlPointSpace, std::map<std::string,int> &newControlPoints, const int phase_id, instrumentedData &allData){
    int n = controlPointSpace.size();
    printSimplex(allData);

    // Compute new point P2 = beta*Ph + (1-beta)*centroid


    p2Phase = allData.phases.size()-1;
    P2.resize(n);
    for(int i=0; i<n; i++){
        P2[i] = ( beta*worst[i] + (1.0-beta)*centroid[i] );
    }

    for(int i=0; i<P2.size(); i++){
        CkPrintf("P2 aka P** dimension %d is %f\n", i, P2[i]);
    }


    // Evaluate point P** by storing new configuration in newControlPoints, and by transitioning to "reflecting" state
    int v = 0;
    for(std::map<std::string, std::pair<int,int> >::iterator cpsIter=controlPointSpace.begin(); cpsIter != controlPointSpace.end(); ++cpsIter){
        const std::string &name = cpsIter->first;
        const std::pair<int,int> &bounds = cpsIter->second;
        const int lb = bounds.first;
        const int ub = bounds.second;
        newControlPoints[name] = keepInRange((int)P2[v],lb,ub);
        CkPrintf("Simplex Tuning: v=%d worstPhase=%d Contracting %s -> %f (ought to be %f )\n", (int)v, (int)worstPhase, (char*)name.c_str(), (double)newControlPoints[name], (double)P[v]);
        v++;
    }
    
    
    // Transition to "contracting" state
    simplexState = contracting;
    CkPrintf("Simplex Tuning: Switched to state: contracting\n");

}




void simplexScheme::computeCentroidBestWorst(std::map<std::string, std::pair<int,int> > & controlPointSpace, std::map<std::string,int> &newControlPoints, const int phase_id, instrumentedData &allData){
    int n = controlPointSpace.size();

    // Find worst performing point in the simplex
    worstPhase = -1;
    worstTime = -1.0;
    bestPhase = 10000000;
    bestTime = 10000000;
    for(std::set<int>::iterator iter = simplexIndices.begin(); iter != simplexIndices.end(); ++iter){
        double t = allData.phases[*iter]->medianTime();
        if(t > worstTime){
            worstTime = t;
            worstPhase = *iter;
        }
        if(t < bestTime){
            bestTime = t;
            bestPhase = *iter;
        }
    }
    CkAssert(worstTime != -1.0 && worstPhase != -1 && bestTime != 10000000 && bestPhase != 10000000);

    best = pointCoords(allData, bestPhase);
    CkAssert(best.size() == n);

    worst = pointCoords(allData, worstPhase);
    CkAssert(worst.size() == n);

    // Calculate centroid of the remaining points in the simplex
    centroid.resize(n);
    for(int i=0; i<n; i++){
        centroid[i] = 0.0;
    }

    int numPts = 0;

    for(std::set<int>::iterator iter = simplexIndices.begin(); iter != simplexIndices.end(); ++iter){
        if(*iter != worstPhase){
            numPts ++;
            // Accumulate into the result vector
            int c = 0;
            for(std::map<std::string,int>::iterator citer = allData.phases[*iter]->controlPoints.begin(); citer != allData.phases[*iter]->controlPoints.end(); ++citer){
                centroid[c] += citer->second;
                c++;
            }

        }
    }

    // Now divide the sums by the number of points.
    for(int v = 0; v<centroid.size(); v++) {
        centroid[v] /= (double)numPts;
    }

    CkAssert(centroid.size() == n);

    for(int i=0; i<centroid.size(); i++){
        CkPrintf("Centroid dimension %d is %f\n", i, centroid[i]);
    }


}



std::vector<double> simplexScheme::pointCoords(instrumentedData &allData, int i){
    std::vector<double> result;
    for(std::map<std::string,int>::iterator citer = allData.phases[i]->controlPoints.begin(); citer != allData.phases[i]->controlPoints.end(); ++citer){
        result.push_back((double)citer->second);
    }
    return result;
}




void ControlPointWriteOutputToDisk(){
    CkAssert(CkMyPe() == 0);
    controlPointManagerProxy.ckLocalBranch()->writeOutputToDisk();
}



#include "ControlPoints.def.h"

#endif

---