/*
	Small neuron network
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <misc.h>
#include <sys/time.h>
#include <tlmplib.h>
#include <netdb.h>
#include <boost/coroutine2/all.hpp>
#include <string>
#include <vector>
#include "moron.h"

using namespace std;

static long long getnow()
{
	struct timeval tv;
	gettimeofday(&tv,NULL);
	return tv.tv_sec*1000000L+tv.tv_usec;
}

static double random01 ()
{
	double r = random();
	return r / RAND_MAX;
}
static  double sigmoid(double z)
{
	return 1.0/(1.0+expl(-z));
}

// Derivative of the sigmoid function.
static double  sigmoid_prime(double z)
{
    return sigmoid(z)*(1-sigmoid(z));
}

class NEURON{
	vector<double> weights;
	double bias;
public:
	NEURON (int nb_input){
		for (int i=0; i<nb_input; i++){
			weights.push_back(random01());
		}
		bias = random01();
	}
	double compute (const vector<double> &inputs){
		double ret = 0;
		if (weights.size() != inputs.size()){
			tlmp_error ("NEURON::compute, wrong input size\n");
		}else{
			for (unsigned i=0; i<inputs.size(); i++){
				ret += weights[i]*inputs[i];
			}
			ret += bias;
			ret = sigmoid(ret);
		}
		return ret;
	}
	unsigned getnbinputs() const {
		return weights.size();
	}
	void save (FILE *fout) const{
		fprintf (fout,"b %lf\n",bias);
		for (auto w:weights) fprintf (fout,"w %lf\n",w);
	}
	int load (FILE *fin){
		int ret = 0;
		char line[100];
		if (fgets(line,sizeof(line),fin)!=NULL && line[0] == 'b'){
			bias = atof(line+2);
			for (auto &w:weights){
				if (fgets(line,sizeof(line),fin)==NULL){
					tlmp_error ("NEURON::load, end of file while reading\n");
					ret = -1;
					break;
				}else if (line[0] != 'w'){
					tlmp_error ("NEURON::load, invalid line %s",line);
					ret = -1;
					break;
				}else{
					w = atof(line+2);
				}
			}
		}else{
			tlmp_error ("NEURON: Expect b line: %s",line);
			ret = -1;
		}
		return ret;
	}
	void dump(){
		printf ("weights:");
		for (auto w:weights) printf (" %lf",w);
		printf ("\n");
		printf ("bias   : %lf\n",bias);
	}
	double getbias() const{
		return bias;
	}
	const vector<double> &getweights() const{
		return weights;
	}
};

class LAYER{
	vector<NEURON> neurons;
public:
	LAYER(int nbneuron, int nbinputs){
		for (int i=0; i<nbneuron; i++){
			neurons.push_back(NEURON(nbinputs));
		}
	}
	vector<double> compute (const vector<double> &inputs){
		vector<double> ret;
		for (auto n:neurons){
			ret.push_back(n.compute(inputs));
		}
		return ret;
	}
	unsigned getnbneurons() const {
		return neurons.size();
	}
	unsigned getnbinputs() const {
		return neurons[0].getnbinputs();
	}
	void save (FILE *fout) const{
		for (auto n:neurons) n.save(fout);
	}
	int load (FILE *fin){
		int ret = 0;
		for (auto n:neurons){
			ret = n.load(fin);
			if (ret != 0) break;
		}
		return ret;
	}
	vector<double> getbiases() const{
		vector<double> ret;
		for (auto n:neurons) ret.push_back(n.getbias());
		return ret;
	}
	vector<double> getweights() const{
		vector<double> ret;
		for (auto n:neurons) for (auto w:n.getweights()) ret.push_back(w);
		return ret;
	}
		
};

class NET{
	vector<LAYER> layers;
public:
	NET (int nbinputs, const vector<int> &topo){
		for (auto n:topo){
			layers.push_back(LAYER(n,nbinputs));
			nbinputs = n;
		}
	}
	NET (int nbinputs, int argc, char *argv[]){
		for (int i=0; i<argc; i++){
			int n = atoi(argv[i]);
			layers.push_back(LAYER(n,nbinputs));
			nbinputs = n;
		}
	}
	NET (const char *fname){
		FILE *fin = fopen (fname,"r");
		if (fin == NULL){
			tlmp_error ("Can't open file %s (%s)\n",fname,strerror(errno));
		}else{
			char line[1000];
			if (fgets(line,sizeof(line),fin)==NULL){
				tlmp_error ("NET::NET empty file, expect netmoron as first line\n");
			}else if (strcmp(line,"netmoron\n")!=0){
				tlmp_error ("NET::NET invalid file, expect netmoron as first line\n");
			}else{
				int ret = 0;
				long pos = 0;	
				while (fgets(line,sizeof(line),fin)!=NULL){
					if (strncmp(line,"layer ",6)==0){
						vector<string> words;
						int n = str_splitline(line+6,' ',words);
						if (n != 2){
							tlmp_error ("NET::NET(fname): Invalid layer line %s",line);
							ret = -1;
							break;
						}else{
							int nbneuron = atoi(words[0].c_str());
							int inputs = atoi(words[1].c_str());
							//printf ("nbneuron=%d inputs=%d\n",nbneuron,inputs);
							layers.push_back (LAYER(nbneuron,inputs));
						}
						pos = ftell(fin);
					}else{
						break;
					}
				}
				if (ret == 0){
					fseek(fin,pos,SEEK_SET);
					for (auto l:layers){
						ret = l.load(fin);
						if (ret == -1) break;
					}
				}
				if (ret == -1) layers.clear();
			}
			fclose (fin);
		}
	}
	bool is_loaded() const {
		return layers.size() > 0;
	}
	vector<double> compute (vector<double> vals){
		for (auto l:layers){
			vals = l.compute(vals);
		}
		return vals;
	}
	int save(const char *fname) const {
		int ret = -1;
		FILE *fout = fopen (fname,"w");
		if (fout == NULL){
			tlmp_error ("Can't open file %s for writing (%s)\n",fname,strerror(errno));
		}else{
			fprintf (fout,"netmoron\n");
			for (auto l:layers){
				fprintf (fout,"layer %u %u\n",l.getnbneurons(),l.getnbinputs());
			}
			for (auto l:layers) l.save(fout);
			fclose (fout);
			ret = 0;
		}
		return ret;
	}
	vector<double> getbiases() const{
		vector<double> ret;
		for (auto l:layers) for (auto b:l.getbiases()) ret.push_back(b);
		return ret;
	}
	vector<double> getweights() const{
		vector<double> ret;
		for (auto l:layers) for (auto w:l.getweights()) ret.push_back(w);
		return ret;
	}
};

struct NIST_IMAGE{
	vector<unsigned char> pixels;
	size_t size() const{
		return pixels.size();
	}
	void clear(){
		pixels.clear();
	}
	void add (unsigned char pixel){
		pixels.push_back(pixel);
	}
	vector<double> getinputs() const {
		vector<double> ret;
		for (auto p:pixels) ret.push_back(p/256.0);
		return ret;
	}
};

static int moron_readnistimage(const char *image_file, unsigned &size, vector<NIST_IMAGE> &images)
{
	int ret = -1;
	string tmp;
	tmp = string("gunzip <") + image_file;
	FILE *fin = popen (tmp.c_str(),"r");
	if (fin == NULL){
		tlmp_error ("Can't open file %s\n",image_file);
	}else{
		struct{
			unsigned magic;
			unsigned nb;
			unsigned rows;
			unsigned cols;
		} header;
		if(fread(&header,1,16,fin)!=16){
			tlmp_error ("can't read nist image header (aborting)\n");
			exit (-1);
		}
		header.magic  = ntohl(header.magic);
		header.nb     = ntohl(header.nb);
		header.rows   = ntohl(header.rows);
		header.cols   = ntohl(header.cols);
		size = header.rows*header.cols;
		printf ("magic=%u nb=%u rows=%u cols=%u size=%u\n",header.magic,header.nb,header.rows,header.cols,size);
		int n;
		unsigned char buf[size];
		while((n=fread(buf,1,size,fin))==(int)size){
			NIST_IMAGE current;
			for (int i=0; i<n; i++){
				current.add(buf[i]);
			}
			images.push_back(current);
		}
		if (n != 0) tlmp_error ("nist image is partial n=%d < %u\n",n,size);
		pclose (fin);
		printf ("nb images = %lu\n",images.size());
		ret = (int)images.size();
	}
	return ret;
}
static int moron_readnistlabels(const char *label_file, vector<unsigned char> &labels)
{
	int ret = -1;
	string tmp;
	tmp = string("gunzip <") + label_file;
	FILE *fin = popen (tmp.c_str(),"r");
	if (fin == NULL){
		tlmp_error ("Can't open file %s\n",label_file);
	}else{
		struct{
			unsigned magic;
			unsigned nb;
		} header;
		if(fread(&header,1,8,fin)!=8){
			tlmp_error ("can't read nist label header (aborting)\n");
			exit (-1);
		}
		header.magic = ntohl(header.magic);
		header.nb    = ntohl(header.nb);
		printf ("magic=%u nb=%u\n",header.magic,header.nb);
		int n;
		unsigned char buf[1000];
		while((n=fread(buf,1,1000,fin))>0){
			for (int i=0; i<n; i++){
				labels.push_back(buf[i]);
			}
		}
		pclose (fin);
		printf ("nb labels = %lu\n",labels.size());
		ret = (int)labels.size();
	}
	return ret;
}
#if 0
static vector<double> np_zeros(const vector<double> &v)
{
	vector<double> ret;
	for (unsigne i=0; i<v.size(); i++) ret.push_back(0.0);
	return ret;
}

static void cost_derivative(output_activations, y)
{
        // Return the vector of partial derivatives \partial C_x
        // partial a for the output activations."""
        return (output_activations-y)
}


/*
        Return a tuple ``(nabla_b, nabla_w)`` representing the
        gradient for the cost function C_x.  ``nabla_b`` and
        ``nabla_w`` are layer-by-layer lists of numpy arrays, similar
        to ``self.biases`` and ``self.weights``.
*/
static void backprop(
	const vector <double> &biases,
	const vector <double> &weights,
	const vector <double> &answer,
	const vector <double> &result,
	vector<double> &nabla_b,
	vector<double> &nabla_w)
{
	nabla_b = np_zeros(biases);
	nabla_w = np_zeros(weights);
        # feedforward
        auto activation = answer;
	vector< typeof activation> activations = {answer}; // [x] # list to store all the activations, layer by layer
        vector<vector<double> > zs;			//  = [] # list to store all the z vectors, layer by layer
        for (auto i:zip(biases,weights)){
            auto z = dot(i.second,activation)+i.first;
            zs.push_back(z)
            activation = sigmoid(z)
            activations.push_back(activation)
	}
        // backward pass
        auto delta = self.cost_derivative(activations[-1], y)
            * sigmoid_prime(zs[-1]);
        nabla_b[-1] = delta
        nabla_w[-1] = np.dot(delta, activations[-2].transpose())
	/*
	# Note that the variable l in the loop below is used a little
	# differently to the notation in Chapter 2 of the book.  Here,
	# l = 1 means the last layer of neurons, l = 2 is the
	# second-last layer, and so on.  It's a renumbering of the
	# scheme in the book, used here to take advantage of the fact
	# that Python can use negative indices in lists.
	*/
        for l in xrange(2, self.num_layers):
            z = zs[-l]
            sp = sigmoid_prime(z)
            delta = np.dot(self.weights[-l+1].transpose(), delta) * sp
            nabla_b[-l] = delta
            nabla_w[-l] = np.dot(delta, activations[-l-1].transpose())
        return (nabla_b, nabla_w)
}

