#ifndef NNLAYER2_INC_CPP
#define NNLAYER2_INC_CPP


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

#include <vector>
#include <algorithm>


typedef float Real;


template<typename T>
bool twrite(const T &x, FILE *f)
  { return fwrite(&x, sizeof(x), 1, f); }
template<typename T>
bool tread(const T &x, FILE *f)
  { return fwrite(&x, sizeof(x), 1, f); }


struct Link;
struct Neuron {
  Real v, d;
  inline Neuron(): v(), d() { }
};

struct Link {
  int nprev, lnext;
  Real w;
  inline Link(): nprev(), lnext(), w((Real)rand()/(Real)RAND_MAX*2 - 1) { }
};


class Layer {
public:
  Layer *prev, *next;

  int size, lsize;
  Neuron *neurons;
  Link *links;
  int lfirst;

  explicit Layer(Layer *prev = nullptr, int size = 0, int lsize = 0):
    prev(), next(), size(), lsize(), neurons(), links(), lfirst()
  {
    while(prev && prev->next) prev = prev->next;
    if (prev) prev->next = this;
    this->prev = prev;
    init(size, lsize);
  }

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

  bool init(int size, int lsize = 0) {
    clear();
    if (size <= 0) return false;
    if (lsize < 0) lsize = 0;

    if (prev ? lsize <= 0 : lsize) return false;
    if (size) neurons = new Neuron[size];
    if (lsize) links = new Link[lsize*size];
    this->size = size;
    this->lsize = lsize;
    return true;
  }

  void clear() {
    if (neurons) delete[] neurons;
    if (links) delete[] links;
    this->size = lsize = 0;
    neurons = nullptr;
    links = nullptr;
    lfirst = 0;
  }

  void prepareBackLinks() {
    if (!links || lsize <= 0) return;

    lfirst = 0;
    int lend = size*lsize;
    for(Link *il = links, *e = il + lend; il < e; ++il) {
      assert(prev && prev->neurons && il->nprev >= 0 && il->nprev < prev->size);
      il->lnext = il - links + 1;
    }

    while(true) {
      bool done = true;
      for(int *il = &lfirst; links[*il].lnext != lend; il = &links[*il].lnext) {
        int a = *il, b = links[a].lnext;
        if (links[a].nprev > links[b].nprev) {
          links[a].lnext = links[b].lnext;
          links[b].lnext = a;
          *il = b;
          done = false;
        }
      }
      if (done) break;
    }

    #ifndef NDEBUG
    if (lsize) {
      int next = 0;
      for(int il = lfirst; il != lend; il = links[il].lnext) {
        assert(links[il].nprev == next || links[il].nprev == next-1);
        if (links[il].nprev == next) ++next;
      }
      assert(next = prev->size);
    }
    #endif
  }

  bool save(FILE *f) const {
    if (!twrite(size, f) || !twrite(lsize, f))
      return false;
    for(Link *il = links, *e = il + lsize*size; il < e; ++il)
      if (!twrite(il->nprev, f) || !twrite(il->w, f))
          return false;
    return true;
  }

  bool load(FILE *f) {
    clear();
    int size = 0, lsize = 0;
    if (!tread(size, f) || !tread(lsize, f) || !init(size, lsize))
      return false;
    for(Link *il = links, *e = il + lsize*size; il < e; ++il)
      if (!tread(il->nprev, f) || !tread(il->w, f) || !prev || !prev->neurons || il->nprev < 0 || il->nprev >= prev->size)
          return false;
    prepareBackLinks();
    return true;
  }

  bool saveAll(FILE *f) const
    { return save(f) && (!next || next->save(f)); }
  bool loadAll(FILE *f)
    { return load(f) && (!next || next->load(f)); }

  bool saveAll(const char *filename) const {
    assert(!prev);
    FILE *f = fopen(filename, "wb");
    if (!f)
      { printf("cannot open file for write: %s\n", filename); return false; }
    int count = totalLayers();
    if (!twrite(count, f) || !saveAll(f))
      { printf("cannot save to file: %s\n", filename); fclose(f); return false; }
    fclose(f);
    return true;
  }

  bool loadAll(const char *filename) {
    assert(!prev);
    FILE *f = fopen(filename, "rb");
    if (!f)
      { printf("cannot open file for read: %s\n", filename); return false; }
    int count = 0;
    if (!tread(count, f) || count <= 0)
      { printf("cannot load from file: %s\n", filename); fclose(f); return false; }
    if (next)
      delete next;
    while(--count)
      new Layer(this);
    if (!loadAll(f))
      { printf("cannot load from file: %s\n", filename); fclose(f); return false; }
    fclose(f);
    return true;
  }

  int totalLayers() const
    { int t = 0; for(const Layer *l = this; l; l = l->next) ++t; return t; }
  int totalNeurons() const
    { int t = 0; for(const Layer *l = this; l; l = l->next) t += l->size; return t; }
  int totalLinks() const
    { int t = 0; for(const Layer *l = this; l; l = l->next) t += l->size*l->lsize; return t; }
  inline size_t totalMemSize() const
    { return totalNeurons()*sizeof(Neuron) + totalLinks()*sizeof(Link); }

  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; }

  void pass(int nb, int ne) {
    assert(prev);
    int lsize = this->lsize;
    Neuron *pn = prev->neurons;
    Link *il = links + nb*lsize;
    for(Neuron *in = neurons + nb, *e = neurons + ne; in < e; ++in) {
      double s = 0;
      for(Link *e = il + lsize; il < e; ++il)
        s += pn[il->nprev].v * il->w;

      // exp sigmoid
      double ss = 1/(1 + exp(-s));
      in->v = ss;
      in->d = ss * (1-ss);

      // 1/(x+1) sigmoid
      //double ss = 1/(1+fabs(s));
      //double ss2 = ss*0.5;
      //in->v = s > 0 ? 1 - ss2 : ss2;
      //in->d = ss2 * ss;
    }
  }

  template<bool hasPrev>
  void backpassTpl(int lb, int le) {
    assert(hasPrev == (bool)prev);
    if (lb == le) return;
    Link *links = this->links;
    Neuron *neurons = this->neurons;
    Neuron *pneurons = prev->neurons;

    double s = 0;
    int ipn = links[lb].nprev;
    for(int il = lb; il != le; ) {
      Link &l = links[il];
      if (hasPrev)
        if (ipn != l.nprev)
          { pneurons[ipn].d *= s; ipn = l.nprev; }
      Neuron &n = neurons[ il/lsize ];
      Real d = n.d;
      if (hasPrev) s += d * l.w;
      l.w += d * n.v;
      il = l.lnext;
    }
    if (hasPrev)
        pneurons[ipn].d *= s;

    assert(le == size*lsize || ipn < links[le].nprev);
  }

  void backpass(int lb, int le) {
    if (prev) backpassTpl<true>(lb, le); else backpassTpl<false>(lb, le);
  }

  void passAll() {
    if (prev) pass(0, size);
    if (next) next->passAll();
  }

  void backpassAll(Real k) {
    if (next) backpass(lfirst, size*lsize);
    if (prev) prev->passAll();
  }
};


#endif