#define cimg_debug 2

#include "CImg.h"
#include <ctype.h>
#include <valarray>
#include <list>
#include <assert.h>
#include <math.h>
#define ACTG char actg[256]; memset(actg, 255, 256); \
	actg[(int)'a']=(char)0; actg[(int)'c']=(char)1; actg[(int)'t']=(char)2; actg[(int)'g']=(char)3; \
	actg[(int)'A']=(char)0; actg[(int)'C']=(char)1; actg[(int)'T']=(char)2; actg[(int)'G']=(char)3; 

using namespace cimg_library;
using namespace std;

class CGenome {
public:
        CGenome(char* filename, int len=7) {
		m_nLen = len;
                ACTG; char ch; unsigned int pos, fileLen, read, i;
                FILE* file = fopen(filename, "rb"); if (!file) { printf("Can't open %s\n", (filename)?filename:"<null>"); }
                fseek(file, 0, SEEK_END); fileLen = ftell(file); assert(fileLen);
                CImg<unsigned char> imgdata = CImg<unsigned char>(fileLen,1); fseek(file, 0, SEEK_SET);
                read = (unsigned int)fread(imgdata.data, 1, fileLen, file); assert(read == fileLen);
                
		pos = 0; i=0; 
		if (imgdata[0] == '>') { // should not be empty, checked in extract.pl
			while (imgdata[i++] != '\n') { if (i>fileLen) return; }
		}

		int width = 1<<m_nLen;
		m_src = CImg<>(width,width).fill(0);
		list<unsigned char> data;
                while (i<fileLen) { 
                        ch = imgdata[i++];
                        if(ch=='>') {
				fasta2asc(data); // process data here
				data.clear();
				while (imgdata[i++] != '\n') { if (i>fileLen) return; }
			}
			char none = (char)255;
                        if (actg[(unsigned char)ch]!=none) data.push_back(ch);
                }
		fasta2asc(data); // process data here
        }
	CImg<float> m_src;
private:
	unsigned int m_nLen;
	void fasta2asc (list<unsigned char>& data) {
		ACTG;	
		char none = (char)255;
		unsigned int i=0, X = 0, Y = 0, width = 1<<m_nLen;
		list<unsigned char>::iterator itr = data.begin();
		while (itr!=data.end()) {
                	char z = actg[*itr];
                	if (z!=none) { 
				X*=2; X+=z%2; X%=width; Y*=2; Y+=z/2; Y%=width;
				if (i>=m_nLen)
					m_src(X,Y)=m_src(X,Y)+1; 
			}
			i++; itr++;
        	}
	}
};

void gen2asc(char* filename, int len) {
	char buf[1024];

	CGenome gen(filename, len);
	*buf = 0; assert(strlen(filename)<1000); 
	sprintf(buf, "%s.%d.asc", filename, len); gen.m_src.save(buf);
	std::printf("%s\n", buf);
}

void usage(char* filename) {
	std::fprintf(stderr, "Usage: %s [-n] <fasta files>\n\tn=1..9\n\n", filename);
	exit(-1);
}

int main(int argc,char **argv) {
	int len = 9, from = 1;
	if (argc < 2)
		usage(argv[0]);

	if (argv[1][0] == '-') {
		len = atoi(argv[1]+1);
		if (len<1)
			usage(argv[0]);
		from = 2;
	}
	
	for (int i=from; i < argc; i++) {
		gen2asc(argv[i], len);
	}

	return 0;
}