autocorrelate.cpp

Go to the documentation of this file.

/* 
* autocorrelate.cpp - Autocorrelation function implementation for Atom Probe
* Copyright (C) 2015  Daniel Haley
* 
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 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, see <http://www.gnu.org/licenses/>.
*/

#include <iostream>
#include <fstream>
#include <limits>

#include "atomprobe/atomprobe.h"

#ifdef _OPENMP
#include <omp.h>
#endif
using namespace AtomProbe;
using namespace std;

#ifdef DEBUG
#include <cassert>
#define ASSERT(f) {assert(f);}
#else
#define ASSERT(f)  {}
#endif

const unsigned int PROGRESS_BAR_SIZE=50;

ProgressBar *treeProgressBar=0;
unsigned int treeProgressValue=0;

bool callback()
{
    if(!treeProgressBar)
        return true;

#ifdef _OPENMP
    if(!omp_get_thread_num())
#endif
    
        treeProgressBar->update(treeProgressValue);

    return true;
}

//Input must be in vector[bin][ions] form
template<class T>
void dumpAutocorrelation(const vector<vector<T> > &autoCorrelateFunc, 
            const RangeFile &rangeFile, float binStep, std::ostream &strm )
{
    strm << "Per-Species Autocorrelation function:" << endl;
    strm << "r " << "\t";
    for(unsigned int uj=0;uj<rangeFile.getNumIons();uj++)
        strm << rangeFile.getName(uj) << "\t" ;
    strm << endl;

    if(!autoCorrelateFunc.size())
        return;

    ASSERT(rangeFile.getNumIons() == autoCorrelateFunc[0].size());

    for(unsigned int ui=0;ui<autoCorrelateFunc.size();ui++)
    {
        strm << binStep*ui << "\t";

        for(unsigned int uj=0;uj<autoCorrelateFunc[ui].size();uj++)
        {
            strm <<autoCorrelateFunc[ui][uj] << "\t";
        }
        strm << endl;
    }
    
}

