LoadPDFgetNFile.cpp

/*WIKI*


Loads a PDFgetN file such as that saved by PDFgetN.
  
PDFgetN generates many types of files for the final result and communication among different executables at internal steps. This algorithm is designed to recognize and load all the ASCII based files created by PDFgetN.  

The file types include 
  * .ain
  * .braw
  * .bsmo
  * .sq
  * .sqa
  * .gr


*WIKI*/
#include "MantidDataHandling/LoadPDFgetNFile.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidAPI/LoadAlgorithmFactory.h"

#include <fstream>
#include <boost/algorithm/string.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/trim.hpp>

using namespace Mantid;
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace std;
using namespace boost;

// FIXME  Add label and unit to output workspace
// FIXME  Consider to output multiple workspaces if there are multiple column data (X, Y1, E1, Y2, E2)

namespace Mantid
{
namespace DataHandling
{

  DECLARE_ALGORITHM(LoadPDFgetNFile)

  // register the algorithm into loadalgorithm factory
  DECLARE_LOADALGORITHM(LoadPDFgetNFile)

  //----------------------------------------------------------------------------------------------
  /** Constructor
   */
  LoadPDFgetNFile::LoadPDFgetNFile()
  {
  }
    
  //----------------------------------------------------------------------------------------------
  /** Destructor
   */
  LoadPDFgetNFile::~LoadPDFgetNFile()
  {
  }

  //----------------------------------------------------------------------------------------------
  /** Init documentation
    */
  void LoadPDFgetNFile::initDocs()
  {
    this->setWikiSummary("Loads a PDFgetN file, which is generated by PDFgetN.");
    this->setOptionalMessage("Types of PDFgetN data files include .sqa, .sq, .gr, and etc.");
  }

  //----------------------------------------------------------------------------------------------
  /** This method does a quick file type check by checking the first 100 bytes of the file
   *  @param filePath- path of the file including name.
   *  @param nread :: no.of bytes read
   *  @param header :: The first 100 bytes of the file as a union
   *  @return true if the given file is of type which can be loaded by this algorithm
   */
  bool LoadPDFgetNFile::quickFileCheck(const std::string& filePath, size_t nread, const file_header& header)
  {
    UNUSED_ARG(nread);
    UNUSED_ARG(header);

    // check the file extension
    std::string extn = extension(filePath);
    bool bascii;
    if (extn.compare("sq"))
      bascii = true;
    else if (extn.compare("sqa"))
      bascii = true;
    else if (extn.compare("sqb"))
      bascii = true;
    else if (extn.compare("gr"))
      bascii = true;
    else if (extn.compare("ain"))
      bascii = true;
    else if (extn.compare("braw"))
      bascii = true;
    else if (extn.compare("bsmo"))
      bascii = true;
    else
      bascii = false;

    /* check the bit of header
    bool is_ascii(true);
    for (size_t i = 0; i < nread; i++)
    {
      if (!isascii(header.full_hdr[i]))
        is_ascii = false;
    }
     return (is_ascii || bascii);
    */

    return (bascii);
  }

  /** checks the file by opening it and reading few lines
   *  @param filePath :: name of the file including its path
   *  @return an integer value how much this algorithm can load the file
   */
  int LoadPDFgetNFile::fileCheck(const std::string& filePath)
  {
    std::ifstream file(filePath.c_str());
    if (!file)
    {
      g_log.error("Unable to open file: " + filePath);
      throw Exception::FileError("Unable to open file: ", filePath);
    }
    std::string str;
    getline(file, str);//workspace title first line
    while (!file.eof())
    {
      getline(file, str);
      if (startsWith(str, "#L"))
      {
        return 80;
      }
    }

    return 0;
  }


  //----------------------------------------------------------------------------------------------
  /** Define input
    */
  void LoadPDFgetNFile::init()
  {
    std::vector<std::string> exts;
    exts.push_back(".sq");
    exts.push_back(".sqa");
    exts.push_back(".sqb");
    exts.push_back(".gr");
    exts.push_back(".ain");
    exts.push_back(".braw");
    exts.push_back(".bsmo");

    auto fileproperty = new FileProperty("Filename", "", FileProperty::Load, exts, Kernel::Direction::Input);
    this->declareProperty(fileproperty, "The input filename of the stored data");

    // auto wsproperty = new WorkspaceProperty<Workspace2D>("OutputWorkspace", "Anonymous", Kernel::Direction::Output);
    // this->declareProperty(wsproperty, "Name of output workspace. ");
    declareProperty(new API::WorkspaceProperty<>("OutputWorkspace", "", Kernel::Direction::Output),
                    "Workspace name to load into.");
  }

