CompuCell/DOXYGEN/core/CompuCell3D/Potts3D/Potts3D.cpp

Go to the documentation of this file.

/*************************************************************************
*    CompuCell - A software framework for multimodel simulations of     *
* biocomplexity problems Copyright (C) 2003 University of Notre Dame,   *
*                             Indiana                                   *
*                                                                       *
* This program is free software; IF YOU AGREE TO CITE USE OF CompuCell  *
*  IN ALL RELATED RESEARCH PUBLICATIONS according to the terms of the   *
*  CompuCell GNU General Public License RIDER you can redistribute it   *
* and/or modify it under the terms of the GNU General Public License as *
*  published by the Free Software Foundation; either version 2 of the   *
*         License, or (at your option) any later version.               *
*                                                                       *
* This program is distributed in the hope that it will be useful, but `  *
*      WITHOUT ANY WARRANTY; without even the implied warranty of       *
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU    *
*             General Public License for more details.                  *
*                                                                       *
*  You should have received a copy of the GNU General Public License    *
*     along with this program; if not, write to the Free Software       *
*      Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.        *
*************************************************************************/


#include "Cell.h"
#include "CellTypeMotilityData.h"
#include "DefaultAcceptanceFunction.h"
#include "AcceptanceFunction.h"
#include "EnergyFunction.h"
#include "CellGChangeWatcher.h"
#include "Stepper.h"
#include "FixedStepper.h"
#include "AttributeAdder.h"
#include <CompuCell3D/Automaton/Automaton.h>

#include <CompuCell3D/Field3D/WatchableField3D.h>

#include <CompuCell3D/Potts3D/TypeTransition.h>
#include "EnergyFunctionCalculator.h"
#include "EnergyFunctionCalculatorStatistics.h"

#include <CompuCell3D/Simulator.h>
//#include <CompuCell3D/plugins/Volume/VolumePlugin.h>
//#include <CompuCell3D/plugins/Surface/SurfacePlugin.h>

#include <BasicUtils/BasicRandomNumberGenerator.h>
#include <BasicUtils/BasicPluginInfo.h>


#include <PublicUtilities/StringUtils.h>
#include <sstream>

#define EXP_STL
#include "Potts3D.h"

using namespace CompuCell3D;
using namespace std;

Potts3D::Potts3D() :
cellFieldG(0),
attrAdder(0),
acceptanceFunction(&defaultAcceptanceFunction), 
energy(0),
depth(1.0),
recentlyCreatedCellId(1),
debugOutputFrequency(1),
sim(0),
automaton(0)
{
        neighbors.assign(100,Point3D());//statically allocated this buffer maybe will come with something better later
        frozenTypeVec.assign(0,0);
        energyCalculator = new EnergyFunctionCalculator();
        energyCalculator->setPotts(this);
        typeTransition=new TypeTransition();
        metropolisFcnPtr=&Potts3D::metropolisFast;
        cellInventory.setPotts3DPtr(this);
}

Potts3D::Potts3D(const Dim3D dim) :
cellFieldG(0),
attrAdder(0),
acceptanceFunction(&defaultAcceptanceFunction), 
energy(0) ,
depth(1.0),
recentlyCreatedCellId(1),
debugOutputFrequency(1),
sim(0),
automaton(0)
{
        neighbors.assign(100,Point3D());//statically allocated this buffer maybe will come with something better later
        frozenTypeVec.assign(0,0);

        createCellField(dim);
        energyCalculator = new EnergyFunctionCalculator();
        energyCalculator->setPotts(this); 
        typeTransition=new TypeTransition();
        metropolisFcnPtr=&Potts3D::metropolisFast;
        cellInventory.setPotts3DPtr(this);

}


Potts3D::~Potts3D() {
        if (cellFieldG) delete cellFieldG;
        if (energyCalculator) delete energyCalculator; energyCalculator=0;
        if (typeTransition) delete typeTransition; typeTransition=0;
        //   if (attrAdder) delete attrAdder; attrAdder=0;
}

