QGpCompatibility/CompatInversionReport.h

Go to the documentation of this file.

/***************************************************************************
**
**  This file is part of QGpCompatibility.
**
**  This file may be distributed and/or modified under the terms of the
**  GNU General Public License version 2 or 3 as published by the Free
**  Software Foundation and appearing in the file LICENSE.GPL included
**  in the packaging of this file.
**
**  This file 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/>.
**
**  See http://www.geopsy.org for more information.
**
**  Created : 2003-05-22
**  Authors:
**    Marc Wathelet
**    Marc Wathelet (ULg, Liège, Belgium)
**    Marc Wathelet (LGIT, Grenoble, France)
**
***************************************************************************/

#ifndef INVERSIONREPORT_H
#define INVERSIONREPORT_H

#include <QGpCoreWave.h>

#include "QGpCompatibilityDLLExport.h"

namespace QGpCompatibility {

class CompatMultiModalFrequency;
class CompatEllipticity;
class CompatEllipticityData;
class CompatMultiModalCurves;
class CompatDispersion;
class CompatDispersionData;

class QGPCOMPATIBILITY_EXPORT CompatInversionReport
{
  //-----------------------------------------------------------------------------------------------------------------------------
  //
  //  1st part: general
  //
  //-----------------------------------------------------------------------------------------------------------------------------
public:
  // Constructor: open file
  CompatInversionReport(bool isWrite, QString reportFile, QString naInFile=QString::null, int nModels=0);
  // Destructor: close filecd ../
  virtual ~CompatInversionReport();
  // Static function to test if report is ready to be opened
  static bool isReady(QString reportFile);
  // Test if file is an inversion report
  static bool isRightFormat(QString reportFile);
  bool isValid() {return _valid;}
  // Test whether this class is of type className
  virtual bool isA(const char * className);
  // Name of this class
  const char * className() {return "CompatInversionReport";}
  // Test format flag at beginning of file
  bool isRightFormat();
  // Checks if two report have the same parameters
  bool isSameParameter(const CompatInversionReport * otherModels) const;
  
  void writeHeader();
  void setNFromNaIn(QString naInfile);
  // Block types to allow block testing
  enum BlockType {Parameters, Omegas, DispersionGoal, ModelDispersion, AutocorrGoal,
                  EllipticityGoal, FKMaxGoal, Information, Unknown, RefraGoal};
  //
  enum paramID {H, Vs, Vp, Rho, Depth, ExpGrad, ModeStdDev, NAParam};
  // Return the block type of block n�2
  BlockType reportType();
  void addParameter(paramID which,int layer);
  void addEllipticityGoal(CompatEllipticityData * ell);
  void modifyEllipticityGoal(CompatEllipticityData * ell);
  // Modify the cost value for model i 
  void setCostValue(int imodel,double costValue);
  
  bool startWritingBlock(int index,BlockType wantedType);
  bool startWritingBlock(BlockType type);
  void endWritingBlock();
  bool startReadingBlock(int index,BlockType wantedType);
  /* Flush the buffers (usefull to call when reading is blocked)
    Introduced to solve a bug noticed 2003-02-10 */
  void flush();
  // Change the flushed state
  void setFlushed(bool s) {_flushed=s;}
  // Load model (and its cost function) given by its index and makes it the current one
  bool loadModel(int modelIndex=-1);
  void loadNAModel(int modelIndex=-1);
  // Load next model (and its cost function) given by its index and makes it the current one
  bool loadNextModel();
  // Get only the cost value from a model
  double loadCostValue(int modelIndex);
  // load the measured ellipticity curve into a new data class (may be null if there's no ellipticity)
  CompatEllipticityData * loadEllipticityGoal();
  // Returns the name of a parameter given by its index (don't forget to delete char)
  QString parameterName(int paramIndex);
  // Returns the value of parameter given by its index in the current model
  double getParamValue(int paramIndex);
  paramID getParamType(int paramIndex) {return _paramAddresses[paramIndex].which;}
  // Returns the current model
  Seismic1DModel *  currentModel();
  // Returns the current NA model
  float *  currentNAModel();
  // Returns the value of the cost function for the current model
  double currentCost() {return _currentCost;}
  // Returns the value of the ID for the current model
  int getID() {return _currentID;}
  /*
     Returns the curve of best fit point versus num models from index startIndex
     startIndex is modified and will contain the index of the first model not checked
     by getBestFit
  */
  QVector<Point> * getBestFits(QVector<Point>& curve,
                                       int& startIndex,double& minCostValue,
                                       int& generatedModels);
  /* Just returns the model index of best model
    */
  int getBestFit();
  void getNAParam(int& itmax, int& nsi, int& ns, int& nr);
  
  int currentID() {return _currentID;}
  int dimension() {return _nd;}
  int countModels();
  void maxOmegasCount(int& dispCount, int& ellCount);
  int currentLoadingIndex() {return _currentIndex;}
  
