shithub: util

Download patch

ref: 02a55f299a829e48e10d1591ea449a225e94f1db
parent: e20feacc2e81c47958683946d2ed81c669559047
author: eli <eli@singularity>
date: Sun Jan 5 15:53:24 EST 2025 

ann refactor


--- a/ann.c
+++ /dev/null
@@ -1,880 +1,0 @@
-#include <u.h>
-#include <libc.h>
-#include <ctype.h>
-
-#define RAND_MAX 0xFFFF
-#define fmax(a,b) (a > b? a: b)
-
-typedef struct Ann Ann;
-typedef struct Layer Layer;
-typedef struct Neuron Neuron;
-typedef struct Weights Weights;
-
-struct Ann {
-	int n;
-	double rate;
-	Layer **layers;
-	Weights **weights;
-	Weights **deltas;
-	void *user;
-	void *internal;
-};
-
-struct Layer {
-	int n;
-	Neuron **neurons;
-};
-
-struct Neuron {
-	double (*activation)(Neuron*);
-	double (*gradient)(Neuron*);
-	double steepness;
-	double value;
-	double sum;
-	void *user;
-	void *internal;
-};
-
-struct Weights {
-	int inputs;
-	int outputs;
-	double **values;
-};
-
-double activation_sigmoid(Neuron*);
-double gradient_sigmoid(Neuron*);
-double activation_tanh(Neuron*);
-double gradient_tanh(Neuron*);
-double activation_leaky_relu(Neuron*);
-double gradient_leaky_relu(Neuron*);
-double activation_piece(Neuron*);
-double gradient_piece(Neuron*);
-
-#define ACTIVATION activation_leaky_relu
-#define GRADIENT gradient_leaky_relu
-
-Ann *anncreate(int, ...);
-Ann *anncreatev(int, int*);
-Layer *layercreate(int, double(*)(Neuron*), double(*)(Neuron*));
-Neuron *neuroninit(Neuron*, double (*)(Neuron*), double (*)(Neuron*), double);
-Neuron *neuroncreate(double (*)(Neuron*), double (*)(Neuron*), double);
-Weights *weightsinitrand(Weights*);
-Weights *weightsinitrandscale(Weights*, double);
-Weights *weightsinitdouble(Weights*, double);
-Weights *weightsinitdoubles(Weights*, double*);
-Weights *weightscreate(int, int, int);
-double *annrun(Ann*, double*);
-double anntrain(Ann*, double*, double*);
-
-typedef struct Adam Adam;
-
-struct Adam {
-	double rate;
-	double beta1;
-	double beta2;
-	Weights **first;
-	Weights **second;
-	double epsilon;
-	int timestep;
-};
-
-double anntrain_adam(Ann*, double*, double*);
-double anntrain_adamax(Ann*, double*, double*);
-
-double
-activation_sigmoid(Neuron *in)
-{
-	return 1.0/(1.0+exp(-in->sum));
-}
-
-double
-gradient_sigmoid(Neuron *in)
-{
-	double y = in->value;
-	return y * (1.0 - y);
-}
-
-double
-activation_tanh(Neuron *in)
-{
-	return tanh(in->sum);
-}
-
-double
-gradient_tanh(Neuron *in)
-{
-	return 1.0 - in->value*in->value;
-}
-
-double
-activation_leaky_relu(Neuron *in)
-{
-	if (in->sum > 0)
-		return in->sum;
-	return in->sum * 0.01;
-}
-
-double
-gradient_leaky_relu(Neuron *in)
-{
-	if (in->sum > 0)
-		return 1.0;
-	return 0.01;
-}
-
-double
-activation_piece(Neuron *in)
-{
-	if (in->sum < -0.5)
-		return 0.0;
-	else if (in->sum > 0.5)
-		return 1.0;
-	return (in->sum + 0.5);
-}
-
-double
-gradient_piece(Neuron *in)
-{
-	if (in->sum > -0.5 && in->sum < 0.5)
-		return 1.0;
-	return 0.01;
-}
-
-Weights*
-weightsinitdoubles(Weights *in, double *init)
-{
-	int i, o;
-
-	for (i = 0; i <= in->inputs; i++)
-		for (o = 0; o < in->outputs; o++)
-			in->values[i][o] = init[o];
-
-	return in;
-}
-
-Weights*
-weightsinitdouble(Weights *in, double init)
-{
-	int i, o;
-
-	for (i = 0; i <= in->inputs; i++)
-		for (o = 0; o < in->outputs; o++)
-			in->values[i][o] = init;
-
-	return in;
-}
-
-Weights*
-weightsinitrandscale(Weights *in, double scale)
-{
-	int i, o;
-
-	srand(time(0));
-	for (i = 0; i <= in->inputs; i++)
-		for (o = 0; o < in->outputs; o++)
-			in->values[i][o] = (((double)rand()/RAND_MAX) - 0.5) * scale;
-
-	return in;
-}
-
-Weights*
-weightsinitrand(Weights *in)
-{
-	weightsinitrandscale(in, 2.0);
-	return in;
-}
-
-Neuron*
-neuroninit(Neuron *in, double (*activation)(Neuron*), double (*gradient)(Neuron*), double steepness)
-{
-	in->activation = activation;
-	in->gradient = gradient;
-	in->steepness = steepness;
-	in->value = 1.0;
-	in->sum = 0;
-	return in;
-}
-
-Neuron*
-neuroncreate(double (*activation)(Neuron*), double (*gradient)(Neuron*), double steepness)
-{
-	Neuron *ret = calloc(1, sizeof(Neuron));
-	neuroninit(ret, activation, gradient, steepness);
-	return ret;
-}
-
-Layer*
-layercreate(int num_neurons, double(*activation)(Neuron*), double(*gradient)(Neuron*))
-{
-	Layer *ret = calloc(1, sizeof(Layer));
-	int i;
-
-	ret->n = num_neurons;
-	ret->neurons = calloc(num_neurons+1, sizeof(Neuron*));
-	for (i = 0; i <= ret->n; i++) {
-		ret->neurons[i] = neuroncreate(activation, gradient, 1.0);
-	}
-	return ret;
-}
-
-Weights*
-weightscreate(int inputs, int outputs, int initialize)
-{
-	int i;
-	Weights *ret = calloc(1, sizeof(Weights));
-	ret->inputs = inputs;
-	ret->outputs = outputs;
-	ret->values = calloc(inputs+1, sizeof(double*));
-	for (i = 0; i <= inputs; i++)
-		ret->values[i] = calloc(outputs, sizeof(double));
-	if (initialize)
-		weightsinitrand(ret);
-	return ret;
-}
-
-Ann*
-anncreate(int num_layers, ...)
-{
-	Ann *ret = calloc(1, sizeof(Ann));
-	va_list args;
-	int arg;
-	int i;
-
-	va_start(args, num_layers);
-	ret->n = num_layers;
-	ret->rate = 0.7;
-	ret->layers = calloc(num_layers, sizeof(Layer*));
-	ret->weights = calloc(num_layers-1, sizeof(Weights*));
-	ret->deltas = calloc(num_layers-1, sizeof(Weights*));
-
-	for (i = 0; i < num_layers; i++) {
-		arg = va_arg(args, int);
-		if (arg < 0 || arg > 1000000)
-			arg = 0;
-		ret->layers[i] = layercreate(arg, ACTIVATION, GRADIENT);
-		if (i > 0) {
-			ret->weights[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 1);
-			ret->deltas[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 0);
-		}
-	}
-
-	va_end(args);
-
-	return ret;
-}
-
-Ann*
-anncreatev(int num_layers, int *layers)
-{
-	Ann *ret = calloc(1, sizeof(Ann));
-	int arg;
-	int i;
-
-	ret->n = num_layers;
-	ret->rate = 0.7;
-	ret->layers = calloc(num_layers, sizeof(Layer*));
-	ret->weights = calloc(num_layers-1, sizeof(Weights*));
-	ret->deltas = calloc(num_layers-1, sizeof(Weights*));
-
-	for (i = 0; i < num_layers; i++) {
-		arg = layers[i];
-		if (arg < 0 || arg > 1000000)
-			arg = 0;
-		if (i < (num_layers-1))
-			ret->layers[i] = layercreate(arg, ACTIVATION, GRADIENT);
-		else
-			ret->layers[i] = layercreate(arg, activation_sigmoid, gradient_sigmoid);
-		if (i > 0) {
-			ret->weights[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 1);
-			ret->deltas[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 0);
-		}
-	}
-
-	return ret;
-}
-
-double*
-annrun(Ann *ann, double *input)
-{
-	int l, i, o;
-	int outputs = ann->layers[ann->n - 1]->n;
-	double *ret = calloc(outputs, sizeof(double));
-	Neuron *O;
-	double sum;
-
-	for (i = 0; i < ann->layers[0]->n; i++)
-		ann->layers[0]->neurons[i]->value = input[i];
-
-	for (l = 1; l < ann->n; l++) {
-		for (o = 0; o < ann->layers[l]->n; o++) {
-			O = ann->layers[l]->neurons[o];
-			O->sum = ann->weights[l-1]->values[ann->weights[l-1]->inputs][o]; // bias
-			sum = O->sum;
-			#pragma omp parallel for reduction (+:sum)
-			for (i = 0; i < ann->layers[l-1]->n; i++)
-				sum += ann->layers[l-1]->neurons[i]->value * ann->weights[l-1]->values[i][o];
-			if (sum < -300.0)
-				sum = -300.0;
-			else if (sum > 300.0)
-				sum = 300.0;
-			O->sum = sum;
-			O->value = O->activation(O);
-		}
-	}
-
-	for (o = 0; o < outputs; o++)
-		ret[o] = ann->layers[ann->n - 1]->neurons[o]->value;
-
-	return ret;
-}
-
-double
-anntrain(Ann *ann, double *inputs, double *outputs)
-{
-	double *error = annrun(ann, inputs);
-	double ret = 0.0;
-	int noutputs = ann->layers[ann->n-1]->n;
-	double acc, sum;
-	int o, i, w, n;
-	Neuron *O, *I;
-	Weights *W, *D, *D2;
-
-	for (o = 0; o < noutputs; o++) {
-		// error = outputs[o] - result
-		error[o] -= outputs[o];
-		error[o] = -error[o];
-		ret += pow(error[o], 2.0) * 0.5;
-		if (error[o] < -.9999999)
-			error[o] = -17.0;
-		else if (error[o] > .9999999)
-			error[o] = 17.0;
-		else
-			error[o] = log((1.0 + error[o]) / (1.0 - error[o]));
-	}
-	D = ann->deltas[ann->n-2];
-	weightsinitdoubles(D, error);
-	for (i = 0; i < (ann->n-2); i++) {
-		D = ann->deltas[i];
-		weightsinitdouble(D, 1.0);
-	}
-
-	// backpropagate MSE
-	D2 = ann->deltas[ann->n-2];
-	for (w = ann->n-2; w >= 0; w--) {
-		D = ann->deltas[w];
-
-		for (o = 0; o < ann->layers[w+1]->n; o++) {
-			O = ann->layers[w+1]->neurons[o];
-			acc = O->gradient(O) * O->steepness;
-			sum = 1.0;
-			if (D2 != D) {
-				W = ann->weights[w + 1];
-				sum = 0.0;
-				#pragma omp parallel for reduction (+:sum)
-				for (n = 0; n < D2->outputs; n++)
-					sum += D2->values[o][n] * W->values[o][n];
-			}
-			for (i = 0; i <= ann->layers[w]->n; i++) {
-			 	D->values[i][o] *= acc * sum;
-			}
-		}
-
-		D2 = D;
-	}
-
-	// update weights
-	for (w = 0; w < ann->n-1; w++) {
-		W = ann->weights[w];
-		D = ann->deltas[w];
-
-		for (i = 0; i <= W->inputs; i++) {
-			I = ann->layers[w]->neurons[i];
-			for (o = 0; o < W->outputs; o++) {
-				W->values[i][o] += D->values[i][o] * ann->rate * I->value;
-			}
-		}
-	}
-
-	free(error);
-	return ret;
-}
-
-Ann*
-adaminit(Ann *ann)
-{
-	int i;
-	Adam *I = calloc(1, sizeof(Adam));
-
-	I->rate = 0.001;
-	I->beta1 = 0.9;
-	I->beta2 = 0.999;
-	I->epsilon = 10e-8;
-	I->timestep = 0;
-	I->first = calloc(ann->n-1, sizeof(Weights*));
-	I->second = calloc(ann->n-1, sizeof(Weights*));
-
-	for (i = 0; i < (ann->n-1); i++) {
-		I->first[i] = weightscreate(ann->layers[i]->n, ann->layers[i+1]->n, 0);
-		I->second[i] = weightscreate(ann->layers[i]->n, ann->layers[i+1]->n, 0);
-	}
-
-	ann->internal = I;
-
-	return ann;
-}
-
-double
-anntrain_adam(Ann *ann, double *inputs, double *outputs)
-{
-	double *error = annrun(ann, inputs);
-	double ret = 0.0;
-	int noutputs = ann->layers[ann->n-1]->n;
-	double acc, sum, m, v;
-	int o, i, w, n;
-	Neuron *O, *I;
-	Weights *W, *D, *D2, *M, *V;
-	Adam *annI;
-
-	if (ann->internal == 0)
-		adaminit(ann);
-	annI = ann->internal;
-	annI->timestep++;
-
-	for (o = 0; o < noutputs; o++) {
-		// error = outputs[o] - result
-		error[o] -= outputs[o];
-		error[o] = -error[o];
-		ret += pow(error[o], 2.0) * 0.5;
-	}
-	D = ann->deltas[ann->n-2];
-	weightsinitdoubles(D, error);
-	for (i = 0; i < (ann->n-2); i++) {
-		D = ann->deltas[i];
-		weightsinitdouble(D, 1.0);
-	}
-
-	// backpropagate MSE
-	D2 = ann->deltas[ann->n-2];
-	for (w = ann->n-2; w >= 0; w--) {
-		D = ann->deltas[w];
-		M = annI->first[w];
-		V = annI->second[w];
-
-		for (o = 0; o < ann->layers[w+1]->n; o++) {
-			O = ann->layers[w+1]->neurons[o];
-			acc = O->gradient(O) * O->steepness;
-			sum = 1.0;
-			if (D2 != D) {
-				W = ann->weights[w+1];
-				sum = 0.0;
-				#pragma omp parallel for reduction (+:sum)
-				for (n = 0; n < D2->outputs; n++)
-					sum += D2->values[o][n] * W->values[o][n];
-			}
-			for (i = 0; i <= ann->layers[w]->n; i++) {
-				I = ann->layers[w]->neurons[i];
-			 	D->values[i][o] *= acc * sum;
-				M->values[i][o] *= annI->beta1;
-				M->values[i][o] += (1.0 - annI->beta1) * D->values[i][o] * I->value;
-				V->values[i][o] *= annI->beta2;
-				V->values[i][o] += (1.0 - annI->beta2) * D->values[i][o] * D->values[i][o] * I->value * I->value;
-			}
-		}
-
-		D2 = D;
-	}
-
-	// update weights
-	for (w = 0; w < ann->n-1; w++) {
-		W = ann->weights[w];
-		M = annI->first[w];
-		V = annI->second[w];
-
-		for (i = 0; i <= W->inputs; i++) {
-			for (o = 0; o < W->outputs; o++) {
-				m = M->values[i][o] / (annI->timestep < 100? (1.0 - pow(annI->beta1, annI->timestep)): 1.0);
-				v = V->values[i][o] / (annI->timestep < 10000? (1.0 - pow(annI->beta2, annI->timestep)): 1.0);
-				W->values[i][o] += (m / (sqrt(v) + annI->epsilon)) * annI->rate;
-			}
-		}
-	}
-
-	free(error);
-	return ret;
-}
-
-double
-anntrain_adamax(Ann *ann, double *inputs, double *outputs)
-{
-	double *error = annrun(ann, inputs);
-	double ret = 0.0;
-	int noutputs = ann->layers[ann->n-1]->n;
-	double acc, sum, m, v;
-	int o, i, w, n;
-	Neuron *O, *I;
-	Weights *W, *D, *D2, *M, *V;
-	Adam *annI;
-
-	if (ann->internal == 0)
-		adaminit(ann);
-	annI = ann->internal;
-	annI->rate = 0.002;
-	annI->timestep++;
-
-	for (o = 0; o < noutputs; o++) {
-		// error = outputs[o] - result
-		error[o] -= outputs[o];
-		error[o] = -error[o];
-		ret += pow(error[o], 2.0) * 0.5;
-	}
-	D = ann->deltas[ann->n-2];
-	weightsinitdoubles(D, error);
-	for (i = 0; i < (ann->n-2); i++) {
-		D = ann->deltas[i];
-		weightsinitdouble(D, 1.0);
-	}
-
-	// backpropagate MSE
-	D2 = ann->deltas[ann->n-2];
-	for (w = ann->n-2; w >= 0; w--) {
-		D = ann->deltas[w];
-		M = annI->first[w];
-		V = annI->second[w];
-
-		for (o = 0; o < ann->layers[w+1]->n; o++) {
-			O = ann->layers[w+1]->neurons[o];
-			acc = O->gradient(O) * O->steepness;
-			sum = 1.0;
-			if (D2 != D) {
-				W = ann->weights[w+1];
-				sum = 0.0;
-				#pragma omp parallel for reduction (+:sum)
-				for (n = 0; n < D2->outputs; n++)
-					sum += D2->values[o][n] * W->values[o][n];
-			}
-			for (i = 0; i <= ann->layers[w]->n; i++) {
-				I = ann->layers[w]->neurons[i];
-			 	D->values[i][o] *= acc * sum;
-				M->values[i][o] *= annI->beta1;
-				M->values[i][o] += (1.0 - annI->beta1) * D->values[i][o] * I->value;
-				V->values[i][o] = fmax(V->values[i][o] * annI->beta2, fabs(D->values[i][o] * I->value));
-			}
-		}
-
-		D2 = D;
-	}
-
-	// update weights
-	for (w = 0; w < ann->n-1; w++) {
-		W = ann->weights[w];
-		M = annI->first[w];
-		V = annI->second[w];
-
-		for (i = 0; i <= W->inputs; i++) {
-			for (o = 0; o < W->outputs; o++) {
-				m = M->values[i][o];
-				v = V->values[i][o];
-				W->values[i][o] += (annI->rate/(1.0 - (annI->timestep < 100? pow(annI->beta1, annI->timestep): 0.0))) * (m/v);
-			}
-		}
-	}
-
-	free(error);
-	return ret;
-}
-
-int
-saveann(char *filename, Ann *ann)
-{
-	Weights *W;
-	int i, j, k;
-	int fd = create(filename, OWRITE, 0644);
-	if (fd < 0) {
-		perror("create");
-		return -1;
-	}
-
-	fprint(fd, "%d\n", ann->n);
-	for (i = 0; i < ann->n; i++)
-		fprint(fd, "%d\n", ann->layers[i]->n);
-
-	for (i = 0; i < (ann->n - 1); i++) {
-		W = ann->weights[i];
-		for (j = 0; j <= W->inputs; j++)
-			for (k = 0; k < W->outputs; k++)
-				fprint(fd, "%f\n", W->values[j][k]);
-	}
-
-	close(fd);
-	return 0;
-}
-
-char *retline = nil;
-
-char*
-readline(int fd)
-{
-	static int length;
-	static int offset;
-	char *epos;
-	int r;
-
-	if (retline == nil) {
-		length = 0x1000;
-		offset = 0;
-		retline = malloc(length);
-		retline[offset] = '\0';
-	}
-
-	r = strlen(retline);
-	if (r > 0) {
-		r++;
-		memmove(retline, &retline[r], length - r);
-		length -= r;
-		offset = strlen(retline);
-	}
-
-	if (length < 0x1000) {
-		retline = realloc(retline, length + 0x1000);
-		length += 0x1000;
-	}
-
-	do {
-		epos = strchr(retline, '\n');
-		if (epos != nil) {
-			*epos = '\0';
-			return retline;
-		}
-
-		r = read(fd, &retline[offset], length - offset - 1);
-		if (r < 0) {
-			perror("read");
-			exits("read");
-		}
-		if (r > 0) {
-			offset += r;
-			retline[offset] = '\0';
-			length += r;
-			retline = realloc(retline, length);
-		}
-	} while(r != 0);
-
-	return nil;
-}
-
-char*
-sreadline(int fd)
-{
-	char *ret = readline(fd);
-	if (ret == nil) {
-		fprint(2, "error: early end of file\n");
-		exits("eof");
-	}
-	return ret;
-}
-
-Ann*
-loadann(char *filename)
-{
-	Weights *W;
-	Ann *ann = nil;
-	char *buf;
-	int i, j, k;
-	int num_layers, *layers;
-	int fd = open(filename, OREAD);
-	if (fd < 0) {
-		perror("open");
-		return ann;
-	}
-
-	buf = sreadline(fd);
-	num_layers = atoi(buf);
-
-	if (num_layers < 2) {
-		fprint(2, "num_layers < 2\n");
-		return ann;
-	}
-
-	layers = calloc(num_layers, sizeof(int));
-	for (i = 0; i < num_layers; i++) {
-		buf = sreadline(fd);
-		layers[i] = atoi(buf);
-	}
-
-	ann = anncreatev(num_layers, layers);
-	for (i = 0; i < (ann->n - 1); i++) {
-		W = ann->weights[i];
-		for (j = 0; j <= W->inputs; j++)
-			for (k = 0; k < W->outputs; k++) {
-				buf = sreadline(fd);
-				W->values[j][k] = atof(buf);
-			}
-	}
-
-	free(retline);
-	retline = nil;
-
-	return ann;
-}
-
-void
-usage(char **argv)
-{
-	fprint(2, "usage: %s [-train] filename [num_layers num_input_layer ... num_output_layer]\n", argv[0]);
-	exits("usage");
-}
-
-void
-main(int argc, char **argv)
-{
-	Ann *ann;
-	char *filename;
-	int train;
-	Dir *dir;
-	int num_layers = 0;
-	int *layers = nil;
-	int i;
-	char *line;
-	double *input;
-	double *output = nil;
-	double *runoutput;
-	int ninput;
-	int noutput;
-	int offset;
-	double f;
-	int trainline;
-	int nline;
-
-	train = 0;
-
-	if (argc < 2)
-		usage(argv);
-
-	filename = argv[1];
-
-	if (argv[1][0] == '-' && argv[1][1] == 't') {
-		if (argc < 3)
-			usage(argv);
-
-		train = 1;
-		filename = argv[2];
-	}
-
-	ann = nil;
-	dir = dirstat(filename);
-	if (dir != nil) {
-		free(dir);
-		ann = loadann(filename);
-		if (ann == nil)
-			exits("loadann");
-	}
-
-	if (argc >= (train + 3)) {
-		num_layers = atoi(argv[train + 2]);
-
-		if (num_layers < 2 || argc != (train + 3 + num_layers))
-			usage(argv);
-
-		layers = calloc(num_layers, sizeof(int));
-
-		for (i = 0; i < num_layers; i++)
-			layers[i] = atoi(argv[train + 3 + i]);
-	}
-
-	if (num_layers > 0) {
-		if (ann != nil) {
-			if (ann->n != num_layers) {
-				fprint(2, "num_layers: %d != %d\n", ann->n, num_layers);
-				exits("num_layers");
-			}
-
-			for (i = 0; i < num_layers; i++) {
-				if (layers[i] != ann->layers[i]->n) {
-					fprint(2, "num_layer_%d: %d != %d\n", i, layers[i], ann->layers[i]->n);
-					exits("num_layer");
-				}
-			}
-		} else {
-			ann = anncreatev(num_layers, layers);
-			if (ann == nil)
-				exits("anncreatev");
-		}
-	}
-
-	if (ann == nil) {
-		fprint(2, "file not found: %s\n", filename);
-		exits("file not found");
-	}
-
-	ninput = ann->layers[0]->n;
-	noutput = ann->layers[ann->n - 1]->n;
-	input = calloc(ninput, sizeof(double));
-	if (train == 1)
-		output = calloc(noutput, sizeof(double));
-
-	trainline = 0;
-	nline = ninput;
-
-	do {
-		int i = 0;
-		while ((line = readline(0)) != nil) {
-			do {
-				if (strlen(line) == 0)
-					break;
-				while(isspace(*line))
-					line++;
-				if (strlen(line) == 0)
-					break;
-				offset = 0;
-				while (isdigit(line[offset]) || line[offset] == '.' || line[offset] == '-')
-					offset++;
-				if (!isspace(line[offset]) && line[offset] != '\0') {
-					fprint(2, "input error: %s\n", line);
-					exits("input");
-				}
-				f = atof(line);
-				if (trainline == 0) {
-					input[i] = f;
-					i++;
-				} else {
-					output[i] = f;
-					i++;
-				}
-				line = &line[offset];
-				while(isspace(*line))
-					line++;
-			} while(i < nline && strlen(line) > 0);
-
-			if (i == nline) {
-				if (trainline == 0) {
-					runoutput = annrun(ann, input);
-					for (i = 0; i < noutput; i++)
-/*						if (runoutput[i] == 0.0)
-							print("0%c", (i == (noutput-1))? '\n': ' ');
-						else if (runoutput[i] == 1.0)
-							print("1%c", (i == (noutput-1))? '\n': ' ');
-						else */
-							print("%f%c", runoutput[i], (i == (noutput-1))? '\n': ' ');
-					free(runoutput);
-				}
-
-				if (train == 1) {
-					if (trainline == 0) {
-						trainline = 1;
-						nline = noutput;
-					} else {
-						anntrain(ann, input, output);
-						trainline = 0;
-						nline = ninput;
-					}
-				}
-				i = 0;
-			}
-		}
-	} while(line != nil);
-
-	if (saveann(filename, ann) != 0)
-		exits("saveann");
-
-	exits(nil);
-}
--- /dev/null
+++ b/ann/ann.c
@@ -1,0 +1,513 @@
+#include <u.h>
+#include <libc.h>
+#include <ctype.h>
+
+#define RAND_MAX 0xFFFF
+#define fmax(a,b) (a > b? a: b)
+
+#include "ann.h"
+
+#define ACTIVATION activation_leaky_relu
+#define GRADIENT gradient_leaky_relu
+
+double
+activation_sigmoid(Neuron *in)
+{
+	return 1.0/(1.0+exp(-in->sum));
+}
+
+double
+gradient_sigmoid(Neuron *in)
+{
+	double y = in->value;
+	return y * (1.0 - y);
+}
+
+double
+activation_tanh(Neuron *in)
+{
+	return tanh(in->sum);
+}
+
+double
+gradient_tanh(Neuron *in)
+{
+	return 1.0 - in->value*in->value;
+}
+
+double
+activation_leaky_relu(Neuron *in)
+{
+	if (in->sum > 0)
+		return in->sum;
+	return in->sum * 0.01;
+}
+
+double
+gradient_leaky_relu(Neuron *in)
+{
+	if (in->sum > 0)
+		return 1.0;
+	return 0.01;
+}
+
+double
+activation_piece(Neuron *in)
+{
+	if (in->sum < -0.5)
+		return 0.0;
+	else if (in->sum > 0.5)
+		return 1.0;
+	return (in->sum + 0.5);
+}
+
+double
+gradient_piece(Neuron *in)
+{
+	if (in->sum > -0.5 && in->sum < 0.5)
+		return 1.0;
+	return 0.01;
+}
+
+Weights*
+weightsinitdoubles(Weights *in, double *init)
+{
+	int i, o;
+
+	for (i = 0; i <= in->inputs; i++)
+		for (o = 0; o < in->outputs; o++)
+			in->values[i][o] = init[o];
+
+	return in;
+}
+
+Weights*
+weightsinitdouble(Weights *in, double init)
+{
+	int i, o;
+
+	for (i = 0; i <= in->inputs; i++)
+		for (o = 0; o < in->outputs; o++)
+			in->values[i][o] = init;
+
+	return in;
+}
+
+Weights*
+weightsinitrandscale(Weights *in, double scale)
+{
+	int i, o;
+
+	srand(time(0));
+	for (i = 0; i <= in->inputs; i++)
+		for (o = 0; o < in->outputs; o++)
+			in->values[i][o] = (((double)rand()/RAND_MAX) - 0.5) * scale;
+
+	return in;
+}
+
+Weights*
+weightsinitrand(Weights *in)
+{
+	weightsinitrandscale(in, 2.0);
+	return in;
+}
+
+Neuron*
+neuroninit(Neuron *in, double (*activation)(Neuron*), double (*gradient)(Neuron*), double steepness)
+{
+	in->activation = activation;
+	in->gradient = gradient;
+	in->steepness = steepness;
+	in->value = 1.0;
+	in->sum = 0;
+	return in;
+}
+
+Neuron*
+neuroncreate(double (*activation)(Neuron*), double (*gradient)(Neuron*), double steepness)
+{
+	Neuron *ret = calloc(1, sizeof(Neuron));
+	neuroninit(ret, activation, gradient, steepness);
+	return ret;
+}
+
+Layer*
+layercreate(int num_neurons, double(*activation)(Neuron*), double(*gradient)(Neuron*))
+{
+	Layer *ret = calloc(1, sizeof(Layer));
+	int i;
+
+	ret->n = num_neurons;
+	ret->neurons = calloc(num_neurons+1, sizeof(Neuron*));
+	for (i = 0; i <= ret->n; i++) {
+		ret->neurons[i] = neuroncreate(activation, gradient, 1.0);
+	}
+	return ret;
+}
+
+Weights*
+weightscreate(int inputs, int outputs, int initialize)
+{
+	int i;
+	Weights *ret = calloc(1, sizeof(Weights));
+	ret->inputs = inputs;
+	ret->outputs = outputs;
+	ret->values = calloc(inputs+1, sizeof(double*));
+	for (i = 0; i <= inputs; i++)
+		ret->values[i] = calloc(outputs, sizeof(double));
+	if (initialize)
+		weightsinitrand(ret);
+	return ret;
+}
+
+Ann*
+anncreate(int num_layers, ...)
+{
+	Ann *ret = calloc(1, sizeof(Ann));
+	va_list args;
+	int arg;
+	int i;
+
+	va_start(args, num_layers);
+	ret->n = num_layers;
+	ret->rate = 0.7;
+	ret->layers = calloc(num_layers, sizeof(Layer*));
+	ret->weights = calloc(num_layers-1, sizeof(Weights*));
+	ret->deltas = calloc(num_layers-1, sizeof(Weights*));
+
+	for (i = 0; i < num_layers; i++) {
+		arg = va_arg(args, int);
+		if (arg < 0 || arg > 1000000)
+			arg = 0;
+		ret->layers[i] = layercreate(arg, ACTIVATION, GRADIENT);
+		if (i > 0) {
+			ret->weights[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 1);
+			ret->deltas[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 0);
+		}
+	}
+
+	va_end(args);
+
+	return ret;
+}
+
+Ann*
+anncreatev(int num_layers, int *layers)
+{
+	Ann *ret = calloc(1, sizeof(Ann));
+	int arg;
+	int i;
+
+	ret->n = num_layers;
+	ret->rate = 0.7;
+	ret->layers = calloc(num_layers, sizeof(Layer*));
+	ret->weights = calloc(num_layers-1, sizeof(Weights*));
+	ret->deltas = calloc(num_layers-1, sizeof(Weights*));
+
+	for (i = 0; i < num_layers; i++) {
+		arg = layers[i];
+		if (arg < 0 || arg > 1000000)
+			arg = 0;
+		if (i < (num_layers-1))
+			ret->layers[i] = layercreate(arg, ACTIVATION, GRADIENT);
+		else
+			ret->layers[i] = layercreate(arg, activation_sigmoid, gradient_sigmoid);
+		if (i > 0) {
+			ret->weights[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 1);
+			ret->deltas[i-1] = weightscreate(ret->layers[i-1]->n, ret->layers[i]->n, 0);
+		}
+	}
+
+	return ret;
+}
+
+double*
+annrun(Ann *ann, double *input)
+{
+	int l, i, o;
+	int outputs = ann->layers[ann->n - 1]->n;
+	double *ret = calloc(outputs, sizeof(double));
+	Neuron *O;
+	double sum;
+
+	for (i = 0; i < ann->layers[0]->n; i++)
+		ann->layers[0]->neurons[i]->value = input[i];
+
+	for (l = 1; l < ann->n; l++) {
+		for (o = 0; o < ann->layers[l]->n; o++) {
+			O = ann->layers[l]->neurons[o];
+			O->sum = ann->weights[l-1]->values[ann->weights[l-1]->inputs][o]; // bias
+			sum = O->sum;
+			#pragma omp parallel for reduction (+:sum)
+			for (i = 0; i < ann->layers[l-1]->n; i++)
+				sum += ann->layers[l-1]->neurons[i]->value * ann->weights[l-1]->values[i][o];
+			if (sum < -300.0)
+				sum = -300.0;
+			else if (sum > 300.0)
+				sum = 300.0;
+			O->sum = sum;
+			O->value = O->activation(O);
+		}
+	}
+
+	for (o = 0; o < outputs; o++)
+		ret[o] = ann->layers[ann->n - 1]->neurons[o]->value;
+
+	return ret;
+}
+
+double
+anntrain(Ann *ann, double *inputs, double *outputs)
+{
+	double *error = annrun(ann, inputs);
+	double ret = 0.0;
+	int noutputs = ann->layers[ann->n-1]->n;
+	double acc, sum;
+	int o, i, w, n;
+	Neuron *O, *I;
+	Weights *W, *D, *D2;
+
+	for (o = 0; o < noutputs; o++) {
+		// error = outputs[o] - result
+		error[o] -= outputs[o];
+		error[o] = -error[o];
+		ret += pow(error[o], 2.0) * 0.5;
+		if (error[o] < -.9999999)
+			error[o] = -17.0;
+		else if (error[o] > .9999999)
+			error[o] = 17.0;
+		else
+			error[o] = log((1.0 + error[o]) / (1.0 - error[o]));
+	}
+	D = ann->deltas[ann->n-2];
+	weightsinitdoubles(D, error);
+	for (i = 0; i < (ann->n-2); i++) {
+		D = ann->deltas[i];
+		weightsinitdouble(D, 1.0);
+	}
+
+	// backpropagate MSE
+	D2 = ann->deltas[ann->n-2];
+	for (w = ann->n-2; w >= 0; w--) {
+		D = ann->deltas[w];
+
+		for (o = 0; o < ann->layers[w+1]->n; o++) {
+			O = ann->layers[w+1]->neurons[o];
+			acc = O->gradient(O) * O->steepness;
+			sum = 1.0;
+			if (D2 != D) {
+				W = ann->weights[w + 1];
+				sum = 0.0;
+				#pragma omp parallel for reduction (+:sum)
+				for (n = 0; n < D2->outputs; n++)
+					sum += D2->values[o][n] * W->values[o][n];
+			}
+			for (i = 0; i <= ann->layers[w]->n; i++) {
+			 	D->values[i][o] *= acc * sum;
+			}
+		}
+
+		D2 = D;
+	}
+
+	// update weights
+	for (w = 0; w < ann->n-1; w++) {
+		W = ann->weights[w];
+		D = ann->deltas[w];
+
+		for (i = 0; i <= W->inputs; i++) {
+			I = ann->layers[w]->neurons[i];
+			for (o = 0; o < W->outputs; o++) {
+				W->values[i][o] += D->values[i][o] * ann->rate * I->value;
+			}
+		}
+	}
+
+	free(error);
+	return ret;
+}
+
+Ann*
+adaminit(Ann *ann)
+{
+	int i;
+	Adam *I = calloc(1, sizeof(Adam));
+
+	I->rate = 0.001;
+	I->beta1 = 0.9;
+	I->beta2 = 0.999;
+	I->epsilon = 10e-8;
+	I->timestep = 0;
+	I->first = calloc(ann->n-1, sizeof(Weights*));
+	I->second = calloc(ann->n-1, sizeof(Weights*));
+
+	for (i = 0; i < (ann->n-1); i++) {
+		I->first[i] = weightscreate(ann->layers[i]->n, ann->layers[i+1]->n, 0);
+		I->second[i] = weightscreate(ann->layers[i]->n, ann->layers[i+1]->n, 0);
+	}
+
+	ann->internal = I;
+
+	return ann;
+}
+
+double
+anntrain_adam(Ann *ann, double *inputs, double *outputs)
+{
+	double *error = annrun(ann, inputs);
+	double ret = 0.0;
+	int noutputs = ann->layers[ann->n-1]->n;
+	double acc, sum, m, v;
+	int o, i, w, n;
+	Neuron *O, *I;
+	Weights *W, *D, *D2, *M, *V;
+	Adam *annI;
+
+	if (ann->internal == 0)
+		adaminit(ann);
+	annI = ann->internal;
+	annI->timestep++;
+
+	for (o = 0; o < noutputs; o++) {
+		// error = outputs[o] - result
+		error[o] -= outputs[o];
+		error[o] = -error[o];
+		ret += pow(error[o], 2.0) * 0.5;
+	}
+	D = ann->deltas[ann->n-2];
+	weightsinitdoubles(D, error);
+	for (i = 0; i < (ann->n-2); i++) {
+		D = ann->deltas[i];
+		weightsinitdouble(D, 1.0);
+	}
+
+	// backpropagate MSE
+	D2 = ann->deltas[ann->n-2];
+	for (w = ann->n-2; w >= 0; w--) {
+		D = ann->deltas[w];
+		M = annI->first[w];
+		V = annI->second[w];
+
+		for (o = 0; o < ann->layers[w+1]->n; o++) {
+			O = ann->layers[w+1]->neurons[o];
+			acc = O->gradient(O) * O->steepness;
+			sum = 1.0;
+			if (D2 != D) {
+				W = ann->weights[w+1];
+				sum = 0.0;
+				#pragma omp parallel for reduction (+:sum)
+				for (n = 0; n < D2->outputs; n++)
+					sum += D2->values[o][n] * W->values[o][n];
+			}
+			for (i = 0; i <= ann->layers[w]->n; i++) {
+				I = ann->layers[w]->neurons[i];
+			 	D->values[i][o] *= acc * sum;
+				M->values[i][o] *= annI->beta1;
+				M->values[i][o] += (1.0 - annI->beta1) * D->values[i][o] * I->value;
+				V->values[i][o] *= annI->beta2;
+				V->values[i][o] += (1.0 - annI->beta2) * D->values[i][o] * D->values[i][o] * I->value * I->value;
+			}
+		}
+
+		D2 = D;
+	}
+
+	// update weights
+	for (w = 0; w < ann->n-1; w++) {
+		W = ann->weights[w];
+		M = annI->first[w];
+		V = annI->second[w];
+
+		for (i = 0; i <= W->inputs; i++) {
+			for (o = 0; o < W->outputs; o++) {
+				m = M->values[i][o] / (annI->timestep < 100? (1.0 - pow(annI->beta1, annI->timestep)): 1.0);
+				v = V->values[i][o] / (annI->timestep < 10000? (1.0 - pow(annI->beta2, annI->timestep)): 1.0);
+				W->values[i][o] += (m / (sqrt(v) + annI->epsilon)) * annI->rate;
+			}
+		}
+	}
+
+	free(error);
+	return ret;
+}
+
+double
+anntrain_adamax(Ann *ann, double *inputs, double *outputs)
+{
+	double *error = annrun(ann, inputs);
+	double ret = 0.0;
+	int noutputs = ann->layers[ann->n-1]->n;
+	double acc, sum, m, v;
+	int o, i, w, n;
+	Neuron *O, *I;
+	Weights *W, *D, *D2, *M, *V;
+	Adam *annI;
+
+	if (ann->internal == 0)
+		adaminit(ann);
+	annI = ann->internal;
+	annI->rate = 0.002;
+	annI->timestep++;
+
+	for (o = 0; o < noutputs; o++) {
+		// error = outputs[o] - result
+		error[o] -= outputs[o];
+		error[o] = -error[o];
+		ret += pow(error[o], 2.0) * 0.5;
+	}
+	D = ann->deltas[ann->n-2];
+	weightsinitdoubles(D, error);
+	for (i = 0; i < (ann->n-2); i++) {
+		D = ann->deltas[i];
+		weightsinitdouble(D, 1.0);
+	}
+
+	// backpropagate MSE
+	D2 = ann->deltas[ann->n-2];
+	for (w = ann->n-2; w >= 0; w--) {
+		D = ann->deltas[w];
+		M = annI->first[w];
+		V = annI->second[w];
+
+		for (o = 0; o < ann->layers[w+1]->n; o++) {
+			O = ann->layers[w+1]->neurons[o];
+			acc = O->gradient(O) * O->steepness;
+			sum = 1.0;
+			if (D2 != D) {
+				W = ann->weights[w+1];
+				sum = 0.0;
+				#pragma omp parallel for reduction (+:sum)
+				for (n = 0; n < D2->outputs; n++)
+					sum += D2->values[o][n] * W->values[o][n];
+			}
+			for (i = 0; i <= ann->layers[w]->n; i++) {
+				I = ann->layers[w]->neurons[i];
+			 	D->values[i][o] *= acc * sum;
+				M->values[i][o] *= annI->beta1;
+				M->values[i][o] += (1.0 - annI->beta1) * D->values[i][o] * I->value;
+				V->values[i][o] = fmax(V->values[i][o] * annI->beta2, fabs(D->values[i][o] * I->value));
+			}
+		}
+
+		D2 = D;
+	}
+
+	// update weights
+	for (w = 0; w < ann->n-1; w++) {
+		W = ann->weights[w];
+		M = annI->first[w];
+		V = annI->second[w];
+
+		for (i = 0; i <= W->inputs; i++) {
+			for (o = 0; o < W->outputs; o++) {
+				m = M->values[i][o];
+				v = V->values[i][o];
+				W->values[i][o] += (annI->rate/(1.0 - (annI->timestep < 100? pow(annI->beta1, annI->timestep): 0.0))) * (m/v);
+			}
+		}
+	}
+
+	free(error);
+	return ret;
+}
--- /dev/null
+++ b/ann/ann.h
@@ -1,0 +1,77 @@
+typedef struct Ann Ann;
+typedef struct Layer Layer;
+typedef struct Neuron Neuron;
+typedef struct Weights Weights;
+
+struct Ann {
+	int n;
+	double rate;
+	Layer **layers;
+	Weights **weights;
+	Weights **deltas;
+	void *user;
+	void *internal;
+};
+
+struct Layer {
+	int n;
+	Neuron **neurons;
+};
+
+struct Neuron {
+	double (*activation)(Neuron*);
+	double (*gradient)(Neuron*);
+	double steepness;
+	double value;
+	double sum;
+	void *user;
+	void *internal;
+};
+
+struct Weights {
+	int inputs;
+	int outputs;
+	double **values;
+};
+
+double activation_sigmoid(Neuron*);
+double gradient_sigmoid(Neuron*);
+double activation_tanh(Neuron*);
+double gradient_tanh(Neuron*);
+double activation_leaky_relu(Neuron*);
+double gradient_leaky_relu(Neuron*);
+double activation_piece(Neuron*);
+double gradient_piece(Neuron*);
+
+Ann *anncreate(int, ...);
+Ann *anncreatev(int, int*);
+Layer *layercreate(int, double(*)(Neuron*), double(*)(Neuron*));
+Neuron *neuroninit(Neuron*, double (*)(Neuron*), double (*)(Neuron*), double);
+Neuron *neuroncreate(double (*)(Neuron*), double (*)(Neuron*), double);
+Weights *weightsinitrand(Weights*);
+Weights *weightsinitrandscale(Weights*, double);
+Weights *weightsinitdouble(Weights*, double);
+Weights *weightsinitdoubles(Weights*, double*);
+Weights *weightscreate(int, int, int);
+double *annrun(Ann*, double*);
+double anntrain(Ann*, double*, double*);
+
+typedef struct Adam Adam;
+
+struct Adam {
+	double rate;
+	double beta1;
+	double beta2;
+	Weights **first;
+	Weights **second;
+	double epsilon;
+	int timestep;
+};
+
+double anntrain_adam(Ann*, double*, double*);
+double anntrain_adamax(Ann*, double*, double*);
+
+int annsave(char *filename, Ann *ann);
+Ann* annload(char *filename);
+char *readline(int fd);
+char *sreadline(int fd);
--- /dev/null
+++ b/ann/annio.c
@@ -1,0 +1,137 @@
+#include <u.h>
+#include <libc.h>
+#include "ann.h"
+
+int
+annsave(char *filename, Ann *ann)
+{
+	Weights *W;
+	int i, j, k;
+	int fd = create(filename, OWRITE, 0644);
+	if (fd < 0) {
+		perror("create");
+		return -1;
+	}
+
+	fprint(fd, "%d\n", ann->n);
+	for (i = 0; i < ann->n; i++)
+		fprint(fd, "%d\n", ann->layers[i]->n);
+
+	for (i = 0; i < (ann->n - 1); i++) {
+		W = ann->weights[i];
+		for (j = 0; j <= W->inputs; j++)
+			for (k = 0; k < W->outputs; k++)
+				fprint(fd, "%f\n", W->values[j][k]);
+	}
+
+	close(fd);
+	return 0;
+}
+
+char *retline = nil;
+
+char*
+readline(int fd)
+{
+	static int length;
+	static int offset;
+	char *epos;
+	int r;
+
+	if (retline == nil) {
+		length = 0x1000;
+		offset = 0;
+		retline = malloc(length);
+		retline[offset] = '\0';
+	}
+
+	r = strlen(retline);
+	if (r > 0) {
+		r++;
+		memmove(retline, &retline[r], length - r);
+		length -= r;
+		offset = strlen(retline);
+	}
+
+	if (length < 0x1000) {
+		retline = realloc(retline, length + 0x1000);
+		length += 0x1000;
+	}
+
+	do {
+		epos = strchr(retline, '\n');
+		if (epos != nil) {
+			*epos = '\0';
+			return retline;
+		}
+
+		r = read(fd, &retline[offset], length - offset - 1);
+		if (r < 0) {
+			perror("read");
+			exits("read");
+		}
+		if (r > 0) {
+			offset += r;
+			retline[offset] = '\0';
+			length += r;
+			retline = realloc(retline, length);
+		}
+	} while(r != 0);
+
+	return nil;
+}
+
+char*
+sreadline(int fd)
+{
+	char *ret = readline(fd);
+	if (ret == nil) {
+		fprint(2, "error: early end of file\n");
+		exits("eof");
+	}
+	return ret;
+}
+
+Ann*
+annload(char *filename)
+{
+	Weights *W;
+	Ann *ann = nil;
+	char *buf;
+	int i, j, k;
+	int num_layers, *layers;
+	int fd = open(filename, OREAD);
+	if (fd < 0) {
+		perror("open");
+		return ann;
+	}
+
+	buf = sreadline(fd);
+	num_layers = atoi(buf);
+
+	if (num_layers < 2) {
+		fprint(2, "num_layers < 2\n");
+		return ann;
+	}
+
+	layers = calloc(num_layers, sizeof(int));
+	for (i = 0; i < num_layers; i++) {
+		buf = sreadline(fd);
+		layers[i] = atoi(buf);
+	}
+
+	ann = anncreatev(num_layers, layers);
+	for (i = 0; i < (ann->n - 1); i++) {
+		W = ann->weights[i];
+		for (j = 0; j <= W->inputs; j++)
+			for (k = 0; k < W->outputs; k++) {
+				buf = sreadline(fd);
+				W->values[j][k] = atof(buf);
+			}
+	}
+
+	free(retline);
+	retline = nil;
+
+	return ann;
+}
--- /dev/null
+++ b/ann/main.c
@@ -1,0 +1,166 @@
+#include <u.h>
+#include <libc.h>
+#include <ctype.h>
+#include "ann.h"
+
+void
+usage(char **argv)
+{
+	fprint(2, "usage: %s [-train] filename [num_layers num_input_layer ... num_output_layer]\n", argv[0]);
+	exits("usage");
+}
+
+void
+main(int argc, char **argv)
+{
+	Ann *ann;
+	char *filename;
+	int train;
+	Dir *dir;
+	int num_layers = 0;
+	int *layers = nil;
+	int i;
+	char *line;
+	double *input;
+	double *output = nil;
+	double *runoutput;
+	int ninput;
+	int noutput;
+	int offset;
+	double f;
+	int trainline;
+	int nline;
+
+	train = 0;
+
+	if (argc < 2)
+		usage(argv);
+
+	filename = argv[1];
+
+	if (argv[1][0] == '-' && argv[1][1] == 't') {
+		if (argc < 3)
+			usage(argv);
+
+		train = 1;
+		filename = argv[2];
+	}
+
+	ann = nil;
+	dir = dirstat(filename);
+	if (dir != nil) {
+		free(dir);
+		ann = annload(filename);
+		if (ann == nil)
+			exits("load");
+	}
+
+	if (argc >= (train + 3)) {
+		num_layers = atoi(argv[train + 2]);
+
+		if (num_layers < 2 || argc != (train + 3 + num_layers))
+			usage(argv);
+
+		layers = calloc(num_layers, sizeof(int));
+
+		for (i = 0; i < num_layers; i++)
+			layers[i] = atoi(argv[train + 3 + i]);
+	}
+
+	if (num_layers > 0) {
+		if (ann != nil) {
+			if (ann->n != num_layers) {
+				fprint(2, "num_layers: %d != %d\n", ann->n, num_layers);
+				exits("num_layers");
+			}
+
+			for (i = 0; i < num_layers; i++) {
+				if (layers[i] != ann->layers[i]->n) {
+					fprint(2, "num_layer_%d: %d != %d\n", i, layers[i], ann->layers[i]->n);
+					exits("num_layer");
+				}
+			}
+		} else {
+			ann = anncreatev(num_layers, layers);
+			if (ann == nil)
+				exits("anncreatev");
+		}
+	}
+
+	if (ann == nil) {
+		fprint(2, "file not found: %s\n", filename);
+		exits("file not found");
+	}
+
+	ninput = ann->layers[0]->n;
+	noutput = ann->layers[ann->n - 1]->n;
+	input = calloc(ninput, sizeof(double));
+	if (train == 1)
+		output = calloc(noutput, sizeof(double));
+
+	trainline = 0;
+	nline = ninput;
+
+	do {
+		int i = 0;
+		while ((line = readline(0)) != nil) {
+			do {
+				if (strlen(line) == 0)
+					break;
+				while(isspace(*line))
+					line++;
+				if (strlen(line) == 0)
+					break;
+				offset = 0;
+				while (isdigit(line[offset]) || line[offset] == '.' || line[offset] == '-')
+					offset++;
+				if (!isspace(line[offset]) && line[offset] != '\0') {
+					fprint(2, "input error: %s\n", line);
+					exits("input");
+				}
+				f = atof(line);
+				if (trainline == 0) {
+					input[i] = f;
+					i++;
+				} else {
+					output[i] = f;
+					i++;
+				}
+				line = &line[offset];
+				while(isspace(*line))
+					line++;
+			} while(i < nline && strlen(line) > 0);
+
+			if (i == nline) {
+				if (trainline == 0) {
+					runoutput = annrun(ann, input);
+					for (i = 0; i < noutput; i++)
+/*						if (runoutput[i] == 0.0)
+							print("0%c", (i == (noutput-1))? '\n': ' ');
+						else if (runoutput[i] == 1.0)
+							print("1%c", (i == (noutput-1))? '\n': ' ');
+						else */
+							print("%f%c", runoutput[i], (i == (noutput-1))? '\n': ' ');
+					free(runoutput);
+				}
+
+				if (train == 1) {
+					if (trainline == 0) {
+						trainline = 1;
+						nline = noutput;
+					} else {
+						anntrain(ann, input, output);
+						trainline = 0;
+						nline = ninput;
+					}
+				}
+				i = 0;
+			}
+		}
+	} while(line != nil);
+
+	if (annsave(filename, ann) != 0)
+		exits("save");
+
+	exits(nil);
+}
--- /dev/null
+++ b/ann/mkfile
@@ -1,0 +1,18 @@
+</$objtype/mkfile
+
+install:V: all
+	cp $O.ann $home/bin/$objtype/ann
+
+all:V: $O.ann
+
+$O.ann: main.$O ann.a$O
+	$LD -o $O.ann main.$O ann.a$O
+
+%.$O: %.c
+	$CC $CFLAGS $stem.c
+
+ann.a$O: ann.$O annio.$O
+	ar vu ann.a$O ann.$O annio.$O
+
+clean:V:
+	rm -f *.[$OS] [$OS].* *.a[$OS]