libfbm.hpp

#ifndef __libfbm_hpp__
#define __libfbm_hpp__

#include <vector>
#include <string>

#include <stdint.h>

#define LIBFBM_MAX_DIM 8

#define LIBFBM_VERSION 0.3

namespace libfbm
{
  class GaussianGenerator
  {
  public:
    virtual ~GaussianGenerator();

    virtual void setSeed(uint32_t seed) = 0;
    virtual void setSeed(const unsigned char *seed_array, size_t len) = 0;

    virtual double getDouble() = 0;
    virtual void getDouble(double *array, size_t len) = 0;
  };

  #define _LIBFBM_GG_POOL_SIZE (312*4)

  class SFMTGaussianGenerator : public GaussianGenerator
  {
  public:
    SFMTGaussianGenerator(uint32_t seed = 0);
    SFMTGaussianGenerator(const unsigned char *seed_array, size_t len);
    ~SFMTGaussianGenerator();

    void setSeed(uint32_t seed);
    void setSeed(const unsigned char *seed_array, size_t len);

    double getDouble();
    void getDouble(double *array, size_t len);

  private:
    void *sfmt;
    uint64_t ipool[_LIBFBM_GG_POOL_SIZE];
    double opool[_LIBFBM_GG_POOL_SIZE];
    size_t ipindex;
    size_t opindex;

    void fill(double *array, size_t len);
  };

  class Seeder
  {
  public:
    Seeder(const char *path);
    ~Seeder();

    void seed(uint32_t index);

    inline size_t size() const { return sizeof(buf); }
    inline operator const unsigned char *() const { return (const unsigned char *)buf; }

  private:
    unsigned char buf[256];
    FILE *fp;

    Seeder(const Seeder&);
    const Seeder& operator=(const Seeder&);
  };

  class zvec
  {
  public:
    inline zvec(size_t dim) : clen(dim) { }

    inline int& operator[](size_t i) { return c[i]; }
    inline const int& operator[](size_t i) const { return c[i]; }

    inline size_t size() const { return clen; }

    inline bool increment(int base)
    {
      for ( size_t i = 0; i < size(); i++ )
        {
          if ( c[i] < base - 1 )
            {
              for ( size_t j = 0; j < i; j++ )
                  c[j] = 0;
              c[i]++;
              return true;
            }
        }
      return false;
    }

    inline bool increment(const zvec& dim)
    {
      for ( size_t i = 0; i < size(); i++ )
        {
          if ( c[i] < dim[i] - 1 )
            {
              for ( size_t j = 0; j < i; j++ )
                  c[j] = 0;
              c[i]++;
              return true;
            }
        }
      return false;
    }

    inline bool increment_but(const zvec& dim, size_t hold)
    {
      for ( size_t i = 0; i < size(); i++ )
        {
          if ( i == hold )
              continue;
          if ( c[i] < dim[i] - 1 )
            {
              for ( size_t j = 0; j < i; j++ )
                  c[j] = 0;
              c[i]++;
              return true;
            }
        }
      return false;
    }

    inline size_t index(size_t base) const
    {
      size_t ret = 0;
      size_t mul = 1;
      for ( size_t i = 0; i < size(); i++ )
        {
          ret += c[i] * mul;
          mul *= base;
        }
      return ret;
    }

    inline size_t index(const zvec& dim) const
    {
      size_t ret = 0;
      size_t mul = 1;
      for ( size_t i = 0; i < size(); i++ )
        {
          ret += c[i] * mul;
          mul *= dim[i];
        }
      return ret;
    }

    inline void zero()
    {
      for ( size_t i = 0; i < size(); i++ )
          c[i] = 0;
    }

    friend inline zvec operator+(const zvec& l, const zvec& r)
    {
      zvec ret(l.size());
      for ( size_t i = 0; i < l.size(); i++ )
          ret[i] = l[i] + r[i];
      return ret;
    }

    friend inline zvec operator-(const zvec& l, const zvec& r)
    {
      zvec ret(l.size());
      for ( size_t i = 0; i < l.size(); i++ )
          ret[i] = l[i] - r[i];
      return ret;
    }