void Potts3D::createEnergyFunction(std::string _energyFunctionType){
        if(_energyFunctionType=="Statistics"){
                if (energyCalculator) delete energyCalculator; energyCalculator=0;
                energyCalculator=new EnergyFunctionCalculatorStatistics();
                energyCalculator->setPotts(this);
                //initialize Statistics Output Energy Finction Here
                return;
        }
        else{
                //default is not to reassign energy finction calculator
                return;
        }
}

LatticeType Potts3D::getLatticeType(){
        return BoundaryStrategy::getInstance()->getLatticeType();
}

void Potts3D::setDepth(double _depth){
        //this function has to be called after initializing bondary strategy and after creating cellFieldG
        //By default Boundary Strategy will precalculate neighbors up to certain depth (4.0). However if user requests more
        //depth additional calculations will be requested here
        depth=_depth;
        float maxDistance=BoundaryStrategy::getInstance()->getMaxDistance();
        if(maxDistance<depth){
                //in such situation user requests depth that is greater than default maxDistance 
                BoundaryStrategy::getInstance()->prepareNeighborLists(depth);
        }
        Dim3D dim=cellFieldG->getDim();
        minCoordinates=Point3D(0,0,0);
        maxCoordinates=Point3D(dim.x,dim.y,dim.z);

        //will calculate necesary space for neighbor storage
        //this may not work for lattice types other than square

        //   Dim3D dim=cellFieldG->getDim();
        // 
        //   minCoordinates=Point3D(0,0,0);
        //   maxCoordinates=Point3D(dim.x,dim.y,dim.z);
        // 
        // 
        //   Point3D middlePt(dim.x/2,dim.y/2,dim.z/2);
        //   
        //    unsigned int token = 0;
        //    int numNeighbors = 0;
        //    double distance;
        //    Point3D testPt;
        //    token =0;
        //    distance=0;
        //    
        //    
        //       while(true){
        //       testPt = cellFieldG->getNeighbor(middlePt, token, distance);
        //       if (distance > depth) break;
        //    
        //          numNeighbors++;
        //       }
        //    //empty previously allocated container for neighbors
        //    neighbors.clear();
        //    neighbors.assign(numNeighbors+1,Point3D());

        maxNeighborIndex=BoundaryStrategy::getInstance()->getMaxNeighborIndexFromDepth(depth);
        cerr<<"\t\t\t\t\t setDepth  maxNeighborIndex="<<maxNeighborIndex<<endl;
        neighbors.clear();
        neighbors.assign(maxNeighborIndex+1,Point3D());

}

void Potts3D::setNeighborOrder(unsigned int _neighborOrder){
        BoundaryStrategy::getInstance()->prepareNeighborListsBasedOnNeighborOrder(_neighborOrder);
        maxNeighborIndex=BoundaryStrategy::getInstance()->getMaxNeighborIndexFromNeighborOrder(_neighborOrder);
        cerr<<"\t\t\t\t\t setNeighborOrder  maxNeighborIndex="<<maxNeighborIndex<<endl;
        Dim3D dim=cellFieldG->getDim();
        minCoordinates=Point3D(0,0,0);
        maxCoordinates=Point3D(dim.x,dim.y,dim.z);

        neighbors.clear();
        neighbors.assign(maxNeighborIndex+1,Point3D());

}


void Potts3D::createCellField(const Dim3D dim) {

        ASSERT_OR_THROW("createCellField() cell field G already created!", !cellFieldG);
        cellFieldG = new WatchableField3D<CellG *>(dim, 0); //added


}

void Potts3D::registerAttributeAdder(AttributeAdder * _attrAdder){
        attrAdder=_attrAdder;
}

// void Potts3D::registerTypeChangeWatcher(TypeChangeWatcher * _typeChangeWatcher){
//    typeTransition->registerTypeChangeWatcher(_typeChangeWatcher);
// }

void Potts3D::registerAutomaton(Automaton* autom) {automaton = autom;}
Automaton* Potts3D::getAutomaton() {return automaton;}