  QString reportName();
  QString reportFileName();
  
protected:
  // Boolean set to true is correclty opened report
  bool _valid;
  // Report Stream
  QDataStream _s;
  // Number of models really inside (all modes)
  int _trueN;
  // Index of current model just loaded (read mode)
  int _currentIndex;
  // Number of blocks, number of models=_n-HEADER_BLOCKS_COUNT
  int _n;
  // Number of parameters
  int _nd;
  void setMaxOmegasCount(int dispCount, int ellCount);
  int _maxDispOmegaCount;
  int _maxEllOmegaCount;
  // Parameters table
  struct paramAddress
  {
    paramID which;
    int layer;
  };
  paramAddress * _paramAddresses;
  // Last Layered Model loaded from stream
  Seismic1DModel * _currentModel;
  // Last cost value loaded from stream
  double _currentCost;
  // Last NA model loaded from stream (version 4)
  float * _currentNAModel;
  // Last ellipticity object loaded from stream
  CompatEllipticity * _currentEllipticity;
  // Temp offset in file (used by startWritingBlock(), endWritingBlock(), startReadingBlock())
  qint64 _ptrStart;
  // Offset table loaded on opening the file to speed up and avoid swap
  qint64 * _offsetTable;
  // Loads the list of parameters in memory
  bool loadParamAddresses();
  // Greatest unique ID, ready for next model (writing)
  // ID of currently load model (reading)
  // Usefull when models are sorted in file, or filtered to keep track of initial
  // number in file (when generated)
  int _currentID;
  // Name of the report
  QString _fileName;
  // initialize current model object
  bool initCurrentModel();
  // initialize current dispersion object
  bool initCurrentEllipticity();
  /* When reading this file while another process is also accessing it in writing mode
    it is sometimes usefull to flush the reading buffer but it is time consuming.
    When has been flushed it is useless to flush again and again while there's no more 
    model added. So this flag is set to true when flushed and returns to false when new models
    are added. Inside this only flush() set it to true. Other objetcts are free
    to set it as they want. */
  bool _flushed;
  // Number of blocks before the inversion results (blocks DispersionModel), set by constructor
  int _headerBlocksCount;
  // Current version
  int _version;
  //  Header block indexes (version dependent)
  qint64 _blockParam, _blockOmegas, _blockGoal, _blockEll;
  void writeNAModel(float * naModel);
  //-----------------------------------------------------------------------------------------------------------------------------
  //
  //  2nd part: specific to dispersion curve
  //
  //-----------------------------------------------------------------------------------------------------------------------------
public:
  // Checks if two report have the same omegas
  bool isSameOmegas(const CompatInversionReport * otherModels) const;
  // Checks if two report have the same goal dispersion
  bool isSameGoalDispersion(CompatInversionReport * otherModels);
  void addOmegas(CompatMultiModalFrequency * omegas);
  void addDispersionGoal(CompatDispersionData * rms);
  /*
     Re-writing over the existing goal curves, extrem caution may be observe to
     the size of the new CompatRMS: it must be exactly the same structure as the
     one added with addDispersionGoal(), only some changes in the values are tolerated
  */
  void modifyDispersionGoal(CompatDispersionData * rms);
  // Add all block in one time
  void addModel(const Seismic1DModel * model, float * naModel,
                const CompatDispersion * disp,
                double cost, int modelID=-1);
  // Add all block in one time
  void addModel(const Seismic1DModel * model, float * naModel,
                const CompatDispersion * disp,
                const CompatEllipticity * ell,
                double cost, int modelID=-1);
  // load the model, the cost and the dispersion curve
  bool loadDispersion(int modelIndex);
  // Same as above except that it loads the dispersion curve in a given dispersion object
  bool loadDispersion(CompatDispersion * disp, CompatEllipticity * ell,int modelIndex);
  // load the measured dispersion curve into a new data class
  CompatDispersionData * loadDispersionGoal();
  // Returns the dispersion of current model
  CompatDispersion * currentDispersion();
  // Returns the ellipticity of current model (may be null if there's no ellipticity)
  CompatEllipticity * currentEllipticity();
  int omegasCount();
  int modesCount();
  // Compute statistical parameter on calculated dispersion curves
  void statDispersion(int ipoint,int imode,int start_model,int end_model,
                      double& mean, double& stddev);
protected:
  // Last dispersion object loaded from stream
  CompatMultiModalCurves * _currentDispersion;
  // initialize current dispersion object
  virtual bool initCurrentDispersion();
};

} // namespace QGpCompatibility

#include "CompatDispersion.h"

namespace QGpCompatibility {

inline CompatDispersion * CompatInversionReport::currentDispersion()
{
  TRACE;
  if(!_currentDispersion) initCurrentDispersion();
  return (CompatDispersion *)_currentDispersion;
}

inline CompatEllipticity * CompatInversionReport::currentEllipticity()
{
  TRACE;
  if(!_currentDispersion) initCurrentDispersion();
  return _currentEllipticity;
}

inline int CompatInversionReport::omegasCount()
{
  TRACE;
  if(!_currentDispersion) initCurrentDispersion();
  return _currentDispersion->omegasCount();
}

inline int CompatInversionReport::modesCount()
{
  TRACE;
  if(!_currentDispersion) initCurrentDispersion();
  return _currentDispersion->modesCount();  
}

inline bool CompatInversionReport::isSameOmegas(const CompatInversionReport * otherModels) const
{
  TRACE;
  return _currentDispersion->isSameOmegas(otherModels->_currentDispersion);
}

inline bool CompatInversionReport::startReadingBlock(int index,BlockType wantedType)
{
  TRACE;
  if(index>=_trueN) {
    fprintf(stderr,"Model index out of range, max=%i, index=%i\n",_trueN-1,index);
    return false;
  }
  _s.device()->seek(_offsetTable[index]);
  int blockType=0;
  _s >> blockType;
  return (blockType==wantedType);
}

inline void CompatInversionReport::writeNAModel(float * naModel)
{
  TRACE;
  if(naModel)
    _s.writeRawData((char *)naModel,_nd*sizeof(float));
  else {
    double * tmp=new double [_nd];
    _s.writeRawData((char *)tmp,_nd*sizeof(float));
    delete [] tmp;
  }
}

} // namespace QGpCompatibility

#endif // COMPATINVERSIONREPORT_H