    friend inline zvec operator*(int f, const zvec& v)
    {
      zvec ret(v.size());
      for ( size_t i = 0; i < v.size(); i++ )
          ret[i] = f * v[i];
      return ret;
    }

    friend inline zvec operator*(const zvec& v, int f)
    {
      zvec ret(v.size());
      for ( size_t i = 0; i < v.size(); i++ )
          ret[i] = f * v[i];
      return ret;
    }

  private:
    int c[LIBFBM_MAX_DIM];
    size_t clen;
  };

  class _Cache
  {
  public:
    _Cache(const std::string& cacheDir, const std::string& cacheName, size_t cacheSize);
    ~_Cache();

    bool load();
    void store(const std::vector<double>& data);

    size_t size() const { return cacheSize; }

    bool inMemory() const { return !data.empty(); }

  private:
    std::vector<double> data;
    FILE *fp;
    std::string cacheDir;
    std::string cacheName;
    size_t cacheSize;

    _Cache(const _Cache& copy);
    const _Cache& operator=(const _Cache& copy);

    friend class _CacheReader;
  };

  class _CacheReader
  {
  public:
    _CacheReader(_Cache& cache, size_t bufferSize);

    inline double operator[](size_t i)
    {
      if ( !buf.empty() && (i < bufpos || i - bufpos >= buf.size()) )
          load(i);
      return data[i - bufpos];
    }

  private:
    _Cache& cache;
    std::vector<double> buf;
    size_t bufpos;
    double *data;

    void load(size_t i);
  };

  class FFT
  {
  public:
    enum edir_t {
      forward = 0,
      backward = 1
    };

    FFT();
    ~FFT();

    void init(const libfbm::zvec& dim, double *data, edir_t direction);

    void execute();

  private:
    void *plan;
  };
  class Field;

  class SGPContext
  {
  public:
    SGPContext(const zvec& fieldDim, const zvec& userDim, const std::string& cacheName);
    virtual ~SGPContext();

    virtual double cov(const zvec& p) = 0; 
    inline const zvec& getDim() const { return userDim; }
    inline const zvec& getFieldDim() const { return fieldDim; }

    void setCacheDir(const std::string& cacheDir);

    inline size_t badEigenCount() const { return bec; }

    bool initCache(bool forceRecalc = false);

  protected:
    virtual void postProcess(Field& field, GaussianGenerator& rng);

    inline void setScaleResult(double f) { scaleResult = f; }


  private:
    zvec fieldDim, userDim;
    _Cache *cache;
    std::string cacheDir;
    std::string cacheName;
    size_t bec;
    double scaleResult;

    friend class Field;

    SGPContext(const SGPContext&);
    const SGPContext& operator=(const SGPContext&);
  };

  class FGNContext : public SGPContext
  {
  public:
    FGNContext(double H, const zvec& dim);

    double cov(const zvec& zvec);

  private:
    double H;
  };

  class FWSContext : public FGNContext
  {
  public:
    FWSContext(double H, const zvec& dim) : FGNContext(H, dim) { }

  protected:
    void postProcess(Field& field, GaussianGenerator& rng);
  };

  class AbstractField
  {
  public:
    inline AbstractField() { }
    virtual ~AbstractField();

    virtual void generate() = 0;
    virtual void generate(GaussianGenerator& rng) = 0;
    virtual void clear() = 0;
    virtual const zvec& getDim() const  = 0;

    inline double operator()(const zvec& p) const
    {
      size_t index = (size_t)2 * p[0];
      for ( size_t i = 1; i < p.size(); i++ )
          index += muls[i] * p[i];
      return datap[index];
    }
    inline double operator()(size_t x) const
    { return datap[2 * x]; }
    inline double operator()(size_t x, size_t y) const
    { return datap[2 * x + muls[1] * y]; }
    inline double operator()(size_t x, size_t y, size_t z) const
    { return datap[2 * x + muls[1] * y + muls[2] * z]; }
    inline double operator()(size_t x, size_t y, size_t z, size_t u) const
    { return datap[2 * x + muls[1] * y + muls[2] * z + muls[3] * u]; }
    inline double operator()(size_t x, size_t y, size_t z, size_t u, size_t v) const
    { return datap[2 * x + muls[1] * y + muls[2] * z + muls[3] * u + muls[4] * v]; }