void Potts3D::registerEnergyFunction(EnergyFunction *function) {

        energyCalculator->registerEnergyFunctionWithName(function,function->toString());
        //    sim->registerSteerableObject(function);

}

void Potts3D::registerEnergyFunctionWithName(EnergyFunction *_function,std::string _functionName){

        energyCalculator->registerEnergyFunctionWithName(_function,_functionName);
        //   sim->registerSteerableObject(_function);

}


void Potts3D::unregisterEnergyFunction(std::string _functionName) {

        energyCalculator->unregisterEnergyFunction(_functionName);
        //    sim->unregisterSteerableObject(_functionName);
        return;

}

double Potts3D::getEnergy() {return energy;}

void Potts3D::setAcceptanceFunctionByName(std::string _acceptanceFunctionName){
        if(_acceptanceFunctionName=="FirstOrderExpansion"){
                acceptanceFunction=&firstOrderExpansionAcceptanceFunction;
                //          cerr<<"setting FirstOrderExpansion"<<endl;
        }else {
                acceptanceFunction=&defaultAcceptanceFunction;
        }

}

void Potts3D::registerAcceptanceFunction(AcceptanceFunction *function) {
        ASSERT_OR_THROW("registerAcceptanceFunction() function cannot be NULL!",
                function);

        acceptanceFunction = function;
}


void Potts3D::registerCellGChangeWatcher(CellGChangeWatcher *_watcher) {
        ASSERT_OR_THROW("registerBCGChangeWatcher() _watcher cannot be NULL!",_watcher);

        cellFieldG->addChangeWatcher(_watcher); 
        //    sim->registerSteerableObject(_watcher);
}


void Potts3D::registerClassAccessor(BasicClassAccessorBase *_accessor) {
        ASSERT_OR_THROW("registerClassAccessor() _accessor cannot be NULL!", _accessor);

        cellFactoryGroup.registerClass(_accessor);
}


void Potts3D::registerStepper(Stepper *stepper) {
        ASSERT_OR_THROW("registerStepper() stepper cannot be NULL!", stepper);

        steppers.push_back(stepper);
}
void Potts3D::registerFixedStepper(FixedStepper *fixedStepper) {
        ASSERT_OR_THROW("registerStepper() stepper cannot be NULL!", fixedStepper);

        fixedSteppers.push_back(fixedStepper);
}


CellG * Potts3D::createCellG(const Point3D pt) {
        ASSERT_OR_THROW("createCell() cellFieldG Point out of range!", cellFieldG->isValid(pt));

        CellG * cell = new CellG();
        cell->extraAttribPtr=cellFactoryGroup.create();
        cell->id=recentlyCreatedCellId;
        ++recentlyCreatedCellId;
        cellInventory.addToInventory(cell);
        cellFieldG->set(pt, cell);
        if(attrAdder){
                attrAdder->addAttribute(cell);
        }
        return cell;   
}


void Potts3D::destroyCellG(CellG *cell,bool _removeFromInventory) {
        if(cell->extraAttribPtr){
                cellFactoryGroup.destroy(cell->extraAttribPtr);
                cell->extraAttribPtr=0;
        }
        if(cell->pyAttrib && attrAdder){
                attrAdder->destroyAttribute(cell);
        }
        //had to introduce these two cases because in the Cell inventory destructor we deallocate memory of pointers stored int the set
        //Since this is done during interation over set changing pointer (cell=0) or removing anything from set might corrupt container or invalidate iterators
        if(_removeFromInventory){
                cellInventory.removeFromInventory(cell);
                delete cell;
                cell=0;
        }else{
                delete cell;
        }

}


double Potts3D::totalEnergy() {
        double energy = 0;
        Dim3D dim = cellFieldG->getDim();

        Point3D pt;
        for (pt.z = 0; pt.z < dim.z; pt.z++)
                for (pt.y = 0; pt.y < dim.y; pt.y++)
                        for (pt.x = 0; pt.x < dim.x; pt.x++)
                                for (unsigned int i = 0; i < energyFunctions.size(); i++)
                                        energy += energyFunctions[i]->localEnergy(pt);

        return energy;

}


