diff --git a/simple/neural/nn-heli.cpp b/simple/neural/nn-heli.cpp new file mode 100644 index 0000000..e79c897 --- /dev/null +++ b/simple/neural/nn-heli.cpp @@ -0,0 +1,167 @@ + +#include +#include +#include + +#include + +#include "nnlayer.inc.cpp" + + +Layer *nl; +Framebuffer fb, fbMin; +Animation fbAnim, fbMinAnim; + +int wasPressed; +double prevX, prevY; + + + +void prepareImage() { + int w, h; + unsigned char *pixels = NULL; + + saveState(); + target(fb); + imageFromViewport(&w, &h, &pixels); + restoreState(); + if (!pixels) return; + + int x0 = w, y0 = h, x1 = 0, y1 = 0; + for(int y = 0; y < h; ++y) { + for(int x = 0; x < w; ++x) { + if (imageGetPixel(w, h, pixels, x, y) != 0x000000ff) { + if (x0 > x) x0 = x; + if (x1 < x) x1 = x; + if (y0 > y) y0 = y; + if (y1 < y) y1 = y; + } + } + } + free(pixels); + pixels = NULL; + + if (x1 < x0 || y1 < y0) return; + + int fw = framebufferGetWidth(fbMin); + int fh = framebufferGetHeight(fbMin); + + double wx = x1 - x0 + 1; + double wy = y1 - y0 + 1; + double s = (fw - 4)/(double)(wx > wy ? wx : wy); + double cx = (x0 + x1)/2.0; + double cy = (y0 + y1)/2.0; + + double xx = fw/2 - s*cx; + double yy = fh/2 - s*cy; + double ww = s*w; + double hh = s*h; + + saveState(); + target(fbMin); + noStroke(); + rectTextured(fbAnim, xx, yy, ww, hh); + imageFromViewport(&w, &h, &pixels); + restoreState(); + + if (!pixels) return; + double *p = nl->a; + for(int y = 0; y < h; ++y) + for(int x = 0; x < w; ++x) + *p++ = colorGetValue(imageGetPixel(w, h, pixels, x, y)); +} + + +void init() { + background(COLOR_BLACK); + stroke(COLOR_WHITE); + fb = createFramebufferEx(512, 512, NULL, FALSE, FALSE, TRUE); + fbMin = createFramebufferEx(28, 28, NULL, FALSE, FALSE, TRUE); + fbAnim = createAnimationFromFramebuffer(fb); + fbMinAnim = createAnimationFromFramebuffer(fbMin); + + saveState(); + target(fb); + clear(); + target(fbMin); + clear(); + restoreState(); + + nl = new Layer(nullptr, 784); + new LayerAdd(*nl, 784); + new LayerSimple(*nl, 256); + new LayerSimple(*nl, 128); + new LayerSimple(*nl, 64); + new LayerSimple(*nl, 64); + new LayerSimple(*nl, 10); + + //nl = new Layer(nullptr, 22*22); + //new LayerConvolution(*nl, 10, 10, 4); + //new LayerConvolution(*nl, 4, 4, 16); + //new LayerConvolution(*nl, 1, 1, 64); + //new LayerPlain(*nl, 64); + //new LayerPlain(*nl, 32); + //new LayerPlain(*nl, 10); + + nl->load("data/weights.bin"); +} + + +void draw() { + saveState(); + + if (mouseDown("left")) { + double x = mouseX(), y = mouseY(); + if (!wasPressed) prevX = x, prevY = y; + + saveState(); + strokeWidth(32); + target(fb); + line(prevX, prevY, x, y); + restoreState(); + + prevX = x, prevY = y; + wasPressed = TRUE; + } else { + wasPressed = FALSE; + } + + if (keyWentDown("space")) { + prepareImage(); + nl->pass(); + saveState(); + target(fb); + clear(); + restoreState(); + } + + noStroke(); + rectTextured(fbAnim, 0, 0, 512, 512); + + stroke(COLOR_WHITE); + rectTextured(fbMinAnim, 16, 16, 22, 22); + + noFill(); + + Layer &nlb = nl->back(); + textSize(8); + int res = 0; + for(int i = 0; i < 10; ++i) { + if (nlb.a[i] > nlb.a[res]) res = i; + textf(16, 90+8*i, "%d: %lf", i, nlb.a[i]); + } + textSize(16); + textf(16, 60, "%d", res); + + restoreState(); +} + + +int main(int largc, char **largv) { + windowSetVariableFrameRate(); + windowSetInit(&init); + windowSetDraw(&draw); + windowRun(); + return 0; +} + diff --git a/simple/neural/nn-trainer.cpp b/simple/neural/nn-trainer.cpp index 6ecc350..2465ed2 100644 --- a/simple/neural/nn-trainer.cpp +++ b/simple/neural/nn-trainer.cpp @@ -2,35 +2,62 @@ #include #include +#include "nnlayer.inc.cpp" +#include "nnlayer.conv.inc.cpp" +#include "nnlayer.lnk.inc.cpp" #include "nntrain.inc.cpp" int main() { srand(time(NULL)); + const char *filename = "data/output/weights.bin"; // 28x28 + //const char *filename = "data/output/weights22.bin"; + printf("load training data\n"); - Trainer t; + Trainer t(0.25); if (!t.loadSymbolMap("data/symbols-data.bin", 784, 10)) return 1; + //if (!t.loadSymbolMap("data/symbols22-data.bin", 22*22, 10)) return 1; printf("create neural network\n"); - double tr = 0.1; - Layer l(nullptr, 784, tr); - new LayerPlain(l, 256, tr); - new LayerPlain(l, 128, tr); - new LayerPlain(l, 64, tr); - new LayerPlain(l, 64, tr); - new LayerPlain(l, 10, tr); + //Layer l(nullptr, 784); + //new LayerSimple(l, 128); + //new LayerSimple(l, 64); + //new LayerSimple(l, 10); + + //Layer l(nullptr, 784); + //new LayerConvolution(l, 12, 12, 4); + //new LayerConvolution(l, 4, 4, 16); + ////new LayerConvolution(l, 1, 1, 128); + ////new LayerSimple(l, 64); + //new LayerSimple(l, 10); + + //Layer l(nullptr, 784); + //new LayerConvolution(l, 22, 22, 1); + //new LayerConvolution(l, 14, 14, 1); + //new LayerConvolution(l, 4, 4, 1); + //new LayerSimple(l, 10); + + Layer l(nullptr, 784); + new LayerLinkConvolution(l, 28, 28, 1, 22, 22, 1, 60, 1, 3); + new LayerLinkConvolution(l, 22, 22, 1, 14, 14, 1, 100, 1, 4); + new LayerLinkConvolution(l, 14, 14, 1, 4, 4, 1, 140, 1, 5); + new LayerSimple(l, 10); + + printf(" neurons: %d, links %d\n", l.totalSize(), l.totalLinks()); - printf("try load previously saved network\n"); - l.load("data/weights.bin"); + //printf("try load previously saved network\n"); + //l.load(filename); printf("train\n"); - double k = pow(0.5, 0.125); - for(double q = k; q > 0.05; q *= k) - t.train(l, 10, 100, q); + t.trainSimple(l, 1000, 0.5, 10000); + //for(double q = 0.5; q > 0.05; q *= 0.5) + // t.trainBlock(l, 50, 20, q, 10000); + //t.trainBlock(l, 50, 20, 0.5, 10000); + //t.trainBlock(l, 10, 100, 0.5, 1); printf("save neural network weights\n"); - if (!l.save("data/output/weights.bin")) return 1; + if (!l.save(filename)) return 1; return 0; } diff --git a/simple/neural/nnlayer.conv.inc.cpp b/simple/neural/nnlayer.conv.inc.cpp new file mode 100644 index 0000000..ff792a5 --- /dev/null +++ b/simple/neural/nnlayer.conv.inc.cpp @@ -0,0 +1,221 @@ +#ifndef NNLAYER_CONV_INC_CPP +#define NNLAYER_CONV_INC_CPP + + +#include "nnlayer.inc.cpp" + + +template +class LayerConvolution: public Layer { +public: + enum { W = Size, WW = W*W, D = Step, P = Padding, P2 = P*2 }; + + int sx, sy, sz; + int psx, psy, psz; + + LayerConvolution(Layer &prev, int sx, int sy, int sz): + Layer(&prev, sx*sy*sz), + sx(sx), sy(sy), sz(sz), + psx((sx-P2-1)*D + W), psy((sy-P2-1)*D + W), psz(this->prev->size/(psx*psy)) + { + assert(sx > 0 && sy > 0 && sz > 0); + assert(psx > 0 && psy > 0 && psz > 0); + assert(psx*psy*psz == this->prev->size); + links = wsize = WW*psz*sz*(sy-P2)*(sx-P2); + w = new double[wsize]; + double k = RaLU ? 1.0/(WW*psz*sz) : 1; + for(double *iw = w, *e = iw + wsize; iw < e; ++iw) + *iw = (rand()/(double)RAND_MAX*2 - 1)*k; + } + + Layer& pass() override { + const int asy = sx*(sy-P2); + const int asx = sx-P2; + const int adz = sx*P2; + const int wdz = WW*psz; + const int pady = D*(psx-asx); + const int paddz = psx*(psy - W); + const int paddy = psx - W; + double *pa = prev->a; + + double *ia = a + P*sx + P; + double *iw = w; + double *ipa = pa; + for(double *e = ia + asy*sz; ia < e; ia += adz, ipa = pa) { + for(double *e = ia + asy; ia < e; ia += P2, ipa += pady) { + for(double *e = ia + asx; ia < e; ++ia, ipa += D) { + double s = 0; + double *iipa = ipa; + for(double *ew = iw + wdz; iw < ew; iipa += paddz) + for(int yy = 0; yy < W; ++yy, iipa += paddy) + for(int xx = 0; xx < W; ++xx, ++iw, ++iipa) + s += *iipa * *iw; + *ia = RaLU ? (s < -1 ? -1 : s > 1 ? 1 : s) : 1/(1 + exp(-s)); // RaLU : sigmoid + } + } + } + return next ? next->pass() : *this; + } + + template + Layer& backpassT(double trainRatio) { + const int asy = sx*(sy-P2); + const int asx = sx-P2; + const int adz = sx*P2; + const int wdz = WW*psz; + const int pady = D*(psx-asx); + const int paddz = psx*(psy - W); + const int paddy = psx - W; + double *pa = prev->a; + double *pda = prev->da; + + double *ia = a + P*sx + P; + double *iw = w; + double *ida = da + P*sx + P; + double *ipa = pa; + double *ipda = pda; + if (Deep) memset(pda, 0, sizeof(*pda)*prev->size); + for(double *e = ia + asy*sz; ia < e; ia += adz, ida += adz, ipa = pa, ipda = pda) { + for(double *e = ia + asy; ia < e; ia += P2, ida += P2, ipa += pady, ipda += pady) { + for(double *e = ia + asx; ia < e; ++ia, ++ida, ipa += D, ipda += D) { + double ds; + if (RaLU) { + double a = *ia; + if (a == -1 || a == 1) { iw += wdz; continue; } + ds = *ida; // RaLU derivation * *ida + } else { + double a = *ia; + ds = a * (1-a) * *ida; // sigmoid derivation * *ida + } + double dst = ds*trainRatio; + + double *iipa = ipa; + double *iipda = ipda; + for(double *ew = iw + wdz; iw < ew; iipa += paddz, iipda += paddz) + for(int yy = 0; yy < W; ++yy, iipa += paddy, iipda += paddy) + for(int xx = 0; xx < W; ++xx, ++iw, ++iipa, ++iipda) { + if (Deep) *iipda += ds * *iw; + *iw += dst * *iipa; + } + } + } + } + + return prev->backpass(trainRatio); + } + + Layer& backpass(double trainRatio) override + { return prev->prev ? backpassT(trainRatio) : backpassT(trainRatio); } +}; + + + +template +class LayerConvolutionShared: public Layer { +public: + enum { W = Size, WW = W*W, D = Step, P = Padding, P2 = P*2 }; + + double *dw; + int sx, sy, sz; + int psx, psy, psz; + + LayerConvolutionShared(Layer &prev, int sx, int sy, int sz): + Layer(&prev, sx*sy*sz), + sx(sx), sy(sy), sz(sz), + psx((sx-P2-1)*D + W), psy((sy-P2-1)*D + W), psz(this->prev->size/(psx*psy)) + { + assert(sx > 0 && sy > 0 && sz > 0); + assert(psx > 0 && psy > 0 && psz > 0); + assert(psx*psy*psz == this->prev->size); + wsize = WW*psz*sz; + links = wsize*(sy-P2)*(sx-P2); + w = new double[wsize + WW*psz]; + dw = w + wsize; + for(double *iw = w, *e = iw + wsize; iw < e; ++iw) + *iw = (rand()/(double)RAND_MAX*2 - 1)*1; + } + + Layer& pass() override { + const int asy = sx*(sy-P2); + const int asx = sx-P2; + const int adz = sx*P2; + const int wdz = WW*psz; + const int pady = D*(psx-asx); + const int paddz = psx*(psy - W); + const int paddy = psx - W; + double *pa = prev->a; + + double *ia = a + P*sx + P; + double *iw = w; + double *ipa = pa; + for(double *e = ia + asy*sz; ia < e; ia += adz, iw += wdz, ipa = pa) { + double *eew = iw + wdz; + for(double *e = ia + asy; ia < e; ia += P2, ipa += pady) { + for(double *e = ia + asx; ia < e; ++ia, ipa += D) { + double s = 0; + double *iipa = ipa; + for(double *iiw = iw; iiw < eew; iipa += paddz) + for(int yy = 0; yy < W; ++yy, iipa += paddy) + for(int xx = 0; xx < W; ++xx, ++iiw, ++iipa) + s += *iipa * *iiw; + *ia = 1/(1 + exp(-s)); // sigmoid + } + } + } + return next ? next->pass() : *this; + } + + template + Layer& backpassT(double trainRatio) { + const int asy = sx*(sy-P2); + const int asx = sx-P2; + const int adz = sx*P2; + const int wdz = WW*psz; + const int pady = D*(psx-asx); + const int paddz = psx*(psy - W); + const int paddy = psx - W; + double *dw = this->dw; + double *edw = dw + wdz; + double *pa = prev->a; + double *pda = prev->da; + + double *ia = a + P*sx + P; + double *iw = w; + double *ida = da + P*sx + P; + double *ipa = pa; + double *ipda = pda; + if (Deep) memset(pda, 0, sizeof(*pda)*prev->size); + for(double *e = ia + asy*sz; ia < e; ia += adz, ida += adz, ipa = pa, ipda = pda) { + memset(dw, 0, sizeof(*dw) * wdz); + + for(double *e = ia + asy; ia < e; ia += P2, ida += P2, ipa += pady, ipda += pady) { + for(double *e = ia + asx; ia < e; ++ia, ++ida, ipa += D, ipda += D) { + const double a = *ia; + const double ds = a * (1-a) * *ida; // sigmoid derivation * *ida + const double dst = ds*trainRatio; + + double *iiw = iw; + double *iipa = ipa; + double *iipda = ipda; + for(double *idw = dw; idw < edw; iipa += paddz, iipda += paddz) + for(int yy = 0; yy < W; ++yy, iipa += paddy, iipda += paddy) + for(int xx = 0; xx < W; ++xx, ++iiw, ++idw, ++iipa, ++iipda) { + if (Deep) *iipda += ds * *iiw; + *idw += dst * *iipa; + } + } + } + + for(double *idw = dw; idw < edw; ++iw, ++idw) + *iw += *idw; + } + + return prev->backpass(trainRatio); + } + + Layer& backpass(double trainRatio) override + { return prev->prev ? backpassT(trainRatio) : backpassT(trainRatio); } +}; + + +#endif diff --git a/simple/neural/nnlayer.inc.cpp b/simple/neural/nnlayer.inc.cpp index a697dc5..a3cef7e 100644 --- a/simple/neural/nnlayer.inc.cpp +++ b/simple/neural/nnlayer.inc.cpp @@ -13,41 +13,45 @@ class Layer { public: Layer *prev, *next; - int size; - double trainRatio; - double *a, *d, *da; + int size, wsize, links; + double *a, *da, *w; - Layer(Layer *prev, int size, double trainRatio = 0): - prev(), next(), size(size), trainRatio(trainRatio) + Layer(Layer *prev, int size): + prev(), next(), size(size), wsize(), links(), w() { assert(size > 0); - a = new double[size*3]; - d = a + size; - da = d + size; - memset(a, 0, sizeof(*a)*size*3); + a = new double[size*2]; + da = a + size; + memset(a, 0, sizeof(*a)*size*2); 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; + if (w) delete[] w; } - 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; } + virtual Layer& backpass(double trainRatio) + { return prev ? prev->backpass(trainRatio) : *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; } + + inline int totalSize() const + { int c = 0; for(const Layer *l = this; l; l = l->next) c += l->size; return c; } + inline int totalLinks() const + { int c = 0; for(const Layer *l = this; l; l = l->next) c += l->links; return c; } + + bool toStream(FILE *f) + { return (!w || fwrite(w, sizeof(double)*wsize, 1, f)) && (!next || next->toStream(f)); } + bool fromStream(FILE *f) + { return (!w || fread (w, sizeof(double)*wsize, 1, f)) && (!next || next->fromStream(f)); } bool save(const char *filename) { assert(!prev); @@ -67,187 +71,88 @@ public: return 1; } - double trainPass(double *x, double *y, double qmin) { + double trainPass(double ratio, double *x, double *y, double qmin) { assert(!prev); - double *fa = a; - a = x; + memcpy(a, x, sizeof(*a)*size); Layer &b = pass(); double qmax = 0; for(double *pa = b.a, *pda = b.da, *e = pa + b.size; pa < e; ++pa, ++pda, ++y) { + assert(*pa == *pa); double d = *y - *pa; *pda = d; + //qmax += d*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)); } + //qmax = sqrt(qmax/b.size); + if (qmax > qmin) + b.backpass(ratio); - 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; - } + //if (qmax < 1e-6) { + // printf("strange:\n"); + // y -= b.size; + // for(double *pa = b.a, *e = pa + b.size; pa < e; ++pa, ++y) + // printf("%f - %f = %f\n", *y, *pa, *y - *pa); + //} - 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(); + return qmax; } }; -class LayerConvolution: public Layer { +class LayerSimple: 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)) + LayerSimple(Layer &prev, int size): + Layer(&prev, size) { - 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; + links = wsize = size * this->prev->size; + w = new double[wsize]; + double k = 1.0/this->prev->size; + for(double *iw = w, *e = iw + wsize; iw < e; ++iw) + *iw = (rand()/(double)RAND_MAX*2 - 1)*k; } - ~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 - } - } + double *pa = prev->a, *ee = pa + prev->size; + double *iw = w; + for(double *ia = a, *e = ia + size; ia < e; ++ia) { + double s = 0; + for(double *ipa = pa; ipa < ee; ++ipa, ++iw) + s += *ipa * *iw; + *ia = 1/(1 + exp(-s)); // sigmoid + //*ia = s > 0 ? s : 0; // RaLU } 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; + Layer& backpass(double trainRatio) override { + double *pa = prev->a, *pda = prev->da, *ee = pa + prev->size; + double *iw = w; 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; } - } + memset(pda, 0, sizeof(*prev->da) * prev->size); + for(double *ia = a, *ida = da, *e = ia + size; ia < e; ++ia, ++ida) { + double a = *ia; + double ds = a * (1-a) * *ida; // sigmoid derivation * ida + //if (!*ia) continue; // RaLU derivation is zero + //double ds = *ida; + double dst = ds*trainRatio; + for(double *ipa = pa, *ipda = pda; ipa < ee; ++ipa, ++ipda, ++iw) { + *ipda += ds * *iw; + *iw += dst * *ipa; } - 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; + for(double *ia = a, *ida = da, *e = ia + size; ia < e; ++ia, ++ida) { + double a = *ia; + double dst = a * (1-a) * *ida * trainRatio; // sigmoid derivation * ida * trainRatio + //if (!*ia) continue; // RaLU derivation is zero + //double dst = *ida * trainRatio; + for(double *ipa = pa; ipa < ee; ++ipa, ++iw) + *iw += dst * *ipa; } } - return prev->backpass(); + return prev->backpass(trainRatio); } }; diff --git a/simple/neural/nnlayer.lnk.inc.cpp b/simple/neural/nnlayer.lnk.inc.cpp new file mode 100644 index 0000000..3efc334 --- /dev/null +++ b/simple/neural/nnlayer.lnk.inc.cpp @@ -0,0 +1,142 @@ +#ifndef NNLAYER_LNK_INC_CPP +#define NNLAYER_LNK_INC_CPP + + +#include + +#include "nnlayer.inc.cpp" + + +class LayerLink: public Layer { +public: + int lsize; + double **wa, **wda; + + LayerLink(Layer &prev, int size, int lsize): + Layer(&prev, size), lsize(lsize) + { + assert(lsize > 0); + links = wsize = size*lsize; + w = new double[wsize]; + wa = new double*[wsize*2]; + wda = wa + wsize; + memset(wa, 0, sizeof(*wa)*wsize*2); + for(double *iw = w, *e = iw + wsize; iw < e; ++iw) + *iw = rand()/(double)RAND_MAX*2 - 1; + } + + ~LayerLink() + { delete[] wa; } + + Layer& pass() override { + double *ia = a; + double *iw = w; + double **iwa = wa; + for(double *e = ia + size; ia < e; ++ia) { + double s = 0; + for(double *e = iw + lsize; iw < e; ++iw, ++iwa) { + assert(*iwa); + s += *iw * **iwa; + } + *ia = 1/(1 + exp(-s)); // sigmoid + } + return next ? next->pass() : *this; + } + + template + Layer& backpassT(double trainRatio) { + double *ia = a; + double *ida = da; + double *iw = w; + double **iwa = wa; + double **iwda = wda; + if (Deep) memset(prev->da, 0, sizeof(*prev->da)*prev->size); + for(double *e = ia + size; ia < e; ++ia, ++ida) { + double a = *ia; + double ds = a * (1-a) * *ida; // sigmoid derivation * *ida + double dst = ds*trainRatio; + for(double *e = iw + lsize; iw < e; ++iw, ++iwa, ++iwda) { + assert(*iwa && *iwda); + if (Deep) **iwda += ds * *iw; + *iw += dst * **iwa; + } + } + return prev->backpass(trainRatio); + } + + Layer& backpass(double trainRatio) override + { return prev->prev ? backpassT(trainRatio) : backpassT(trainRatio); } +}; + + + +class LayerLinkConvolution: public LayerLink { +public: + LayerLinkConvolution(Layer &prev, int psx, int psy, int psz, int sx, int sy, int sz, int lsize, int step = 1, int pad = 0): + LayerLink(prev, sx*sy*sz, lsize*psz) + { + assert(psx > 0 && psy > 0 && psz > 0); + assert(sx > 0 && sy > 0 && sz > 0); + assert(step > 0 && pad >= 0); + assert(psx*psy*psz == this->prev->size); + assert(pad + (sx-1)*step < psx); + assert(pad + (sy-1)*step < psy); + + int hs = (int)ceil(sqrt(lsize)) + 1; + int s = hs*2 + 1; + struct Point { + int x, y, r; + inline bool operator< (const Point &p) const { return r < p.r; } + } *points = new Point[s*s], *p = points; + + for(int y = -hs; y <= hs; ++y) { + for(int x = -hs; x <= hs; ++x, ++p) { + int sector = x >= 0 && y > 0 ? 1 + : y >= 0 && x < 0 ? 2 + : x <= 0 && y < 0 ? 3 + : y <= 0 && x > 0 ? 4 + : 0; + int subsector = sector % 2 ? abs(x) >= abs(y) : abs(y) >= abs(x); + p->x = x; + p->y = y; + p->r = (x*x + y*y)*10 + sector*2 + subsector; + } + } + std::sort(points, points + s*s); + + int *order = new int[lsize]; + for(int z = 0; z < sz; ++z) { + for(int y = 0; y < sy; ++y) { + for(int x = 0; x < sx; ++x) { + p = points; + int *io = order; + for(int i = 0; i < lsize; ++i, ++io) { + int xx, yy; + do { + xx = pad + x*step + p->x; + yy = pad + y*step + p->y; + ++p; + } while(xx < 0 || yy < 0 || xx >= psx || yy >= psy); + *io = yy*psx + xx; + } + std::sort(order, order + lsize); + + double **iwa = &wa [ (z*sx*sy + y*sx + x)*this->lsize ]; + double **iwda = &wda[ (z*sx*sy + y*sx + x)*this->lsize ]; + for(int pz = 0; pz < psz; ++pz) { + for(int i = 0; i < lsize; ++i, ++iwa, ++iwda) { + *iwa = &this->prev->a [pz*psx*psy + order[i]]; + *iwda = &this->prev->da[pz*psx*psy + order[i]]; + } + } + } + } + } + + delete[] order; + delete[] points; + } +}; + + +#endif diff --git a/simple/neural/nntrain.inc.cpp b/simple/neural/nntrain.inc.cpp index 1dcc54b..5432828 100644 --- a/simple/neural/nntrain.inc.cpp +++ b/simple/neural/nntrain.inc.cpp @@ -7,20 +7,23 @@ class Trainer { public: + double trainRatio; int sizeX, sizeY, count; double *x, *y; - Trainer(): sizeX(), sizeY(), count(), x(), y() { } + explicit Trainer(double trainRatio = 0.5): trainRatio(trainRatio), sizeX(), sizeY(), count(), x(), y() { } - Trainer(int sizeX, int sizeY, int count): Trainer() + Trainer(double trainRatio, int sizeX, int sizeY, int count): + Trainer(trainRatio) { init(sizeX, sizeY, count); } ~Trainer() - { if (count) delete[] x; } + { deinit(); } void init(int sizeX, int sizeY, int count) { + deinit(); assert(sizeX > 0); assert(sizeY > 0); assert(count > 0); @@ -40,12 +43,67 @@ public: } - double train(Layer &l, int successCount, int blockSize, double qmin) { + double trainSimple(Layer &l, int successCount, double qmin, int reportStep) { + assert(count); + assert(!l.prev); + assert(sizeX == l.size); + assert(sizeY == l.back().size); + assert(successCount > 0); + assert(qmin > 0); + + printf("training: %d, %lf\n", successCount, qmin); + double *res = new double[successCount]; + double *rp = res, *re = res + successCount; + double rsum = 0; + int rcount = 0; + memset(res, 0, sizeof(*res)*successCount); + + int success = 0, total = 0, nextReport = reportStep; + double avg = 0; + for(int i = 0; i < 1000000000; ++i) { + int index = rand() % count; + + double target = avg*0.5; + double q = 0; + for(int i = 0; i < 10; ++i) { + double qq = l.trainPass(trainRatio, x + sizeX*index, y + sizeY*index, target); + if (q < qq) q = qq; + ++total; + if (qq <= target) break; + break; + } + + + rcount += (q > qmin) - (*rp > qmin); + rsum += q - *rp; + *rp++ = q; + if (rp == re) rp = res; + + int cnt = i+1 < successCount ? i+1 : successCount; + avg = rsum/cnt; + + if (q > qmin) success = 0; else ++success; + + if (total >= nextReport || success >= successCount) { + printf(" iterations: %d, error rate: %lf, avg res: %lf\n", total, rcount/(double)cnt, avg); + nextReport = total + reportStep; + } + if (success >= successCount) break; + } + + delete[] res; + printf("done\n"); + return rsum/successCount; + } + + + double trainBlock(Layer &l, int successCount, int blockSize, double qmin, int reportStep = 0) { assert(count); assert(!l.prev); assert(sizeX == l.size); assert(sizeY == l.back().size); assert(blockSize > 0 && qmin > 0); + assert(reportStep >= 0); printf("training: %d, %lf\n", blockSize, qmin); double **blockXY = new double*[blockSize*2]; @@ -53,10 +111,10 @@ public: double qmin3 = qmin2*0.9; int success = 0; - int total = 0; + int total = 0, nextReport = reportStep; int repeats, blockRepeats; double qmax0, qsum0, qmax, qsum; - for(int i = 0; i < 10000; ++i) { + for(int i = 0; i < 1000000; ++i) { for(int i = 0; i < blockSize; ++i) { int index = rand() % count; blockXY[i*2 + 0] = x + sizeX*index; @@ -71,7 +129,7 @@ public: for(int i = 0; i < blockSize; ++i, xy += 2) { double q0 = 0; for(int i = 0; i < 100; ++i) { - double q = l.trainPass(xy[0], xy[1], qmin3); + double q = l.trainPass(trainRatio, xy[0], xy[1], qmin3); if (!i) q0 = q; ++repeats; if (q < qmin3) break; @@ -85,14 +143,18 @@ public: } total += repeats; - printf(" blocks %d (samples: %d, total: %d, repeats: %3d (%lf)): %lf -> %lf, %lf -> %lf\n", - i+1, (i+1)*blockSize, total, blockRepeats-1, repeats/(double)(blockRepeats*blockSize) - 1, qmax0, qmax, qsum0/blockSize, qsum/blockSize); - if (qmax0 > qmin) success = 0; else - if (++success == successCount) break; + if (qmax0 > qmin) success = 0; else ++success; + + if (total >= nextReport || success >= successCount) { + nextReport = total + reportStep; + printf(" blocks %d (samples: %d, total: %d, repeats: %3d (%lf)): %lf -> %lf, %lf -> %lf\n", + i+1, (i+1)*blockSize, total, blockRepeats-1, repeats/(double)(blockRepeats*blockSize) - 1, qmax0, qmax, qsum0/blockSize, qsum/blockSize); + } + if (success >= successCount) break; } - free(blockXY); + delete[] blockXY; printf("done\n"); return qmax0; } @@ -116,7 +178,7 @@ public: unsigned char *data = new unsigned char[testSize*count]; memset(data, 0, testSize*count); if (!fread(data, testSize*count, 1, f)) - return printf("cannot read from file '%s'\n", filename), free(data), fclose(f), false; + return printf("cannot read from file '%s'\n", filename), delete[] data, fclose(f), false; fclose(f); @@ -124,7 +186,8 @@ public: const unsigned char *pd = data; const double delta = 0; double *ey = y + sizeY*count; - for(double *py = y; py < ey; ++py) *py = delta; + for(double *py = y; py < ey; ++py) + *py = delta; for(double *px = x, *py = y; py < ey; py += sizeY) { for(double *ex = px + sizeX; px < ex; ++px, ++pd) *px = *pd/255.0; diff --git a/simple/neural/test-matrix.c b/simple/neural/test-matrix.c new file mode 100644 index 0000000..7f27171 --- /dev/null +++ b/simple/neural/test-matrix.c @@ -0,0 +1,103 @@ + +#include +#include + +#include + + +#define WIDTH 64 +#define HEIGHT 48 +#define COUNT 256 +#define CELL 16 + + +typedef struct { + int x, y, r; +} Point; + + +int matrix[HEIGHT][WIDTH]; +Point order[ (HEIGHT*2 - 1)*(WIDTH*2 - 1) ]; + + +void initOrder() { + for(int y = 1-HEIGHT; y < HEIGHT; ++y) { + for(int x = 1-WIDTH; x < WIDTH; ++x) { + Point *p = &order[(y + HEIGHT - 1)*(WIDTH*2 - 1) + x + WIDTH - 1]; + p->x = x; + p->y = y; + + int sector = x >= 0 && y > 0 ? 1 + : y >= 0 && x < 0 ? 2 + : x <= 0 && y < 0 ? 3 + : y <= 0 && x > 0 ? 4 + : 0; + int subsector = sector % 2 ? abs(x) >= abs(y) : abs(y) >= abs(x); + + p->r = (x*x + y*y)*10 + sector*2 + subsector; + } + } + + int size = sizeof(order)/sizeof(*order); + Point p; + while(1) { + int found = 0; + for(int i = 1; i < size; ++i) { + Point *p1 = &order[i], *p0 = p1 - 1; + if (p0->r > p1->r) { + memcpy(&p, p0, sizeof(p)); + memcpy(p0, p1, sizeof(p)); + memcpy(p1, &p, sizeof(p)); + found = 1; + } + } + if (!found) break; + } +} + + +void prepareMatrix(int x, int y) { + if (x < 0) x = 0; + if (x >= WIDTH) x = WIDTH-1; + if (y < 0) y = 0; + if (y >= HEIGHT) y = HEIGHT-1; + memset(matrix, 0, sizeof(matrix)); + Point *p = order; + for(int i = 0; i < COUNT; ++i) { + int xx, yy; + do { + xx = x + p->x; + yy = y + p->y; + ++p; + } while(xx < 0 || yy < 0 || xx >= WIDTH || yy >= HEIGHT); + matrix[yy][xx] = 1; + } +} + + +void init() { + initOrder(); + windowSetSize(CELL*WIDTH, CELL*HEIGHT); +} + + +void draw() { + noStroke(); + prepareMatrix(mouseX()/CELL, mouseY()/CELL); + for(int y = 0; y < HEIGHT; ++y) { + for(int x = 0; x < WIDTH; ++x) { + fill( matrix[y][x] ? COLOR_BLUE : COLOR_LIGHTBLUE ); + rect(CELL*(x + 0.1), CELL*(y + 0.1), CELL*0.8, CELL*0.8); + } + } +} + + +int main() { + windowSetVariableFrameRate(TRUE); + windowSetInit(&init); + windowSetDraw(&draw); + windowRun(); + return 0; +} +