    inline double& at(const zvec& p)
    {
      size_t index = 2 * p[0];
      for ( size_t i = 1; i < p.size(); i++ )
          index += muls[i] * p[i];
      return datap[index];
    }

    inline const double& at(const zvec& p) const
    {
      size_t index = 2 * p[0];
      for ( size_t i = 1; i < p.size(); i++ )
          index += muls[i] * p[i];
      return datap[index];
    }

    inline const size_t *getStrides() const { return muls; }

  protected:
    std::vector<double> Z;
    double *datap;
    size_t muls[LIBFBM_MAX_DIM];

  private:
    AbstractField(const AbstractField&);
    const AbstractField& operator=(const AbstractField&);
  };

  class Field : public AbstractField
  {
  public:
    Field(SGPContext& context, bool allowCorrelated = false);
    ~Field();

    void setBufferSize(size_t bufferSize);

    void generate();

    void generate(GaussianGenerator& rng);

    void clear();

    inline double operator()(const zvec& p) const
    {
      size_t index = (size_t)2 * p[0];
      for ( size_t i = 1; i < p.size(); i++ )
          index += muls[i] * p[i];
      return datap[index];
    }
    inline double operator()(size_t x) const
    { return datap[2 * x]; }
    inline double operator()(size_t x, size_t y) const
    { return datap[2 * x + muls[1] * y]; }
    inline double operator()(size_t x, size_t y, size_t z) const
    { return datap[2 * x + muls[1] * y + muls[2] * z]; }
    inline double operator()(size_t x, size_t y, size_t z, size_t u) const
    { return datap[2 * x + muls[1] * y + muls[2] * z + muls[3] * u]; }
    inline double operator()(size_t x, size_t y, size_t z, size_t u, size_t v) const
    { return datap[2 * x + muls[1] * y + muls[2] * z + muls[3] * u + muls[4] * v]; }

    void integrate();

    inline const zvec& getDim() const { return context.getDim(); }

    inline double& at(const zvec& p)
    {
      size_t index = 2 * p[0];
      for ( size_t i = 1; i < p.size(); i++ )
          index += muls[i] * p[i];
      return datap[index];
    }

    inline const double& at(const zvec& p) const
    {
      size_t index = 2 * p[0];
      for ( size_t i = 1; i < p.size(); i++ )
          index += muls[i] * p[i];
      return datap[index];
    }

    inline const size_t *getStrides() const { return muls; }

  private:
    FFT fft;
    SGPContext& context;
    bool allowCorrelated;
    size_t bufferSize;
    zvec psel;

    Field(const Field&);
    const Field& operator=(const Field&);
  };

  class SGPTester
  {
  public:
    void test(SGPContext& context, size_t printInterval = 1);
  };

  class FBMSteinContext : public SGPContext
  {
  public:
    FBMSteinContext(double H, size_t dim, size_t size, double Rhint = -1, bool mapDim = true);

    double cov(const zvec& zvec);

  protected:
    void postProcess(Field& field, GaussianGenerator& rng);

  public:
    void disablePostProcessing();

    double getR() const { return R; }

    static size_t userSize2FieldSize(size_t size, double R);
    static size_t fieldSize2UserSize(size_t size, double R);
    static double getRForH(double H, size_t dim, double Rhint = -1);

  private:
    double H;
    double R;
    double beta, c0, c2;
    double Rscale;
    bool pp;
  };

  class PowerLawContext : public SGPContext
  {
  public:
    PowerLawContext(double H, const zvec& dim, double Var = 2.0);

    double cov(const zvec& dim);

  private:
    double H, Var;
  };

  class PLFPSField : public AbstractField
  {
  public:
    PLFPSField(const zvec& dim, double H);
    ~PLFPSField();

    void generate();
    void generate(GaussianGenerator& rng);

    void clear();

    const zvec& getDim() const { return dim; }

  private:
    zvec dim;
    FFT fft;
    double H;
    std::vector<double> cache;
    size_t len;

    bool initCache();
  };

  extern SFMTGaussianGenerator gaussianRandomGenerator;
};

#endif // !__libfbm_hpp__