double Potts3D::changeEnergy(Point3D pt, const CellG *newCell,
                                                                          const CellG *oldCell) {
                                                                                  double change = 0;
                                                                                  for (unsigned int i = 0; i < energyFunctions.size(); i++)
                                                                                          change += energyFunctions[i]->changeEnergy(pt, newCell, oldCell);

                                                                                  return change;
}


unsigned int Potts3D::metropolis(const unsigned int steps, const double temp) {
        return (this->*metropolisFcnPtr)(steps,temp);
}

unsigned int Potts3D::metropolisList(const unsigned int steps, const double temp) {
        ASSERT_OR_THROW("Potts3D: cell field G not initialized", cellFieldG);


        flips = 0;
        attemptedEC=0;
        BasicRandomNumberGenerator *rand = BasicRandomNumberGenerator::getInstance();
        Dim3D dim = cellFieldG->getDim();
        ASSERT_OR_THROW("Potts3D: You must supply an acceptance function!",
                acceptanceFunction);

        //   numberOfAttempts=steps;
        numberOfAttempts=(int)(maxCoordinates.x-minCoordinates.x)*(maxCoordinates.y-minCoordinates.y)*(maxCoordinates.z-minCoordinates.z)*sim->getFlip2DimRatio();

        BoundaryStrategy * boundaryStrategy = BoundaryStrategy::getInstance();
        for (unsigned int i = 0; i < numberOfAttempts; i++) {

                currentAttempt=i;
                //run fixed steppers - they are executed regardless whether spin flip take place or not . Note, regular steopers are executed only after spin flip attepmts takes place 
                for (unsigned int j = 0; j < fixedSteppers.size(); j++)
                        fixedSteppers[j]->step();



                Point3D pt;

                // Pick a random point
                pt.x = rand->getInteger(minCoordinates.x, maxCoordinates.x - 1);
                pt.y = rand->getInteger(minCoordinates.y, maxCoordinates.y - 1);
                pt.z = rand->getInteger(minCoordinates.z, maxCoordinates.z - 1);


                //     pt.x = rand->getInteger(0, dim.x - 1);
                //     pt.y = rand->getInteger(0, dim.y - 1);
                //     pt.z = rand->getInteger(0, dim.z - 1);

                CellG *cell = cellFieldG->get(pt);

                if(sizeFrozenTypeVec && cell){
                        if(checkIfFrozen(cell->type))
                                continue;
                }

                unsigned int token = 0;
                int numNeighbors = 0;
                double distance;

                token =0;
                distance=0;
                while(true){

                        neighbors[numNeighbors] = cellFieldG->getNeighbor(pt, token, distance);
                        if (distance > depth) break;

                        if (cellFieldG->get(neighbors[numNeighbors]) != cell)
                                numNeighbors++;

                }

                // If no boundry neighbors were found start over
                if (!numNeighbors) continue;

                // Pick a random neighbor
                Point3D changePixel = neighbors[rand->getInteger(0, numNeighbors - 1)];




                if(sizeFrozenTypeVec && cellFieldG->get(changePixel)){
                        if(checkIfFrozen(cellFieldG->get(changePixel)->type))
                                continue;
                }
                ++attemptedEC;  

                flipNeighbor=pt;


                // Calculate change in energy

                double change = energyCalculator->changeEnergy(changePixel, cell, cellFieldG->get(changePixel),i); 

                // Acceptance based on probability
                double motility=0.0;
                if(cellTypeMotilityVec.size()){
                        unsigned int newCellTypeId=(cell ? (unsigned int)cell->type :0);
                        unsigned int oldCellTypeId=(cellFieldG->get(changePixel)? (unsigned int)cellFieldG->get(changePixel)->type :0);
                        if(newCellTypeId && oldCellTypeId)
                                motility=(cellTypeMotilityVec[newCellTypeId]<cellTypeMotilityVec[oldCellTypeId] ? cellTypeMotilityVec[newCellTypeId]:cellTypeMotilityVec[oldCellTypeId]);
                        else if(newCellTypeId){
                                motility=cellTypeMotilityVec[newCellTypeId];
                        }else if (oldCellTypeId){
                                motility=cellTypeMotilityVec[oldCellTypeId];
                        }else{//should never get here
                                motility=0;
                        }
                }else{
                        motility=temp;
                }
                double prob = acceptanceFunction->accept(motility, change);



                if (prob >= 1 || rand->getRatio() < prob) {
                        // Accept the change
                        energy += change;

                        cellFieldG->set(changePixel, cell);
                        flips++;
                        energyCalculator->setLastFlipAccepted(true);
                }else{
                        energyCalculator->setLastFlipAccepted(false);
                }


                // Run steppers
                for (unsigned int j = 0; j < steppers.size(); j++)
                        steppers[j]->step();
        }
        unsigned int currentStep=sim->getStep();
        if(! (currentStep % debugOutputFrequency) ){
                cerr<<"Number of Attempted Energy Calculations="<<attemptedEC<<endl;
        }
        return flips;
}

