public class ClientGlobals {
	
	public static boolean useVmatch=true;
	public static String vmatchCommand;
	
	public static final String VERSION = "2.0";
	public static final String AUTHOR_EMAIL = "contact.cube@univie.ac.at";
	public static final String WEB_ADRESS = "http://cube.univie.ac.at/gepard";
	
	public static String SETTINGS_DIR = System.getProperty("user.home") + 
		System.getProperty("file.separator") + ".gepard/"; 
	
	public static final int CONTROLPANEL_WIDTH = 260;
	public static final int MINWORDLEN = 1;
	public static final int MAXWORDLEN = 32767;
	public static final int MINWINDOWSIZE = 0;
	public static final int MAXWINDOWSIZE = 10000;
	public static final int MINZOOM = 0;
	public static final int MAXZOOM = 1000000;
	public static final int MAXNAMELEN_TITLE = 50;
	
	
	
	public static final int MINVMATCHLENGTH = 50000;

	public static final float ZOOMINFACTOR =0.5f;
	public static final float ZOOMOUTFACTOR =2f;

	public static final int BIGKEYBOARDSTEP = 25;
	
	public static String APPNAME = "Gepard";
	
	public static final String PATH_MATRICES = "/resources/matrices/";
	public static final String FILE_MATRICES = "matrices.xml";
	public static final String AUTO_NUCL_MATRIX = "edna.mat";
	public static final String AUTO_NUCL_MATRIX_NAME = "DNA";
	public static final String AUTO_PROT_MATRIX = "blosum62.mat";
	public static final String AUTO_PROT_MATRIX_NAME = "BLOSUM62";

	public static final String FILE_HELP = "/resources/help/help.xml";
	public static final String FILE_USERMSG = "/resources/help/usermsg.txt";
	public static final String FILE_REFERENCE = "/resources/help/reference.txt";

	public static final String CONFIG_FILE = "config.xml";
	
	public static final int LOCALCLICK_SHOWALIGN = 0;
	public static final int LOCALCLICK_EXPORT = 1;
	public static final int LOCALCLICK_NOTHING = 2;
	
	public static final String NOGENE_STRING = "-";
	
	public static final String RES_APPICON = "/resources/images/gepard.gif";
	
	public static final String GENE_NAME_HOR_PREFIX  = "Hor : ";
	public static final String GENE_NAME_VERT_PREFIX = "Vert: ";
	

	public static String cutString(String str, int maxlen) {
		
		if (maxlen <= 3)
			return "";
		else if (str.length() > maxlen)
			return str.substring(0, maxlen - 3) + "...";
		else
			return str;
	}

	public static SubstitutionMatrix getEDNA() throws IOException, InvalidSubMatFileException {
//		return SubstitutionMatrix.debugTrivialMatrix();
		return SubstitutionMatrix.loadFromResource(ClientGlobals.PATH_MATRICES + "edna.mat");
	}

	public static void errMessage(String msg) {
		JOptionPane.showMessageDialog( 
				null, msg, "Error", JOptionPane.ERROR_MESSAGE );
	}
	
	public static void warnMessage(String msg) {
		JOptionPane.showMessageDialog( 
				null, msg, "Warning", JOptionPane.WARNING_MESSAGE );
	}
	
	public static void infoMessage(String msg) {
		JOptionPane.showMessageDialog( 
				null, msg, "Information", JOptionPane.INFORMATION_MESSAGE );
	}
	
	public static void copyTextToClipboard(String text) {
		Clipboard systemClipboard =
			Toolkit.getDefaultToolkit().getSystemClipboard();
		// create transferrable text object and set clipboard text
		Transferable transferableText =	new StringSelection(text);
		systemClipboard.setContents(transferableText,null);
	}
	
	
	private static NumberFormat FMT_FILESIZE = new DecimalFormat("###0.000", new DecimalFormatSymbols(Locale.ENGLISH));
	
	public static String convenientFileSize(long fileSize) {
		if (fileSize >= 1073741824) 
			return FMT_FILESIZE.format((double)fileSize / (double)1073741824) + " gb";
		else if (fileSize >= 1048576) 
			return FMT_FILESIZE.format((double)fileSize / (double)1048576) + " mb";
		else if (fileSize >= 1024) 
			return FMT_FILESIZE.format((double)fileSize / (double)1024) + " kb";
		else
			return fileSize + " b";
		
	}
	
	private static NumberFormat FMT_NUMBERS = new DecimalFormat("###,###,###,###", new DecimalFormatSymbols(Locale.ENGLISH));

	public static String formatNumber(int number) {
		return FMT_NUMBERS.format(number);
	}
	
	
	public static void unexpectedError(Exception e, Controller ctrl) {
		
		// handles any unexpected errors
		// asks the user to send an error report to the developer
		
		String guidump = ctrl.getGUIDump();
		String msg;
		msg = "An unexpected error occured!\nStack trace & GUI details copied to clipboard.";
		JOptionPane.showMessageDialog(null, msg, "Error", JOptionPane.ERROR_MESSAGE);
		
		String report = stack2string(e) + "\n\n\n\n" + guidump;
		StringSelection stringSelection = new StringSelection (report);
		Clipboard clpbrd = Toolkit.getDefaultToolkit ().getSystemClipboard ();
		clpbrd.setContents (stringSelection, null);
	
	}
	
	
	public static String stack2string(Exception e) {
		try {
			StringWriter sw = new StringWriter();
			PrintWriter pw = new PrintWriter(sw);
			e.printStackTrace(pw);
			return "------\r\n" + sw.toString() + "------\r\n";
		} catch (Exception e2) {
			return "bad stack2string";
		}
	}
	
	
	public static String extractFilename(String filename) {
		return filename.substring(filename.lastIndexOf(File.separator)+1);
	}
	
	public static String extractDirectory(String filename) {
		return filename.substring(0,filename.lastIndexOf(File.separator)+1);
	}
	
	public static String createFilenameFromName(String name, int maxlen) {
		
		name = name.replaceAll(" ", "_");
		if (name.length() > maxlen)
			name = name.substring(0, maxlen);
		
		return name + ".fa";
	}
	
