#ifndef NNLAYER_INC_CPP
#define NNLAYER_INC_CPP


#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>



class Layer {
public:
  Layer *prev, *next;
  int size;
  double trainRatio;
  double *a, *d, *da;

  Layer(Layer *prev, int size, double trainRatio = 0):
    prev(), next(), size(size), trainRatio(trainRatio)
  {
    assert(size > 0);
    a = new double[size*3];
    d = a + size;
    da = d + size;
    memset(a, 0, sizeof(*a)*size*3);
    if (prev) (this->prev = &prev->back())->next = this;
  }

  inline Layer& front()
    { Layer *l = this; while(l->prev) l = l->prev; return *l; }
  inline Layer& back()
    { Layer *l = this; while(l->next) l = l->next; return *l; }

  virtual ~Layer() {
    if (next) delete next;
    if (prev) prev->next = nullptr;
    delete[] a;
  }

  virtual bool toStream(FILE *f)
    { return next ? next->toStream(f) : true; }
  virtual bool fromStream(FILE *f)
    { return next ? next->fromStream(f) : true; }

  virtual Layer& pass()
    { return next ? next->pass() : *this; }
  virtual Layer& backpass()
    { return prev ? prev->backpass() : *this; }

  bool save(const char *filename) {
    assert(!prev);
    FILE *f = fopen(filename, "wb");
    if (!f) return printf("cannot open file '%s' for write\n", filename), false;
    if (!toStream(f)) return printf("cannot write to file '%s'\n", filename), fclose(f), false;
    fclose(f);
    return 1;
  }

  bool load(const char *filename) {
    assert(!prev);
    FILE *f = fopen(filename, "rb");
    if (!f) return printf("cannot open file '%s' for read\n", filename), false;
    if (!fromStream(f)) return printf("cannot read from file '%s'\n", filename), fclose(f), false;
    fclose(f);
    return 1;
  }

  double trainPass(double *x, double *y, double qmin) {
    assert(!prev);
    double *fa = a;
    a = x;
    Layer &b = pass();

    double qmax = 0;
    for(double *pa = b.a, *pda = b.da, *e = pa + b.size; pa < e; ++pa, ++pda, ++y) {
      double d = *y - *pa;
      *pda = d;
      double q = fabs(d);
      if (qmax < q) qmax = q;
    }
    if (qmax > qmin) b.backpass();

    a = fa;
    return qmax;
  }
};


class LayerPlain: public Layer {
public:
  double *b, *w;

  LayerPlain(Layer &prev, int size, double trainRatio = 0):
    Layer(&prev, size, trainRatio)
  {
    b = new double[size*(1 + prev.size)];
    w = b + size;
    for(double *p = b, *e = p + size*(1 + prev.size); p < e; ++p)
      *p = rand()/(double)RAND_MAX*2 - 1;
  }

  ~LayerPlain()
    { delete[] b; }

  bool toStream(FILE *f) override
    { return fwrite(b, sizeof(double)*size*(1 + prev->size), 1, f) && (!next || next->toStream(f)); }
  bool fromStream(FILE *f)
    { return fread (b, sizeof(double)*size*(1 + prev->size), 1, f) && (!next || next->toStream(f)); }

  Layer& pass() override {
    double *prevA = prev->a, *ee = prevA + prev->size;
    double *pw = w;
    for(double *pa = a, *pb = b, *pd = d, *e = pa + size; pa < e; ++pa, ++pb, ++pd) {
      double s = *pb;
      for(double *ppa = prevA; ppa < ee; ++ppa, ++pw)
        s += *ppa * *pw;
      double ex = exp(-s);
      double ex1 = ex + 1;
      double ex2 = 1/ex1;
      *pa = ex2;        // sigmoid
      *pd = ex*ex2*ex2; // sigmoid derivation
    }
    return next ? next->pass() : *this;
  }

  Layer& backpass() override {
    double tr = trainRatio;
    double *prevA = prev->a, *prevDa = prev->da, *ee = prevA + prev->size;
    double *pw = w;
    if (prev->prev) {
      memset(prevDa, 0, sizeof(*prev->da) * prev->size);
      for(double *pb = b, *pd = d, *pda = da, *e = pb + size; pb < e; ++pb, ++pd, ++pda) {
        double ds = *pd * *pda;
        double dst = ds*tr;
        *pb += dst;
        for(double *ppa = prevA, *ppda = prevDa; ppa < ee; ++ppa, ++ppda, ++pw) {
          *ppda += ds * *pw;
          *pw += dst * *ppa;
        }
      }
    } else {
      for(double *pb = b, *pd = d, *pda = da, *e = pb + size; pb < e; ++pb, ++pd, ++pda) {
        double dst = *pd * *pda * tr;
        *pb += dst;
        for(double *ppa = prevA; ppa < ee; ++ppa, ++pw)
          *pw += dst * *ppa;
      }
    }
    return prev->backpass();
  }
};