  //----------------------------------------------------------------------------------------------
  /** Main executor
    */
  void LoadPDFgetNFile::exec()
  {
    // 1. Parse input file
    std::string inpfilename = getProperty("Filename");

    parseDataFile(inpfilename);

    // 2. Generate output workspace
    generateDataWorkspace();

    setProperty("OutputWorkspace", outWS);

    return;
  }
  
  //----------------------------------------------------------------------------------------------
  /** Parse PDFgetN data file to
    * 1. a 2D vector for column data
    * 2. a 1D string vector for column name
    */
  void LoadPDFgetNFile::parseDataFile(std::string filename)
  {
    // 1. Open file
    std::ifstream ifile;
    ifile.open((filename.c_str()));

    if (!ifile.is_open())
    {
      stringstream errmsg;
      errmsg << "Unable to open file " << filename << ".  Quit!";
      g_log.error() << errmsg.str() << std::endl;
      throw std::runtime_error(errmsg.str());
    }
    else
    {
      g_log.notice() << "Open PDFgetN File " << filename << std::endl;
    }

    // 2. Parse
    bool readdata = false;

    char line[256];
    while(ifile.getline(line, 256))
    {
      string sline(line);

      if (!readdata && startsWith(sline, "#L"))
      {
        // a) Find header line for the data segment in the file
        parseColumnNameLine(sline);
        readdata = true;

        // Set up the data structure
        size_t numcols = mColumnNames.size();
        for (size_t i = 0; i < numcols; ++i)
        {
          std::vector<double> tempvec;
          mData.push_back(tempvec);
        }

      }
      else if (readdata)
      {
        // b) Parse data
        parseDataLine(sline);
      }
      else
      {
        // c) Do nothing
        ;
      }
    } // ENDWHILE

    if (!readdata)
    {
      stringstream errmsg;
      errmsg << "Unable to find a line staring with #L as the indicator of data segment. ";
      g_log.error() << errmsg.str() << std::endl;
      throw std::runtime_error(errmsg.str());
    }

    return;
  }

  //----------------------------------------------------------------------------------------------
  /** Check whether the line starts with some specific character
    */
  bool LoadPDFgetNFile::startsWith(std::string s, std::string header)
  {
    bool answer = true;

    if (s.size() < header.size())
    {
      answer = false;
    }
    else
    {
      size_t numchars = header.size();
      for (size_t i = 0; i < numchars; ++i)
      {
        char c0 = s[i];
        char c1 = header[i];
        if (c0 != c1)
        {
          answer = false;
          break;
        }
      }
    }

    return answer;
  }

  /** Parse column name line staring with \#L
    */
  void LoadPDFgetNFile::parseColumnNameLine(std::string line)
  {
    // 1. Split
    vector<string> terms;
    boost::split(terms, line, is_any_of(" \t\n"), token_compress_on);

    // 2. Validate
    if (terms.empty())
    {
      throw std::runtime_error("There is nothing in the input line!");
    }

    string header = terms[0];
    if (header.compare("#L") != 0)
    {
      stringstream errmsg;
      errmsg << "Expecting header as #L.  Input line has header as " << header
             << ". Unable to proceed. ";
      g_log.error() << errmsg.str() << endl;
      throw std::runtime_error(errmsg.str());
    }

    // 3. Parse
    size_t numcols = terms.size()-1;
    stringstream msgss;
    msgss << "Column Names: ";
    for (size_t i = 0; i < numcols; ++i)
    {
      this->mColumnNames.push_back(terms[i+1]);
      msgss << setw(-3) << i << ": " << setw(-10) << mColumnNames[i];
    }
    g_log.information() << msgss.str() << endl;

    return;
  }

  /** Parse data line
    */
  void LoadPDFgetNFile::parseDataLine(string line)
  {
    // 1. Trim (stripg) and Split
    boost::trim(line);
    vector<string> terms;
    boost::split(terms, line, is_any_of(" \t\n"), token_compress_on);

    // 2. Validate
    size_t numcols = mData.size();
    if (line[0] == '#')
    {
      // Comment/information line to indicate the start of another section of data
      return;
    }
    else if (terms.size() != numcols)
    {
      // Data line with incorrect number of columns
      stringstream warnss;
      warnss <<  "Line (" << line << ") has incorrect number of columns other than " << numcols << " as expected. ";
      g_log.warning(warnss.str());
      return;
    }

    // 3. Parse
    for (size_t i = 0; i < numcols; ++i)
    {
      string temps = terms[i];
      double tempvalue;
      if (temps.compare("NaN") == 0)
      {
        // FIXME:  Need to discuss with Peter about how to treat NaN value
        //tempvalue = DBL_MAX-1.0;
        tempvalue = 0.0;
      }
      else if (temps.compare("-NaN") == 0)
      {
        //tempvalue = -DBL_MAX+1.0;
        // FIXME:  Need to discuss with Peter about how to treat NaN value
        tempvalue = 0.0;
      }
      else
      {
        tempvalue = atof(temps.c_str());
      }

      mData[i].push_back(tempvalue);
    }

    return;
  }