int main(int argc, char *argv[])
{
    if(!(argc == 6  || argc ==7) )
    {
        cerr << "USAGE : " << argv[0] << " POSFILE RANGEFILE Radius binStep Sampling [outputfile]" << endl;
        return 1;
    }

    std::string posFilename,rangeFilename,  outputFile;
    float sampleFactor,scanRadius,binStep;

    posFilename=argv[1];
    rangeFilename=argv[2];
    if(argc ==7)
        outputFile=argv[6];

    if(stream_cast(sampleFactor,argv[5]))
    {
        cerr << "unable to understand sampling factor. Aborting." << endl;
        return 1;
    }
    if(stream_cast(scanRadius,argv[3]))
    {
        cerr << "unable to understand scan radius. Aborting." << endl;
        return 1;
    }   
    
    if(stream_cast(binStep,argv[4]))
    {
        cerr << "unable to understand bin step size. Aborting." << endl;
        return 1;
    }   

    if(sampleFactor < 0.0f || sampleFactor > 1.0f)
    {
        cerr << "sampling factor must lie in range [0,1]"<<endl ;
        return 1;
    }

    if(scanRadius <=0.0f)
    {
        cerr << "Scan radius must be positive" << endl;
        return 1;
    }

    if(binStep <=0.0f)
    {
        cerr << "Bin step size must be positive" << endl;
        return 1;
    }

    std::ofstream of;
    if(!outputFile.empty())
    {
        of.open(outputFile.c_str());
        if(!of)
        {
            cerr << "Unable to open outputfile :" << outputFile << endl;
            return 1;
        }
    }
    //Load ion data
    //--
    unsigned int errCode;
    vector<IonHit> ions;
    cerr << "Loading pos data from " <<posFilename << "..." << endl;
    errCode=loadPosFile(ions,posFilename.c_str());  
    
    if(errCode)
    {
        cerr << "Error loading posfile :" << endl;  
        cerr << getPosFileErrString(errCode) << endl;
        return 1;
    }

    cerr << "\tLoaded :" << ions.size() << " ions." << endl;
    //--

    //load range data
    //--
    RangeFile rangeFile;
    if(!rangeFile.openGuessFormat(rangeFilename.c_str()))
    {
        cerr << "Error loading rangefile :" <<
             rangeFile.getErrString() << endl;
        return 1;
    }   
    //--

    //Apply the rangefile
    //--
    cerr << "Ranging...";
    rangeFile.range(ions);
    cerr << "Done" << endl;
    cerr << "\tRanged:" << ions.size() << " ions." << endl;

    if(ions.empty())
    {
        cerr << "No ions after ranging. Aborting." << endl; 
        return 1;
    }
    //--


    //Sample the data
    //--
    if(sampleFactor < 1.0f)
    {
        {
        cerr << "Sampling...";
        vector<IonHit> sampled;
        sampleIons(ions,sampleFactor,sampled);
        sampled.swap(ions);
        }
        cerr << "Done" << endl;
        cerr << "\tSampled:" << ions.size() << " ions." << endl;

        if(ions.empty())
        {
            cerr << "No ions after sampling. Aborting." << endl;    
            return 1;
        }
    }
    //--
    
            
    K3DTreeExact tree;
    BoundCube b;
    //Initialise the KD tree
    //--
    cerr << "Building search structures....";   
    treeProgressBar= new ProgressBar();
    treeProgressBar->setLength(PROGRESS_BAR_SIZE);
    treeProgressBar->init();

    tree.setCallback(callback);
    tree.setProgressPointer(&treeProgressValue);
    tree.resetPts(ions,false);
    tree.build();
    delete treeProgressBar;
    treeProgressBar=0;
    
    tree.getBoundCube(b);

    if(b.isFlat())
    {
        cerr << "Data structure not OK. Aborting" << endl;
        cerr << "Does the input data form a volume?" << endl;
        return 1;
    }
    //--


    unsigned int binCount = scanRadius/binStep;

    if(!binCount)
    {
        cerr << "Search radius :" << scanRadius << " is too small, needs to be bigger than the bin size:" << binStep << endl;
        return 1;
    }

    //First dimension is bin count, second is composition

    cerr << "Computing autocorrelation..." ;
    ProgressBar progressBar;
    progressBar.setLength(PROGRESS_BAR_SIZE);
    progressBar.init(); 


    unsigned int numCounted =0;

    vector<vector<unsigned int> > overallCompositionTable;
    //Zero final composition table
    //--
    overallCompositionTable.resize(binCount);   
    for(unsigned int ui=0;ui<overallCompositionTable.size();ui++)
        overallCompositionTable[ui].resize(rangeFile.getNumIons(),0);
    //--

    

    //Perform composition scan to find atom counts as a function of radius and atom type
    //--
#pragma omp parallel for
    for(unsigned int ui=0;ui<ions.size();ui++)
    {
        vector<size_t> pointIdx;
        vector<vector<size_t> > compositionTable;

        //Zero composition table
        compositionTable.resize(binCount);  
        for(unsigned int uj=0;uj<compositionTable.size();uj++)
            compositionTable[uj].resize(rangeFile.getNumIons(),0);

        //Find all points in radius
        tree.findUntaggedInRadius(ions[ui].getPos(),b,scanRadius, pointIdx);
        numCounted++;
#ifdef _OPENMP
        if(!omp_get_thread_num())
#endif
            progressBar.update((float)numCounted/ions.size()*100.0f);

    

        if(pointIdx.empty())
            continue;

        //Compute local composition frequency table. 
        // This is essentially the building blocks for an RDF
        for(size_t uj=0;uj<pointIdx.size();uj++)
        {   
            unsigned int ionId;
            float mass,radius;
            //target ion
            IonHit p;   
            p = ions[tree.getOrigIndex(pointIdx[uj])];
            mass=p.getMassToCharge();
            ionId = rangeFile.getIonID(mass);
            ASSERT(ionId < rangeFile.getNumIons());

            //compute distance between ions
            radius = sqrt(p.getPos().sqrDist(ions[ui].getPos()));

            //Do not record points that are at the same location
            if(radius <=sqrt(std::numeric_limits<float>::epsilon()))
                continue;

            unsigned int ionBin;
            ionBin = (unsigned int) (radius/binStep);
            if(ionBin == binCount && ionBin)
                ionBin--;
        
            compositionTable[ionBin][ionId]++;
        }

        #pragma omp critical
        {
            for(size_t uj=0;uj<compositionTable.size();uj++)
            {
                for(size_t uk=0;uk<compositionTable[uj].size();uk++)
                {
                    overallCompositionTable[uj][uk]+=compositionTable[uj][uk];
                }
            }
        }

    }
    //--

    //force finalisation  of progress bar
    progressBar.update(100);

    //Convert the raw count data to composition frequency
    //First entry is bin, second is ion ID
    vector<vector<float> > compPct;

    compPct.resize(binCount);
    for(unsigned int ui=0;ui<binCount;ui++)
        computeComposition(overallCompositionTable[ui],compPct[ui]);


    //switch from vector[bin][iontype] to vector[iontype][bin]
    AtomProbe::transposeVector(compPct);

    //Compute mean composition, per ion 
    vector<float> meanComp;
    meanComp.resize(rangeFile.getNumIons());    

    //loop over ion types
    for(unsigned int ui=0;ui<compPct.size();ui++)
    {
        double res;
        res=0;
        //loop over radius
        for(unsigned int uj=0;uj<compPct[ui].size();uj++)
            res+=compPct[ui][uj];
        
        res/=compPct[ui].size();
        meanComp[ui]=res;
    }   

    //Zero-centre the concentration data
    //loop over ions
    for(unsigned int ui=0;ui<compPct.size();ui++)
    {
        //loop over radius
        for(unsigned int uj=0;uj<compPct[ui].size();uj++)
            compPct[ui][uj]-=meanComp[ui];
    }

    vector<vector<float> > &deviationComp=compPct;

    //Compute the auto correlation function, per species
    //---

    vector<vector<float> > autoCorrelateFunc;
    vector<float> denominator;
    autoCorrelateFunc.resize(binCount);
    //Compute denominator
    unsigned int numIons=rangeFile.getNumIons();
    denominator.resize(numIons);
    for(unsigned int ui=0;ui<numIons;ui++)
    {
        double normConstant;
        normConstant=0;
        for(unsigned int uj=0;uj<deviationComp[ui].size();uj++)
            normConstant+=deviationComp[ui][uj]*deviationComp[ui][uj];

        //double->float conversion
        denominator[ui] = normConstant;
    }

    //Compute numerator. Note that in "Atom Probe Field Ion Microscopy"
    // ISBN : 0198513879 ; pp 322,  Equation 5.64 does not really handle
    // boundary conditions, but instead reduces the autocorrelation lag maximum
    autoCorrelateFunc.resize(numIons);
    for(unsigned int ui=0;ui<numIons;ui++)
    {
        autoCorrelateFunc[ui].resize(binCount);
        // Ki is the index to k in equation 5.64
        for(unsigned int ki=0; ki<binCount; ki++)
        {
            double numValue;
            numValue=0;
            for(unsigned int uj=0;uj<binCount-ki; uj++)
            {
                numValue+=deviationComp[ui][uj]*deviationComp[ui][uj+ki];
            }

            if(denominator[ui] > std::numeric_limits<float>::epsilon())
                autoCorrelateFunc[ui][ki] = numValue / denominator[ui];
            else
                autoCorrelateFunc[ui][ki] = 0;
            
        }

    
    }   

    //Autocorrelation vector in form vector[ion][bin],
    // so transpose to vector[bin][ion]
    transposeVector(autoCorrelateFunc);

    if(outputFile.empty())
        dumpAutocorrelation(autoCorrelateFunc,rangeFile,binStep,cerr);
    else
    {
        cerr << "Wrote output to:" << outputFile << endl;
        dumpAutocorrelation(autoCorrelateFunc,rangeFile,binStep,of);
        std::ofstream debugOf;
        std::string s;
        s= outputFile + "-debugrdf";
        debugOf.open(s.c_str());
        dumpAutocorrelation(overallCompositionTable,rangeFile,binStep,debugOf);
    }
        

}