class LayerConvolution: public Layer {
public:
  enum { W = 4, D = W-2, WW = W*W };

  double *b, *w, *dw;
  int sx, sy, sz;
  int psx, psy, psz;

  LayerConvolution(Layer &prev, int sx, int sy, int sz, double trainRatio = 0):
    Layer(&prev, sx*sy*sz, trainRatio),
    sx(sx), sy(sy), sz(sz),
    psx(sx*2 + D), psy(sy*2 + D), psz(prev.size/(psx*psy))
  {
    assert(sx > 0 && sy > 0 && sz > 0);
    assert(psx > 0 && psy > 0 && psz > 0);
    assert(psx*psy*psz == prev.size);
    b = new double[size + WW*psz*(sz+1)];
    w = b + size;
    dw = w + WW*psz*sz;
    for(double *p = b, *e = p + (size + WW*psz*sz); p < e; ++p)
      *p = rand()/(double)RAND_MAX*2 - 1;
  }

  ~LayerConvolution()
    { delete[] b; }

  bool toStream(FILE *f) override
    { return fwrite(b, sizeof(double)*(size + WW*psz*sz), 1, f) && (!next || next->toStream(f)); }
  bool fromStream(FILE *f)
    { return fread (b, sizeof(double)*(size + WW*psz*sz), 1, f) && (!next || next->toStream(f)); }

  Layer& pass() override {
    double *pa = prev->a;
    int sxy = sx*sy, stepW = psz*WW, stepPA = psx*2, stepIPA = psx - W;
    for(double *ia = a, *ib = b, *id = d, *iw = w, *ea = ia + size; ia < ea; iw += stepW) {
      double *eew = iw + stepW;
      for(double *ipa = pa, *ea = ia + sxy; ia < ea; ipa += stepPA) {
        for(double *ea = ia + sx; ia < ea; ++ia, ++ib, ++id, ipa += 2) {
          double s = *ib;
          for(double *iiw = iw, *iipa = ipa; iiw < eew; iipa = ipa)
            for(int yy = 0; yy < W; ++yy, iipa += stepIPA)
              for(int xx = 0; xx < W; ++xx, ++iipa, ++iiw)
                s += *iipa * *iiw;
          double ex = exp(-s);
          double ex1 = ex + 1;
          double ex2 = 1/ex1;
          *ia = ex2;        // sigmoid
          *id = ex*ex2*ex2; // sigmoid derivation
        }
      }
    }
    return next ? next->pass() : *this;
  }

  Layer& backpass() override {
    double tr = trainRatio;
    int sxy = sx*sy, stepW = psz*WW, stepPA = psx*2, stepIPA = psx - W;
    double *dw = this->dw, *edw = dw + stepW, *pa = prev->a, *pda = prev->da;
    if (prev->prev) {
      memset(pda, 0, sizeof(*pda) * prev->size);
      for(double *ib = b, *id = d, *ida = da, *iw = w, *eda = ida + size; ida < eda; iw += stepW) {
        double *eew = iw + stepW;
        memset(dw, 0, sizeof(*dw) * stepW);
        for(double *ipa = pa, *ipda = pda, *eda = ida + sxy; ida < eda; ipa += stepPA, ipda += stepPA) {
          for(double *eda = ida + sx; ida < eda; ++ib, ++id, ++ida, ipa += 2) {
            double ds = *id * *ida;
            double dst = ds*tr;
            *ib += dst;
            for(double *iiw = iw, *idw = dw, *iipa = ipa, *iipda = ipda; idw < edw; iipa = ipa, iipda = ipda)
              for(int yy = 0; yy < W; ++yy, iipa += stepIPA, iipda += stepIPA)
                for(int xx = 0; xx < W; ++xx, ++iipa, ++iipda, ++iiw, ++idw)
                  { *iipda += ds * *iiw; *idw += dst * *iipa; }
          }
        }
        for(double *idw = dw; iw < eew; ++iw, ++idw)
          *iw += *idw;
      }
    } else {
      for(double *ib = b, *id = d, *ida = da, *iw = w, *eda = ida + size; ida < eda; iw += stepW) {
        memset(dw, 0, sizeof(*dw) * stepW);
        for(double *ipa = pa, *eda = ida + sxy; ida < eda; ipa += stepPA) {
          for(double *eda = ida + sx; ida < eda; ++ib, ++id, ++ida, ipa += 2) {
            double dst = *id * *ida * tr;
            *ib += dst;
            for(double *idw = dw, *iipa = ipa; idw < edw; iipa = ipa)
              for(int yy = 0; yy < W; ++yy, iipa += stepIPA)
                for(int xx = 0; xx < W; ++xx, ++iipa, ++idw)
                  *idw += dst * *iipa;
          }
        }
        for(double *idw = dw; idw < edw; ++iw, ++idw)
          *iw += *idw;
      }
    }
    return prev->backpass();
  }
};


#endif