unsigned int Potts3D::metropolisFast(const unsigned int steps, const double temp) {
        ASSERT_OR_THROW("Potts3D: cell field G not initialized", cellFieldG);


        flips = 0;
        attemptedEC=0;
        BasicRandomNumberGenerator *rand = BasicRandomNumberGenerator::getInstance();
        Dim3D dim = cellFieldG->getDim();
        ASSERT_OR_THROW("Potts3D: You must supply an acceptance function!",
                acceptanceFunction);
        //cerr<<"steps="<<steps<<" temp="<<temp<<" acceptanceFunction="<<acceptanceFunction<<endl;

        numberOfAttempts=(int)(maxCoordinates.x-minCoordinates.x)*(maxCoordinates.y-minCoordinates.y)*(maxCoordinates.z-minCoordinates.z)*sim->getFlip2DimRatio();
        //   numberOfAttempts=steps;


        cerr<<"FAST numberOfAttempts="<<numberOfAttempts<<endl;
        BoundaryStrategy * boundaryStrategy = BoundaryStrategy::getInstance();
        for (unsigned int i = 0; i < numberOfAttempts; i++) {

                currentAttempt=i;
                //run fixed steppers - they are executed regardless whether spin flip take place or not . Note, regular steopers are executed only after spin flip attepmts takes place 
                for (unsigned int j = 0; j < fixedSteppers.size(); j++)
                        fixedSteppers[j]->step();



                Point3D pt;

                // Pick a random point
                pt.x = rand->getInteger(minCoordinates.x, maxCoordinates.x - 1);
                pt.y = rand->getInteger(minCoordinates.y, maxCoordinates.y - 1);
                pt.z = rand->getInteger(minCoordinates.z, maxCoordinates.z - 1);
                //    cerr<<"pt flip="<<pt<<endl;
                //     pt.x = rand->getInteger(0, dim.x - 1);
                //     pt.y = rand->getInteger(0, dim.y - 1);
                //     pt.z = rand->getInteger(0, dim.z - 1);
                //     cout<<"pt="<<pt<<" rand->getSeed()="<<rand->getSeed()<<endl;


                CellG *cell = cellFieldG->get(pt);

                if(sizeFrozenTypeVec && cell){
                        if(checkIfFrozen(cell->type))
                                continue;
                }

                unsigned int directIdx=rand->getInteger(0, maxNeighborIndex);
                //      cerr<<"directIdx="<<directIdx<<endl;

                Neighbor n=boundaryStrategy->getNeighborDirect(pt,directIdx);
                //      cerr<<"n.pt"<<n.pt<<" n.distance="<<n.distance<<endl;

                if(!n.distance){
                        //if distance is 0 then the neighbor returned is invalid
                        continue;
                }
                Point3D changePixel = n.pt;
                //       cerr<<"pt="<<pt<<" n.pt="<<n.pt<<" difference="<<pt-n.pt<<endl;
                //check if changePixel refers to different cell. 

                if (cellFieldG->get(changePixel) == cell){
                        continue;//skip the rest of the loop if change pixel points to the same cell as pt
                }else{
                        ;

                }

                //     cout<<"pt="<<pt<<" changePixel="<<changePixel<<endl;

                if(sizeFrozenTypeVec && cellFieldG->get(changePixel)){
                        if(checkIfFrozen(cellFieldG->get(changePixel)->type))
                                continue;
                }
                ++attemptedEC;  

                flipNeighbor=pt;

                // Calculate change in energy
                //cerr<<"steps="<<steps<<" temp="<<temp<<" acceptanceFunction="<<acceptanceFunction<<endl;

                double change = energyCalculator->changeEnergy(changePixel, cell, cellFieldG->get(changePixel),i); 

                //cerr<<"This is change: "<<change<<endl;       

                //cerr<<"steps="<<steps<<" temp="<<temp<<" acceptanceFunction="<<acceptanceFunction<<endl;

                // Acceptance based on probability
                double motility=0.0;
                if(cellTypeMotilityVec.size()){
                        unsigned int newCellTypeId=(cell ? (unsigned int)cell->type :0);
                        unsigned int oldCellTypeId=(cellFieldG->get(changePixel)? (unsigned int)cellFieldG->get(changePixel)->type :0);
                        if(newCellTypeId && oldCellTypeId)
                                motility=(cellTypeMotilityVec[newCellTypeId]<cellTypeMotilityVec[oldCellTypeId] ? cellTypeMotilityVec[newCellTypeId]:cellTypeMotilityVec[oldCellTypeId]);
                        else if(newCellTypeId){
                                motility=cellTypeMotilityVec[newCellTypeId];
                        }else if (oldCellTypeId){
                                motility=cellTypeMotilityVec[oldCellTypeId];
                        }else{//should never get here
                                motility=0;
                        }
                }else{
                        motility=temp;
                }

                //cerr<<"newCell="<<cell<<" oldCell="<<cellFieldG->get(changePixel)<<" motility="<<motility<<endl;

                double prob = acceptanceFunction->accept(motility, change);
                //     cerr<<"prob="<<prob<<endl;


                //     cout<<"change E="<<change<<" prob="<<prob<<endl;
                if (prob >= 1 || rand->getRatio() < prob) {
                        // Accept the change
                        //        cout<<"accepted flip "<< change<<endl;
                        energy += change;
                        //       cerr<<"energy="<<energy<<endl;
                        cellFieldG->set(changePixel, cell);
                        flips++;
                        energyCalculator->setLastFlipAccepted(true);
                }else{
                        energyCalculator->setLastFlipAccepted(false);
                }


                // Run steppers
                for (unsigned int j = 0; j < steppers.size(); j++)
                        steppers[j]->step();

                //     exit(0);
        }
        unsigned int currentStep=sim->getStep();
        if(! (currentStep % debugOutputFrequency) ){
                cerr<<"Number of Attempted Energy Calculations="<<attemptedEC<<endl;
        }
        //    exit(0);
        return flips;

}