static void update_mini_batch(
	NET &net,
	const vector<vector<double> > &answers,
	const vector<vector<double> > &results,
	int eta)
{
	auto biases  = net.getbiases();
	auto weights = net.getweights();
	auto nabla_b = np_zeros(biases);
	auto nabla_w = np_zeros(weights);
	for (unsigned i=0; i<answer.size(); i++){
		auto &answer = answers[i];
		auto &result = results[i];
		//delta_nabla_b, delta_nabla_w = self.backprop(x, y)
		backprop(answer,result);
		for (unsigned j=0; j<nabla_b.size(); j++) nabla_b[j] += delta_nabla_b[j];
		for (unsigned j=0; j<nabla_w.size(); j++) nabla_w[j] += delta_nabla_w[j];
	}
	self.weights = [w-(eta/len(mini_batch))*nw
                        for w, nw in zip(self.weights, nabla_w)]
        self.biases = [b-(eta/len(mini_batch))*nb
                       for b, nb in zip(self.biases, nabla_b)]
}
#endif

<mod>
int main (int argc, char *argv[])
{
	glocal int ret = -1;
	glocal int nbinputs = 10;
	glocal const char *nist_images = NULL;
	glocal const char *nist_labels = NULL;
	glocal bool test = false;
	glocal vector<int> topo;
	glocal const char *netfile = NULL;
	glocal const char *textfile = NULL;
	glocal.ret = <call tlmpprogram>(argc,argv);
	<f init>
		setproginfo ("moron","0.0"
			,"Minimalist neuron network to learn (me learing, not the network)\n"
			 "\n"
			 "moron n1 n2 n3...\n"
			 "n1 = number of neuron in first layer\n"
			 "n2 = number of neuron in second\n"
			);
		setarg ('i',"nbinputs","Number of input data",glocal.nbinputs,false);
		setarg ('I',"nist_images","nist images file",glocal.nist_images,false);
		setarg ('T',"nist_labels","nist labels file",glocal.nist_labels,false);
		setarg (' ',"topo","Network topology",glocal.topo,false);
		setarg (' ',"netfile","Parameter from a previous learn session",glocal.netfile,false);
		setarg (' ',"test","Various test",glocal.test,false);
		setarg (' ',"testfile","Read a text file",glocal.textfile,false);
	</f>
	<f main_noarg>
		int ret = -1;
		if (glocal.nist_images != NULL){
			vector<NIST_IMAGE> images;
			unsigned size=0;
			moron_readnistimage(glocal.nist_images,size,images);
			vector<unsigned char> labels;
			if (glocal.nist_labels != NULL){
				moron_readnistlabels(glocal.nist_labels,labels);
			}
			if (glocal.netfile != NULL){
				NET net(glocal.netfile);
				if (net.is_loaded()){
					int nbmatch = 0;
					for (auto img:zip(images,labels)){
						vector<double> vals = img.first.getinputs();
						vals = net.compute(vals);
						auto ret = maxval(vals);
						if (ret.first == img.second) nbmatch++;
					}
					printf ("nbmatch=%d\n",nbmatch);
				}
			}else{
				NET net(size,glocal.topo);
				int nbmatch = 0;
				for (auto img:zip(images,labels)){
					vector<double> vals = img.first.getinputs();
					vals = net.compute(vals);
					auto ret = maxval(vals);
					if (ret.first == img.second) nbmatch++;
				}
				printf ("nbmatch=%d\n",nbmatch);
				net.save ("/tmp/n.net");
			}
			ret = 0;		
		}else if (glocal.test){
			const int nbdouble=64;
			typedef double v4si __attribute__ ((vector_size (nbdouble*sizeof(double))));
			v4si v1,v2,res;
			for (int i=0; i<nbdouble; i++){
				v1[i] = 100+i;
				v2[i] = 100+i;
			}
			res = v1 * v2 + 1;
			for (int i=0; i<nbdouble; i++){
				printf ("%lf %lf %lf %lf\n",v1[i],v2[i],res[i],v1[i]*v2[i]);
			}
			const int loop=10000000;
			long long start = getnow();
			double total1 = 0;
			for (int i=0; i<loop; i++){
				res = v1*v2 + i;
				total1 += res[0];
			}
			long long end1 = getnow();
			double total2 = 0;
			for (int i=0; i<loop; i++){
				for (int j=0; j<nbdouble; j++){
					res[j] = v1[j]*v2[j] + i;
				}
				total2 += res[0];
			}
			long long end2 = getnow();
			printf ("%Ld %Ld  %lf %lf\n",end1-start,end2-end1,total1,total2);
			ret = 0;
		}else if (glocal.textfile){
			FILE *fin = fopen (glocal.textfile,"r");
			if (fin != NULL){
				int n = 0;
				for (auto x:readlines(fin)){
					n++;
					if (n > 10) break;
					printf ("     %s",x.c_str());
				}
				fclose(fin);
			}
		}else{
			NET net (glocal.nbinputs,glocal.topo);
			vector<double> inputs;
			for (int i=0; i<glocal.nbinputs; i++) inputs.push_back(i);
			vector<double> vals = inputs;
			vals = net.compute(vals);
			printf ("vals: ");
			for (auto v:vals) printf (" %lf",v);
			printf ("\n");
			ret = 0;
		}
		return ret;
	</f>
	</call>
	return glocal.ret;
}
</mod>