diff --git a/TREES_code/simulation.cpp b/TREES_code/simulation.cpp
index 304a669..9ff6acb 100644
--- a/TREES_code/simulation.cpp
+++ b/TREES_code/simulation.cpp
@@ -45,7 +45,8 @@ Simulation::Simulation(std::string& parameterFileName){
     verbose = pf.getBool("verbose");
     tmax = pf.getPositiveLong("t_max");
     sampleInterval = pf.getPositiveLong("sample_interval");
-
+    checkpointInterval = pf.getLong("checkpoint_interval"); // checkPoints==0 turns this feature off
+    keepOldCheckpoints = pf.getBool("keep_old_checkpoints");
     std::string seedS = pf.getStringLower("seed");
     if (seedS[0]=='r') {
         seed = randomize();
@@ -59,8 +60,6 @@ Simulation::Simulation(std::string& parameterFileName){
     
     pf.close();
 
-    sampleList.clear();
-
     // save copy of parameterFile:
     std::ifstream f(parameterFileName);
     if (f) {
@@ -77,14 +76,26 @@ Simulation::~Simulation() {
     }
 }
 
-void Simulation::runAndSave(int replicate) {
+void Simulation::initialize(int replicate) {
     thePop->initialize();
+    sampleList.clear();
+    checkpointList.clear();
+    makeFileName(replicate);
     // Sample generation 0:
     sampleList.push_back(thePop->makeSample());
+}
+
+void Simulation::makeFileName(int replicate) {
+    std::ostringstream fileName;//(simName, std::ios_base::app);
+    fileName << simName << "_results_" << replicate << ".sim";
+    resultsFileName = fileName.str();
+}
+
+void Simulation::runFromGeneration(timeType gen) {
     double t1 = getNow();
     double prevSec = 0;
     
-    for(timeType t=1; t<=tmax; ++t) {
+    for(timeType t=gen+1; t<=tmax; ++t) {
         //std::cout << "make gen " << t << '\n';
         thePop->makeNextGeneration();
         if (thePop->size()==0) {
@@ -101,26 +112,116 @@ void Simulation::runAndSave(int replicate) {
                 prevSec = nowSec;
             }
         }
+        if (checkpointInterval>0 && t%checkpointInterval==0) {
+            unsigned seed = rand();
+            if (keepOldCheckpoints) {
+                checkpointList.push_back(thePop->makeCheckpoint(seed));
+            } else {
+                checkpointList.assign(1,thePop->makeCheckpoint(seed));
+            }
+            randseed(seed);
+            save();
+        }
+    }
+    if (checkpointInterval==0 || tmax%checkpointInterval != 0) {
+        save();
     }
-    save(replicate);
 }
 
-void Simulation::save(int replicate) {
-    int fileVersion = 1;
-    std::ostringstream fileName;//(simName, std::ios_base::app);
-    fileName << simName << "_results_" << replicate << ".sim";
-    std::cout << "Saving " << fileName.str() << '\n';
-    Simfile simfile(fileName.str());
-    simfile.write(fileVersion);
-    simfile.write((int)seed);
+
+void Simulation::resume(std::string resultsFile, timeType resumeGen, std::string newResultsFile) {
+    // read all samples and checkpoints:
+    load(resultsFile);
+    if (resumeGen==-1) {
+        resumeGen = checkpointList.back()->getGeneration();
+    } else {
+        // Find the right checkpoint (they may not be saved at equal intervals):
+        int cpi = 0;
+        timeType gen = checkpointList[cpi]->getGeneration();
+        while (gen<resumeGen && cpi<checkpointList.size()-1) {
+            gen = checkpointList[++cpi]->getGeneration();
+        }
+        if (gen!=resumeGen) {
+            std::cout << "Error. Can't find checkpoint for generation " << resumeGen << " in file " << resultsFile << '\n';
+            exit(0);
+        }
+        // truncate lists:
+        checkpointList.resize(cpi+1);
+    }
+    int si = 0;
+    timeType gen = sampleList[si]->getGeneration();
+    while (gen<resumeGen && si<sampleList.size()-1) {
+        gen = sampleList[++si]->getGeneration();
+    }
+    if (gen>resumeGen) {
+        --si;
+    }
+    sampleList.resize(si+1);
+    thePop->resumeFromCheckpoint(*checkpointList.back());
+    randseed(checkpointList.back()->seed);
+    resultsFileName = newResultsFile;
+    runFromGeneration(resumeGen);
+}
+
+
+void Simulation::save() {
+    int fileVersion = 2; // As of v. 1.1
+    std::cout << "Saving " << resultsFileName << '\n';
+    oSimfile simfile(resultsFileName);
+    simfile.write<int>(fileVersion);
+    simfile.write<int>((int)seed);
     simfile.writeString(parameterFileCopy.str());
-    simfile.write(thePop->getGenetics().getLoci());
+    simfile.write<int>(thePop->getGenetics().getLoci());
     // save all samples:
-    simfile.write((int)sampleList.size());
-    for(int si=0; si<sampleList.size(); ++si) {
-        simfile << *sampleList.at(si);
+    simfile.write<int>((int)sampleList.size());
+    for(Sample* s : sampleList) {
+        thePop->writeSample(*s,simfile);
+    }
+    if (thePop->isTrackingGenes()) {
+        thePop->getGenetics().writeGeneLists(simfile);
+    }
+    if (checkpointInterval>0) {
+        simfile.write<int>((int)checkpointList.size());
+        for(Checkpoint* cp : checkpointList) {
+            thePop->writeCheckpoint(*cp, simfile);
+        }
+    }
+}
+
+// Load samples and checkpoints from results file. Don't change simulation parameters.
+void Simulation::load(std::string& resultsFileName) {
+
+    iSimfile resFile(resultsFileName);
+    int fileVersion = resFile.read<int>();
+    if( fileVersion!=2 ) {
+        std::cout << "Wrong file version of results file " << resultsFileName << '\n';
+        exit(0);
+    }
+    // read seed (not used):
+    resFile.read<int>();
+    // read (old) parameter-file copy (not used):
+    resFile.readString();
+    // read loci (not used):
+    resFile.read<int>();
+    
+    int sampleCount = resFile.read<int>();
+    sampleList.clear();
+    sampleList.reserve(sampleCount);
+    for (int si=0; si<sampleCount; ++si) {
+        sampleList.push_back(thePop->readSample(resFile));
     }
     if (thePop->isTrackingGenes()) {
-        thePop->getGenetics().saveGeneLists(simfile);
+        thePop->getGenetics().readGeneLists(resFile);
     }
+
+    if (checkpointInterval>0) {
+        int count = resFile.read<int>();
+        checkpointList.clear();
+        checkpointList.reserve(count);
+        for (int cpi=0; cpi<count; ++cpi) {
+            checkpointList.push_back(thePop->readCheckpoint(resFile));
+        }
+    }
+
 }
+
diff --git a/TREES_code/simulation.h b/TREES_code/simulation.h
index ce02851..ba1f1d6 100644
--- a/TREES_code/simulation.h
+++ b/TREES_code/simulation.h
@@ -36,6 +36,7 @@
 // and simulation output (samples).
 // It runs the simulation through a simple for-loop over time and
 // eventually saves all parameters and population samples to a 'sim'-file.
+// Checkpoints, if activated, are saved to the same file.
 ///////////////////////////////////////////////
 class Simulation {
 protected:
@@ -43,15 +44,22 @@ class Simulation {
     bool verbose;
     timeType tmax;
     timeType sampleInterval;
+    timeType checkpointInterval; // how often to save checkpoints
+    bool keepOldCheckpoints;
     std::string  simName;
+    std::string  resultsFileName;
     std::ostringstream  parameterFileCopy;
     Population* thePop;
     std::vector<Sample*> sampleList;
-    void save(int replicate);
-
+    std::vector<Checkpoint*> checkpointList;
+    void load(std::string& resultsFile);
+    void save();
+    void makeFileName(int replicate);
 public:
     Simulation(std::string& parameterFileName);
     ~Simulation();
-    void runAndSave(int iterations);
+    void initialize(int replicate);
+    void runFromGeneration(timeType gen);
+    void resume(std::string resultsFile, timeType gen, std::string newResultsFile);
 };
 #endif /* defined(__Species__simulation__) */
diff --git a/TREES_code/space.cpp b/TREES_code/space.cpp
index 83d634b..d70222e 100644
--- a/TREES_code/space.cpp
+++ b/TREES_code/space.cpp
@@ -48,6 +48,8 @@ void NullSpace::nextGeneration() {}
 void NullSpace::addChild(int mom, int dad) {}
 void NullSpace::compactData(std::vector<bool>& alive) {}
 void NullSpace::addToSample(Sample& theSample) {}
+void NullSpace::readToSample(Sample& theSample, iSimfile& isf, int n) {};
+int NullSpace::resumeFromCheckpoint(Checkpoint& cp, int dataIndex) {return dataIndex; }
 double NullSpace::getPosition(int individual, int dim) { return 0.0; }
 double NullSpace::getDist2(int ind1, int ind2) {return 0.0;}
 
@@ -135,12 +137,12 @@ double DiscreteSpaceImp::getDist2(int ind1, int ind2) { // squared distance
 }
 
 int DiscreteSpaceImp::popSizeInPatch(int linearPatch) {
-    if ((int)pop.getAge()!=generationLastCount) {
+    if (pop.getAge()!=generationLastCount) {
         patchPopSizes.assign(getPatchCount(), 0);
         for (int i=0; i<pop.size(); ++i) {
             patchPopSizes.at(getLinearPatch(i)) += 1;
         }
-        generationLastCount = (int)pop.getAge();
+        generationLastCount = pop.getAge();
     }
     return patchPopSizes.at(linearPatch);
 }
@@ -149,15 +151,7 @@ int DiscreteSpaceImp::popSizeInPatch(int linearPatch) {
 void DiscreteSpaceImp::initialize(int n0) {
     // Put everyone in initialPatch
     patches.assign(n0*Dims, initialPatch);
-    if (Dims==1) {
-        linearPatches = patches;
-    } else {
-        int p = 0;
-        for (int d=0; d<Dims; ++d) {
-            p = length*p + initialPatch;
-        }
-        linearPatches.assign(patches.size(), p);
-    }
+    assignLinearPatches();
 }
 
 void DiscreteSpaceImp::disperse() {
@@ -203,9 +197,11 @@ void DiscreteSpaceImp::disperse() {
                     default:
                         break;
                 }
+                linearPatches[i] = calcLinearPatch(i);
             }
-        }
-    }
+        } // for (int i=0; i<pop.size(); ++i) {
+        //assignLinearPatches();
+    } // if (length>1)
 }
 
 int DiscreteSpaceImp::treatBoundaries(int pos, int individual) {
@@ -245,16 +241,24 @@ void DiscreteSpaceImp::prepareNewGeneration(int size) {
 
 void DiscreteSpaceImp::nextGeneration() {
     patches = newPatches;
+    assignLinearPatches();
+}
+
+int DiscreteSpaceImp::calcLinearPatch(int individual) {
+    int p = 0;
+    for (int d=0; d<Dims; ++d) {
+        p = length*p + getPatch(individual, d);
+    }
+    return p;
+}
+
+void DiscreteSpaceImp::assignLinearPatches() {
     if (Dims==1) {
         linearPatches = patches;
     } else {
         linearPatches.assign(patches.size(), 0);
         for (int i=0; i<patches.size(); ++i) {
-            int p = 0;
-            for (int d=0; d<Dims; ++d) {
-                p = length*p + getPatch(i, d);
-            }
-            linearPatches[i] = p;
+            linearPatches[i] = calcLinearPatch(i);
         }
     }
 }
@@ -283,9 +287,19 @@ void DiscreteSpaceImp::compactData(std::vector<bool>& alive) {
 
 void DiscreteSpaceImp::addToSample(Sample& theSample) {
     // add sample spatial position:
-    theSample.addData(new IntData(patches));
+    theSample.addData(new XData<int>(patches));
+}
+void DiscreteSpaceImp::readToSample(Sample& theSample, iSimfile& isf, int n) {
+    theSample.addData(new XData<int>(isf,n));
 }
 
+int DiscreteSpaceImp::resumeFromCheckpoint(Checkpoint &cp, int dataIndex) {
+    patches = (dynamic_cast<XData<int>&>(cp.getData(dataIndex++))).getData();
+    patchPopSizes.assign(getPatchCount(), 0);
+    generationLastCount = -1;
+    assignLinearPatches();
+    return dataIndex;
+}
 
 ContinuousSpace::ContinuousSpace(Population& p, ParameterFile& pf) :
 Space(p) {
@@ -322,7 +336,7 @@ ContinuousSpace::~ContinuousSpace() {
 
 bool ContinuousSpace::isDiscrete() { return false; }
 
-double& ContinuousSpace::getCoord(int indiviudal, int dim) {
+ContinuousSpace::positionType& ContinuousSpace::getCoord(int indiviudal, int dim) {
     return pos[indiviudal*Dims + dim];
 }
 
@@ -335,9 +349,9 @@ double ContinuousSpace::getDist2(int ind1, int ind2) { // squared distance
         double dx = getCoord(ind1, 0)-getCoord(ind2, 0);
         return dx*dx;
     } else {
-        double dist2=0.0;
-        double *p1 = &getCoord(ind1, 0);
-        double *p2 = &getCoord(ind2, 0);
+        double dist2 = 0.0;
+        positionType *p1 = &getCoord(ind1, 0);
+        positionType *p2 = &getCoord(ind2, 0);
         for (int d=0; d<Dims; ++d) {
             double dx = getDistance(*p1, *p2);
             dist2 += dx*dx;
@@ -363,7 +377,7 @@ void ContinuousSpace::initialize(int n0) {
     pos.assign(n0*Dims, initialPosition);
 }
 
-double ContinuousSpace::treatBoundaries(double pos, int individual) {
+ContinuousSpace::positionType ContinuousSpace::treatBoundaries(ContinuousSpace::positionType pos, int individual) {
     if (pos>maxPos || pos<0) {
         switch (boundaryCondition) {
             case Reflective:
@@ -453,6 +467,15 @@ void ContinuousSpace::compactData(std::vector<bool>& alive) {
 
 void ContinuousSpace::addToSample(Sample& theSample) {
     // add sample spatial position:
-    theSample.addData(new FloatData(pos));
+    theSample.addData(new XData<float>(pos));
+}
+
+void ContinuousSpace::readToSample(Sample& theSample, iSimfile& isf, int n) {
+    theSample.addData(new XData<float>(isf,n*Dims));
+}
+
+int ContinuousSpace::resumeFromCheckpoint(Checkpoint &cp, int dataIndex) {
+    pos = (dynamic_cast<XData<float>&>(cp.getData(dataIndex++))).getData();
+    return dataIndex;
 }
 
diff --git a/TREES_code/space.h b/TREES_code/space.h
index f199f71..ac07302 100644
--- a/TREES_code/space.h
+++ b/TREES_code/space.h
@@ -55,6 +55,8 @@ class Space {
     virtual void addChild(int mom, int dad)=0;
     virtual void compactData(std::vector<bool>& alive)=0;
     virtual void addToSample(Sample& theSample)=0;
+    virtual void readToSample(Sample& theSample, iSimfile& isf, int n)=0;
+    virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex)=0;
     virtual double getPosition(int individual, int dimension)=0;
     virtual double getDist2(int ind1, int ind2)=0; // squared distance
     virtual bool isDiscrete()=0;
@@ -77,6 +79,8 @@ class DiscreteSpace : public Space {
     virtual void addChild(int mom, int dad)=0;
     virtual void compactData(std::vector<bool>& alive)=0;
     virtual void addToSample(Sample& theSample)=0;
+    virtual void readToSample(Sample& theSample, iSimfile& isf, int n)=0;
+    virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex)=0;
     virtual double getPosition(int individual, int dim)=0;
     virtual double getDist2(int ind1, int ind2)=0; // squared distance
     virtual bool isDiscrete();
@@ -89,7 +93,7 @@ class DiscreteSpaceImp : public DiscreteSpace {
     std::vector<int> newPatches;
     std::vector<int> patchPopSizes;
     std::vector<int> linearPatches;
-    int generationLastCount;
+    timeType generationLastCount;
     int length; // size of space = length^Dims
     int initialPatch; // starting patch for all individuals (in all dimensions)
     Trait* PDisp;
@@ -114,6 +118,8 @@ class DiscreteSpaceImp : public DiscreteSpace {
             return dx;
         }
     }
+    int calcLinearPatch(int individual);
+    void assignLinearPatches();
 public:
     DiscreteSpaceImp(Population& p, ParameterFile& pf);
     virtual ~DiscreteSpaceImp();
@@ -130,6 +136,8 @@ class DiscreteSpaceImp : public DiscreteSpace {
     virtual void addChild(int mom, int dad);
     virtual void compactData(std::vector<bool>& alive);
     virtual void addToSample(Sample& theSample);
+    virtual void readToSample(Sample& theSample, iSimfile& isf, int n);
+    virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex);
     virtual double getPosition(int individual, int dim);
     virtual double getDist2(int ind1, int ind2); // squared distance
 };
@@ -153,6 +161,8 @@ class NullSpace : public DiscreteSpace {
     virtual void addChild(int mom, int dad);
     virtual void compactData(std::vector<bool>& alive);
     virtual void addToSample(Sample& theSample);
+    virtual void readToSample(Sample& theSample, iSimfile& isf, int n);
+    virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex);
     virtual double getPosition(int individual, int dimension);
     virtual double getDist2(int ind1, int ind2); // squared distance
 };
@@ -160,18 +170,19 @@ class NullSpace : public DiscreteSpace {
 
 class ContinuousSpace : public Space {
 protected:
-    std::vector<double> pos;
-    std::vector<double> newPos;
+    typedef float positionType;
+    std::vector<positionType> pos;
+    std::vector<positionType> newPos;
     double initialPosition;  // starting position for all individuals (in all dimensions)
     Trait* PDisp;
-    double dispdist; // mean dispersal distance
-    double maxPos; // size in each dimension
+    positionType dispdist; // mean dispersal distance
+    positionType maxPos; // size in each dimension
     enum Boundary {Circular, Reflective, Absorbing};
     Boundary boundaryCondition;
-    double treatBoundaries(double pos, int individual);
+    positionType treatBoundaries(positionType pos, int individual);
     // inline for speed, may be negative:
-    inline double getDistance(double pos1, double pos2) {
-        double dx = pos1-pos2;
+    inline positionType getDistance(positionType pos1, positionType pos2) {
+        positionType dx = pos1-pos2;
         if (boundaryCondition==Circular) {
             if (dx>maxPos/2) {
                 return maxPos-dx;
@@ -187,7 +198,7 @@ class ContinuousSpace : public Space {
     ContinuousSpace(Population& p, ParameterFile& pf);
     virtual ~ContinuousSpace();
     virtual double getPosition(int individual, int dim);
-    double& getCoord(int individual, int dim); // position with reference
+    positionType& getCoord(int individual, int dim); // position with reference
     virtual void initialize(int n0);
     virtual void disperse();
     virtual void prepareNewGeneration(int size);
@@ -195,6 +206,8 @@ class ContinuousSpace : public Space {
     virtual void addChild(int mom, int dad);
     virtual void compactData(std::vector<bool>& alive);
     virtual void addToSample(Sample& theSample);
+    virtual void readToSample(Sample& theSample, iSimfile& isf, int n);
+    virtual int resumeFromCheckpoint(Checkpoint& cp, int dataIndex);
     virtual double getDist2(int ind1, int ind2); // squared distance
     virtual bool isDiscrete();
 };