unsigned int Potts3D::metropolisBoundaryWalker(const unsigned int steps, const double temp) {
        ASSERT_OR_THROW("Potts3D: cell field G not initialized", cellFieldG);


        flips = 0;
        attemptedEC=0;
        BasicRandomNumberGenerator *rand = BasicRandomNumberGenerator::getInstance();
        Dim3D dim = cellFieldG->getDim();
        ASSERT_OR_THROW("Potts3D: You must supply an acceptance function!",
                acceptanceFunction);
        //cerr<<"steps="<<steps<<" temp="<<temp<<" acceptanceFunction="<<acceptanceFunction<<endl;

        //numberOfAttempts=(int)(maxCoordinates.x-minCoordinates.x)*(maxCoordinates.y-minCoordinates.y)*(maxCoordinates.z-minCoordinates.z)*sim->getFlip2DimRatio();
        //   numberOfAttempts=steps;
        //Number of attempts is equal to number of boundary pixels
        numberOfAttempts=boundaryPixelSet.size();

        cerr<<"numberOfAttempts="<<numberOfAttempts<<endl;
        long boundaryPointIndex;
        long counter=0;
        set<Point3D>::iterator sitr;
        vector<Point3D> boundaryPointVector;
        boundaryPointVector.assign(boundaryPixelSet.begin(),boundaryPixelSet.end());

        BoundaryStrategy * boundaryStrategy = BoundaryStrategy::getInstance();
        for (unsigned int i = 0; i < numberOfAttempts; i++) {

                currentAttempt=i;
                //run fixed steppers - they are executed regardless whether spin flip take place or not . Note, regular steopers are executed only after spin flip attepmts takes place 
                for (unsigned int j = 0; j < fixedSteppers.size(); j++)
                        fixedSteppers[j]->step();



                Point3D pt;

                // Pick a random integer
                boundaryPointIndex=rand->getInteger(minCoordinates.x, boundaryPointVector.size()-1);     
                pt=boundaryPointVector[boundaryPointIndex];

                //    boundaryPointIndex=rand->getInteger(minCoordinates.x, boundaryPixelSet.size()-1);//we use most up-to-date set size
                //  counter=0;
                //for (set<Point3D>::iterator sitr=boundaryPixelSet.begin() ; sitr!=boundaryPixelSet.end(); ++sitr){
                // if(counter==boundaryPointIndex){
                //      pt=*sitr;
                //      break;
                // }
                // ++counter;
                //}

                //    cerr<<"pt flip="<<pt<<endl;
                //     pt.x = rand->getInteger(0, dim.x - 1);
                //     pt.y = rand->getInteger(0, dim.y - 1);
                //     pt.z = rand->getInteger(0, dim.z - 1);
                //     cout<<"pt="<<pt<<" rand->getSeed()="<<rand->getSeed()<<endl;


                CellG *cell = cellFieldG->get(pt);

                if(sizeFrozenTypeVec && cell){
                        if(checkIfFrozen(cell->type))
                                continue;
                }

                unsigned int directIdx=rand->getInteger(0, maxNeighborIndex);
                //      cerr<<"directIdx="<<directIdx<<endl;

                Neighbor n=boundaryStrategy->getNeighborDirect(pt,directIdx);
                //      cerr<<"n.pt"<<n.pt<<" n.distance="<<n.distance<<endl;

                if(!n.distance){
                        //if distance is 0 then the neighbor returned is invalid
                        continue;
                }
                Point3D changePixel = n.pt;
                //       cerr<<"pt="<<pt<<" n.pt="<<n.pt<<" difference="<<pt-n.pt<<endl;
                //check if changePixel refers to different cell. 

                if (cellFieldG->get(changePixel) == cell){
                        continue;//skip the rest of the loop if change pixel points to the same cell as pt
                }else{
                        ;

                }

                //     cout<<"pt="<<pt<<" changePixel="<<changePixel<<endl;

                if(sizeFrozenTypeVec && cellFieldG->get(changePixel)){
                        if(checkIfFrozen(cellFieldG->get(changePixel)->type))
                                continue;
                }
                ++attemptedEC;  

                flipNeighbor=pt;

                // Calculate change in energy
                //cerr<<"steps="<<steps<<" temp="<<temp<<" acceptanceFunction="<<acceptanceFunction<<endl;

                double change = energyCalculator->changeEnergy(changePixel, cell, cellFieldG->get(changePixel),i); 

                //cerr<<"This is change: "<<change<<endl;       

                //cerr<<"steps="<<steps<<" temp="<<temp<<" acceptanceFunction="<<acceptanceFunction<<endl;

                // Acceptance based on probability
                double motility=0.0;
                if(cellTypeMotilityVec.size()){
                        unsigned int newCellTypeId=(cell ? (unsigned int)cell->type :0);
                        unsigned int oldCellTypeId=(cellFieldG->get(changePixel)? (unsigned int)cellFieldG->get(changePixel)->type :0);
                        if(newCellTypeId && oldCellTypeId)
                                motility=(cellTypeMotilityVec[newCellTypeId]<cellTypeMotilityVec[oldCellTypeId] ? cellTypeMotilityVec[newCellTypeId]:cellTypeMotilityVec[oldCellTypeId]);
                        else if(newCellTypeId){
                                motility=cellTypeMotilityVec[newCellTypeId];
                        }else if (oldCellTypeId){
                                motility=cellTypeMotilityVec[oldCellTypeId];
                        }else{//should never get here
                                motility=0;
                        }
                }else{
                        motility=temp;
                }
                double prob = acceptanceFunction->accept(motility, change);
                //     cerr<<"prob="<<prob<<endl;


                //     cout<<"change E="<<change<<" prob="<<prob<<endl;
                if (prob >= 1 || rand->getRatio() < prob) {
                        // Accept the change
                        //        cout<<"accepted flip "<< change<<endl;
                        energy += change;
                        //       cerr<<"energy="<<energy<<endl;
                        cellFieldG->set(changePixel, cell);
                        flips++;
                        energyCalculator->setLastFlipAccepted(true);
                }else{
                        energyCalculator->setLastFlipAccepted(false);
                }


                // Run steppers
                for (unsigned int j = 0; j < steppers.size(); j++)
                        steppers[j]->step();

                //     exit(0);
        }
        unsigned int currentStep=sim->getStep();
        if(! (currentStep % debugOutputFrequency) ){
                cerr<<"Number of Attempted Energy Calculations="<<attemptedEC<<endl;
        }
        //    exit(0);
        return flips;

}