	public static String insertBeforeExtension(String filename, String insert) {
		int lastDot = filename.lastIndexOf('.');
		int lastSep = filename.lastIndexOf(File.separator);

		if (lastDot > lastSep && lastDot >= 0) 
			return filename.substring(0, lastDot) + insert+ filename.substring(lastDot);
		else
			return filename + insert;
	}
	
} public interface DrawTarget {

	void setPreferredSize(Dimension dimension);

	int getWidth();

	int getHeight();
	

} public class Plotter {

	// member variables
	// int dpwidth, dpheight; // dotplot dimensions
	DotMatrix dm;

	// selection coordinates
	int selX1 = -1, selY1 = -1, selX2 = -1, selY2 = -1;

	public static final int XOFF = 85;
	public static final int YOFF = 110;
	private static final int GC_RATIO_OFF = 160;
	private static final Font INFO_FONT = new Font("Monospaced", Font.PLAIN, 12);
	private static final Font GENE_FONT = new Font("Monospaced", Font.PLAIN, 11);

	private static final int NAME_MIN_HOR_MARGIN = 10;
	private static final int NAME_VERT_MARGIN = 30;

	// crosshair
	private int crossX = -1, crossY = -1;
	private int crossPosX, crossPosY;

	// gene tooltip
	private int genex = -1, geney = -1;

	private int minPlotWidth = 0;

	private BufferedImage offscreen; // " "

	private String infoPlotTitle;
	private String infoZoom;
	private String infoWordLen;
	private String infoWindowSize;
	private String infoSubmat;
	private String infoProgram;
	private String infoGCRatio1;
	private String infoGCRatio2;

	private String seq1Name, seq2Name;

	private int imgWidth, imgHeight;

	private DrawTarget dp;

	private String genehor;
	private String genevert;

	// private Color greyscale[]; // greyscale helper array

	public Plotter(DotMatrix dm, DrawTarget dp) {

		this.dm = dm;
		this.dp = dp;

		// create backbuffer
		offscreen = new BufferedImage(dm.getWidth(), dm.getHeight(), BufferedImage.TYPE_INT_RGB);

		// create information strings (except plot title)
		infoZoom = "Zoom: " + dm.getParameterSet().zoom + " : 1";
		infoWordLen = "Word length: " + dm.getParameterSet().wordLength;
		infoWindowSize = "Window size: " + dm.getParameterSet().windowSize;
		infoSubmat = "Matrix: " + dm.getParameterSet().submatName;
		infoProgram = "Program: " + ClientGlobals.APPNAME + " (" + ClientGlobals.VERSION + ")";

		NumberFormat fmt = new DecimalFormat("#####0.0000", new DecimalFormatSymbols(Locale.ENGLISH));
		infoGCRatio1 = "GC ratio seq1: " + fmt.format(dm.getGCratio1());
		infoGCRatio2 = "GC ratio seq2: " + fmt.format(dm.getGCratio2());

		// calculate the minimum width of the dp, need for scrollpane & picture
		// export
		// minimum width is determined using text width of seq vs. seq string,
		// the
		// GC ratio string and the program name string
		FontMetrics fontMetrics = offscreen.getGraphics().getFontMetrics(INFO_FONT);
		int infoGCRatioWidth = XOFF + GC_RATIO_OFF + fontMetrics.stringWidth(infoGCRatio1);
		int infoProgramWidth = XOFF + GC_RATIO_OFF + fontMetrics.stringWidth(infoProgram);

		// store sequence names
		int hormargin = Math.max(Math.max(NAME_MIN_HOR_MARGIN, fontMetrics.stringWidth(dm.getParameterSet().seq1Start + "") * 2),
				fontMetrics.stringWidth(dm.getParameterSet().seq1Stop + "") * 2);

		seq1Name = ClientGlobals.cutString(dm.getSeq1Name(), (int) ((dm.getWidth() - hormargin) / fontMetrics.charWidth(' ')));
		seq2Name = ClientGlobals.cutString(dm.getSeq2Name(), (int) ((dm.getHeight() - NAME_VERT_MARGIN) / fontMetrics.charWidth(' ')));

		// calc minimum (=> maximum of the two widths above
		minPlotWidth = Math.max(infoGCRatioWidth, infoProgramWidth);

		// calc image dimensions
		imgWidth = (minPlotWidth > XOFF + dm.getWidth()) ? minPlotWidth + 5 : XOFF + dm.getWidth() + 5;
		imgHeight = YOFF + dm.getHeight() + 5;

		// now we need to create the plot title with corretly cut strings
		String name1 = dm.getSeq1Name();
		String name2 = dm.getSeq2Name();
		infoPlotTitle = "";// seq1Name + " vs. " + seq2Name;
		// step 1: calc max amount of chars for plot title
		int charsfortitle = ((imgWidth - XOFF) / fontMetrics.charWidth(' ')) - 5;
		int halfchars = charsfortitle / 2;
		// step 2 calc max and min seq name length
		int maxnamelen = Math.max(name1.length(), name2.length());
		int minnamelen = Math.min(name1.length(), name2.length());
		// step 3: distribute characters
		// if both names are shorter than half of the available characters:
		// don't cut
		if (maxnamelen <= halfchars)
			infoPlotTitle = name1 + " vs. " + name2;
		// if both names are longer than the half... cut both!
		else if (minnamelen > halfchars)
			infoPlotTitle = ClientGlobals.cutString(name1, halfchars) + " vs. " + ClientGlobals.cutString(name2, halfchars);
		// if one seq is shorter and one is longer than the half: completely
		// print shorter sequence
		// and leave rest to the longer one
		else {
			if (name1.length() < name2.length())
				// name1 is the short one
				infoPlotTitle = name1 + " vs. " + ClientGlobals.cutString(name2, charsfortitle - name1.length());
			else
				// name2 is the short one
				infoPlotTitle = ClientGlobals.cutString(name1, charsfortitle - name2.length()) + " vs. " + name2;

		}

		dp.setPreferredSize(new Dimension(imgWidth, imgHeight));

	}

	// sets the current crosshair position
	// recalculates all dotplot data
	public void reCalc(float colorlower, float colorupper, float greyscalestart) {

		// store dot matrix data
		float[][] matrix = dm.getDotMatrix();
		

		int colorDot = 0;
		int r = 0, g = 0, b = 0; // final color values

		for (int i = 0; i < dm.getWidth(); i++) {
			for (int j = 0; j < dm.getHeight(); j++) {

				if (matrix[i][j] >= colorlower) {
					float dotscore = matrix[i][j];
					if (dotscore > colorupper)
						dotscore = colorupper;

					float ratio = greyscalestart + ((float) (dotscore - colorlower) / (float) (colorupper - colorlower)) * (1f - greyscalestart);
					// calculate grey value
					colorDot = (int) ((1f - ratio) * 255f);
				} else
					colorDot = 255;
				r = colorDot;
				g = colorDot;
				b = colorDot;

				offscreen.setRGB(i, j, r * 0x10000 + g * 0x100 + b); // is this
																		// still
																		// needed?
			}
		}
	}

	public void setSelection(int x1, int y1, int x2, int y2) {
		// correct values if outside of dotplot
		if (x1 < XOFF)
			x1 = XOFF;
		if (x2 < XOFF)
			x2 = XOFF;
		if (y1 < YOFF)
			y1 = YOFF;
		if (y2 < YOFF)
			y2 = YOFF;

		if (x1 >= XOFF + dm.getWidth())
			x1 = XOFF + dm.getWidth() - 1;
		if (x2 >= XOFF + dm.getWidth())
			x2 = XOFF + dm.getWidth() - 1;
		if (y1 >= YOFF + dm.getHeight())
			y1 = YOFF + dm.getHeight() - 1;
		if (y2 >= YOFF + dm.getHeight())
			y2 = YOFF + dm.getHeight() - 1;

		selX1 = x1;
		selY1 = y1;
		selX2 = x2;
		selY2 = y2;
	}

	public void setCrossHair(int x, int y, int posx, int posy) {
		// correct values if outside plot
		if (x < XOFF)
			x = XOFF;
		if (y < YOFF)
			y = YOFF;
		if (x >= XOFF + dm.getWidth())
			x = XOFF + dm.getWidth() - 1;
		if (y >= YOFF + dm.getHeight())
			y = YOFF + dm.getHeight() - 1;
		// store values
		crossX = x;
		crossY = y;
		crossPosX = posx;
		crossPosY = posy;
	}

	public void moveCrossHair(int xmove, int ymove) {
		crossPosX += xmove;
		crossPosY += ymove;
		// correct values if necessary
		/*
		 * if (crossPosX<0) crossPosX = 0; if (crossPosX>=dm.getSeq1Length())
		 * crossPosX = dm.getSeq1Length(); if (crossPosY<0) crossPosY = 0; if
		 * (crossPosY>=dm.getSeq2Length()) crossPosY = dm.getSeq2Length();
		 */

	}

	public void removeCrossHair() {
		crossX = -1;
		crossY = -1;
	}

	public void removeSelection() {
		selX1 = -1;
		selY1 = -1;
	}

	public void redraw(Graphics2D g) {

		// whiten (image dimensions but at least draw panel dimensions
		int whitenWidth = Math.max(imgWidth, dp.getWidth());
		int whitenHeight = Math.max(imgHeight, dp.getHeight());
		g.setColor(Color.WHITE);
		g.fillRect(0, 0, whitenWidth, whitenHeight);
		// draw frame
		g.setColor(Color.BLACK);
		g.drawRect(XOFF - 1, YOFF - 1, dm.getWidth() + 1, dm.getHeight() + 1);

		// draw plot information
		g.setFont(INFO_FONT);

		FontMetrics fontMetrics = g.getFontMetrics(g.getFont());

		// draw first sequence name
		g.drawString(seq1Name, XOFF + (dm.getWidth() / 2) - (fontMetrics.stringWidth(seq1Name) / 2), YOFF - 7);
		// draw second sequence name, rotated
		Font oldFont = g.getFont();
		Font rotated = getRotatedFont(oldFont);
		g.setFont(rotated);
		g.drawString(seq2Name, XOFF - 7, YOFF + (dm.getHeight() / 2) + (fontMetrics.stringWidth(seq2Name) / 2));
		g.setFont(oldFont);

		g.drawString(infoPlotTitle, XOFF, 13);
		g.drawString(infoZoom, XOFF, 28);
		g.drawString(infoWordLen, XOFF, 43);
		g.drawString(infoWindowSize, XOFF, 58);
		g.drawString(infoSubmat, XOFF, 73);
		// draw GC ratio (if nucleotide dotplot)
		if (this.dm.isNucleotideMatrix()) {
			g.drawString(infoGCRatio1, XOFF + GC_RATIO_OFF, 43);
			g.drawString(infoGCRatio2, XOFF + GC_RATIO_OFF, 58);
		}
		g.drawString(infoProgram, XOFF + GC_RATIO_OFF, 73);

		// draw start/end coordinates + lines

		g.drawLine(XOFF, YOFF - 10, XOFF, YOFF);
		g.drawString(dm.getParameterSet().seq1Start + "", XOFF - 3, YOFF - 17);
		g.drawLine(XOFF + dm.getWidth() - 1, YOFF - 10, XOFF + dm.getWidth() - 1, YOFF);
		g.drawString(dm.getParameterSet().seq1Stop + "", XOFF + dm.getWidth() - fontMetrics.stringWidth(dm.getParameterSet().seq1Stop + "") + 3, YOFF - 17);
		g.drawLine(XOFF - 10, YOFF, XOFF, YOFF);
		g.drawString(dm.getParameterSet().seq2Start + "", XOFF - 17 - fontMetrics.stringWidth(dm.getParameterSet().seq2Start + ""), YOFF + 4);
		g.drawLine(XOFF - 10, YOFF + dm.getHeight() - 1, XOFF, YOFF + dm.getHeight() - 1);
		g.drawString(dm.getParameterSet().seq2Stop + "", XOFF - 17 - fontMetrics.stringWidth(dm.getParameterSet().seq2Stop + ""), YOFF + dm.getHeight() + 2);
		// draw plot
		g.drawImage(offscreen, XOFF, YOFF, null);
		// draw selection
		if (selX1 != -1) {
			int x = 0, y = 0, width = 0, height = 0;
			if (selX1 < selX2) {
				x = selX1;
				width = selX2 - selX1;
			} else {
				x = selX2;
				width = selX1 - selX2;
			}
			if (selY1 < selY2) {
				y = selY1;
				height = selY2 - selY1;
			} else {
				y = selY2;
				height = selY1 - selY2;
			}

			g.setColor(Color.BLUE);
			g.drawRect(x, y, width, height);
		}

		if (crossX > -1) {
			// draw crosshair
			g.setColor(Color.BLUE);
			g.drawLine(crossX, YOFF, crossX, crossY - 1);
			g.drawLine(crossX, crossY + 1, crossX, dm.getHeight() + YOFF);
			g.drawLine(XOFF, crossY, crossX - 1, crossY);
			g.drawLine(crossX + 1, crossY, dm.getWidth() + XOFF, crossY);
			// draw position
			g.drawString(crossPosX + ", " + crossPosY, crossX + 5, crossY - 5);
		}

		// draw gene tool tip
		if (genex > -1) {

			// get gene name widths
			FontMetrics genemetrics = offscreen.getGraphics().getFontMetrics(GENE_FONT);
			int texthorwidth = genemetrics.stringWidth(genehor);
			int textvertwidth = genemetrics.stringWidth(genevert);
			int textwidth = (texthorwidth > textvertwidth) ? texthorwidth : textvertwidth;
			// and height
			int textheight = genemetrics.getHeight();
			// and the strange vertical shifting
			int yshift = genemetrics.getAscent();

			final int GENE_YOFF = 20;
			final int GENE_XOFF = 0;

			// draw gene name box
			g.setFont(GENE_FONT);
			g.setColor(Color.BLACK);
			g.drawRect(genex + GENE_XOFF, geney + GENE_YOFF, textwidth + 2, textheight * 2 + 2);
			g.setColor(Color.WHITE);
			g.fillRect(genex + GENE_XOFF + 1, geney + GENE_YOFF + 1, textwidth + 2 - 1, textheight * 2 + 2 - 1);

			// draw gene names
			g.setColor(Color.BLACK);
			g.drawString(genehor, genex + GENE_XOFF + 1, geney + GENE_YOFF + yshift + 1);
			g.drawString(genevert, genex + GENE_XOFF + 1, geney + GENE_YOFF + textheight + yshift + 1);

		}

	}

	private Font getRotatedFont(Font original) {
		AffineTransform fontAT = new AffineTransform();
		fontAT.rotate(Math.PI / 2 * 3);
		return original.deriveFont(fontAT);
	}

	public BufferedImage getFullImage() {

		// create new image which is large enough to hold the total plot
		BufferedImage fullImage = new BufferedImage(imgWidth, imgHeight, BufferedImage.TYPE_INT_RGB);
		// create graphics context
		Graphics2D g = (Graphics2D) fullImage.getGraphics();
		// draw plot into this image & return
		redraw(g);

		return fullImage;
	}

	public void setGeneToolTip(int x, int y, String genehor, String genevert) {
		genex = x;
		geney = y;
		this.genehor = ClientGlobals.GENE_NAME_HOR_PREFIX + genehor;
		this.genevert = ClientGlobals.GENE_NAME_VERT_PREFIX + genevert;
	}

	public String getGeneNameInformation() {
		return this.genehor + "\n" + this.genevert;
	}

	public void noToolTips() {
		genex = -1;
	}

} public interface AbortionChecker {
	public boolean dotplotAborted();
} public class DotMatrix implements Serializable {

	private static final long serialVersionUID = 5086505082934770263L;

	// members
	private float[][] dotmatrix = null;

	private int width = 0, height = 0;

	private float maxdotscore, mindotscore;
	private float avgdotscore;
	
	private float GCratio1=0, GCratio2;

	private int seq1Length, seq2Length;

	private ParameterSet params;

	private String seq1Name, seq2Name;

	private SuffixArray sa;

	private SubstitutionMatrix submat;
	
	private boolean bFuncats=false;
	
	private boolean nuclmatrix;

	// simple accessor methods
	public float[][] getDotMatrix() {
		return dotmatrix;
	}

	public boolean isNucleotideMatrix() {
		return nuclmatrix;
	}


	public boolean usingFuncats() {
		return bFuncats;
	}
	
	public float getGCratio1() {
		return GCratio1;
	}
	
	public float getGCratio2() {
		return GCratio2;
	}
	
	public SuffixArray getSuffixArray() {
		return sa;
	}

	public String getSeq1Name() {
		return seq1Name;
	}

	public String getSeq2Name() {
		return seq2Name;
	}

	// public float[][]getNoiseMatrix(){return noisematrix;}

	public ParameterSet getParameterSet() {
		return params;
	}

	// public int getRatio(){return ratio;}
	public int getSeq1Length() {
		return seq1Length;
	}

	public int getSeq2Length() {
		return seq2Length;
	}

	// public int getWordLength(){return wordlen;}
	// public int getWindowSize(){return windowsize;}
	public int getWidth() {
		return width;
	}

	public int getHeight() {
		return height;
	}

	public float getMaxDotScore() {
		return maxdotscore;
	}

	public float getMinDotScore() {
		return mindotscore;
	}
	
	public float getAvgDotScore() {
		return avgdotscore;
	}

	// additional parameters method
	public DotMatrix(byte[] s1, byte[] s2, String seq1n, String seq2n,
			SuffixArray sa, ParameterSet params, SubstitutionMatrix submat, 
			DotMatrixCallback callback, AbortionChecker ac, boolean SAforSecondSeq ) {
		
		
		// store sequences suffix array, word length, ratio & char distribution
		// flag
		// seq1 = s1;
		// seq2 = s2;
		seq1Length = s1.length;
		seq2Length = s2.length;
		seq1Name = seq1n;
		seq2Name = seq2n;
		this.sa = sa;
		this.submat = submat;
		// wordlen = params.wordLength;
		// ratio = params.ratio;
		// windowsize = params.windowSize;
		this.params = params;
		//int x = params.windowSize;
		// if coordinates of one sequence are both 0 -> full sequence
		if ((params.seq1Start == 0) && (params.seq1Stop == 0))
			params.seq1Stop = s1.length - 1;
		if ((params.seq2Start == 0) && (params.seq2Stop == 0))
			params.seq2Stop = s2.length - 1;
		// calculate dimensions
		width = ((params.seq1Stop - params.seq1Start + 1) / params.zoom);
		height = ((params.seq2Stop - params.seq2Start + 1) / params.zoom);

		nuclmatrix = submat.isNucleotideMatrix();
		// compute GC ratio
		if (nuclmatrix)
			calcGCRatios(s1,s2,params);

		// calculate dot matrix
		if (params.wordLength != 0)
			calcDotMatrix(sa, s1, s2, callback, ac, SAforSecondSeq );
		else 
			calcDotMatrixWindowMode(sa, s1, s2, callback, ac);
		
	}
	
	// matrix calculation
	private void calcDotMatrix(SuffixArray sa, byte[] seq1, byte[] seq2, DotMatrixCallback callback, AbortionChecker ac, boolean SAforSecondSeq ) {

		// determine callback step
		int callbackStep = (callback != null) ? callback.tellCallbackStep(params.wordLength, params.windowSize ) : 2147483647;
		// create dot matrix array
		dotmatrix = new float[width][height];
		
		// set correct sequence order
		byte[] nonSAsequence;
		int from, to;
		int SAseqStart, SAseqStop;
		int nonSAseqStart, nonSAseqStop;
		if (!SAforSecondSeq) {	// SA for first sequence
			nonSAsequence = seq2;
			from = params.seq2Start;
			to = params.seq2Stop - params.wordLength + 1;
			SAseqStart = params.seq1Start;
			SAseqStop  = params.seq1Stop;
			nonSAseqStart = params.seq2Start;
			nonSAseqStop  = params.seq2Stop;
		} else {				// SA for second sequence
			nonSAsequence = seq1;
			from = params.seq1Start;
			to = params.seq1Stop - params.wordLength + 1;
			SAseqStart = params.seq2Start;
			SAseqStop  = params.seq2Stop;
			nonSAseqStart = params.seq1Start;
			nonSAseqStop  = params.seq1Stop;
		}
		
		int[] hits = null;
		
		// set pseudodelimiter (for pseudocontigs)   ... don't dare to remove this anymore, old relic from remote functionality
		byte pseudodelimiter;
		boolean nuclmatrix = submat.isNucleotideMatrix();
		if (nuclmatrix)
			pseudodelimiter = submat.map( (byte)'Z');
		else
			pseudodelimiter=-99;
		
		// create identity and complementary maps
		byte[] idmap = new byte[128];
		for (byte i=0; i<128 && i >= 0; i++) idmap[i] = i;
		// use complementary map if no protein sequences are used
		byte[] complmap ;
		if (submat.isNucleotideMatrix())
			complmap = Sequence.getComplementaryMap();
		else 
			complmap = idmap;
		
		// iterate through non-SA sequence
		for (int i = from; i < to; i++) {
			// only do any searching if first letter != the pseudodelimiter Z
			// also ignore N's and X's in nucleotide matrices
			if ( (nonSAsequence[i] != pseudodelimiter )&&
					!(  nuclmatrix && (     nonSAsequence[i] ==  5         ) ) ) {

			//	if (submat.reverseMap(nonSAsequence[i]) == 'N') System.exit(1);

				// search current word of the non-SA sequence in the suffix array (in both directions)
				int[] forwardhits = sa.search(nonSAsequence, i, params.wordLength, false, idmap);

				int[] reversehits;
				if (i >= params.wordLength-1) {
					reversehits = sa.search(nonSAsequence, i, params.wordLength, true, complmap);
				} else
					reversehits = new int[0];
				
				// iterate through both hit arrays
				hits=forwardhits;
				boolean firsthits=true;
				int j=0;
				
				
				//for (int j = 0; j < hits.length; j++) {
				while (true) {
					// check indices
					
					if (j == hits.length)
						 if (firsthits) {
							 if (reversehits.length > 0) {
								 firsthits = false;
								 hits = reversehits;
								 j=0;
							 } else
								 break;
						 } else
							 break;
					
					// check if hit is in bounds
					if ((hits[j] >= SAseqStart)	&& (hits[j] <= SAseqStop)) {
						// transform coordinates
						float tx = (float) (hits[j] - SAseqStart) / params.zoom;
						float ty = (float) (i - nonSAseqStart) / params.zoom;
						
						float value = 0;
					//	// use windowsize?
					//	if (params.windowSize == 0) // no -> score 1
							value = 1;
					//	else
					//		value = getWindowScore(seq1, seq2, hits[j], i, params.windowSize);

						int x1=0,x2=0,y1=0,y2=0;
					 
						// set proper coordinates (switch if SA for second sequence)
						if (!SAforSecondSeq) {
							// determine the 4 involved dot matrix spots
							x1 = (int) tx;
							x2 = ((int) tx) + 1;
							y1 = (int) ty;
							y2 = ((int) ty) + 1;		
					
						} else {
							// determine the 4 involved dot matrix spots
							x1 = (int) ty;
							x2 = ((int) ty) + 1;
							y1 = (int) tx;
							y2 = ((int) tx) + 1;		
						}
						// calc fractions
						float fracx2 = tx - (float) ((int) tx);
						float fracx1 = 1 - fracx2;
						float fracy2 = ty - (float) ((int) ty);
						float fracy1 = 1 - fracy2;
						
						// distribute values to array
						if ((x1 < width) && (y1 < height)) 
							dotmatrix[x1][y1] += fracx1 * fracy1 * value;
						if ((x2 < width) && (y1 < height))
							dotmatrix[x2][y1] += fracx2 * fracy1 * value;
						if ((x1 < width) && (y2 < height)) 
							dotmatrix[x1][y2] += fracx1 * fracy2 * value;
						if ((x2 < width) && (y2 < height)) 
							dotmatrix[x2][y2] += fracx2 * fracy2 * value;
						
						
					}
					
					// increase index and shift hit array if necessary
					j++;
				
				}
				
				// check for callback
				if (i % callbackStep == 0) {
					callback.dotmatrixCalcStatus(
						// from:  params.seq2Start
						// to  : params.seq2Stop - params.wordLength + 1
						(float)(i-nonSAseqStart) / (float)(nonSAseqStop - params.wordLength + 1-nonSAseqStart) * 100
					);
					// abort?
					if (ac.dotplotAborted()) {
						callback.tellAborted();
						return;
					}
				}
			
			}
			
		}			
		// now determine highest and lowest value
		calcMaxMinAvg(false);
		
	}
	/*
	// matrix calculation (wordlength == 0)
	private void calcDotMatrixWordlenZeronew(SuffixArray sa, byte[] seq1,
			byte[] seq2, DotMatrixCallback callback, AbortionChecker ac) {
		
		// determine callback step
//		int callbackStep = (callback != null) ? callback.tellCallbackStep(params.wordLength, params.windowSize ) : 2147483647;
	
		// create dot matrix array
		dotmatrix = new float[width][height];
		
		
		// create identity and complementary maps
		byte[] idmap = new byte[128];
		for (byte i=0; i<128 && i >= 0; i++) idmap[i] = i;
		// use complementary map if no protein sequences are used
		byte[] complmap ;
		if (submat.isNucleotideMatrix())
			complmap = Sequence.getComplementaryMap();
		else 
			complmap = idmap;
		
		int winsize = params.windowSize;
		
		long progress=0;
		int width = params.seq1Stop-params.seq1Start+1;
		int height = params.seq2Stop-params.seq2Start+1;
		long progresstotal = ((width*height) - (width+height-1)) * 2;
		
		
		// *** RUN 1, forward matches
		
		// iterate over left and upper border
		int i=params.seq1Start,j=params.seq2Stop;
		
		while (i<=params.seq1Stop) {
			
			int x = i;
			int y = j;
			
			// calculate initial window bounds
			int backward = (winsize/2);
			if (x-backward < 0)
				backward -= (backward-x);
			if (y-backward < 0)
				backward -= (backward-y);
			int forward = (winsize/2);
			
			float score=0;
			// calculate initial window
			for (int k=-backward; (k<=forward) && (x + k<seq1.length) && (y + k<seq2.length) ; k++) {
				// calculate score
				score += submat.getScore(seq1[x + k], seq2[y + k]);
			}
			
			// write score to dotmatrix
			distributeScore(x, y, score);
			
			// now extend the diagonal as long as possible
			x+=1; 
			y+=1;
			while ((x<seq1.length) && (y<seq2.length)) {
				// check if the upper left end of the current diagonal was in range
				if ((x-backward-1 >= 0) && (y-backward-1 >= 0))  
					// subtract upper left end
					score -= submat.getScore(seq1[x-backward-1], seq2[y-backward-1]);
				// check if the lower right end of the current diagonal is in range
				if ((x+forward < seq1.length) && (y+forward<seq2.length)) 
					// add new lower right end
					score += submat.getScore(seq1[x], seq2[y]);
				// write score to dotmatrix
				distributeScore(x, y, score);
				// extend further
				x++; y++;
				
				progress++;
			}
			
			// check for callback
			callback.dotmatrixCalcStatus(
					(float)progress/(float)progresstotal*100f
			);
			// abort?
			if (ac.dotplotAborted()) {
				callback.tellAborted();
				return;
			}
			
			// increase variables
			if (j > params.seq2Start) j--;
			else i++;
		}
		
		// *** RUN 2, reverse matches
		
		// iterate over upper and right border
		i=params.seq1Stop;
		j=params.seq2Stop;
		
		while (i>=params.seq1Start) {
			
			int x = i;
			int y = j;
			
			
			// calculate initial window bounds
			int backward = (winsize/2);
			if (x+backward > seq1.length)
				backward = seq1.length - x -1;
			if (y-backward < 0)
				backward -= (backward-y);
			int forward = (winsize/2);
			
			
			float score=0;
			// calculate initial window
			for (int k=-backward; (k<=forward) && (x - k>=0) && (y + k<seq2.length) ; k++) {
				// calculate score
				score += submat.getScore(seq1[x - k], complmap[seq2[y + k]]);
			}
			
			// write score to dotmatrix
			distributeScore(x, y, score);
			
			// now extend the diagonal as long as possible
			x-=1; 
			y+=1;
			while ((x>=0) && (y<seq2.length)) {
				// check if the upper left end of the current diagonal was in range
				if ((x+backward+1 < seq1.length) && (y-backward-1 >= 0))  
					// subtract upper left end
					score -= submat.getScore(seq1[x+backward+1], complmap[seq2[y-backward-1]]);
				// check if the lower right end of the current diagonal is in range
				if ((x-forward >= 0) && (y+forward<seq2.length)) 
					// add new lower right end
					score += submat.getScore(seq1[x], complmap[seq2[y]]);
				// write score to dotmatrix
				distributeScore(x, y, score);
				// extend further
				x--; y++;
				
				progress++;
			}
			
			// check for callback
			callback.dotmatrixCalcStatus(
					(float)progress/(float)progresstotal*100f
			);
			// abort?
			if (ac.dotplotAborted()) {
				callback.tellAborted();
				return;
			}
		
			// increase variables
			if (j > params.seq2Start) j--;
			else i--;
			
		}
		

		// now determine highest and lowest value
		calcMaxAndMin();
		
	}
	*/
	/*
	private void distributeScore(int x, int y, float score) {
		// transform coordinates
		float tx = (float) (x - params.seq1Start) / params.zoom;
		float ty = (float) (y - params.seq2Start) / params.zoom;
		// determine the 4 involved dot matrix spots
		int x1 = (int) tx;
		int x2 = ((int) tx) + 1;
		int y1 = (int) ty;
		int y2 = ((int) ty) + 1;
		// calc fractions
		float fracx2 = tx - (float) ((int) tx);
		float fracx1 = 1 - fracx2;
		float fracy2 = ty - (float) ((int) ty);
		float fracy1 = 1 - fracy2;
		// distribute value to array
		if ((x1 < width) && (y1 < height))
			dotmatrix[x1][y1] += fracx1 * fracy1 * score;
		if ((x2 < width) && (y1 < height))
			dotmatrix[x2][y1] += fracx2 * fracy1 * score;
		if ((x1 < width) && (y2 < height))
			dotmatrix[x1][y2] += fracx1 * fracy2 * score;
		if ((x2 < width) && (y2 < height))
			dotmatrix[x2][y2] += fracx2 * fracy2 * score;
	}
*/

	private void calcMaxMinAvg(boolean windowmode) {
		
		// <CAP> calc cap
		float cap;
		if (windowmode)
			cap = Float.POSITIVE_INFINITY;
		else
			cap = params.zoom  / 4f;
		
		
		// find minimum and maximum dot score
		maxdotscore = Float.NEGATIVE_INFINITY;
		mindotscore = Float.POSITIVE_INFINITY;
		// iterate over whole matrix
		for (int i = 0; i < width; i++) {
			for (int j = 0; j < height; j++) {
				if (dotmatrix[i][j] > cap)
					dotmatrix[i][j] = cap;
				if (dotmatrix[i][j] > maxdotscore)
					maxdotscore = dotmatrix[i][j];
				if (dotmatrix[i][j] < mindotscore)
					mindotscore = dotmatrix[i][j];
				avgdotscore += dotmatrix[i][j];
			}
		}
		avgdotscore /= (width * height);
	}

	private void calcGCRatios(byte[] s1, byte[] s2, ParameterSet params) {
		
		// iterate through first sequence
		int GCcount=0;
		for (int i=params.seq1Start ; i<=params.seq1Stop ; i++) {
			if (s1[i] == 3 || s1[i] == 4)
				GCcount++;
		}
		
		// store value
		GCratio1 = (float)GCcount / (float)(params.seq1Stop - params.seq1Start + 1);
		
		// iterate through second sequence
		GCcount=0;
		for (int i=params.seq2Start; i<=params.seq2Stop; i++) {
			if (s2[i] == 3 || s2[i] == 4)
				GCcount++;
		}
		
		// store value
		GCratio2 = (float)GCcount / (float)(params.seq2Stop - params.seq2Start + 1);
	}
	
	
	// constructor for existing dotplot data
	// previously used to re-construct remote dotplot, not used anymore
	public DotMatrix(float[][] matrix, float maxDotScore, float minDotScore,
			String seq1Name, int seq1Length, float seq1GCRatio,
			String seq2Name, int seq2Length, float seq2GCRatio,
			ParameterSet params) {
		
		// just copy values
		this.dotmatrix = matrix;
		this.maxdotscore = maxDotScore;
		this.mindotscore = minDotScore;
		this.seq1Name = seq1Name;
		this.seq1Length = seq1Length;
		this.GCratio1 = seq1GCRatio;
		this.seq2Name = seq2Name;
		this.seq2Length = seq2Length;
		this.GCratio2 = seq2GCRatio;
		this.params = params;
		this.nuclmatrix = true;
		
		// set additional values
		width = matrix.length;
		height = matrix[0].length;
	}
	

	public CompressedDotMatrix getCompressedMatrix() {
		byte[] intmatrix = new byte[width*height];
		float range = maxdotscore - mindotscore;
		int newmaxscore=Integer.MIN_VALUE, newminscore=Integer.MAX_VALUE;
		// iterate row-by-row
		for (int j=0; j<height; j++) {
			for (int i=0; i<width; i++) {
				int newval = Math.round ((dotmatrix[i][j] - mindotscore) / range * 255f );
				if (newval < newminscore) newminscore = newval;
				if (newval > newmaxscore) newmaxscore = newval;
				intmatrix[j*width+i] = (byte)newval;
			}
		}
		
		return new CompressedDotMatrix(intmatrix, newmaxscore, newminscore);
	}
	
	
	////////////////////////////////////////////////////////////////////
	
	// matrix calculation (wordlength == 0)
	private void calcDotMatrixWindowMode(SuffixArray sa, byte[] seq1,
			byte[] seq2, DotMatrixCallback callback, AbortionChecker ac) {
		
		// determine callback step
		int callbackStep = (callback != null) ? callback.tellCallbackStep(params.wordLength, params.windowSize ) : 2147483647;
	
		// create dot matrix array
		dotmatrix = new float[width][height];
		
		// create identity and complementary maps
		byte[] idmap = new byte[128];
		for (byte i=0; i<128 && i >= 0; i++) idmap[i] = i;
		// use complementary map if no protein sequences are used
		byte[] complmap ;
		if (submat.isNucleotideMatrix())
			complmap = Sequence.getComplementaryMap();
		else 
			complmap = idmap;
		// iterate through sequences
		//for (int i = params.seq2Start; i < params.seq2Stop - params.wordLength	+ 1; i++) {
		//	for (int j = params.seq1Start; j < params.seq1Stop- params.wordLength + 1; j++) {
		for (int i = params.seq2Start; i <= params.seq2Stop; i++) {
			for (int j = params.seq1Start; j <= params.seq1Stop; j++) {
				// transform coordinates
				float tx = (float) (j - params.seq1Start) / params.zoom;
				float ty = (float) (i - params.seq2Start) / params.zoom;

				float value = 0;
				// get window score
				value = getWindowScore(seq1, seq2, j, i, params.windowSize)
						+ getWindowScoreBackward(seq1, seq2, j, i, params.windowSize, complmap);
				
				// determine the 4 involved dot matrix spots
				int x1 = (int) tx;
				int x2 = ((int) tx) + 1;
				int y1 = (int) ty;
				int y2 = ((int) ty) + 1;

				// calc fractions
				float fracx2 = tx - (float) ((int) tx);
				float fracx1 = 1 - fracx2;
				float fracy2 = ty - (float) ((int) ty);
				float fracy1 = 1 - fracy2;

				// distribute value to array
				if ((x1 < width) && (y1 < height))
					dotmatrix[x1][y1] += fracx1 * fracy1 * value;
				if ((x2 < width) && (y1 < height))
					dotmatrix[x2][y1] += fracx2 * fracy1 * value;
				if ((x1 < width) && (y2 < height))
					dotmatrix[x1][y2] += fracx1 * fracy2 * value;
				if ((x2 < width) && (y2 < height))
					dotmatrix[x2][y2] += fracx2 * fracy2 * value;

			}
			
			// check for callback
			if (i % callbackStep == 0) {
				callback.dotmatrixCalcStatus(
					// from:  params.seq2Start
					// to  : params.seq2Stop - params.wordLength + 1
					(float)(i-params.seq2Start) / (float)(params.seq2Stop - params.wordLength + 1-params.seq2Start) * 100
				);
				// abort?
				if (ac.dotplotAborted()) {
					callback.tellAborted();
					return;
				}
			}			
		}

		// now determine highest and lowest value
		calcMaxMinAvg(true);
		
	}
	
	
	
//	 calculate the score of a given sequence window
	// (used for wordlength==0 calculations)
	private float getWindowScore(byte[] seq1, byte[] seq2, int x, int y,int windowsize) {
		
		float score = 0;
		int consecutives=0;
		
		int backforward = windowsize/2;
		
		if (nuclmatrix) {
			for (int i = -backforward; (i <= backforward) ; i++) {
				if ((x+i>0) && (y+i>0) && (x + i < seq1.length) && (y + i < seq2.length)) {
					if (seq1[x + i] == seq2[y + i]) {
						score += 1;
						score += consecutives;
						if (consecutives < 7) consecutives++;
					}
					else {
						//score--;
						consecutives=0;
					}
				}
			}
		} else {
			// protein dotplot
			for (int i = -backforward; (i <= backforward) ; i++) {
				if ((x+i>0) && (y+i>0) && (x + i < seq1.length) && (y + i < seq2.length)) {
					// add score determined from substitution matrix
					score += submat.getScore(seq1[x + i], seq2[y + i]);
					
//					System.out.println(submat.getScore(seq1[x + i], seq2[y + i]));
				}
				
			}
		}
		
		
		return score;
	}
	
	
	// calculate the backward score of a given  sequence window
	// (used for wordlength==0 calculations)
	private float getWindowScoreBackward(byte[] seq1, byte[] seq2, int x, int y,
			int windowsize, byte[] map) {
		
		float score = 0;
		int consecutives=0;
		int backforward = windowsize/2;
		
		if (nuclmatrix) {
			for (int i = -backforward; (i <= backforward) ; i++) {
				if ((x-i>0) && (y+i>0) && (x - i < seq1.length) && (y + i < seq2.length)) {
					if (seq1[x - i] == map[seq2[y + i]]) {
						score += 1;
						score += consecutives;
						if (consecutives < 7) consecutives++;
					}
					else {
						//score--;
						consecutives=0;
					}
				}
			}
		}
		
		return score;

	/*	for (int i = 0; (i < windowsize) && (x - i >= 0)
				&& (y + i < seq2.length); i++) {
			// add score determined from substitution matrix
			score += submat.getScore(seq1[x - i], map[seq2[y + i]]);
		}
		return score;*/
	}

	
} public interface DotMatrixCallback {
	public void dotmatrixCalcStatus(final float percent);
	public int tellCallbackStep(int wordlen, int windowsize);
	public void tellAborted();
} public class FASTAReader {
	
	private static final float NUCLEOTIDE_THRESHOLD = 0.95f;
	private static final int CHARACTERS_FOR_AUTOMATRIX = 250;

	public static Sequence readFile(String file, SubstitutionMatrix substmat)
			throws FileNotFoundException, IOException, InvalidFASTAFileException {
		String name = null;
		byte[] sequence = null;
		
		boolean invalidchars=false;
		
		int ACGTcount = 0;
		boolean[] isACGT = getACGTChecker();
		int charcount=0;
		
		// get lowercase->uppercase mapper
		byte[] LCToUC = getUppercaseMapper();
		
		// get file length
		File f = new File(file);
		int filelen = (int) f.length();
		if (filelen<0) 
			throw new InvalidFASTAFileException("Files larger than 2.1GB are not supported");
		if (filelen==0) 
			throw new InvalidFASTAFileException("Empty file");

		// get input stream
		DataInputStream in = new DataInputStream(new FileInputStream(file));
		// read full contents into byte array
		byte[] contents = new byte[filelen];
		in.readFully(contents, 0, filelen);
		in.close();
		// check for status line
		if (contents.length == 0 || contents[0] != '>')
			// invalid format
			throw new InvalidFASTAFileException("No status line found");

		// search first line break, store name
		int i;
		for (i = 0; contents[i] != 10 && contents[i] != 13; i++)	;
		// one more?
		if (contents[i+1] == 10) i++;
		name = new String(contents, 1, i - 1).trim() ;
		// now parse rest of file

		// first run: count valid characters, overall characters and ACGT amount
		int validcount = 0;
		for (int j = i; j < contents.length; j++) {
			
			// new sequence?
			while (j<contents.length && contents[j] == '>') {
				// ignore status line
				while (j<contents.length && contents[j] != 10 && contents[j] != 13) j++;
				j++;
				// another LF here?
				if (j<contents.length && contents[j] == 10) j++;
				// 
			}
			
			if (j>=contents.length)
				break;
			
					
			// valid character?
			if (contents[j]<0)
				throw new InvalidFASTAFileException("Invalid character in FASTA file: " + (char)contents[j]);
	 
			// convert character to uppercase
			contents[j] = LCToUC[contents[j]]; 
			// valid character for current marix?
			if (substmat.map(contents[j]) > 0)
				validcount++;
			// any character except 13 and 10?
			if (substmat.map(contents[j]) >= 0)
				charcount++;
			// A, C, G or T?
			if (isACGT[contents[j]]) 
				ACGTcount++;

		}
		
		boolean isnucleotide = ((float)ACGTcount/(float)charcount) >= NUCLEOTIDE_THRESHOLD;
		
		// prepare multi-fasta
		List<String> allNames = new ArrayList<String>();
		List<Integer> allStarts = new ArrayList<Integer>();
		// add first one
		allNames.add(name);
		allStarts.add(0);
		
		// second run: create and fill array
		sequence = new byte[validcount];
		int k = 0;
		for (int j = i; j < contents.length; j++) {
			
			// new sequence?
			while (j<contents.length && contents[j] == '>') {
				int statusStart=j;
				// ignore status line
				while (j<contents.length && contents[j] != 10 && contents[j] != 13) j++;
				j++;
				allNames.add(new String(contents,statusStart+1,j-statusStart-2));
				// another LF here?
				if (j<contents.length && contents[j] == 10) j++;
				allStarts.add(k);
			}
			if (j>=contents.length)
				break;

			// add to sequence or output error
			byte mapval = substmat.map(contents[j]);
			// map U to T
			if (isnucleotide && mapval == 16)
				mapval = 2;
			// map unknown chars to N
			if (isnucleotide && mapval > 5)
				mapval = 5;

			if (mapval > 0)
				sequence[k++] = mapval;	
			else if (mapval == 0) 
				invalidchars = true;
		}
		
		// create and return LocalSequence object
		Sequence seq = null;
		if (allNames.size()>1) {
			name = "Multi FASTA (" + ClientGlobals.extractFilename(file) + ")";
			seq = new Sequence(sequence, name, isnucleotide, invalidchars);
			seq.setMulti(allNames, allStarts);
		} else
			seq = new Sequence(sequence, name, isnucleotide, invalidchars);
		
		return seq; 
	}
	
	private static byte[] getUppercaseMapper() {
		byte[] ret = new byte[256];
		Arrays.fill(ret,(byte)10);	// convert all unsupported character to LF=10
									// because LF is ignored anyway
		// write upper case letter A-Z
		for (byte i=65; i<=90; i++)
			ret[i] = i;
		// write lower case letters a-z and convert to uppercase 
		for (byte i=97; i<=122; i++)
			ret[i] = (byte)(i-32);
		
		return ret;
		
	}
	
	private static boolean[] getACGTChecker() {
		boolean[] ret = new boolean[256];
		ret['A'] = true;
		ret['C'] = true;
		ret['G'] = true;
		ret['T'] = true;
		ret['N'] = true;
		ret['X'] = true;
		return ret;
	}
	
	public static boolean isNucleotideFile(String file) throws IOException, InvalidFASTAFileException {
		
		// get lowercase->uppercase mapper && ACGT checker
		byte[] LCToUC = getUppercaseMapper();
		boolean[] isACGT = getACGTChecker();
		
		// get input stream
		DataInputStream in = new DataInputStream(new FileInputStream(file));
		// read first CHARACTERS_FOR_AUTOMATRIX bytes
		int readbytes = (int)Math.min(CHARACTERS_FOR_AUTOMATRIX, new File(file).length());
		byte[] contents = new byte[readbytes];
		in.readFully(contents, 0, readbytes);
		in.close();
		
		// check for status line
		if (contents.length == 0 || contents[0] != '>')
			// invalid format
			throw new InvalidFASTAFileException("No status line found");
		
		int start=0;
		
		// skip status line      
		while (start < contents.length && contents[start] != 10 && contents[start] != 13) start++;
		// did we find no line break?
		if (start == contents.length) throw new InvalidFASTAFileException("End of FASTA status line could not be identified");
		// one more?
		if (contents[start+1] == 10) start++;
		// and the last one to reach the new line
		start++;
		
		
		
		int totalCount=0, acgtCount=0;
		
		// count ACGT and total characters
		for (int i=start; i<contents.length; i++) {
		
			while (i<contents.length && contents[i] == '>') {
				// ignore status line
				while (i<contents.length && contents[i] != 10 && contents[i] != 13) i++;
				i++;
				// another LF here?
				if (i<contents.length && contents[i] == 10) i++;
			}
			if (i>=contents.length)
				break;
			
		
			if (contents[i] != 13 && contents[i] != 10) {
				if (contents[i]<0)
					throw new InvalidFASTAFileException("Invalid character in FASTA file: " + (char)contents[i]);
				if (isACGT[LCToUC[contents[i]]])
					acgtCount++;
				totalCount++;
			}
		}
		
		// check ACGT/total ratio
		if  ( ((double)acgtCount / (double)totalCount ) >= NUCLEOTIDE_THRESHOLD )
			return true;
		else
			return false;
		
	}
	
	public static int countNuclSeqLength(String file) throws IOException {
		
		int charcount=0;
		
		BufferedReader reader = new BufferedReader(new FileReader(file));
		
		// ignore first line -> status
		//reader.readLine();
		
		String line;
		while ((line = reader.readLine()) != null) {
			// ignore status lines
			if (!line.startsWith(">"))
				charcount += line.length();
		}
		
		reader.close();
		
		return charcount;
		
	}
	
	public static void main(String[] args) throws FileNotFoundException, IOException, InvalidFASTAFileException, InvalidSubMatFileException {
				
		Sequence s = readFile("/home/jan/work/workspace/Gepard/stuff/bad.fa", ClientGlobals.getEDNA());
		System.out.println(s.isMulti());
//		
//		s = readFile("/home/jan/work/workspace/Gepard/stuff/ecoli_part_1.fa", ClientGlobals.getEDNA());
//		System.out.println(s.isMulti());
		
		System.out.println(s.hasInvalidChars());
		System.out.println(s.getLength());
		System.out.println(s.likelyNucleotides());
//		System.out.println(isNucleotideFile("/home/jan/work/workspace/Gepard/stuff/bad.fa"));
	}
	
} public class InvalidFASTAFileException extends Exception {

	private static final long serialVersionUID = 1234756066690198181L;

	public InvalidFASTAFileException(String msg) {
		super(msg);
	}
} public class InvalidSubMatFileException extends Exception {

	private static final long serialVersionUID = 6116954616817540722L;

	public InvalidSubMatFileException(String msg) {
		super(msg);
	}
} public class Needleexport {
	public static void export(Sequence sequence1, Sequence sequence2,int tX,int tY, String path1, int window, boolean complementary, SubstitutionMatrix submat) {
		int backward = window / 2;
		int forward = window / 2;
		byte[] seq1 = sequence1.getSequenceData();
		byte[] seq2 = sequence2.getSequenceData();
		// empty position byte arrays
		char [] strpos1 = new char [window];
		char [] strpos2 = new char [window];
		Arrays.fill(strpos1, ' '); 
		Arrays.fill(strpos2, ' ');
		
		byte[] complmap = Sequence.getComplementaryMap();
		char[] strseq1 = new char[window];
		char[] strseq2 = new char[window];
		char [] strsim = new char [window];
		// iterate through range
		int j=0;
		for (int i=-backward; i<=forward&&j<window; i++) {
			byte a=-1, b=-1;
			// get characters from both sequences
			if (tX+i >= 0 && tX+i < seq1.length) 
				a = seq1[tX+i];
			if (!complementary) {
				if (tY+i >= 0 && tY+i < seq2.length)
					b = seq2[tY+i];
			} else {
				if (tY-i >= 0 && tY-i < seq2.length)
					b = complmap[seq2[tY-i]];
			}
			
				// insert characters into strings
			if (a != -1) 	strseq1[j] = submat.reverseMap(a);
			else 			strseq1[j] = ' ';
			
			if (b != -1) 	strseq2[j] = submat.reverseMap(b);
			else 			strseq2[j] = ' ';
			j++;	
		}
		int [][] scoringmatrix = {
				{5,-3,-3,-3},
				{-3,5,-3,-3},
				{-3,-3,5,-3},
				{-3,-3,-3,5}
			};
		if (!submat.isNucleotideMatrix()) {
			scoringmatrix = submat.getMatrix();
		}
		NeedlemanWunsch matrices = NeedlemanWunsch.GetAlignmentMatrix(String.valueOf(strseq1), String.valueOf(strseq2), scoringmatrix, submat);
		List<String> ans = new ArrayList<String>(NeedlemanWunsch.GetAlignment(String.valueOf(strseq1), String.valueOf(strseq2), scoringmatrix, matrices.align_matrix, matrices.matrix_pointer));
		strseq1 = new char[window];
		strseq2 = new char[window];
		char [] finalstr1 = new char[ans.get(0).length()];
		char [] finalstr2 = new char[ans.get(0).length()];
		for (int a = 0; a<ans.get(0).length(); a ++) {
			finalstr1[a] = (char) ans.get(0).charAt(a);
		}
		for (int a = 0; a<ans.get(1).length(); a ++) {
			finalstr2[a] = (char) ans.get(1).charAt(a);
		}
		for (int a = 0; a<strseq1.length&&a<window; a++) {
			if (strseq1[a] == strseq2[a] && strseq1[a] != ' ') {
				strsim[a] = getNuclSimCode(submat.map((byte) strseq1[a]), submat.map((byte) strseq1[a]), submat);;
			}
			else {
				strsim[a] = ' ';
			}
		}
		int FASTA_COLS = 70;
		// open file
		PrintWriter writer;
		try {
			writer = new PrintWriter(new FileWriter(path1));
			String rangeString = ", range: " + ((tX -(window / 2)) + Integer.parseInt(ans.get(2))) + "-"+((tX-(window / 2)) + Integer.parseInt(ans.get(4)));
			writer.println(">" + ClientGlobals.cutString(sequence1.getName(), FASTA_COLS -rangeString.length()-1) + rangeString);
			// iterate over range
			int col=0;
			for (int i=0; i<finalstr1.length; i++) {
				writer.print(finalstr1[i]);
				col++;
				if (col % FASTA_COLS == 0) writer.println();
			}
			if (col % FASTA_COLS != 0) writer.println();
			String rangeString2 = ", range: " + ((tY -(window / 2)) + Integer.parseInt(ans.get(3))) + "-"+((tY -(window / 2)) + Integer.parseInt(ans.get(5)));
			writer.println(">" + ClientGlobals.cutString(sequence2.getName(), FASTA_COLS -rangeString2.length()-1) + rangeString2);
			// iterate over range
			col=0;
			for (int i=0; i<finalstr2.length; i++) {
				writer.print(finalstr2[i]);
				col++;
				if (col % FASTA_COLS == 0) writer.println();
			}
			writer.close();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		Arrays.fill(strseq1, ' '); 
		Arrays.fill(strseq2, ' ');
		Arrays.fill(finalstr1, ' '); 
		Arrays.fill(finalstr2, ' ');
	}

	private static char getNuclSimCode(byte map, byte map2, SubstitutionMatrix submat) {
		if (map != -1 && map2 != -1) {
			if (map==map2)
				return ':';
			else {
				int score = submat.getScore(map,map2);
				if (score > 0)
					return '.';
				else
					return ' ';
			}
		}
		else
			return ' ';
	}
} public class ParameterSet implements Serializable {

	private static final long serialVersionUID = 2991408205215417902L;

	public int zoom;
	public int wordLength;
	public int windowSize;
	public int seq1Start, seq1Stop;
	public int seq2Start, seq2Stop;
	public String submatName;
	
	
	private static final String paramDivider = "&";

	public ParameterSet(int iratio, int iwordlen, int iwinsize, int iseq1start,
			int iseq1stop, int iseq2start, int iseq2stop, String submatName) {
		
		this.zoom = iratio;
		this.wordLength = iwordlen;
		this.windowSize = iwinsize;
		this.seq1Start = iseq1start;
		this.seq1Stop = iseq1stop;
		this.seq2Start = iseq2start;
		this.seq2Stop = iseq2stop;
		this.submatName = submatName;
	}
	
	
	
	public ParameterSet() {
	}

	public ParameterSet getClone() {
		return new ParameterSet(zoom, wordLength, windowSize, seq1Start,
				seq1Stop, seq2Start, seq2Stop, submatName);
	}
	
	public String serialize() {
		return zoom + paramDivider + wordLength + paramDivider + windowSize + paramDivider +
			seq1Start + paramDivider + seq1Stop + paramDivider + seq2Start + 
			paramDivider + seq2Stop + paramDivider + submatName + paramDivider;
	}
	
	public String toString() {
		return "Zoom: " + zoom+ ", wordLength: " + wordLength+ ", windowSize: " + windowSize + ", " +
			"seq1: " + this.seq1Start + "-" + this.seq1Stop + ", " + 
			"seq2: " + this.seq2Start + "-" + this.seq2Stop + "  ";
	}
	
	// constructor which creates a parameter set from serialized data
	public ParameterSet(String serialized) {
		String[] values = serialized.split(paramDivider);

		zoom = Integer.parseInt(values[0]);
		wordLength = Integer.parseInt(values[1]);
		windowSize = Integer.parseInt(values[2]);
		seq1Start = Integer.parseInt(values[3]);
		seq1Stop = Integer.parseInt(values[4]);
		seq2Start = Integer.parseInt(values[5]);
		seq2Stop = Integer.parseInt(values[6]);
		submatName = values[7];
	}

} public class Sequence {
	
	private byte[] sequenceData;
	private String name;
	private boolean likelynucleotides;
	private boolean invalidchars;
	private boolean complementary=false;
	
	private boolean isMulti;
	private String[] multiNames;
	private int[] multiStarts;
	
	public Sequence(byte[] sequence, String name, boolean likelynucleotides, boolean invalidchars) {
		this.name = name;
		this.sequenceData = sequence;
		this.likelynucleotides = likelynucleotides;
		this.invalidchars = invalidchars;
	}
	
	// accessor methods
	public byte[] getSequenceData() { return sequenceData; }
	public String getName() { return ( (complementary) ? "[compl.]" : "").concat(name); }
	public int getLength() { return sequenceData.length; }
	public boolean likelyNucleotides() { return likelynucleotides; }
	public boolean hasInvalidChars() { return invalidchars; }
	
	public void setComplementaryFlag() {
		complementary = true;
	}
	
	
	private static byte[] complement=null;
	public static void complementarizeSequence(Sequence sequence) {
		// initialize table?
		if (complement == null) 
			complement = getComplementaryMap();
		
		
		byte[] seq = sequence.getSequenceData();
		for (int i=0; i<seq.length; i++)
			seq[i] = complement[seq[i]];
		
		
		sequence.setComplementaryFlag();
		
		// 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
		// A  T  G  C  S  W  R  Y  K  M  B  V  H  D  N  U  Z
		// T  A  C  G  S  W  Y  R  M  K  V  B  D  H  N  A  Z
		
		/*
		A         A         T
		C         C         G
		G         G         C
		T/U         T         A
		M         A/C         K
		R         A/G         Y
		W         A/T         W
		S         C/G         S
		Y         C/T         R
		K         G/T         M
		V         A/C/G         B
		H         A/C/T         D
		D         A/G/T         H
		B         C/G/T         V
		X/N         A/C/G/T         X
		.         None         .
		*/
	}
	
	public static byte[] getComplementaryMap() {
		byte[] complement = new byte[128];
		complement[1] = 2;
		complement[2] = 1;
		complement[3] = 4;
		complement[4] = 3;
		complement[5] = 5;
		/*
		complement[6] = 6;
		complement[7] = 8;
		complement[8] = 7;
		complement[9] = 10;
		complement[10] = 9;
		complement[11] = 12;
		complement[12] = 11;
		complement[13] = 14;
		complement[14] = 13;
		complement[15] = 15;
		complement[16] = 1;
		complement[17] = 17;*/
		return complement;
	}

	
	// multi-fasta stuff
	public void setMulti(List<String> names, List<Integer> starts) {
		multiNames = names.toArray(new String[0]);
		multiStarts = new int[starts.size()];
		int i=0;
		for (int v : starts)
			multiStarts[i++]=v;
		isMulti = true;
	}
	
	public boolean isMulti() {
		return isMulti;
	}
	
	public String[] getMultiNames() {
		return multiNames;
	}

	public int[] getMultiStarts() {
		return multiStarts;
	}
} public class Smithexport {
	public static void export(Sequence sequence1, Sequence sequence2,int tX,int tY, String path1, int window, boolean complementary, SubstitutionMatrix submat) throws IOException {
		int backward = window / 2;
		int forward = window / 2;
		byte[] seq1 = sequence1.getSequenceData();
		byte[] seq2 = sequence2.getSequenceData();
		// empty position byte arrays
		char [] strpos1 = new char [window];
		char [] strpos2 = new char [window];
		Arrays.fill(strpos1, ' '); 
		Arrays.fill(strpos2, ' ');
		
		byte[] complmap = Sequence.getComplementaryMap();
		char[] strseq1 = new char[window];
		char[] strseq2 = new char[window];
		char [] strsim = new char [window];
		// iterate through range
		int j=0;
		for (int i=-backward; i<=forward&&j<window; i++) {
			byte a=-1, b=-1;
			// get characters from both sequences
			if (tX+i >= 0 && tX+i < seq1.length) 
				a = seq1[tX+i];
			if (!complementary) {
				if (tY+i >= 0 && tY+i < seq2.length)
					b = seq2[tY+i];
			} else {
				if (tY-i >= 0 && tY-i < seq2.length)
					b = complmap[seq2[tY-i]];
			}
			
				// insert characters into strings
			if (a != -1) 	strseq1[j] = submat.reverseMap(a);
			else 			strseq1[j] = ' ';
			
			if (b != -1) 	strseq2[j] = submat.reverseMap(b);
			else 			strseq2[j] = ' ';
			j++;	
		}
		int[][] scoringmatrix = {
				{5,-3,-3,-3},
				{-3,5,-3,-3},
				{-3,-3,5,-3},
				{-3,-3,-3,5}
		};
		if (!submat.isNucleotideMatrix()) {
			scoringmatrix = submat.getMatrix();
		}
		SmithWaterman matrices = SmithWaterman.GetAlignmentMatrix(String.valueOf(strseq1), String.valueOf(strseq2), scoringmatrix, submat);
		List<String> ans = new ArrayList<String>(SmithWaterman.GetAlignment(String.valueOf(strseq1), String.valueOf(strseq2), scoringmatrix, matrices.align_matrix, matrices.matrix_pointer));
		strseq1 = new char[window];
		strseq2 = new char[window];
		char [] finalstr1 = new char[ans.get(0).length()];
		char [] finalstr2 = new char[ans.get(0).length()];
		for (int a = 0; a<ans.get(0).length(); a ++) {
			finalstr1[a] = (char) ans.get(0).charAt(a);
		}
		for (int a = 0; a<ans.get(1).length(); a ++) {
			finalstr2[a] = (char) ans.get(1).charAt(a);
		}
		for (int a = 0; a<strseq1.length&&a<window; a++) {
			if (strseq1[a] == strseq2[a] && strseq1[a] != ' ') {
				strsim[a] = getNuclSimCode(submat.map((byte) strseq1[a]), submat.map((byte) strseq1[a]), submat);;
			}
			else {
				strsim[a] = ' ';
			}
		}
		int FASTA_COLS = 70;
		// open file
		PrintWriter writer;
			writer = new PrintWriter(new FileWriter(path1));
			String rangeString = ", range: " + ((tX- (window / 2)) + Integer.parseInt(ans.get(2))) + "-"+((tX -(window / 2)) + Integer.parseInt(ans.get(4)));
			writer.println(">" + ClientGlobals.cutString(sequence1.getName(), FASTA_COLS -rangeString.length()-1) + rangeString);
			// iterate over range
			int col=0;
			for (int i=0; i<finalstr1.length; i++) {
				writer.print(finalstr1[i]);
				col++;
				if (col % FASTA_COLS == 0) writer.println();
			}
			if (col % FASTA_COLS != 0) writer.println();
			String rangeString2 = ", range: " + ((tY -(window / 2)) + Integer.parseInt(ans.get(3))) + "-"+((tY -(window / 2)) + Integer.parseInt(ans.get(5)));
			writer.println(">" + ClientGlobals.cutString(sequence2.getName(), FASTA_COLS -rangeString2.length()-1) + rangeString2);
			// iterate over range
			col=0;
			for (int i=0; i<finalstr2.length; i++) {
				writer.print(finalstr2[i]);
				col++;
				if (col % FASTA_COLS == 0) writer.println();
			}
			writer.close();
		

		Arrays.fill(strseq1, ' '); 
		Arrays.fill(strseq2, ' ');
		Arrays.fill(finalstr1, ' '); 
		Arrays.fill(finalstr2, ' ');
	}

	private static char getNuclSimCode(byte map, byte map2, SubstitutionMatrix submat) {
		if (map != -1 && map2 != -1) {
			if (map==map2)
				return ':';
			else {
				int score = submat.getScore(map,map2);
				if (score > 0)
					return '.';
				else
					return ' ';
			}
		}
		else
			return ' ';
	}
} public class SubstitutionMatrix {

	private static final int PROTEIN_THRESHOLD = 20;
	
	private int[][] scoreMatrix;
	private byte[] charMapping;
	private char[] revMapping;
	private boolean bNucleotideMatrix=true;

	private SubstitutionMatrix(int[][] scorematrix, byte[] charMapping, char[] revMapping) {
		this.scoreMatrix = scorematrix;
		this.charMapping = charMapping;
		this.revMapping = revMapping;
	}

	public int getScore(byte a, byte b) {
		return scoreMatrix[a][b];
	}
	
	public int[][] getMatrix() {
		return scoreMatrix;
	}
	
	public byte map(byte input) {
		return charMapping[input];
	}
	
	public char reverseMap(byte input) {
		return revMapping[input];
	}
	
	public int getAlphabetSize() {
		return scoreMatrix[0].length-1;
	}
	
	public boolean isNucleotideMatrix() {
		return bNucleotideMatrix;
	}
	
	public void setNucleotideMatrix(boolean b) {
		bNucleotideMatrix = b;
	}

	public static SubstitutionMatrix loadFromResource(String file)
			throws IOException, InvalidSubMatFileException {
		
		// open file
		BufferedReader br = new BufferedReader(new InputStreamReader(
				SubstitutionMatrix.class.getResourceAsStream(file)));
		
		return doLoad(br, file);
	}
	
	public static SubstitutionMatrix loadFromFile(String file)
			throws IOException, InvalidSubMatFileException {
		
		// open file
		BufferedReader br = new BufferedReader(new FileReader(file));

		return doLoad(br, file);
	}
	
		
	private static SubstitutionMatrix doLoad(BufferedReader br, String file) throws IOException, InvalidSubMatFileException {

		String curLine;
		boolean firstLine = true;
		// create character mapping array
		byte[] revMapping = new byte[256];
		// initialize with -2 -> invalid character
		//Arrays.fill(charMapping,(byte)-2);
		// set CR and LF to ignore -> -1
		revMapping[13] = -1;
		revMapping[10] = -1;
		// backmapping array
		char[] backMapping = new char[256];
		
		int[][] matrix = null;
		int mapPos = 1;
		int linenum=0;
		int reallinenum=0;
		// parse, read line by line
		while ((curLine = br.readLine()) != null) {
			reallinenum++;
			// ignore comments
			if ((curLine.length() > 0 ) && (curLine.charAt(0) != '#')) {
				if (firstLine) {
					// first line containing data -> table heading
					// e.g. A R N D C Q E G H I L K M F P S T W Y V B Z X *
					// create character mapping
					for (int i = 0; i < curLine.length(); i++) {
						if (curLine.charAt(i) != ' ') {
							// throw exception if character already defined
							if (revMapping[curLine.charAt(i)] != 0)
								throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - character defined twice: " + curLine.charAt(i));
							else {
								revMapping[curLine.charAt(i)] = (byte)mapPos++;
								// if this is not the last character a whitespace must follow
								if ( (i<(curLine.length()-1)) && (curLine.charAt(i+1)!=' '))
									throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - invalid expression: " + curLine.substring(i,i+2) );
								
							}
						}						
						
					}
					// create backward mapping 
					for (int i=0; i<256; i++)
						if (revMapping[i] >= 0)
							backMapping[revMapping[i]] = (char)i;
					// create scorematrix array
					matrix = new int[mapPos][mapPos];
					// set flag
					firstLine = false;

				} else {
					// all the other lines contain data seperated by spaces
					// e.g. A 2 -2 0 0 -2 0 0 1 -1 -1 -2 -1 -1 -3 1 1 1 -6 -3 0
					// 0 0 0 -8
					// first token is ignored; the rest is parsed and written
					// into the array
					
					// check if we already have all needed lines
					if (linenum == (mapPos-1))
						throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - too many lines");
					int toknum = 0;
					int poslast = 0;
					boolean lastchar=false;
					boolean foundnonspace = false;
					for (int i = 0; i < curLine.length(); i++) {
						// last character ?
						if (i == curLine.length()-1) {
							if (curLine.charAt(i) != ' ') {
								lastchar=true;
								i++;
							}
						}
						
						if (lastchar || ((curLine.charAt(i) == ' ') && foundnonspace)) {
							if (toknum > 0) {
								// check if there are too many tokens
								if (toknum >= mapPos)
									throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - too many tokens");
								// extract integer value and insert into array
								try {
									matrix[toknum][linenum+1] = Integer.parseInt(curLine.substring(poslast, i).trim());
								} catch (NumberFormatException nfe) {
									throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - invalid integer value: " + curLine.substring(poslast, i).trim());
								}

							}
							else {
								// first token: check if correct character is given (matrix must be symmetric)
								if (revMapping[curLine.trim().charAt(0)] != (linenum+1)) 
									throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - invalid character: " + curLine.trim().charAt(0) + ", exptected: " + backMapping[linenum+1]);
								// a white space must follow now!
								if (curLine.trim().charAt(1) != ' ')
									throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - invalid expression: " + curLine.substring(0,2) );
									

							}

							toknum++;
							poslast = i;
							foundnonspace = false;
						} else if (curLine.charAt(i) != ' ')
							foundnonspace = true;
					}
					// decrease toknum by one if last characters where whitespaces
					if (!foundnonspace)
						toknum--;
					if (toknum > 0) {  // any tokens found or empty line?
						if (toknum < (mapPos-1)) // incomplete line
							throw new InvalidSubMatFileException("Error in " + file + ":" + reallinenum + " - incomplete line");
						else
							linenum++;
					} 
				}
			}
		}
		// check if we are missing lines
		if (linenum < (mapPos-1)) 
			throw new InvalidSubMatFileException("Error in " + file + " - there are lines missing");
		
		
		return new SubstitutionMatrix(matrix,revMapping,backMapping);
	}
	
	public void guessIfSubstitutionMatrix() {
		bNucleotideMatrix = (getAlphabetSize() < PROTEIN_THRESHOLD);
	}
	
} public class SuffixArray {
	
	private int[] suftab;
	private byte[] orgdata;
	private short[] lcp;
	
	private AbortionChecker ac;
	
	boolean isgepardsa=true;
	
	public SuffixArray(Sequence seq, int K, AbortionChecker ac) {
		
		this.ac = ac;
		
		// copy byte array to int array and add padding zeros
		byte[] s = seq.getSequenceData();
		int si[] = new int[s.length+3];
		
		for (int i=0; i<s.length; i++)
			si[i] = s[i];
		
		si[s.length]=si[s.length+1]=si[s.length+2]=0;
//		si[s.length]=si[s.length+1]=si[s.length+2]=K;
			
		// allocate space for suffix array
		suftab = new int[s.length];
		// generate suffix array
		genArray(si,suftab,s.length,K+1);
		
		// check for abortion
		if ((ac != null) && (ac.dotplotAborted())) {
			return;
		}
		
		// generate lcp table
		lcp = genLCP(suftab,s);
		// store original data
		orgdata = s;	
		
		
		
	}
	
	
	private SuffixArray(int[] st, short[] l,  byte[] od, boolean isgepardsa) {
		suftab = st;
		orgdata = od;
		lcp = l;
		this.isgepardsa = isgepardsa;
	}

	// search in 'query' from 'qstart' for 'qlen' bytes in the direction 'backwards?'
	// it uses the 'map' array for identity or complementary character mapping
	public int[] search(byte[] query, int qstart, int qlen, boolean backwards, byte[] map) {
	
		int[] ret = new int[0];
		int r,l,m;
		
		l=0; r=orgdata.length-1;
		m = (l+r)/2;
		
		while (l<=r) {
			m = (l+r)/2;	
			
			// compare current position
			int comp;
//			if (suftab[m] <= orgdata.length-qlen) {
				comp = compare(orgdata,suftab[m],query,qstart,qlen, backwards, map);
//			}
//			else
//				comp = -1; // probably wrong!!!
				
			if (comp == 1)       // orgdata > query at current pos -> search to the left 
				r = m-1;
			else if (comp == -1) // orgdata < query at current pos -> search to the right
				l=m+1;
			else {
				
				
				
				// find left boundary
				l=m;
				while ( (l>0) && (lcp[l-1] >= qlen)) l--;
				// find right boundary
				r=m;
				while ( r<(orgdata.length-1) && (lcp[r] >= qlen)) r++;

				// number of hits:  r-l+1
				// prepare result array
				ret = new int[r-l+1];
				
				// fill result array
				for (int i=0; i<r-l+1; i++) {
					ret[i] = suftab[l+i];
				}
				
				// set r<l to exit the loop
				r=0; l=1;
			} 
		}
			

		return ret;
		
	}
	
	// results: -1 -> s1<s2      0 ->  s1==s2    1 -> s1>s2
	private int compare(byte[] s1, int off1, byte[] s2, int off2, int p, boolean backwards, byte[] map) {
		
		// characters of the second sequence are mapped (to complementary when searching backwards!)
		
		if (!backwards) {
			int i;
			for (i=0; i<p && (i+off1<s1.length) && (i+off2<s2.length) ;i++) {
				if (s1[i+off1] < map[s2[i+off2]])
					return -1;
				else if (s1[i+off1] > map[s2[i+off2]])
					return 1;
			}
			if (i<p) {
				if (isgepardsa) {
					if (i+off1 == s1.length)
						return -1;
					else 
						return 1;
				} else {
					if (i+off1 == s1.length)
						return 1;
					else 
						return -1;
				}
			}
			else
				return 0;
		} else {
			int i;
			for (i=0; i<p && (i+off1<s1.length) && (i+off2<s2.length);i++) {
				if (s1[i+off1] < map[s2[-i+off2]])
					return -1;
				else if (s1[i+off1] > map[s2[-i+off2]])
					return 1;
			}
			if (i<p) {
				if (isgepardsa) {
					if (i+off1 == s1.length)
						return -1;
					else 
						return 1;
				} else {
					if (i+off1 == s1.length)
						return 1;
					else 
						return -1;
				}
			}
			else 
				return 0;
		
		}
		
	/*	if (!backwards) {
			for (int i=0; i<p;i++) {
				if (s1[i+off1] < s2[i+off2])
					return -1;
				else if (s1[i+off1] > s2[i+off2])
					return 1;
			}
			return 0;
		} else {
			for (int i=0; i<p;i++) {
				if (s1[i+off1] < s2[-i+off2])
					return -1;
				else if (s1[i+off1] > s2[-i+off2])
					return 1;
			}
			return 0;
		
		}*/
	}
	

	
	private void radixPass(int[] a, int[] b, int[] s, int offset, int n, int K) {
		
		int[] c = new int[K+1];					// create counter array
		for (int i=0; i<n; i++)	 c[s[a[i]+offset]]++;	// count occurences
		for (int i=0, sum=0; i<=K; i++) {		// calculate first indices in target array
			int t = c[i]; c[i] = sum; sum+=t;	//    for each character from 0..K
		}
		
		for (int i=0; i<n; i++)		// assign new positions in target array (sort)
			b[c[s[a[i]+offset]]++] = a[i];
			
	}	
	
	
	public int[] getRawArray() {
		return suftab;
	}
	
	
	// generate suffix array from input array s[] with key 0..K
	private void genArray(int[] s, int[] SA, int n, int K) {
				
		// check for abortion
		if ((ac != null) && (ac.dotplotAborted())) {
			suftab = null;
			orgdata = null;
			lcp = null;
			return;
		}
		
		// calculate amount of suffixes index mod 3 = 0, = 1, =2
		int n0 = (n+2)/3,	n1 = (n+1)/3,	n2=n/3;
		int n02 = n0 + n2;
		
		
		// generate helper arrays
		int s12a[] = new int[n02+3];
		int s12b[] = new int[n02+3];
	
		// generate indices with mod 3 != 0
		for (int i=0, j=0; i<n+(n0-n1); i++) if (i%3 != 0) s12a[j++] = i;
				
		// radix sort mod 3 != 0 suffixes stored in s12a
		// 3 passes for each triple position
		radixPass(s12a, s12b, s ,2,n02, K);
		radixPass(s12b, s12a, s ,1,n02, K);
		radixPass(s12a, s12b, s ,0,n02, K);
		// (potential) suffix array is now stored in s12b
		
		// find lexographic names
		int name=0;
		int c0=-1,c1=-1,c2=-1;
		for (int i=0; i<n02; i++) {
			// check if current triple is different from predecessor
			if ( (s[s12b[i]] != c0) || (s[s12b[i]+1] != c1) || (s[s12b[i]+2] != c2) ) {
				name++;		// new lexographic name
				c0=s[s12b[i]];	c1=s[s12b[i]+1];	c2=s[s12b[i]+2];	// store last values
			}
			// store lex. names in s12a (mod3=1 -> left half, mod3=2 -> right half)
			if (s12b[i] % 3 == 1)	s12a[s12b[i]/3] = name;			// left half
			else					s12a[s12b[i]/3 + n0] = name;	// right half
		}
		
		// names already unique?
		if (name < n02) { // no -> recurse
			genArray(s12a, s12b, n02, name);
			
			// check for abortion
			if ((ac != null) && (ac.dotplotAborted())) {
				suftab = null;
				orgdata = null;
				lcp = null;
				return;
			}
			
			// create and store unique names
			for (int i=0; i<n02; i++) s12a[s12b[i]] = i+1;
		} else { // yes -> generate suffix array from names
			for (int i=0; i<n02; i++) s12b[s12a[i]-1]=i;
		}
	
		int s0a[] = new int[n0];
		int s0b[] = new int[n0];
		// <DONTKNOW>
		for (int i=0,j=0; i<n02; i++)
			if (s12b[i] < n0) s0a[j++] = 3*s12b[i];
		// </DONTKNOW>
		radixPass(s0a, s0b, s,0, n0, K);
		s0a=null;
		
		// merge arrays now
		
		// t = current pos in SA12 (s12b) array
		// p = current pos in SA0 (s0b) array
		// i = pos of current 12 suffix
		// j = pos of current 0 suffix
		// k = current pos in final suffix array cSA
		
		for (int p=0, t=n0-n1, k=0; k<n; k++) {
			// get real pos of current mod3 =1,=2 suffix
			int i =    (s12b[t]<n0) ? s12b[t] * 3 + 1 : (s12b[t]-n0)*3+2;
			int j = s0b[p];
			
			boolean q=false;
			if (s12b[t] < n0) {
				int a1=s[i],a2=s12a[s12b[t] + n0],b1=s[j],b2=s12a[j/3];
				q=		(a1 < b1 || (a1 == b1 && a2 <= b2));		
			} else {
				int a1=s[i],a2=s[i+1],a3=s12a[s12b[t]-n0+1] ,      b1=s[j],b2=s[j+1],b3=s12a[j/3+n0];
				q = (a1 < b1 ||(a1 == b1 && (a2 < b2 || (a2 == b2 && a3 <= b3))));
			}
			if (q) {
				SA[k] = i; t++;
				if (t == n02) { // only SA0 suffixes left
					for (k++; p<n0; p++,k++) SA[k] = s0b[p];
				}
					
			} else {
				SA[k] = j; p++;
				if (p == n0) { // only SA12 suffixes left
					for (k++; t<n02; t++,k++) SA[k] =    (s12b[t]<n0) ? s12b[t] * 3 + 1 : (s12b[t]-n0)*3+2;
				}
			}
		}
	}

	
	public void saveToFile(String filename) throws IOException {
		// create output stream
		LEDataOutputStream safile =	new LEDataOutputStream(
			new BufferedOutputStream(new GZIPOutputStream(new FileOutputStream(filename))));
		
		// vector to recall big lcps
		Vector<BigLCP> biglcps = new Vector<BigLCP>();
		
		// write type 0 (gepard suffix array)
		safile.write(0);
		
		// write number of array entries
		safile.writeInt(suftab.length);
		// write suffix array
		for (int i=0; i<suftab.length; i++) 
			safile.writeInt(suftab[i]);
		// write lcp table
		for (int i=0; i<lcp.length; i++) {
			// remember big lcp?
			if (lcp[i] < 255) {
				safile.write(unsignedshort2byte(lcp[i]));
			} else {
				biglcps.add(new BigLCP(i, lcp[i]));
				safile.write(unsignedshort2byte((byte)255));
			}
		}
		
		// now write big lcps to file
		for (BigLCP lcp : biglcps) {
			safile.writeInt(lcp.suffix+1);
			safile.writeInt(lcp.lcp);
		}
			
		
		safile.close();
	}
	
	public static SuffixArray loadFromFile(File file, byte[] orgdata) throws FileNotFoundException, IOException {
		
		LEDataInputStream in = new LEDataInputStream(
				new BufferedInputStream(new GZIPInputStream(new FileInputStream(file))));
		
		// read type
		byte type = in.readByte();
		boolean isgepardsa = (type == 0);
		
		// read sequence length
		int seqlength = in.readInt();
		// create arrays
		int[] suftab = new int[seqlength];
		short[] lcp = new short[seqlength-1];
		// read suffix array
		for (int i=0; i<seqlength; i++)
	        suftab[i] = in.readInt();
		// read lcp table
		for (int i=0; i<seqlength-1; i++) 
        	lcp[i] = byte2unsignedshort(in.readByte());
		
		// now load remaining big lcp information
		boolean eof=false;
		while (!eof) {
			try {
				int cursuf = in.readInt();
				int curlcp = in.readInt();
				
				// restrict to maximum positive short value
				if (curlcp > Short.MAX_VALUE) curlcp = Short.MAX_VALUE;

				lcp[cursuf-1] = (short)curlcp;
				
			} catch (EOFException e) {
				eof = true;
			}
		}
        		
		in.close();
		
		
		return new SuffixArray(suftab,lcp,orgdata, isgepardsa);
	       
		   
	}
	
	private static short byte2unsignedshort(byte i) {
		if (i < 0) {
			return (short)(256 + i);
		} else {
			return i;
		}
	}
	
	private static byte unsignedshort2byte(short i) {
		if (i > 127)
			return (byte)(i-256);
		else
			return (byte)i;
	}
	
	private short[] genLCP(int[] suftab, byte[] orgdata) {
		// shortcut
		int len = suftab.length;
		// create lcp array
		short[] ret = new short[len-1];
		// iterate through suffix table
		for (int i=0; i<len-1; i++) {
			// determine LCP of suffix i and i+1
			int  j=0;
			while (j<32768 && (suftab[i]+j < len) && (suftab[i+1]+j <len)
				 && (orgdata[suftab[i]+j] == orgdata[suftab[i+1]+j]) ) {
				j++;
			//	System.out.println(j);
			}
			// store LCP
			ret[i] = (short)j;
			
			
		}
		
		return ret;
	}

	public static String getSAFilename(String sequenceFile) {
		// extract filename
		String filename = extractFilename(sequenceFile);
		// construct suffix array file name using original file size and extension .sa
		return filename + "_" + new File(sequenceFile).length() + ".sa";
	}
	
	private static String extractFilename(final String source) {
		int sep = source.lastIndexOf(System.getProperty("file.separator"));
		if (sep > -1) {
			return source.substring(sep+1);
		} else
			return source;
	}
	
	public void compareSAs(SuffixArray sa) {
		
		int j=0;
		for (int i=0; i<sa.suftab.length; i++) {
			if (this.suftab[i] != sa.suftab[i]) {
				System.out.println("differing at suffix " + i  + ", this: " + this.suftab[i] + ", other: " + sa.suftab[i]);
				j++;
				/*
				for (int k=0; k<=8 && suftab[i]+k < suftab.length; k++)
					System.out.print(mapit[orgdata[suftab[i]+k]]);
				System.out.println();*/
				System.out.println();
				dump("vmatch  ", orgdata, this.suftab[i]);
				dump("gepard  ", orgdata, sa.suftab[i]);
				
			}
			if (j >= 10)
			break;
		}
		/*
		j=0;
		for (int i=0; i<sa.lcp.length; i++) {
			if (this.lcp[i]!= sa.lcp[i]) {
				System.out.println("differing at lcp " + i  + ", this: " + this.lcp[i] + ", other: " + sa.lcp[i]);
			}
		if (j > 100)
				break;
		}*/
	}
	
	private static void dump(String name ,  byte[] orgdata, int pos) {
		char mapit[] = {' ', 'A', 'T', 'G' , 'C', 'N'};
		
		System.out.print(name);
		for (int k=0; k<50 && pos+k < orgdata.length; k++)
			System.out.print(mapit[orgdata[pos+k]]);
		System.out.println();
	}
	
	public void printYourself() {
		char mapit[] = {' ', 'A', 'T', 'G' , 'C'};
		
		for (int m=0; m<suftab.length; m++) {
			System.out.print(m + ": ");
			for (int i=0; i<=8 && suftab[m]+i < suftab.length; i++)
				System.out.print(mapit[orgdata[suftab[m]+i]]);
			System.out.println();
		}
	}
	
	private static class BigLCP {
		
		public int suffix;
		public int lcp;
		
		public BigLCP(int suffix, int lcp) {
			this.suffix = suffix;
			this.lcp = lcp;
		}
		
	}
	
	public void verifyArray() {
		
		byte[] idmap = new byte[128];
		for (byte i=0; i<128 && i >= 0; i++) idmap[i] = i;
		
		for (int i=0; i<suftab.length-1; i++) {
			
			if (compare(orgdata, suftab[i], orgdata, suftab[i+1], 1000, false, idmap) > 0) {
				System.out.println("WAAAAA, bad one: " + suftab[i+1]);
				dump(i+ "\t", orgdata, suftab[i]);
				dump((i+1)+"\t", orgdata, suftab[i+1]);
//				System.exit(1);
			}
				
		//	System.out.println();
			
		}
	}
	
	public void showFirst(int n) {
		for (int i=0; i<n; i++)
			dump("", orgdata, suftab[i]);
	}
	
	public int getLength() {
		return suftab.length;
	}
	
	public void doubles() {
		byte[] check = new byte[suftab.length];
		
		System.out.println("seq: " + suftab.length);
		
		for (int i=0; i<check.length; i++) {
			
			if (suftab[i] >= check.length)
				System.out.println("Was: " + i);
			
			if (check[suftab[i]] == 1) 
				System.err.println("NOOOOOOOOOOOOOO");
			else
				check[suftab[i]] = 1;
		}
	}
	
}
import java.awt.Dimension;
import java.io.File;
import java.io.IOException;
import org.gepard.client.ClientGlobals;
import org.gepard.client.DrawTarget;
import org.gepard.client.Plotter;
import org.gepard.common.AbortionChecker;
import org.gepard.common.DotMatrix;
import org.gepard.common.DotMatrixCallback;
import org.gepard.common.FASTAReader;
import org.gepard.common.InvalidFASTAFileException;
import org.gepard.common.InvalidSubMatFileException;
import org.gepard.common.Needleexport;
import org.gepard.common.ParameterSet;
import org.gepard.common.Sequence;
import org.gepard.common.Smithexport;
import org.gepard.common.SubstitutionMatrix;
import org.gepard.common.SuffixArray;
		Sequence seq2 = null;
		try {
			seq2 = FASTAReader.readFile(sequence2, submat);
		} catch (IOException e) {
			System.err.println("Error: Could not load second sequence file:");
			System.err.println(e.getMessage());
			System.exit(1);
		} catch (InvalidFASTAFileException e) {
			System.err.println("Error: File could not be recognized as FASTA format file:");
			System.err.println(sequence2);
			System.err.println("Problem: " + e.getMessage());
			System.exit(1);
		}
		// complementary
		if (arguments.isSet("secondcomp"))
			Sequence.complementarizeSequence(seq2);

		// dimension and dotplot parameters
		int word = -1, window = -1;
		try {

			// load dotplot calculation parameters
			if (arguments.isSet("word"))
				word = Integer.parseInt(arguments.getValue("word"));
			else if (arguments.isSet("window"))
				word = 0;
			else
				word = DEF_WORD;
			window = (arguments.isSet("window")) ? Integer.parseInt(arguments.getValue("window")) : DEF_WINDOW;
		} catch (NumberFormatException e) {
			System.err.println("Error: Invalid integer number:");
			System.err.println(e.getMessage());
			System.exit(1);
		}

		// display parameters
		// lower
		float lower = -1;
		if (arguments.isSet("lower"))
			lower = Float.parseFloat(arguments.getValue("lower")) / 100f;
		else {
			if (word > 0)
				lower = 0f;
			else
				lower = Float.NaN; // will be set later when we have dotplot
									// scores
		}
		// upper
		float upper = -1;
		if (arguments.isSet("upper")) {
			upper = Float.parseFloat(arguments.getValue("upper")) / 100f;
		} else {
			if (word > 0)
				upper = 0.2f;
			else
				upper = 1.0f;
		}
		// greyscale
		float greyscale = -1;
		if (arguments.isSet("greyscale")) {
			greyscale = Float.parseFloat(arguments.getValue("greyscale")) / 100f;
		} else {
			if (word > 0)
				greyscale = 0f;
			else
				greyscale = 0.5f;
		}

		// get parameter set
		ParameterSet params = new ParameterSet(-1, word, window, -1, -1, -1, -1, ClientGlobals.extractFilename(arguments.getValue("matrix").toLowerCase()));
		// set correct coordinates
		try {
			setCoordinates(params, seq1, seq2, arguments);

			// plot size
			// zoom specified from user?
			if (arguments.isSet("zoom")) {
				// yes
				int zoom = Integer.parseInt(arguments.getValue("zoom"));
				width = XMARGIN + (int) Math.ceil((float) seq1.getLength() / (float) zoom);
				height = YMARGIN + (int) Math.ceil((float) seq2.getLength() / (float) zoom);

				/*
				 * // plot too small? if (width-XMARGIN <= 1 || height-YMARGIN
				 * <= 1) {
				 * 
				 * }
				 */

				params.zoom = zoom;

			} else {
				// get size from user of def values
				width = (arguments.isSet("maxwidth")) ? Integer.parseInt(arguments.getValue("maxwidth")) : DEF_WIDTH;
				height = (arguments.isSet("maxheight")) ? Integer.parseInt(arguments.getValue("maxheight")) : DEF_HEIGHT;

				// calc zoom factor
				params.zoom = calcZoom(width, height, params.seq1Start, params.seq1Stop, params.seq2Start, params.seq2Stop);
			}

			// plot too small?
			if (width - XMARGIN <= 1 || height - YMARGIN <= 1) {
				System.err.println("Error: This plot would become too small.");
				System.err.println("Please use bigger plot size or smaller zoom factor.");
				System.exit(1);
			}

			// invalid values?
			if (params.seq1Start >= params.seq1Stop || params.seq2Start >= params.seq2Stop) {
				System.err.println("Error: To coordinates may not be greater than from coordinates!");
				System.exit(1);
			}
		} catch (NumberFormatException e) {
			System.err.println("Error: Invalid integer number:");
			System.err.println(e.getMessage());
			System.exit(1);
		}

		// get suffix array
		boolean saSecondSeq = false;
		SuffixArray sa = null;