  //----------------------------------------------------------------------------------------------
  void LoadPDFgetNFile::setUnit(Workspace2D_sptr ws)
  {
    // 1. Set X
    string xcolname = mColumnNames[0];

    if (xcolname.compare("Q") == 0)
    {
      string unit = "MomentumTransfer";
      ws->getAxis(0)->setUnit(unit);
    }
    else if (xcolname.compare("r") == 0)
    {
      ws->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
      Unit_sptr unit = ws->getAxis(0)->unit();
      boost::shared_ptr<Units::Label> label = boost::dynamic_pointer_cast<Units::Label>(unit);
      label->setLabel("AtomicDistance", "Angstrom");
    }
    else
    {
      stringstream errss;
      errss << "X axis " << xcolname << " is not supported for unit. " << endl;
      g_log.warning() << errss.str() << endl;
    }

    // 2. Set Y
    string ycolname = mColumnNames[1];
    string ylabel("");
    if (ycolname.compare("G(r)") == 0)
    {
      ylabel = "PDF";
    }
    else if (ycolname.compare("S") == 0)
    {
      ylabel = "S";
    }
    else
    {
      ylabel = "Intensity";
    }
    ws->setYUnitLabel(ylabel);

    return;
  }

  /** Generate output data workspace
    * Assumption.  One data set must contain more than 1 element.
    */
  void LoadPDFgetNFile::generateDataWorkspace()
  {
    // 0. Check
    if (mData.size() == 0)
    {
      throw runtime_error("Data set has not been initialized. Quit!");
    }

    // 1. Figure out direction of X and number of data set
    bool xascend = true;
    if (mData.size() >= 2 && mData[0][1] < mData[0][0])
    {
      xascend = false;
    }

    size_t numsets = 0;
    vector<size_t> numptsvec;
    size_t arraysize = mData[0].size();
    if (arraysize <= 1)
    {
      throw runtime_error("Number of columns in data is less and equal to 1.  It is unphysically too small.");
    }

    double prex = mData[0][0];
    size_t vecsize = 1;
    for (size_t i = 1; i < arraysize; ++i)
    {
      double curx = mData[0][i];
      if ( ((xascend) && (curx < prex)) || ((!xascend) && (curx > prex)) )
      {
        // X in ascending order and hit the end of one set of data
        // X in descending order and hit the end of one set of data
        // Record the current data set information and start the next data set
        numsets += 1;
        numptsvec.push_back(vecsize);
        vecsize = 1;
      }
      else
      {
        // In the middle of a set of data
        ++ vecsize;
      }
      // Loop variable udpate
      prex = curx;
    } // ENDFOR
    // Record the last data set information
    ++ numsets;
    numptsvec.push_back(vecsize);

    bool samesize = true;
    for (size_t i = 0; i < numsets; ++i)
    {
      if (i > 0)
      {
        if (numptsvec[i] != numptsvec[i-1])
        {
          samesize = false;
        }
      }
      g_log.information() << "Set " << i << ":  Number of Points = " << numptsvec[i] << std::endl;
    }
    if (!samesize)
    {
      stringstream errmsg;
      errmsg << "Multiple bank (number of banks = " << numsets << ") have different size of data array.  Unable to handle this situation.";
      g_log.error() << errmsg.str() << std::endl;
      throw std::runtime_error(errmsg.str());
    }
    size_t size = numptsvec[0];

    // 2. Generate workspace2D object and set the unit
    outWS = boost::dynamic_pointer_cast<Workspace2D>(
          API::WorkspaceFactory::Instance().create("Workspace2D", numsets, size, size));

    setUnit(outWS);

    // 3. Set number
    size_t numspec = outWS->getNumberHistograms();
    for (size_t i = 0; i < numspec; ++i)
    {
      MantidVec& X = outWS->dataX(i);
      MantidVec& Y = outWS->dataY(i);
      MantidVec& E = outWS->dataE(i);

      size_t baseindex = i*size;
      for (size_t j = 0; j < size; ++j)
      {
        size_t index;
        if (xascend)
          index = j;
        else
          index = (size-1)-j;

        X[index] = mData[0][baseindex+j];
        Y[index] = mData[1][baseindex+j];
        E[index] = mData[2][baseindex+j];
      }
    }

    return;
  }

} // namespace DataHandling
} // namespace Mantid