void Potts3D::setMetropolisAlgorithm(std::string _algName){

        string algName=_algName;
        changeToLower(algName);

        if(algName=="list"){
                metropolisFcnPtr=&Potts3D::metropolisList;
        }else if (algName=="fast"){
                metropolisFcnPtr=&Potts3D::metropolisFast;
        }else if (algName=="boundarywalker"){
                metropolisFcnPtr=&Potts3D::metropolisBoundaryWalker;
        }else{
                metropolisFcnPtr=&Potts3D::metropolisFast;
        }

}



void Potts3D::setCellTypeMotilityVec(std::vector<float> & _cellTypeMotilityVec){
        cellTypeMotilityVec=_cellTypeMotilityVec;
}
void Potts3D::initializeCellTypeMotility(std::vector<CellTypeMotilityData> & _cellTypeMotilityVector){

        ASSERT_OR_THROW("AUTOMATON IS NOT INITIALIZED",automaton);

        unsigned int typeIdMax=0;
        //finding max typeId
        for(int i =0 ; i < _cellTypeMotilityVector.size() ;++i ){
                //cerr<<"CHECKING "<<_cellTypeMotilityVector[i].typeName<<endl;
                unsigned int id=automaton->getTypeId(_cellTypeMotilityVector[i].typeName);
                if(id>typeIdMax)
                        typeIdMax=id;
        }

        cellTypeMotilityVec.assign(typeIdMax+1,0.0);
        for(int i =0 ; i < _cellTypeMotilityVector.size() ;++i ){
                unsigned int id = automaton->getTypeId(_cellTypeMotilityVector[i].typeName);
                cellTypeMotilityVec[id]=_cellTypeMotilityVector[i].motility;
        }

        //for(int i =0 ; i < _cellTypeMotilityVector.size() ;++i ){
        //      
        //      cerr<<"cellTypeMotilityVec["<<i<<"]="<<cellTypeMotilityVec[i]<<endl;
        //}
        //
        //exit(0);
}


bool Potts3D::checkIfFrozen(unsigned char _type){



        for (unsigned int i = 0 ; i< frozenTypeVec.size(); ++i ){
                if(frozenTypeVec[i]==_type)
                        return true;
        }
        return false;

}
void