package ErrorCorrection;
import java.util.*;
import java.io.*;
import java.math.BigDecimal;

import static ErrorCorrection.Poisson.*;
import static ErrorCorrection.Exponent.*;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.biojava.nbio.core.exceptions.CompoundNotFoundException;

public class DataSet {
	int freqThr;
	int lengthThr;
	int additionalfreqThrMult;
	int maxAllErrorsPerc;
	public ArrayList<Read> reads;
	ArrayList<Read> badReads;
	ArrayList<Read> goodReads;
	ArrayList<Kmer_general> allKmers;
	ArrayList<Haplotype> haplotypes;
	int k;
	String file_name;
        String file_name_short;
	int[] valuesDistribution;
	int[] valuesDistribution1;
        int finderrorsseglen;
        int avProc;
	public DataSet (String addr) throws IOException
	{
		lengthThr = 100;
		additionalfreqThrMult = 10;
		maxAllErrorsPerc = 40;
                freqThr = -1;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = addr;
                finderrorsseglen = k;
		FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
		
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();

		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		
		HashMap<String,Integer> seq = new HashMap<String, Integer>();
                
		
		try
		{
		while (s!= null)
		{
			System.out.println("Reading " + name);
			
			if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
				//	reads.add(new Read(nucl, name));
				
//					System.out.println(name);
				
					if (seq.containsKey(nucl))
						seq.put(nucl, seq.get(nucl)+1);
					else
						seq.put(nucl, 1);
                                        
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
		}
		catch (Exception e)
		{
			e.printStackTrace();
		}
	//	reads.add(new Read(nucl, name));
		
		if (seq.containsKey(nucl))
			seq.put(nucl, seq.get(nucl)+1);
		else
			seq.put(nucl, 1);		
		Iterator it = seq.entrySet().iterator();
		int count = 1;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucleo = (String) me.getKey(); 
			int f = (Integer) me.getValue();
			String nm = "read" + count;
			if (nucleo.length() > lengthThr)
				reads.add(new Read(nucleo,nm,f));
			count++;
		}
		
		fr.close();
	}
        public DataSet (String addr, int lt) throws IOException
	{
		lengthThr = lt; // 50
		additionalfreqThrMult = 10;
		maxAllErrorsPerc = 40;
                freqThr = -1;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = addr;
                finderrorsseglen = k;
		FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);

                String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();


		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();

		HashMap<String,Integer> seq = new HashMap<String, Integer>();

		try
		{
		while (s!= null)
		{
//			System.out.println("Reading " + name);

			if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
				//	reads.add(new Read(nucl, name));

//					System.out.println(name);

					if (seq.containsKey(nucl))
						seq.put(nucl, seq.get(nucl)+1);
					else
						seq.put(nucl, 1);

					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
		}
		catch (Exception e)
		{
			e.printStackTrace();
		}
	//	reads.add(new Read(nucl, name));

		if (seq.containsKey(nucl))
			seq.put(nucl, seq.get(nucl)+1);
		else
			seq.put(nucl, 1);
		Iterator it = seq.entrySet().iterator();
		int count = 1;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucleo = (String) me.getKey();
			int f = (Integer) me.getValue();
			String nm = "read" + count;
			if (nucleo.length() > lengthThr)
				reads.add(new Read(nucleo,nm,f));
			count++;
		}

		fr.close();
	}
        public DataSet (String addr, int minlen, int maxlen) throws IOException
	{
		additionalfreqThrMult = 10;
		maxAllErrorsPerc = 40;
                freqThr = -1;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = addr;
                finderrorsseglen = k;
		FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
		
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();

		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		
		HashMap<String,Integer> seq = new HashMap<String, Integer>();
                
		
		try
		{
		while (s!= null)
		{
//			System.out.println("Reading " + name);
			
			if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
				//	reads.add(new Read(nucl, name));
				
//					System.out.println(name);
				
					if (seq.containsKey(nucl))
						seq.put(nucl, seq.get(nucl)+1);
					else
						seq.put(nucl, 1);
                                        
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
		}
		catch (Exception e)
		{
			e.printStackTrace();
		}
	//	reads.add(new Read(nucl, name));
		
		if (seq.containsKey(nucl))
			seq.put(nucl, seq.get(nucl)+1);
		else
			seq.put(nucl, 1);		
		Iterator it = seq.entrySet().iterator();
		int count = 1;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucleo = (String) me.getKey(); 
			int f = (Integer) me.getValue();
			String nm = "read" + count;
			if ((nucleo.length() >= minlen)&&(nucleo.length() <= maxlen))
				reads.add(new Read(nucleo,nm,f));
			count++;
		}
		
		fr.close();
	}
	public DataSet (String addr, char c) throws IOException
	{
            // keep read name
		lengthThr = 50;
		additionalfreqThrMult = 3;
		maxAllErrorsPerc = 40;
                finderrorsseglen = k;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = addr;
		FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
		
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();
		String nucl = new String();
		String name = new String();
		name = "";
                int count = 1;
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		while (s!= null)
		{
                        if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
					reads.add(new Read(nucl, name));
//                                        System.out.println("Read " + count + " added");
                                        count++;
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
		reads.add(new Read(nucl, name));
		fr.close();
	}
        public DataSet (String addr, char c, int minlen) throws IOException
	{
            // keep read name
		lengthThr = minlen;
		additionalfreqThrMult = 3;
		maxAllErrorsPerc = 40;
                finderrorsseglen = k;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = addr;
		FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
		
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();
		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		while (s!= null)
		{
                        if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
					reads.add(new Read(nucl, name));	
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
		reads.add(new Read(nucl, name));
		fr.close();
	}
	public DataSet(DataSet ds)
	{
		freqThr = ds.freqThr;
		lengthThr = ds.lengthThr;
                this.avProc = ds.avProc;
		additionalfreqThrMult = ds.additionalfreqThrMult;
		maxAllErrorsPerc = ds.maxAllErrorsPerc;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = ds.file_name;
                file_name_short = ds.file_name_short;
		k = ds.k;
                finderrorsseglen =ds.finderrorsseglen;
		
		HashMap<String,Integer> seq = new HashMap<String, Integer>();
                HashMap<String,String> genot = new HashMap<String, String>();
                
                int iread = 0;
		
		for(Read r : ds.reads)
		{
                        iread++;
			String nucl = "";
			for (int i = 0; i < r.getLength(); i++)
				if (r.corrections[i] != null)
				{
					if (r.corrections[i].type == 'R')
						nucl = nucl + r.corrections[i].replacement;
					if (r.corrections[i].type == 'I')
						nucl = nucl + r.nucl.charAt(i) + r.corrections[i].replacement;
				}
				else
					nucl = nucl + r.nucl.charAt(i);
                        String nucl1 = "";
                        for (int i = 0; i < nucl.length(); i++)
                            if (nucl.charAt(i)!='-')
                                nucl1 = nucl1 + nucl.charAt(i);
                        nucl = nucl1;
			if (nucl.length() > lengthThr)
                        {
/*                            Corrector cr = new Corrector();
                            System.out.println("Checking consistency: read" + iread + "//" + ds.reads.size());
                            if (!cr.checkConsistencyKmers(nucl, ds))
                                continue; */
                            
				if (seq.containsKey(nucl))
                                {
					seq.put(nucl, seq.get(nucl)+r.frequency);
                                        genot.put(nucl, r.genotype);
                                }
				else
                                {
					seq.put(nucl, r.frequency);
                                        genot.put(nucl, r.genotype);
                                }
                        }
		}
		Iterator it = seq.entrySet().iterator();
		int count = 1;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucleo = (String) me.getKey(); 
			int f = (Integer) me.getValue();
                        String gen = (String) genot.get(nucleo);
			String nm = "read" + count;
			if (nucleo.length() > lengthThr)
				reads.add(new Read(nucleo,nm,f,gen));
			count++;
		}
	}
        public DataSet(DataSet ds, boolean noGenotyping)
	{
		freqThr = ds.freqThr;
		lengthThr = ds.lengthThr;
                this.avProc = ds.avProc;
		additionalfreqThrMult = ds.additionalfreqThrMult;
		maxAllErrorsPerc = ds.maxAllErrorsPerc;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = ds.file_name;
                file_name_short = ds.file_name_short;
		k = ds.k;
                finderrorsseglen =ds.finderrorsseglen;
		
		HashMap<String,Integer> seq = new HashMap<String, Integer>();
                
                int iread = 0;
		
		for(Read r : ds.reads)
		{
                        iread++;
			String nucl = "";
			for (int i = 0; i < r.getLength(); i++)
				if (r.corrections[i] != null)
				{
					if (r.corrections[i].type == 'R')
						nucl = nucl + r.corrections[i].replacement;
					if (r.corrections[i].type == 'I')
						nucl = nucl + r.nucl.charAt(i) + r.corrections[i].replacement;
				}
				else
					nucl = nucl + r.nucl.charAt(i);
                        String nucl1 = "";
                        for (int i = 0; i < nucl.length(); i++)
                            if (nucl.charAt(i)!='-')
                                nucl1 = nucl1 + nucl.charAt(i);
                        nucl = nucl1;
			if (nucl.length() > lengthThr)
                        {
/*                            Corrector cr = new Corrector();
                            System.out.println("Checking consistency: read" + iread + "//" + ds.reads.size());
                            if (!cr.checkConsistencyKmers(nucl, ds))
                                continue; */
                            
				if (seq.containsKey(nucl))
					seq.put(nucl, seq.get(nucl)+r.frequency);
				else
					seq.put(nucl, r.frequency);	
                        }
		}
		Iterator it = seq.entrySet().iterator();
		int count = 1;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucleo = (String) me.getKey(); 
			int f = (Integer) me.getValue();
			String nm = "read" + count;
			if (nucleo.length() > lengthThr)
				reads.add(new Read(nucleo,nm,f));
			count++;
		}
	}
        public DataSet(DataSet ds, String idmethod)
	{
		freqThr = ds.freqThr;
		lengthThr = ds.lengthThr;
		additionalfreqThrMult = ds.additionalfreqThrMult;
		maxAllErrorsPerc = ds.maxAllErrorsPerc;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
		file_name = ds.file_name;
                file_name_short = ds.file_name_short;
		k = ds.k;
                this.avProc = ds.avProc;
                finderrorsseglen =ds.finderrorsseglen;
                
                String tag = "";
                if (idmethod.equalsIgnoreCase("ET"))
                {
                    StringTokenizer st = new StringTokenizer(ds.reads.get(0).name,"_");
                    tag = st.nextToken() + "_" + st.nextToken();
                }
		
		HashMap<String,Integer> seq = new HashMap<String, Integer>();
		
		for(Read r : ds.reads)
		{
			String nucl = "";
			for (int i = 0; i < r.getLength(); i++)
				if (r.corrections[i] != null)
				{
					if (r.corrections[i].type == 'R')
						nucl = nucl + r.corrections[i].replacement;
					if (r.corrections[i].type == 'I')
						nucl = nucl + r.nucl.charAt(i) + r.corrections[i].replacement;
				}
				else
					nucl = nucl + r.nucl.charAt(i);
                        String nucl1 = "";
                        for (int i = 0; i < nucl.length(); i++)
                            if (nucl.charAt(i)!='-')
                                nucl1 = nucl1 + nucl.charAt(i);
                        nucl = nucl1;
			if (nucl.length() > lengthThr)
				if (seq.containsKey(nucl))
					seq.put(nucl, seq.get(nucl)+r.frequency);
				else
					seq.put(nucl, r.frequency);	
		}
		Iterator it = seq.entrySet().iterator();
		int count = 1;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucleo = (String) me.getKey(); 
			int f = (Integer) me.getValue();
                        String nm = "";
                        if (idmethod.equalsIgnoreCase("ET"))
                            nm = tag + "_" + count + "_" + f;
                        else
                            nm = "read" + count;
			if (nucleo.length() > lengthThr)
				reads.add(new Read(nucleo,nm,f));
			count++;
		}
	}
        public DataSet(String addr, String idmethod) throws FileNotFoundException, IOException
        {            
                int ntokenfreqFromEnd = 1;
                this.file_name = addr;
                
                reads = new ArrayList<Read>();
		haplotypes = new ArrayList<Haplotype>();
                
                FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
		
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();
		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		while (s!= null)
		{
                        if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
//                                        System.out.println(name);
                                        int f = 0;
                                        if (idmethod.equalsIgnoreCase("ET"))
                                        {
                                            StringTokenizer st = new StringTokenizer(name,"_");
                                            int nnn = st.countTokens();
                                            for (int tok = 1; tok <= nnn-ntokenfreqFromEnd; tok++)
                                                st.nextToken();
                                            f = (int) (Double.parseDouble(st.nextToken()));
//                                            System.out.println();
                                        }
                                        
                                        if (idmethod.equalsIgnoreCase("KEC"))
                                        {
                                            StringTokenizer st = new StringTokenizer(name,"|");
                                            st.nextToken();
                                            f = (int) (10000*Double.parseDouble(st.nextToken()));
                                        }
                                        if (idmethod.equalsIgnoreCase("EPLD"))
                                        {
                                            StringTokenizer st = new StringTokenizer(name,"-");
                                            if (st.countTokens() == 2)
                                                f = 1;
                                            if (st.countTokens() > 2)
                                            {
 //                                               System.out.println(st.countTokens());
                                                int ct = st.countTokens();
                                                for (int tok = 1; tok <= ct-1; tok++)
                                                    st.nextToken();
                                                f = (int) (Integer.parseInt(st.nextToken()));
                                            }
                                        }
                                                    
					reads.add(new Read(nucl, name,f));	
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
                int f = 0;
                if (idmethod.equalsIgnoreCase("ET"))
                {
                    StringTokenizer st = new StringTokenizer(name,"_");
                    int nnn = st.countTokens();
                    for (int tok = 1; tok <= nnn-ntokenfreqFromEnd; tok++)
                        st.nextToken();
                    f = (int) (Double.parseDouble(st.nextToken()));
                }                
                if (idmethod.equalsIgnoreCase("KEC"))
                {
                    StringTokenizer st = new StringTokenizer(name,"|");
                    st.nextToken();
                    f = (int) (10000*Double.parseDouble(st.nextToken()));
                }
                if (idmethod.equalsIgnoreCase("EPLD"))
                {
                    StringTokenizer st = new StringTokenizer(name,"-");
                    if (st.countTokens() == 2)
                        f = 1;
                    if (st.countTokens() > 2)
                    {
                        for (int tok = 1; tok <= st.countTokens()-1; tok++)
                            st.nextToken();
                        f = (int) (Integer.parseInt(st.nextToken()));
                    }
                 }
                
		reads.add(new Read(nucl, name,f));
		fr.close();           
        }
        public DataSet(String addr, String idmethod, int ntokenfreqFromEnd) throws FileNotFoundException, IOException
        {            
                reads = new ArrayList<Read>();
		haplotypes = new ArrayList<Haplotype>();
                
                FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
		
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();
		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		while (s!= null)
		{
                        if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
			{
                                        int f = 0;
                                        if (idmethod.equalsIgnoreCase("ET"))
                                        {
                                            StringTokenizer st = new StringTokenizer(name,"_");
                                            int nnn = st.countTokens();
                                            for (int tok = 1; tok <= nnn-ntokenfreqFromEnd; tok++)
                                                st.nextToken();
                                            f = (int) (Double.parseDouble(st.nextToken()));
//                                            System.out.println();
                                        }
                                        
                                        if (idmethod.equalsIgnoreCase("KEC"))
                                        {
                                            StringTokenizer st = new StringTokenizer(name,"|");
                                            st.nextToken();
                                            f = (int) (10000*Double.parseDouble(st.nextToken()));
                                        }
                                                    
					reads.add(new Read(nucl, name,f));	
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
                int f = 0;
                if (idmethod.equalsIgnoreCase("ET"))
                {
                    StringTokenizer st = new StringTokenizer(name,"_");
                    int nnn = st.countTokens();
                    for (int tok = 1; tok <= nnn-ntokenfreqFromEnd; tok++)
                        st.nextToken();
                    f = (int) (Double.parseDouble(st.nextToken()));
                }                
                if (idmethod.equalsIgnoreCase("KEC"))
                {
                    StringTokenizer st = new StringTokenizer(name,"|");
                    st.nextToken();
                    f = (int) (10000*Double.parseDouble(st.nextToken()));
                }
                
		reads.add(new Read(nucl, name,f));
		fr.close();           
        }
        public DataSet(DataSet ds, int i)
        {
                reads = new ArrayList();
                lengthThr = 50;
                additionalfreqThrMult = 10;
                maxAllErrorsPerc = 40;
                this.avProc = ds.avProc;
                freqThr = -1;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
                finderrorsseglen = k;
                for (Read r : ds.reads)
                    this.reads.add(r);
        }
        public DataSet()
        {
                reads = new ArrayList();
                lengthThr = 50;
                additionalfreqThrMult = 10;
                maxAllErrorsPerc = 40;
                freqThr = -1;
		reads = new ArrayList<Read>();
		allKmers = new ArrayList<Kmer_general>();
		haplotypes = new ArrayList<Haplotype>();
                finderrorsseglen = k;
        }
	public void setK(int i)
	{
		k=i;
	}
	void setFreqThr(int i)
	{
		freqThr = i;
	}
	public void setLenThr(int i)
	{
		lengthThr = i;
	}
	void setMaxAllErrorsPerc(int i)
	{
		maxAllErrorsPerc = i;
	}
        public void setFileName(String s)
        {
            this.file_name = s;
        }
        public void setFileNameShort(String s)
        {
            this.file_name_short = s;
        }
        public void setAvProc(int i)
        {
            this.avProc = i;
        }
	public void calculateKMersAndKCounts()
	{
//                allKmers = new  ArrayList<Kmer_general>(); 
		HashMap<String, ArrayList<Kmer_occur>> kmers= new HashMap<String, ArrayList<Kmer_occur>>();
		for (int j = 0; j < reads.size(); j++)
		{
                        System.out.println("Calculate kmers: " + j + "//" + reads.size());
			for (int i = 0; i< reads.get(j).nucl.length()-k+1; i++)
			{
				String s = reads.get(j).nucl.substring(i, i+k);
				if (kmers.containsKey(s))
					((ArrayList<Kmer_occur>) kmers.get(s)).add(new Kmer_occur(reads.get(j),j,i,k));
				else
				{
					ArrayList<Kmer_occur> a = new ArrayList<Kmer_occur>();
					a.add(new Kmer_occur(reads.get(j),j,i,k));
					kmers.put(s, a);
				}
			}
		}
		Iterator it = kmers.entrySet().iterator();
		int maxKCount = 0;
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			String nucl =(String) me.getKey();
			ArrayList<Kmer_occur> readsOfKmer = ((ArrayList<Kmer_occur>) me.getValue());
			int kCount = 0;
			int kCount1 = 0;
//			HashSet<String> readsWhichContKmer = new HashSet<String>();
			for (int i = 0; i < readsOfKmer.size(); i++)
			{
				kCount+= readsOfKmer.get(i).r.frequency;
/*				if (!readsWhichContKmer.contains(readsOfKmer.get(i).r.nucl))
				{
					kCount1+=readsOfKmer.get(i).r.frequency;
					readsWhichContKmer.add(readsOfKmer.get(i).r.nucl);
				}*/
			}
//			int kCount = readsOfKmer.size();
			if (kCount > maxKCount)
				maxKCount = kCount;
			allKmers.add(new Kmer_general(nucl,kCount));
			for (int i = 0; i < readsOfKmer.size(); i++)
			{
				Kmer km = new Kmer(nucl,readsOfKmer.get(i).stPos);
				km.setKCount(kCount);
//				km.setKCount1(kCount1);
				readsOfKmer.get(i).r.addKmer(km);
			}
		}
			
		
		valuesDistribution = new int[maxKCount+1];
		it = kmers.entrySet().iterator();
		while (it.hasNext())
		{
			Map.Entry me = (Map.Entry)it.next();
			ArrayList<Kmer_occur> readsOfKmer = ((ArrayList<Kmer_occur>) me.getValue());
			int kCount = 0;
			for (int i = 0; i < readsOfKmer.size(); i++)
				kCount+= readsOfKmer.get(i).r.frequency;
			valuesDistribution[kCount]++;
		}
		
		valuesDistribution1 = new int[maxKCount+1];
		for (Read r : reads)
			for(Kmer km : r.kmers)
                        {
                            if (km.kCount > maxKCount)
                                DynamicOut.printStep(r.name + " " + km.StartPos);
                            valuesDistribution1[km.kCount]+=r.frequency;
                        }
		DynamicOut.finishSteps();
	}
	
	void clustering() throws IOException
        {
		 int count = 0;
                 int count1 = 0;
		 for (Read r : reads)
		 {
			 count++;
                         count1++;
                         if (r.nucl.length() < this.k)
                         {
                             System.out.println("The value of k is higher than the length of some reads. Choose smaller k!");
                             System.exit(0);
                         }
                         System.out.println("Clustering: read " + count + "/" + reads.size());
			 File f = new File("kmer" + count + "data.txt");
			 FileWriter fw = new FileWriter(f);
//			 fw.write(r.kmers.size() + " " + 2 + "\n");
//			 for (Kmer km : r.kmers)
//				 fw.write(km.StartPos+ " " + km.getCount() + "\n");
			 fw.write(r.kmers.size() + " " + 1 + "\n");
			 for (Kmer km : r.kmers)
				 fw.write(km.getCount() + "\n");
			 fw.close();
			 
			 Runtime run=Runtime.getRuntime(); 
			 Process p=null;
                         if (count1 == 3000)
                         {
                            System.gc();
                            count1 = 0;
                         }
                                                 
                         String exPath = "fams"+File.separator+"fams.exe 0 0 200 " + "kmer" + count + "data ."+ File.separator;
			 p=run.exec(exPath); 
//			 p=run.exec("fams 0 0 200 " + "kmer" + count + "data .\\"); 
			 try {
				p.waitFor();
			} catch (InterruptedException e) {
				e.printStackTrace();
			}
			
			 if (p.exitValue()!= 0)
			 {
				 System.out.println("Error in clustering program");
				 continue;
			 }
			 
			 File fout = new File("out_" + "kmer" + count + "data.txt");
			 
		/*	 while (!fout.exists())
			 {
				 System.out.println("Waiting when process finish");
			 }*/
			 
			 FileReader fr = new FileReader(fout);
			 BufferedReader br = new BufferedReader(fr);
			 
			 String s = br.readLine();
			 int i = 0;
			 while(s!= null)
			 {
		/*		 StringTokenizer st = new StringTokenizer(s, " \t\n");
				 st.nextToken();
				 r.kmers.get(i).setClusterParam(Double.parseDouble(st.nextToken()));*/
				 r.kmers.get(i).setClusterParam(Double.parseDouble(s));
				 i++;
				 s = br.readLine();
			 }
                         br.close();
                         fr.close();
			 
			 f.delete();
                         f = new File("pilot_200_" + "kmer" + count + "data.txt");
                         f.delete();
                         f = new File("modes_" + "kmer" + count + "data.txt");
                         f.delete();
		 }
                 DynamicOut.finishSteps();
	}
        
        void clusteringParallel() throws IOException, InterruptedException, ExecutionException
        {
		 int count = 0;
                 int count1 = 0;
                 List< Future > futuresList = new ArrayList< Future >();
//                 int nrOfProcessors = Runtime.getRuntime().availableProcessors();
                 int nrOfProcessors = this.avProc;
                 ExecutorService eservice = Executors.newFixedThreadPool(nrOfProcessors);
                 
		 for (Read r : reads)
		 {
			 count++;
                         count1++;
                         if (r.nucl.length() < this.k)
                         {
                             System.out.println("The value of k is higher than the length of some reads. Choose smaller k!");
                             System.exit(0);
                         }
                         ClustTask ct = new ClustTask(r,count);
                        futuresList.add(eservice.submit(ct));
			 
		 }
                  Object taskResult;
                  for(Future future:futuresList) 
                  {
                        taskResult = future.get();
                        System.out.println("Clustered: read "+taskResult + "/" + reads.size());
                  }
                 DynamicOut.finishSteps();
	}

        void findFreqThresholdPoisson() throws IOException
        {
            double P = 1;
            // put value of lamba into Poisson class (static for simplicity)
                        int r = 1;
            while ((r < valuesDistribution1.length) && (valuesDistribution1[r] != 0))
                r++;
            System.out.println("r= " + r);
            double[] valuesDistributionS = new double[r-1];
            for (int i = 0; i < r-1; i++)
                valuesDistributionS[i] = valuesDistribution1[i+1];
           
            
            estimateLambdaPoisson(valuesDistributionS, 1);
            int n = 1;
            while (P > 0.01) {
                try {
                    P = 1 - pPoiss(n);
                } catch (Exception ex) {
                    Logger.getLogger(DataSet.class.getName()).log(Level.SEVERE, "", ex);
                }
                n++;
            }
            freqThr = n;

            System.out.println("Threshold =  " + n);
        }

         void findFreqThresholdExponent()
        {
            double P = 1;
            // put value of lamba into Exponent class (static for simplicity)
            int r = 1;
            while (valuesDistribution1[r] != 0)
                r++;
            double[] valuesDistributionS = new double[r-1];
            for (int i = 0; i < r-1; i++)
                valuesDistributionS[i] = valuesDistribution1[i+1];
            estimateLambdaExponent(valuesDistributionS, 1);
            int n = 1;
            while (P > 0.01) {
                try {
                    P = 1 - pExpon(n);
                } catch (Exception ex) {
                    Logger.getLogger(DataSet.class.getName()).log(Level.SEVERE, "", ex);
                }
                n++;
            }
            freqThr = n;

            System.out.println("Threshold =  " + n);
        }
        void findFreqThreshold()
        {

                if (finderrorsseglen == 0)
                {
                    int i = 1;
                    while ((valuesDistribution[i] != 0)&&(i < valuesDistribution.length))
                        i++;
                    if (i >= valuesDistribution.length)
                    {
                        System.out.println("No such number of consequtive zeros. Please specify a smaller number");
                        System.exit(1);
                    }
                    freqThr = i;
                    return;
                }

                double[] variations = new double[valuesDistribution.length+1];
                int varseglen = finderrorsseglen;
                
                if (valuesDistribution.length - varseglen <= 0)
                    {
                        System.out.println("No such number of consequtive zeros. Please specify a smaller number");
                        System.exit(1);
                    }
/*                for (int i = 0; i < valuesDistribution.length; i++)
                    System.out.println(i + " " + valuesDistribution[i]);*/
                for (int i = 1; i < valuesDistribution.length - varseglen; i++)
                {
                    double cnt = 0;
                    for (int j = i+1; j <= i + varseglen; j++)
                        cnt +=Math.abs(valuesDistribution1[j] - valuesDistribution1[j-1]);
                    variations[i] = cnt;
                }
                int i = 1;
                while ((variations[i]!= 0)&&(i < variations.length))
                    i++;
                if (i >= variations.length)
                    {
                        System.out.println("No such number of consequtive zeros. Please specify a smaller number");
                        System.exit(1);
                    }
                while ((variations[i]== 0)&&(i < variations.length))
                    i++;
                if (i >= variations.length)
                    {
                        System.out.println("No such number of consequtive zeros. Please specify a smaller number");
                        System.exit(1);
                    }
                freqThr = i+varseglen-1;
                System.out.println();
        }

	public void PrintReadsStat(String dir) throws IOException
	{
		File f = new File(dir);
		if (!f.exists())
			f.mkdir();
		int count = 1;
		int count1 = 0;
		String fol_sign = "//";
/*		for (Read r : reads)
		{
			System.out.println("Writing results: read " + count + "/" + reads.size()+"\n");	
			count++;
			count1+=r.frequency;
			String name;
			if (badReads != null)
			{
				if (badReads.contains(r))
					name = dir + fol_sign + "kmers" +"(k=" + k + ")" + "(" + r.name+")"+"bad"+"("+r.frequency+")"+ ".txt";
				else
					name = dir + fol_sign + "kmers" +"(k=" + k + ")" + "(" + r.name+")"+"good"+"("+r.frequency+")"+ ".txt";
			}
			else
				name = dir + fol_sign + "kmers" +"(k=" + k + ")" + "(" + r.name+")"+"("+r.frequency+")"+ ".txt";
			FileWriter fw = new FileWriter(name);
//			fw.write(r.kmers.size() + " " + 2 + "\n");
			int maxKcount = 0;
			for (Kmer km : r.kmers)
			{
			//	fw.write(km.nucl + " " + km.StartPos+ " " + km.getCount() + "\n");
				fw.write(km.StartPos+ " " + km.getCount() + " " + r.getNucl(km.StartPos) + " " + km.clusterParam1 + " " + r.NuclType[km.StartPos] + "\n");
				if (km.getCount() > maxKcount)
					maxKcount = km.getCount();
			}
			for (int i = r.getLength()-k+1; i < r.getLength(); i++)
				fw.write(i+ " " + "N/A" + " " + r.getNucl(i) + " " + "N/A" + " " + "N/A" + "\n");
		
			fw.close();
		}
		System.out.println("Total number of reads =" + count1);
		
		FileWriter fw = new FileWriter(dir + fol_sign + "kvalues_distribution" +"(k=" + k + ")" +"("+dir+")"+ ".txt");
		for (int i = 1; i < valuesDistribution.length; i++)
			fw.write(i + " " + valuesDistribution[i] + "\n");
		fw.close();*/



                double[] clust = new double[valuesDistribution.length+1];
                double[] clust1 = new double[valuesDistribution.length+1];
/*		File fclust = new File(file_name + "data.txt");
		FileWriter fwclust = new FileWriter(fclust);
                int z = valuesDistribution1.length-1;
                fwclust.write(z + " " + 1 + "\n");
		for (int i = 1; i < valuesDistribution1.length; i++)
                {
                         int o = valuesDistribution1[i] + i;
			 fwclust.write(o + "\n");
//                         fwclust.write(i + " " + valuesDistribution1[i] + "\n");
                }
		fwclust.close();

		Runtime run=Runtime.getRuntime();
		Process p=null;
		p=run.exec("fams//fams.exe 0 0 200 " + file_name + "data .\\");
		try {
			p.waitFor();
		} catch (InterruptedException e) {
			e.printStackTrace();
		}

		if (p.exitValue()!= 0)
		{
			 System.out.println("Error in clustering program");
		}

		File fout = new File("out_" + file_name + "data.txt");


		FileReader fr = new FileReader(fout);
		BufferedReader br = new BufferedReader(fr);

		String s = br.readLine();
		int iclust = 1;
		while(s!= null)
		{
                    StringTokenizer st = new StringTokenizer(s," ");
                    while (st.hasMoreTokens())
                    {
                        clust[iclust] = Double.parseDouble(st.nextToken());
 //                       clust1[iclust] = Double.parseDouble(st.nextToken());
                    }
                    iclust++;
                    s = br.readLine();
		}

		fclust.delete();*/

                double[] variations = new double[valuesDistribution.length+1];
                int varseglen = 25;
                for (int i = 1; i < valuesDistribution.length - varseglen; i++)
                {
                    double cnt = 0;
                    for (int j = i+1; j <= i + varseglen; j++)
                        cnt +=Math.abs(valuesDistribution1[j] - valuesDistribution1[j-1]);
                    variations[i] = cnt;
                }

//                FileWriter fw1 = new FileWriter(dir + fol_sign + "kvalues_distribution1" +"(k=" + k + ")" +"("+dir+")"+ ".txt");
                FileWriter fw1 = new FileWriter(dir + File.separator + "kvalues_distribution1.txt");
		for (int i = 1; i < valuesDistribution1.length; i++)
			fw1.write(i + " " + valuesDistribution1[i] + " " + clust[i] +" "+ clust1[i] +" " + variations[i] + "\n");
		fw1.close();


/*                FileWriter fw2 = new FileWriter(dir + fol_sign + "kvalues_distribution1_boxes" +"(k=" + k + ")" +"("+dir+")"+ ".txt");
                int count2 = 0;
                int sum = 0;
		for (int i = 1; i < valuesDistribution1.length; i++)
                {
                    if (count2 < 20)
                    {
                        count2++;
                        sum+=valuesDistribution1[i];
                    }
                    else
                    {
			fw2.write(i + " " + sum + "\n");
                        sum = 0;
                        count2 = 0;
                    }
                }
		fw2.close();*/
	}
	void clear()
	{
		int count = 0;
		 for (Read r : reads)
		 {
			 count++;
//			 System.out.println("Deleting files: read " + count + "/" + reads.size()+"\n");	
			 File f = new File("modes_kmer" + count + "data.txt");			
			 f.deleteOnExit();
			 f = new File("pilot_200_kmer" + count + "data.txt");
			 f.deleteOnExit();
			 f = new File("out_kmer" + count + "data.txt");
                         f.delete();
			 int i;
			 f.deleteOnExit();
			 
		 }
	}
	
	void deleteStrangeNucl()
	{
		ArrayList<Read> toDel = new ArrayList<Read>();
		for (Read r : reads)
		{
			if (r.nucl.contains("N"))
				toDel.add(r);
		}
		reads.removeAll(toDel);
			
	}
	void findErrors()
	{
		badReads = new ArrayList<Read>();
		goodReads = new ArrayList<Read>();
		for (Read r : reads)
		{
			boolean bad = false;
			for (Kmer k : r.kmers)
				if (k.kCount <= freqThr)
				{
					badReads.add(r);
					bad = true;
					break;
				}
			if (!bad)
				goodReads.add(r);
		}
		int bad = 0;
		int all = 0;
		for (Read r : badReads)
			bad+=r.frequency;
		for (Read r : reads)
			all+=r.frequency;
		System.out.println("bad = " + bad + "/ all = " + all);
		
/*		HashSet hs = new HashSet();
		for (Read r : goodReads)
			for (Kmer k: r.kmers)
				hs.add(k);
		
		System.out.println(hs.size());*/
		
		for (Read r : badReads)
		{
//			r.errorRegions = new ArrayList<ErrorRegion>();
			int a[] = new int[r.getLength()-k+1];
			double b[] = new double[r.getLength()-k+1];
			for (Kmer k : r.kmers)
			{
				a[k.getStartPos()] = k.getKCount();
				b[k.getStartPos()] = k.clusterParam1;
			}
			int i = 0;
			int l = 0;
			String s = "";
			int start = 0;
			double clustprm = 0;
			while (i <= r.getLength()-k)
			{
				if (a[i] <= freqThr)
				{
					if (l==0)
					{
						start = i;
						clustprm = b[i];
					}
					l++;
					s+=r.getNucl(i);
					r.NuclType[i] = 1;
				}
				else
				{
					if (l > 0)
					{
						if ((Math.abs(b[i] - clustprm) < 0.1) && (a[i] <= freqThr*additionalfreqThrMult)&&(clustprm >= 0))
						{
							l++;
							s+=r.getNucl(i);
							r.NuclType[i] = 1;
						}
						else
						{
							int upd_start = start-1;
							while (upd_start >=0)
							{
								if ((Math.abs(b[upd_start] - clustprm) < 0.1) && (a[upd_start] <= freqThr*additionalfreqThrMult)&&(clustprm >= 0))
								{
									l++;
									s = r.getNucl(upd_start) + s;
									r.NuclType[upd_start] = 1;
									start = upd_start;
									upd_start--;
								}
								else
									break;
							}
							r.addErrorRegion(new ErrorRegion(s,start,i-1,l));
							l = 0;
							s = "";
						}
					}
					else
					{
						l = 0;
						s = "";
					}
				}
				i++;
			}
			if (l > 0)
			{
				int upd_start = start-1;
				while (upd_start >=0)
				{
					if ((Math.abs(b[upd_start] - clustprm) < 0.1) && (a[upd_start] <= freqThr*additionalfreqThrMult)&&(clustprm >= 0))
					{
						l++;
						s = r.getNucl(upd_start) + s;
						r.NuclType[upd_start] = 1;
						start = upd_start;
						upd_start--;
					}
					else
						break;
				}
				r.addErrorRegion(new ErrorRegion(s,start,r.getLength()-k,l));
			}
			
			i = 0;
/*			while (i < r.errorRegions.size()-1)
			{
				if (r.errorRegions.get(i).end >= r.errorRegions.get(i+1).begin)
				{
// slow, can be faster					
					String union = r.nucl.substring(r.errorRegions.get(i).begin, r.errorRegions.get(i+1).end+1);
					int s_new = r.errorRegions.get(i).begin;
					int e_new = r.errorRegions.get(i+1).end;
					ErrorRegion r_new = new ErrorRegion(union,s_new,e_new,e_new-k+1-s_new);
					r.errorRegions.remove(i+1);
					r.errorRegions.remove(i);
					r.errorRegions.add(i, r_new);
				}
				else
					i++;
			}*/
		}
	}
	
	void printErrorsStat(String dir) throws IOException
	{
		File f = new File(dir);
		if (!f.exists())
			f.mkdir();
		String fol_sign = "//";
		String name = dir + fol_sign + "error_stat(k=" + k +").txt";
		FileWriter fw = new FileWriter(name); 
		
// all errors
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
			{
				String type = new String();
				if (er.begin == 0)
					type = "begin";
				if((er.end + k  == r.getLength())&&(er.begin != 0))
					type = "end";
				if((er.end + k  != r.getLength())&&(er.begin != 0))
				{
					if (r.getNucl(er.end) == r.getNucl(er.end+1))
						type = "ins";
					else
						type = "delOrRepl";
				}
				
				String s = new String();
				s = er.length + " " + r.name  + " " + er.begin + " " + type + " " + (double) er.length/r.getLength()+ " " + r.getLength()+ "\n";
				fw.write(s);
			}
		
		fw.close();
		
		// how many errors of different lengths
		fw = new FileWriter(dir+fol_sign +"error_stat_lendistrib_(k=" + k +").txt");
		int ml = 0;
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
				if (er.length > ml)
					ml = er.length;
		int[] a = new int[ml+1];
		for (Read r : badReads)
		{
			for (ErrorRegion er : r.errorRegions)
				a[er.length]++;
		}
		for (int i = 1; i <= ml; i++)
			fw.write(i + " " + a[i]+"\n");
		fw.close();
		
		FileWriter fw1 = new FileWriter(dir+fol_sign +"kerrors_number_repl(k=" + k +").txt");
		FileWriter fw2 = new FileWriter(dir+fol_sign +"(k-2)errors_stat(k=" + k +").txt");
		FileWriter fw3 = new FileWriter(dir+fol_sign +"(k-3)errors_stat(k=" + k +").txt");
		FileWriter fw4 = new FileWriter(dir+fol_sign +"(k-4)errors_stat(k=" + k +").txt");
		FileWriter fwm1 = new FileWriter(dir+fol_sign +"(k-1)errors_stat(k=" + k +").txt");
		fw = new FileWriter(dir+fol_sign +"kerrors_stat(k=" + k +").txt");
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
			{
				if ((er.length == k)&&(er.begin!= 0)&&(er.end + k != r.getLength()))
				{
						String s = r.nucl.substring(er.begin, er.end);
						char c = r.nucl.charAt(er.end);
						char repl = c;
						int maxKcount = 0;
                                                int count = 0;
						for (Kmer_general km : allKmers)
							if ((km.kCount > freqThr)&&(km.nucl.startsWith(s)))
								if (km.nucl.charAt(k-1)!= c)
                                                                {
                                                                        count++;
									if (km.kCount > maxKcount)
									{
										maxKcount = km.kCount;
                                                                                repl = km.nucl.charAt(k-1);
									}
                                                                }
						fw.write(r.name + " " + er.begin + " " + maxKcount + " " + count + " ");
                                                if (repl == r.getNucl(er.end+1))
                                                    fw.write("ins" + "\n");
                                                else
                                                    fw.write("repl" + "\n");
					
				}
				if (er.length == k-1)
				{
//					System.out.println("(k-1)-errors considering: " + r.name);
					fwm1.write(r.name  + " " + er.begin + " ");
					if (r.getNucl(er.end) == r.getNucl(er.end+1))
						fwm1.write(1 + "\n");
					else
						fwm1.write(0 + "\n");
				}
				if (er.length == k-2)
				{
//					System.out.println("(k-2)-errors considering: " + r.name);
					fw2.write(r.name  + " " + er.begin + " ");
					int l = 0;
					int i = er.end+1;
					while((r.getNucl(i) == r.getNucl(er.end))&&(i < r.getLength()))
					{
						l++;
						i++;
					}
					fw2.write(l + " ");
					String s = r.getNucl(er.begin, er.end+1);
					s = s + r.getNucl(er.end+1);
					int bestKcount = 0;
					for (Kmer_general km : allKmers)
						if ((km.kCount > freqThr)&&(km.nucl.contentEquals(s)))
						{
							if (km.kCount > bestKcount)
								bestKcount = km.kCount;
						}
					fw2.write(bestKcount + "\n");
					
				}
				if (er.length == k-3)
				{
//					System.out.println("(k-3)-errors considering: " + r.name);
					fw3.write(r.name  + " " + er.begin + " ");
					int l = 0;
					int q = 0;
					int i = er.end+1;
					int j = er.end + 2;
					while((r.getNucl(i) == r.getNucl(er.end))&&(i < r.getLength()))
					{
						l++;
						i++;
					}
					if (r.getNucl(er.end)!= r.getNucl(er.end + 1))
						while((r.getNucl(j) == r.getNucl(er.end+1))&&(j < r.getLength()))
						{
							q++;
							j++;
						}
					fw3.write(l + " " + q + " ");
					int bestKcount = 0;
					if (q == 1)
					{
						String s = r.getNucl(er.begin, er.end+1+q);
						s = s + r.getNucl(er.end+1);
						
						for (Kmer_general km : allKmers)
							if ((km.kCount > freqThr)&&(km.nucl.contentEquals(s)))
							{
								if (km.kCount > bestKcount)
									bestKcount = km.kCount;
							}
					}
					fw3.write(bestKcount + "\n");
					
				}
				if (er.length == k-4)
				{
//					System.out.println("(k-4)-errors considering: " + r.name);
					fw4.write(r.name  + " " + er.begin + " ");
					int l = 0;
					int q = 0;
					int i = er.end+1;
					int j = er.end + 2;
					while((r.getNucl(i) == r.getNucl(er.end))&&(i < r.getLength()))
					{
						l++;
						i++;
					}
					if (r.getNucl(er.end)!= r.getNucl(er.end + 1))
						while((r.getNucl(j) == r.getNucl(er.end+1))&&(j < r.getLength()))
						{
							q++;
							j++;
						}
					fw4.write(l + " " + q + " ");
					int bestKcount = 0;
					if (q == 2)
					{
						String s = r.getNucl(er.begin, er.end+1+q);
						s = s + r.getNucl(er.end+1);
						
						for (Kmer_general km : allKmers)
							if ((km.kCount > freqThr)&&(km.nucl.contentEquals(s)))
							{
								if (km.kCount > bestKcount)
									bestKcount = km.kCount;
							}
					}
					fw4.write(bestKcount + "\n");
					
				}
			}
		
		System.out.println(allKmers.size());
		int bad = 0;
		int good = 0;
		for (Read r : badReads)
			bad += r.frequency;
		for (Read r : goodReads)
			good += r.frequency;
		System.out.println(bad + "/" + good);
		fw.close();
		fw1.close();
		fw2.close();
		fw3.close();
		fw4.close();
		fwm1.close();
		
		FileWriter fws = new FileWriter(dir+fol_sign +"strange_nucl(k=" + k +").txt");
		for (Read r : reads)
		{
			for (int i = 0; i < r.getLength(); i++)
				if ((r.getNucl(i)!='A')&&(r.getNucl(i)!='T')&&(r.getNucl(i)!='G')&&(r.getNucl(i)!='C'))
					fws.write(r.name + " " + i + "\n");
		}
		fws.close();

/*		int i = 0;
		for (Read r : reads)
		{
			if(r.nucl.contains("GACCAGTTCATCATCATATCCCTAN"))
				i++;
		}
		System.out.println("AAAAA=" + i);*/
		FileWriter fw_er_len_ratio = new FileWriter(dir+fol_sign +"errors_reg_len_ratio(k=" + k +").txt");
		for (Read r : badReads)
		{
			double l = 0;
			int lr = r.getLength();
			for (ErrorRegion er : r.errorRegions)
				l+=er.length;
			double d =  100.0*(double)(l/lr);
			fw_er_len_ratio.write(r.name + " " + l + " " + r.getLength() + " " + d + "\n");
		}
		fw_er_len_ratio.close();
		
	}
	void cutTails()
	{
		for (Read r : badReads)
		{
			System.out.println("Deleting tails: " + r.name +"\n" );
			ArrayList<ErrorRegion> corrected = new ArrayList<ErrorRegion>();
			for (ErrorRegion er : r.errorRegions)
			{
				if (er.begin == 0)
					for (int i = er.begin; i <= er.end; i++)
						r.corrections[i] = new Correction('D');
				
				if ((er.end + k  == r.getLength())&&(er.begin != 0))
					for (int i = er.begin+k-1; i < r.getLength(); i++)
						r.corrections[i] = new Correction('D');
			/*	if (er.length > r.getLength()*(maxErrorlenPerc/100))
					reads.remove(r);*/
			}
		}
	}
	void correctErrors_k()
	{
		
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
				if ((er.length == k)&&(er.begin!= 0)&&(er.end + k != r.getLength()))
				{
					System.out.println("k-errors correction: " + r.name);
					String s = r.nucl.substring(er.begin, er.end);
					char c = r.nucl.charAt(er.end);
					char repl = c;
					int maxKcount = 0;
					for (Kmer_general km : allKmers)
						if ((km.kCount > freqThr)&&(km.nucl.startsWith(s)))
							if (km.nucl.charAt(k-1)!= c)
								if (km.kCount > maxKcount)
								{
									maxKcount = km.kCount;
									repl = km.nucl.charAt(k-1);
								}
					r.corrections[er.end] = new Correction('R',""+repl);
				}
		
		
	}
	void correctErrors_kminus1()
	{
		
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
				if ((er.length == k-1)&&(er.begin!= 0)&&(er.end + k < r.getLength()-1))
				{
					System.out.println("(k-1)-errors correction: " + r.name);
					if (r.getNucl(er.end) == r.getNucl(er.end+1))
						r.corrections[er.end] = new Correction('D');
				}
		
		
	}
	void correctErrors_kminus2()
	{
		
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
				if (er.length == k-2)
				{
					int l = 0;
					int i = er.end+1;
					while((r.getNucl(i) == r.getNucl(er.end))&&(i < r.getLength()))
					{
						l++;
						i++;
					}
					String s = r.getNucl(er.begin, er.end+1);
					s = s + r.getNucl(er.end+1);
					int bestKcount = 0;
					for (Kmer_general km : allKmers)
						if ((km.kCount > freqThr)&&(km.nucl.contentEquals(s)))
						{
							if (km.kCount > bestKcount)
								bestKcount = km.kCount;
						}
					if (l == 2)
						r.corrections[er.end] = new Correction('D');
					if (l==0)
						if (bestKcount > 0)
							r.corrections[er.end+1] = new Correction('I',r.getNucl(er.end+1));
				}
		
		
	}
	void correctErrors_kminus3()
	{
		
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
				if (er.length == k-3)
				{
					int l = 0;
					int q = 0;
					int i = er.end+1;
					int j = er.end + 2;
					while((r.getNucl(i) == r.getNucl(er.end))&&(i < r.getLength()))
					{
						l++;
						i++;
					}
					if (r.getNucl(er.end)!= r.getNucl(er.end + 1))
						while((r.getNucl(j) == r.getNucl(er.end+1))&&(j < r.getLength()))
						{
							q++;
							j++;
						}
					int bestKcount = 0;
					if (q == 1)
					{
						String s = r.getNucl(er.begin, er.end+1+q);
						s = s + r.getNucl(er.end+1);
						
						for (Kmer_general km : allKmers)
							if ((km.kCount > freqThr)&&(km.nucl.contentEquals(s)))
							{
								if (km.kCount > bestKcount)
									bestKcount = km.kCount;
							}
					}
					if (l == 3)
						r.corrections[er.end] = new Correction('D');
					if (q==1)
						if (bestKcount > 0)
							r.corrections[er.end+1+q] = new Correction('I',r.getNucl(er.end+1));
				}
		
		
	}
	void correctErrors_kminus4()
	{		
		for (Read r : badReads)
			for (ErrorRegion er : r.errorRegions)
				if (er.length == k-4)
				{
					int l = 0;
					int q = 0;
					int i = er.end+1;
					int j = er.end + 2;
					while((r.getNucl(i) == r.getNucl(er.end))&&(i < r.getLength()))
					{
						l++;
						i++;
					}
					if (r.getNucl(er.end)!= r.getNucl(er.end + 1))
						while((r.getNucl(j) == r.getNucl(er.end+1))&&(j < r.getLength()))
						{
							q++;
							j++;
						}
					int bestKcount = 0;
					if (q == 2)
					{
						String s = r.getNucl(er.begin, er.end+1+q);
						s = s + r.getNucl(er.end+1);
						
						for (Kmer_general km : allKmers)
							if ((km.kCount > freqThr)&&(km.nucl.contentEquals(s)))
							{
								if (km.kCount > bestKcount)
									bestKcount = km.kCount;
							}
					}
					if (l == 4)
						r.corrections[er.end] = new Correction('D');
					if (q==2)
						if (bestKcount > 0)
							r.corrections[er.end+1+q] = new Correction('I',r.getNucl(er.end+1));
				}	
	}
	void PrintCorrectedReads(String outfile) throws IOException

	{
		FileWriter fw = new FileWriter(outfile);
		for(Read r : reads)
		{
			String s = "";
			for (int i = 0; i < r.getLength(); i++)
				if (r.corrections[i] != null)
				{
					if (r.corrections[i].type == 'R')
						s = s + r.corrections[i].replacement;
					if (r.corrections[i].type == 'I')
						s = s + r.nucl.charAt(i) + r.corrections[i].replacement;
				}
				else
					s = s + r.nucl.charAt(i);
			if (s.length() > lengthThr)
				for (int i = 0; i < r.frequency; i++)
					fw.write(">" + r.name + "_" + i + "\n" + s + "\n");
		}
		fw.close();
	}
	public void PrintHaplotypes(String outfile) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		HaplComp hc = new HaplComp();
		Collections.sort(haplotypes, hc);
		
		double count = 0;
		for (Haplotype h : haplotypes)
			count += h.concentration;
		for (Haplotype h : haplotypes)
		{
			int decimalPlace = 4;
			BigDecimal bd = new BigDecimal(100*h.concentration);
			bd = bd.setScale(decimalPlace,BigDecimal.ROUND_UP);
			double x = bd.doubleValue();
			bd = new BigDecimal(100*h.concentration/count);
			bd = bd.setScale(decimalPlace,BigDecimal.ROUND_UP);
			double y = bd.doubleValue();
			fw.write(">" + h.name + "_" + h.frequency + "_" + x + "_"+ "|" + y + "|" + "\n" + h.nucl + "\n");
		}
		fw.close();
	}
        public void PrintHaplotypes(String outfile, String idmethod) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		HaplComp hc = new HaplComp();
		Collections.sort(haplotypes, hc);
		
		double count = 0;
		for (Haplotype h : haplotypes)
			count += h.concentration;
		for (Haplotype h : haplotypes)
		{
			int decimalPlace = 4;
                        if (!idmethod.equalsIgnoreCase("ET"))
                        {
                            BigDecimal bd = new BigDecimal(100*h.concentration);
                            bd = bd.setScale(decimalPlace,BigDecimal.ROUND_UP);
                            double x = bd.doubleValue();
                            bd = new BigDecimal(100*h.concentration/count);
                            bd = bd.setScale(decimalPlace,BigDecimal.ROUND_UP);
                            double y = bd.doubleValue();
                            fw.write(">" + h.name + "_" + h.frequency + "_" + x + "_"+ "|" + y + "|" + "\n" + h.nucl + "\n");
                        }
                        else
                        {
                            fw.write(">" + h.name + "_" + h.concentration + "\n"+ h.nucl + "\n");
                        }
		}
		fw.close();
	}
        public void PrintHaplotypesWithNameTag(String outfile, String tag) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		HaplComp hc = new HaplComp();
		Collections.sort(haplotypes, hc);
		
		double count = 0;
		for (Haplotype h : haplotypes)
			count += h.concentration;
		for (Haplotype h : haplotypes)
		{
			int decimalPlace = 4;
			BigDecimal bd = new BigDecimal(100*h.concentration);
			bd = bd.setScale(decimalPlace,BigDecimal.ROUND_UP);
			double x = bd.doubleValue();
			bd = new BigDecimal(100*h.concentration/count);
			bd = bd.setScale(decimalPlace,BigDecimal.ROUND_UP);
			double y = bd.doubleValue();
			fw.write(">"+ h.name + "_" + tag + "_" + h.frequency + "_" + x + "_"+ "|" + y + "|" + "\n" + h.nucl + "\n");
		}
		fw.close();
	}
	public void PrintUniqueReads(String outfile) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		for (Read r : reads)
		{
			fw.write(">" + r.name + "_" + r.frequency + "\n" + r.nucl + "\n");		
		}
		fw.close();
	}
        public void PrintUniqueReadsWithTagGenotype(String outfile, String tag) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		for (int i = 0; i < reads.size(); i++)
		{
			fw.write(">" + tag + "_"  + reads.get(i).genotype + "_" + (i+1) + "_" + reads.get(i).frequency + "\n" + reads.get(i).nucl + "\n");		
		}
		fw.close();
	}
        public void PrintUniqueReadsNoFreqTag(String outfile) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		for (Read r : reads)
		{
			fw.write(">" + r.name + "\n" + r.nucl + "\n");		
		}
		fw.close();
	}
	public void findHaplotypes() throws IOException
	{
		haplotypes = new ArrayList<Haplotype>();
		StrComp sc = new StrComp();
		Collections.sort(reads, sc);
		System.out.println("Finding haplotypes");
                int goodhapl = 0;
                int badhapl = 0;
                int lenratio = 90; //90
                int i = 1;
		for (Read r : reads)
		{
                    System.out.println("Hapl: read " + i + "\\" + reads.size());
			boolean toAdd = true;
			int nhap = 0;
			for (Haplotype h : haplotypes)
				if (h.containsRead(r))
				{
					toAdd = false;
					nhap++;
				}
			if (toAdd)
                        {
				haplotypes.add(new Haplotype(r));
                       /*         if (badReads.contains(r))
                                    badhapl++;
                                else
                                    goodhapl++;*/
                        }
			else
				for (Haplotype h : haplotypes)
					if (h.containsRead(r))
						h.concentration += r.frequency/nhap;
		}
                
                double allconc = 0;
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                
		for (Haplotype h : haplotypes)
			h.concentration = h.concentration / allconc;


                ArrayList<Haplotype> shorthapl = new ArrayList<Haplotype>();
                allconc = 0;
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                double avlen = 0;
                for (Haplotype h : haplotypes)
                    avlen += h.getLength()*((double) h.concentration/allconc);
                for (Haplotype h : haplotypes)
                    if (100 * ((double) h.getLength()/avlen) < lenratio)
                        shorthapl.add(h);
                haplotypes.removeAll(shorthapl);
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                for (Haplotype h : haplotypes)
                    h.concentration /= allconc;


                System.out.println("Goodhapl=" + goodhapl);
                System.out.println("Badhapl=" + badhapl);
	}
        public void findHaplotypes(String idmethod) throws IOException
	{
		haplotypes = new ArrayList<Haplotype>();
		StrComp sc = new StrComp();
		Collections.sort(reads, sc);
		System.out.println("Finding haplotypes");
                int goodhapl = 0;
                int badhapl = 0;
                int lenratio = 90; //90
		for (Read r : reads)
		{
			boolean toAdd = true;
			int nhap = 0;
			for (Haplotype h : haplotypes)
				if (h.containsRead(r))
				{
					toAdd = false;
					nhap++;
				}
			if (toAdd)
                        {
				haplotypes.add(new Haplotype(r));
                       /*         if (badReads.contains(r))
                                    badhapl++;
                                else
                                    goodhapl++;*/
                        }
			else
				for (Haplotype h : haplotypes)
					if (h.containsRead(r))
						h.concentration += r.frequency/nhap;
		}
                if (!idmethod.equalsIgnoreCase("ET"))
                    for (Haplotype h : haplotypes)
			h.concentration = h.concentration / getNreads();


                ArrayList<Haplotype> shorthapl = new ArrayList<Haplotype>();
                double allconc = 0;
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                double avlen = 0;
                for (Haplotype h : haplotypes)
                    avlen += h.getLength()*((double) h.concentration/allconc);
                for (Haplotype h : haplotypes)
                    if (100 * ((double) h.getLength()/avlen) < lenratio)
                        shorthapl.add(h);
                haplotypes.removeAll(shorthapl);
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                
                if (!idmethod.equalsIgnoreCase("ET"))
                    for (Haplotype h : haplotypes)
                        h.concentration /= allconc;


                System.out.println("Goodhapl=" + goodhapl);
                System.out.println("Badhapl=" + badhapl);
	}
        public void findHaplotypes(String idmethod, String idsort) throws IOException
	{
		haplotypes = new ArrayList<Haplotype>();
                if (idsort.equalsIgnoreCase("Length"))
                {
                    StrComp sc = new StrComp();
                    Collections.sort(reads, sc);
                }
                
                if (idsort.equalsIgnoreCase("Frequency"))
                {
                    ReadComp sc = new ReadComp();
                    Collections.sort(reads, sc);
                }
                
		System.out.println("Finding haplotypes");
                int goodhapl = 0;
                int badhapl = 0;
                int lenratio = 90; //90
		for (Read r : reads)
		{
			boolean toAdd = true;
			int nhap = 0;
			for (Haplotype h : haplotypes)
				if (h.containsRead(r))
				{
					toAdd = false;
					nhap++;
				}
			if (toAdd)
                        {
				haplotypes.add(new Haplotype(r));
                       /*         if (badReads.contains(r))
                                    badhapl++;
                                else
                                    goodhapl++;*/
                        }
			else
				for (Haplotype h : haplotypes)
					if (h.containsRead(r))
						h.concentration += r.frequency/nhap;
		}
                if (!idmethod.equalsIgnoreCase("ET"))
                    for (Haplotype h : haplotypes)
			h.concentration = h.concentration / getNreads();


                ArrayList<Haplotype> shorthapl = new ArrayList<Haplotype>();
                double allconc = 0;
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                double avlen = 0;
                for (Haplotype h : haplotypes)
                    avlen += h.getLength()*((double) h.concentration/allconc);
                for (Haplotype h : haplotypes)
                    if (100 * ((double) h.getLength()/avlen) < lenratio)
                        shorthapl.add(h);
                haplotypes.removeAll(shorthapl);
                for (Haplotype h : haplotypes)
                    allconc+=h.concentration;
                
                if (!idmethod.equalsIgnoreCase("ET"))
                    for (Haplotype h : haplotypes)
                        h.concentration /= allconc;


                System.out.println("Goodhapl=" + goodhapl);
                System.out.println("Badhapl=" + badhapl);
	}
        public void findHaplotypesAlign(String idmethod, String idsort, int gapop, int gapext) throws IOException, CompoundNotFoundException
	{
		haplotypes = new ArrayList<Haplotype>();
                if (idsort.equalsIgnoreCase("Length"))
                {
                    StrComp sc = new StrComp();
                    Collections.sort(reads, sc);
                }
                
                if (idsort.equalsIgnoreCase("Frequency"))
                {
                    ReadComp sc = new ReadComp();
                    Collections.sort(reads, sc);
                }
                
		System.out.println("Finding haplotypes");
		for (int i = 0; i < reads.size(); i++)
		{
                        Read r = reads.get(i);
                        System.out.println(i + "//" + reads.size());
                    
			boolean toAdd = true;
                        ArrayList<Haplotype> conthapls = new ArrayList();
			for (Haplotype h : haplotypes)
				if (h.containsReadAlign(r,gapop,gapext))
                                    conthapls.add(h);
			if (conthapls.size() == 0)
				haplotypes.add(new Haplotype(r));
			else
				for (Haplotype h : conthapls)
					h.frequency += r.frequency/conthapls.size();
		}

                double allconc = 0;
                for (Haplotype h : haplotypes)
                    allconc+=h.frequency;
                
                    for (Haplotype h : haplotypes)
                    {
                        if (!idmethod.equalsIgnoreCase("ET"))
                            h.concentration = ((double) h.frequency) / allconc;
                        else
                            h.concentration = h.frequency;
                    }
                    

	}
	void compareHaplWithClones(String clonefile) throws IOException
	{
		DataSet cds = new DataSet(clonefile,'c');
		for (Haplotype h : haplotypes)
		{
			String clns = "";
			for (Read r : cds.reads)
				if (h.nucl.contains(r.nucl) || r.nucl.contains(h.nucl))
					clns = clns + r.name;
			if (clns != "")
				h.clone = clns;
		}
	}
	void delLargeErrors()
	{
		for (Read r : badReads)
		{
			double l = 0;
			int lr = r.getLength();
			for (ErrorRegion er : r.errorRegions)
				l+=er.length;
			double d =  100.0*(double)(l/lr);
			if (d > r.getLength()*(maxAllErrorsPerc/100))
				reads.remove(r);
		}
	}
	int getNerrors()
	{
		int count = 0;
		for (Read r : reads)
			count+= r.errorRegions.size();
		return count;
	}
	int maxReadLen()
	{
		int m = 0;
		for (Read r : reads)
			if (r.getLength() > m)
				m = r.getLength();
		return m;
	}
	int minReadLen()
	{
		int m = Integer.MAX_VALUE;
		for (Read r : reads)
			if (r.getLength() < m)
				m = r.getLength();
		return m;
	}
	public int getNreads()
	{
		int n = 0;
		for (Read r : reads)
			n+= r.frequency;
		return n;
	}
	int getNuniquereads()
	{
		return reads.size();
	}
	double AvReadLen()
	{
		double m = 0;
		for (Read r : reads)
			m += r.getLength();
		return m/reads.size();
	}
	int getNUniqueErrorReads()
	{
		return badReads.size();
	}
        void PrintReadsDistrib(String dir) throws IOException
	{
		String fol_sign = "//";

		FileWriter fw = new FileWriter(dir + fol_sign + "kvalues_distribution" +"(k=" + k + ")" +"("+dir+")"+ ".txt");
		for (int i = 1; i < valuesDistribution.length; i++)
			fw.write(i + " " + valuesDistribution[i] + "\n");
		fw.close();

                double[] clust = new double[valuesDistribution.length+1];
                double[] clust1 = new double[valuesDistribution.length+1];
/*		File fclust = new File(file_name + "data.txt");
		FileWriter fwclust = new FileWriter(fclust);
                int z = valuesDistribution1.length-1;
                fwclust.write(z + " " + 1 + "\n");
		for (int i = 1; i < valuesDistribution1.length; i++)
                {
                         int o = valuesDistribution1[i] + i;
			 fwclust.write(o + "\n");
//                         fwclust.write(i + " " + valuesDistribution1[i] + "\n");
                }
		fwclust.close();

		Runtime run=Runtime.getRuntime();
		Process p=null;
		p=run.exec("fams//fams.exe 0 0 200 " + file_name + "data .\\");
		try {
			p.waitFor();
		} catch (InterruptedException e) {
			e.printStackTrace();
		}

		if (p.exitValue()!= 0)
		{
			 System.out.println("Error in clustering program");
		}

		File fout = new File("out_" + file_name + "data.txt");


		FileReader fr = new FileReader(fout);
		BufferedReader br = new BufferedReader(fr);

		String s = br.readLine();
		int iclust = 1;
		while(s!= null)
		{
                    StringTokenizer st = new StringTokenizer(s," ");
                    while (st.hasMoreTokens())
                    {
                        clust[iclust] = Double.parseDouble(st.nextToken());
 //                       clust1[iclust] = Double.parseDouble(st.nextToken());
                    }
                    iclust++;
                    s = br.readLine();
		}

		fclust.delete();*/

                double[] variations = new double[valuesDistribution.length+1];
                int varseglen = 25;
                for (int i = 1; i < valuesDistribution.length - varseglen; i++)
                {
                    double cnt = 0;
                    for (int j = i+1; j <= i + varseglen; j++)
                        cnt +=Math.abs(valuesDistribution1[j] - valuesDistribution1[j-1]);
                    variations[i] = cnt;
                }

                FileWriter fw1 = new FileWriter(dir + fol_sign + "kvalues_distribution1" +"(k=" + k + ")" +"("+dir+")"+ ".txt");
		for (int i = 1; i < valuesDistribution1.length; i++)
			fw1.write(i + " " + valuesDistribution1[i] + " " + clust[i] +" "+ clust1[i] +" " + variations[i] + "\n");
		fw1.close();


/*                FileWriter fw2 = new FileWriter(dir + fol_sign + "kvalues_distribution1_boxes" +"(k=" + k + ")" +"("+dir+")"+ ".txt");
                int count2 = 0;
                int sum = 0;
		for (int i = 1; i < valuesDistribution1.length; i++)
                {
                    if (count2 < 20)
                    {
                        count2++;
                        sum+=valuesDistribution1[i];
                    }
                    else
                    {
			fw2.write(i + " " + sum + "\n");
                        sum = 0;
                        count2 = 0;
                    }
                }
		fw2.close();*/
	}
        void PrintReadsFreqDistrib(String addr) throws IOException
        {
            FileWriter fw = new FileWriter(addr);
            for (Read r:reads)
                fw.write(r.frequency + "\n");
            fw.close();
        }
        void calcDirection(String ref_file, double gapop, double gapext) throws IOException
        {
            DataSet ref = new DataSet(ref_file);
            FileWriter fw = new FileWriter(this.file_name +"_refStat.txt");
            for (Read r : reads)
            {
                fw.write(r.name + " ");
                for (int i = 0; i < 2; i++)
                {
                    String s_ref = ref.reads.get(0).nucl;
                    String q = "Direct";
                    Runtime run=Runtime.getRuntime();
                    Process p=null;
                    FileWriter fw_alin = new FileWriter("allign_input.fas");
                    fw_alin.write(r.name + "\n");
                    fw_alin.write(r.nucl + "\n");
                    fw_alin.write(">reference" + "\n");
                    fw_alin.write(s_ref + "\n");
                    fw_alin.close();
                    String param = "ClustalW2//clustalw2 -INFILE=" +"allign_input.fas" + " -OUTFILE=" + "allign_output.fas";
                    param+= " -OUTPUT=FASTA -DNAMATRIX=IUB -GAPOPEN=" + gapop +" -GAPEXT=" + gapext + " -TYPE=DNA -PWDNAMATRIX=IUB -PWGAPOPEN=" + gapop+ " -PWGAPEXT=" + gapext;

                    p = run.exec(param);
                    BufferedReader stdInput = new BufferedReader(new InputStreamReader(p.getInputStream()));
                    String s1;
                    while ((s1 = stdInput.readLine()) != null) {
    //                  System.out.println(s1);
                    }
                    if (p.exitValue()!= 0)
                    {
                        System.out.println("Error in allgnment program");
                    }
                    DataSet allignment = new DataSet("allign_output.fas",'c');
                    String al_read = "";
                    String al_ref = "";
                    for (Read r1 : allignment.reads)
                        if (!r1.name.contains("reference"))
                            al_read = r1.getNucl();
                        else
                            al_ref = r1.getNucl();
                    int match = 0;
                    for (int j = 0; j < al_read.length(); j++)
                        if (al_read.charAt(j) == al_ref.charAt(j))
                            match++;
                    double score = (double) match/r.getLength();
                    fw.write(q + " " + score + " ");
                    String rev_ref = "";
                    for (int j = s_ref.length() - 1; j >=0; j --)
                    {
                        if (s_ref.charAt(j) == 'A')
                            rev_ref += 'T';
                        if (s_ref.charAt(j) == 'T')
                            rev_ref += 'A';
                        if (s_ref.charAt(j) == 'G')
                            rev_ref += 'C';
                        if (s_ref.charAt(j) == 'C')
                            rev_ref += 'G';
                    }
                    s_ref = rev_ref;
                    q = "Reverse";
                    File f = new File("allign_input.fas");
                    f.delete();
                    f = new File("allign_output.fas");
                    f.delete();
                    f = new File("allign_input.dnd");
                    f.delete();
                }
                fw.write("\n");
            }
            fw.close();
        }
        void convertZoyaData(String addr) throws FileNotFoundException, IOException
        {
                FileReader fr = new FileReader(addr);
		BufferedReader br = new BufferedReader(fr);
                FileWriter fw = new FileWriter(addr + "_converted.fas");
	
		String s = br.readLine();
                while(s.equalsIgnoreCase(""))
                    s = br.readLine();

		String nucl = new String();
		String name = new String();
		name = "";
		int i = 1;
		int j = s.length();
		while ((i < s.length())&&(s.charAt(i)!=' '))
		{
			name+=s.charAt(i);
			i++;
		}
		s = br.readLine();
		try
		{
		while (s!= null)
		{
//			System.out.println("Reading " + name);

			if (s.length() == 0)
			{
				s = br.readLine();
				continue;
			}
			if(s.charAt(0) == '>')
                        {
                                        StringTokenizer st = new StringTokenizer(name,"_");
                                        st.nextToken();
                                        st.nextToken();
                                        st.nextToken();
                                        st.nextToken();
                                        st.nextToken();
                                        int mult = Integer.parseInt(st.nextToken());
                                        for (int u = 0; u < mult; u++)
                                            fw.write(">" + name + "_" + u + "\n" + nucl + "\n");
					nucl = "";
					name = "";
					i = 1;
					while ((i < s.length())&&(s.charAt(i)!=' '))
					{
						name+=s.charAt(i);
						i++;
					}
			}
			else
				nucl+=s.toUpperCase();
			s = br.readLine();
		}
		}
		catch (Exception e)
		{
			e.printStackTrace();
		}
	//	reads.add(new Read(nucl, name));
                StringTokenizer st = new StringTokenizer(name,"_");
                st.nextToken();
                st.nextToken();
                st.nextToken();
                st.nextToken();
                st.nextToken();
                int mult = Integer.parseInt(st.nextToken());
                for (int u = 0; u < mult; u++)
                    fw.write(">" + name + "_" + u + "\n" + nucl + "\n");
		fr.close();
                fw.close();
        }
        void calculateRealErrorsStat(String ref, int homoplen, double gapop, double gapext) throws IOException
        {
            Read clone = (new DataSet(ref)).reads.get(0);
            clone.name = "Clone";
            FileWriter fw = new FileWriter(this.file_name + "_" + ref + "_realerrors.txt");
            fw.write("read freq len averrorsdist repl ins del ");
            for (int i = 1; i <= homoplen; i++)
                fw.write("homop" + i + " ");
            fw.write("allerr nreplat nreplac nreplag nrepltc nrepltg nreplgc" + "\n");
            FileWriter fw_un = new FileWriter(this.file_name + "_" + ref + "_realerrorsun.txt");
/*            fw_un.write("read freq len averrorsdist repl ins del ");
            for (int i = 1; i <= homoplen; i++)
                fw_un.write("homop" + i + " ");
            fw_un.write("allerr nreplat nreplac nreplag nrepltc nrepltg nreplgc" + "\n");*/
            
            
            for (Read r : reads)
            {
//                    DynamicOut.printStep(r.name + "//" + reads.size());
                    System.out.println(r.name + "//" + reads.size());
                    Runtime run=Runtime.getRuntime();
                    Process p=null;
                    String al_in = this.file_name + "_allign_input.fas";
                    String al_out = this.file_name + "_allign_output.fas";
                    FileWriter fw_alin = new FileWriter(al_in);

                    fw_alin.write(">" + r.name + "\n" + r.nucl + "\n");
                    fw_alin.write(">" + clone.name + "\n" + clone.nucl + "\n");
                    fw_alin.close();
                    String param = "ClustalW2//clustalw2 -INFILE=" +al_in + " -OUTFILE=" + al_out;
                    param+= " -OUTPUT=FASTA -DNAMATRIX=IUB -GAPOPEN=" + gapop +" -GAPEXT=" + gapext + " -TYPE=DNA -PWDNAMATRIX=IUB -PWGAPOPEN=" + gapop+ " -PWGAPEXT=" + gapext;
                    p = run.exec(param);
                    BufferedReader stdInput = new BufferedReader(new InputStreamReader(p.getInputStream()));
                    String s;
                    while ((s = stdInput.readLine()) != null)
                    {
//                        System.out.println(s);
                    }
                    if (p.exitValue()!= 0)
                    {
                        System.out.println("Error in allgnment program");
                    }
                    DataSet alignment = new DataSet(al_out,'c');
/*                    File f = new File("allign_input.fas");
                    f.delete();
                    f = new File("allign_output.fas");
                    f.delete();
                    f = new File("allign_input.dnd");
                    f.delete();*/
                    String read_align = "";
                    String clone_align = "";
                    for (Read r1: alignment.reads)
                    {
                        if (r1.name.equalsIgnoreCase(r.name))
                            read_align = r1.nucl;
                        if (r1.name.equalsIgnoreCase(clone.name))
                            clone_align = r1.nucl;
                    }
                    
                    int l = read_align.length();
                    int pos_r = -1;
                    int pos_c = -1;

                    int startcompdiff = 0;
                    int endcompdiff = l-1;
                    String sz = "";
                    if (read_align.charAt(0) == '-')
                        sz = read_align;
                    else
                        sz = clone_align;
                    while (sz.charAt(startcompdiff) == '-')
                        startcompdiff++;
                    if (read_align.charAt(l-1) == '-')
                        sz = read_align;
                    else
                        sz = clone_align;
                    while (sz.charAt(endcompdiff) == '-')
                        endcompdiff--;

                    int errors[] = new int[homoplen+2];
                    int nrepl = 0;
                    int nins = 0;
                    int ndel = 0;
                    int nreplat = 0;
                    int nreplac = 0;
                    int nreplag = 0;
                    int nrepltc = 0;
                    int nrepltg = 0;
                    int nreplgc = 0;
                    ArrayList<Integer> erpos = new ArrayList<Integer>();
                    for (int i = startcompdiff; i <= endcompdiff; i++)
                    {
                        if (read_align.charAt(i) == clone_align.charAt(i))
                            continue;
                        if ((read_align.charAt(i)!= '-')&&(clone_align.charAt(i)!='-'))
                        {
                            nrepl++;
                            int readpos = 0;
                            for (int j = 0; j <= i; j++)
                                if (read_align.charAt(j)!='-')
                                    readpos++;
                            erpos.add(readpos);
                            if ((read_align.charAt(i) == 'A')&&(clone_align.charAt(i) == 'T'))
                                nreplat++;
                            if ((read_align.charAt(i) == 'T')&&(clone_align.charAt(i) == 'A'))
                                nreplat++;
                            if ((read_align.charAt(i) == 'A')&&(clone_align.charAt(i) == 'C'))
                                nreplac++;
                            if ((read_align.charAt(i) == 'C')&&(clone_align.charAt(i) == 'A'))
                                nreplac++;
                            if ((read_align.charAt(i) == 'A')&&(clone_align.charAt(i) == 'G'))
                                nreplag++;
                            if ((read_align.charAt(i) == 'G')&&(clone_align.charAt(i) == 'A'))
                                nreplag++;
                            if ((read_align.charAt(i) == 'T')&&(clone_align.charAt(i) == 'C'))
                                nrepltc++;
                            if ((read_align.charAt(i) == 'C')&&(clone_align.charAt(i) == 'T'))
                                nrepltc++;
                            if ((read_align.charAt(i) == 'T')&&(clone_align.charAt(i) == 'G'))
                                nrepltg++;
                            if ((read_align.charAt(i) == 'G')&&(clone_align.charAt(i) == 'T'))
                                nrepltg++;
                            if ((read_align.charAt(i) == 'G')&&(clone_align.charAt(i) == 'C'))
                                nreplgc++;
                            if ((read_align.charAt(i) == 'C')&&(clone_align.charAt(i) == 'G'))
                                nreplgc++;
                            continue;
                        }
                        char z;
                        if (read_align.charAt(i)!= '-')
                            z = read_align.charAt(i);
                        else
                            z = clone_align.charAt(i);
                        int bhomr = i;
                        int ehomr = i;
                        int bhomc = i;
                        int ehomc = i;
                        while (((read_align.charAt(ehomr)==z)||(read_align.charAt(ehomr) == '-'))&&(ehomr <= endcompdiff))
                            ehomr++;
                        while ((read_align.charAt(bhomr)==z)||(read_align.charAt(bhomr) == '-')&&(bhomr >= startcompdiff))
                            bhomr--;
                        while (((clone_align.charAt(ehomc)==z)||(clone_align.charAt(ehomc) == '-'))&&(ehomc <= endcompdiff))
                            ehomc++;
                        while ((clone_align.charAt(bhomc)==z)||(clone_align.charAt(bhomc) == '-')&&(bhomc >= startcompdiff))
                            bhomc--;
                        int b = Math.max(bhomr+1, bhomc+1);
                        int e = Math.min(ehomr-1, ehomc-1);
                        if (e == b)
                        {
                            if (read_align.charAt(i) == '-')
                                ndel++;
                            else
                                nins++;
                        }
                        else
                        {
                            int len = 0;
                            for (int j = b; j <= e; j++)
                                if (clone_align.charAt(j) == z)
                                    len++;
                            errors[len]++;
                        }
                        int q = b;
                        for (int j = b; j <=e; j++)
                            if (read_align.charAt(q)==z)
                            {
                                q = j;
                                break;
                            }
                        int readpos = 0;
                        for (int j = 0; j <= q; j++)
                            if (read_align.charAt(j)!='-')
                                readpos++;
                        erpos.add(readpos);
//                        int len = Math.min(e-b+1, homoplen);
                        i = e;
                    }
                    int erdist = 0;
                    for (int i = 1; i < erpos.size(); i++)
                        erdist += (erpos.get(i)-erpos.get(i-1)-1);
                    double averrorsdist = -1;
                    if (erpos.size() >= 2)
                        averrorsdist = (double) erdist/(erpos.size()-1);
                    int allerr = nrepl + nins + ndel;
                    for (int i = 1; i <= homoplen; i++)
                        allerr += errors[i];
                    fw.write(r.name + " " + r.frequency + " " +r.getLength() + " " + averrorsdist + " " + nrepl*r.frequency + " " + nins*r.frequency + " " + ndel*r.frequency + " ");
                    for (int i = 1; i <= homoplen; i++)
                        fw.write(errors[i]*r.frequency + " ");
                    fw.write(allerr*r.frequency + " " + nreplat*r.frequency + " " + nreplac*r.frequency + " " + nreplag*r.frequency + " " + nrepltc*r.frequency + " " + nrepltg*r.frequency + " " + nreplgc*r.frequency + "\n");
                    fw_un.write(r.frequency + " " +r.getLength() + " " + averrorsdist + " " + nrepl + " " + nins + " " + ndel + " ");
                    for (int i = 1; i <= homoplen; i++)
                        fw_un.write(errors[i] + " ");
                    fw_un.write(allerr + " " + nreplat+ " " + nreplac+ " " + nreplag+ " " + nrepltc+ " " + nrepltg+ " " + nreplgc + "\n");
                        
            }
            DynamicOut.finishSteps();
            fw.close();
            fw_un.close();
            
        }
        void setFindErrorsSeglen(int i)
        {
            finderrorsseglen = i;
        }
        void fixDirection(double gapop, double gapext) throws IOException
        {
            String s_ref = "";
            for (Read r : reads)
                if (r.getLength() > s_ref.length())
                    s_ref = r.getNucl();
            for (Read r : reads)
            {
                System.out.println("Fixing direction: " + r.name);
                int bestmatch = 0;
                int bestdirect = -1;
                String nucl = r.getNucl();
                for (int i = 0; i <= 1; i++)
                {
                    Runtime run=Runtime.getRuntime();
                    Process p=null;
                    FileWriter fw_alin = new FileWriter("allign_input.fas");
                    fw_alin.write(">" + r.name + "\n");
                    fw_alin.write(nucl + "\n");
                    fw_alin.write(">reference" + "\n");
                    fw_alin.write(s_ref + "\n");
                    fw_alin.close();
                    String param = "ClustalW2//clustalw2 -INFILE=" +"allign_input.fas" + " -OUTFILE=" + "allign_output.fas";
                    param+= " -OUTPUT=FASTA -DNAMATRIX=IUB -GAPOPEN=" + gapop +" -GAPEXT=" + gapext + " -TYPE=DNA -PWDNAMATRIX=IUB -PWGAPOPEN=" + gapop+ " -PWGAPEXT=" + gapext;

                    p = run.exec(param);
                    BufferedReader stdInput = new BufferedReader(new InputStreamReader(p.getInputStream()));
                    String s1;
                    while ((s1 = stdInput.readLine()) != null) {
    //                  System.out.println(s1);
                    }
                    if (p.exitValue()!= 0)
                    {
                        System.out.println("Error in allgnment program");
                    }
                    DataSet allignment = new DataSet("allign_output.fas",'c');
                    String al_read = "";
                    String al_ref = "";
                    for (Read r1 : allignment.reads)
                        if (!r1.name.contains("reference"))
                            al_read = r1.getNucl();
                        else
                            al_ref = r1.getNucl();
                    int match = 0;
                    for (int j = 0; j < al_read.length(); j++)
                        if (al_read.charAt(j) == al_ref.charAt(j))
                            match++;
                    if (match > bestmatch)
                    {
                        bestmatch = match;
                        bestdirect = i;
                    }
                    File f = new File("allign_input.fas");
                    f.delete();
                    f = new File("allign_output.fas");
                    f.delete();
                    f = new File("allign_input.dnd");
                    f.delete();
                    nucl = r.getRevComp().getNucl();
                }
                if (bestdirect == 1)
                    r.nucl = r.getRevComp().getNucl();
            }
        }
        void fixDirection(Read ref, double gapop, double gapext) throws IOException
        {
            String s_ref = ref.getNucl();
            for (Read r : reads)
            {
                System.out.println("Fixing direction: " + r.name);
                int bestmatch = 0;
                int bestdirect = -1;
                String nucl = r.getNucl();
                for (int i = 0; i <= 1; i++)
                {
                    Runtime run=Runtime.getRuntime();
                    Process p=null;
                    FileWriter fw_alin = new FileWriter("allign_input.fas");
                    fw_alin.write(">" + r.name + "\n");
                    fw_alin.write(nucl + "\n");
                    fw_alin.write(">reference" + "\n");
                    fw_alin.write(s_ref + "\n");
                    fw_alin.close();
                    String param = "ClustalW2//clustalw2 -INFILE=" +"allign_input.fas" + " -OUTFILE=" + "allign_output.fas";
                    param+= " -OUTPUT=FASTA -DNAMATRIX=IUB -GAPOPEN=" + gapop +" -GAPEXT=" + gapext + " -TYPE=DNA -PWDNAMATRIX=IUB -PWGAPOPEN=" + gapop+ " -PWGAPEXT=" + gapext;

                    p = run.exec(param);
                    BufferedReader stdInput = new BufferedReader(new InputStreamReader(p.getInputStream()));
                    String s1;
                    while ((s1 = stdInput.readLine()) != null) {
    //                  System.out.println(s1);
                    }
                    if (p.exitValue()!= 0)
                    {
                        System.out.println("Error in allgnment program");
                    }
                    DataSet allignment = new DataSet("allign_output.fas",'c');
                    String al_read = "";
                    String al_ref = "";
                    for (Read r1 : allignment.reads)
                        if (!r1.name.contains("reference"))
                            al_read = r1.getNucl();
                        else
                            al_ref = r1.getNucl();
                    int match = 0;
                    for (int j = 0; j < al_read.length(); j++)
                        if (al_read.charAt(j) == al_ref.charAt(j))
                            match++;
                    if (match > bestmatch)
                    {
                        bestmatch = match;
                        bestdirect = i;
                    }
                    File f = new File("allign_input.fas");
                    f.delete();
                    f = new File("allign_output.fas");
                    f.delete();
                    f = new File("allign_input.dnd");
                    f.delete();
                    nucl = r.getRevComp().getNucl();
                }
                if (bestdirect == 1)
                    r.nucl = r.getRevComp().getNucl();
            }
        }
        void fixDirectionRefParallel(DataSet refs, int gapop, int gapext) throws InterruptedException, ExecutionException
        {
                 List< Future > futuresList = new ArrayList< Future >();
                 int nrOfProcessors = Runtime.getRuntime().availableProcessors();
                 ExecutorService eservice = Executors.newFixedThreadPool(nrOfProcessors);
                 
		 for (Read r : reads)
		 {
                        RevCompTask rt = new RevCompTask(r,refs,gapop,gapext);
                        futuresList.add(eservice.submit(rt));
			 
		 }
                  Object taskResult;
                  int i = 0;
                  for(Future future:futuresList) 
                  {
                        i++;
                        taskResult = future.get();
                        System.out.println("RevComp: read "+i + "/" + reads.size());                            
                  }
        }
        public void fixDirectionGenotypingRefParallel(DataSet refs, int gapop, int gapext) throws InterruptedException, ExecutionException
        {
                 List< Future > futuresList = new ArrayList< Future >();
//                 int nrOfProcessors = Runtime.getRuntime().availableProcessors();
                 int nrOfProcessors = this.avProc;
                 ExecutorService eservice = Executors.newFixedThreadPool(nrOfProcessors);
                 
		 for (Read r : reads)
		 {
                        RevCompGenTask rt = new RevCompGenTask(r,refs,gapop,gapext);
                        futuresList.add(eservice.submit(rt));
			 
		 }
                  Object taskResult;
                  int i = 0;
                  for(Future future:futuresList) 
                  {
                        i++;
                        taskResult = future.get();
                        System.out.println("RevComp: read "+i + "/" + reads.size());                            
                  }
        }
        public void fixDirectionGenotypingRefParallel(DataSet refs, int gapop, int gapext, double distRefPerc) throws InterruptedException, ExecutionException
        {
                 List< Future > futuresList = new ArrayList< Future >();
//                 int nrOfProcessors = Runtime.getRuntime().availableProcessors();
                 int nrOfProcessors = this.avProc;
                 ExecutorService eservice = Executors.newFixedThreadPool(nrOfProcessors);
                 
		 for (Read r : reads)
		 {
                        RevCompGenTask rt = new RevCompGenTask(r,refs,gapop,gapext,distRefPerc);
                        futuresList.add(eservice.submit(rt));
			 
		 }
                  Object taskResult;
                  int i = 0;
                  for(Future future:futuresList) 
                  {
                        i++;
                        taskResult = future.get();
                        System.out.println("RevComp: read "+i + "/" + reads.size());                            
                  }
        }
        void fixDirectionRef(DataSet refs, int gapop, int gapext) throws InterruptedException, ExecutionException, IOException, CompoundNotFoundException
        {
                 int i = 1;
		 for (Read r : reads)
		 {
                    System.out.println("RevComp: read "+i + "/" + reads.size()); 
                    i++;
                     
                    int bestDist = Integer.MAX_VALUE;
                    int toRev = 0;
                    for (Read ref : refs.reads)
                    {
                        int dDir = r.calcEditDistAbsAlignWithGaps(ref, gapop, gapext);
                        int dRev = r.getRevComp().calcEditDistAbsAlignWithGaps(ref, gapop, gapext);
                        bestDist = Math.min(bestDist, dDir);
                        bestDist = Math.min(bestDist, dRev);
                        if (dDir == bestDist)
                            toRev = 0;
                        if (dRev == bestDist)
                            toRev = 1; 
                    }
                    if (toRev == 1)
                        r.RevComp();			 
		 }
        }

        public void PrintReads(String outfile) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		for(Read r : reads)
                {
				for (int i = 0; i < r.frequency; i++)
					fw.write(">" + r.name + "_" + i + "\n" + r.getNucl() + "\n");
		}
		fw.close();
	}
        public void PrintReadsNoCopyNumber(String outfile) throws IOException
	{
		FileWriter fw = new FileWriter(outfile);
		for(Read r : reads)
                {
				for (int i = 0; i < r.frequency; i++)
					fw.write(">" + r.name + "\n" + r.getNucl() + "\n");
		}
		fw.close();
	}
        void fixDirection(String refFile, double gapop, double gapext) throws IOException
        {
            DataSet refds = new DataSet(refFile);
            String s_ref = refds.reads.get(0).getNucl();
            for (Read r : reads)
            {
                System.out.println("Fixing direction: " + r.name + "/" + reads.size());
                int bestmatch = 0;
                int bestdirect = -1;
                String nucl = r.getNucl();
                for (int i = 0; i <= 1; i++)
                {
                    Runtime run=Runtime.getRuntime();
                    Process p=null;
                    FileWriter fw_alin = new FileWriter("allign_input.fas");
                    fw_alin.write(">" + r.name + "\n");
                    fw_alin.write(nucl + "\n");
                    fw_alin.write(">reference" + "\n");
                    fw_alin.write(s_ref + "\n");
                    fw_alin.close();
                    String param = "ClustalW2//clustalw2 -INFILE=" +"allign_input.fas" + " -OUTFILE=" + "allign_output.fas";
                    param+= " -OUTPUT=FASTA -DNAMATRIX=IUB -GAPOPEN=" + gapop +" -GAPEXT=" + gapext + " -TYPE=DNA -PWDNAMATRIX=IUB -PWGAPOPEN=" + gapop+ " -PWGAPEXT=" + gapext;

                    p = run.exec(param);
                    BufferedReader stdInput = new BufferedReader(new InputStreamReader(p.getInputStream()));
                    String s1;
                    while ((s1 = stdInput.readLine()) != null) {
    //                  System.out.println(s1);
                    }
                    if (p.exitValue()!= 0)
                    {
                        System.out.println("Error in allgnment program");
                    }
                    DataSet allignment = new DataSet("allign_output.fas",'c');
                    String al_read = "";
                    String al_ref = "";
                    for (Read r1 : allignment.reads)
                        if (!r1.name.contains("reference"))
                            al_read = r1.getNucl();
                        else
                            al_ref = r1.getNucl();
                    int match = 0;
                    for (int j = 0; j < al_read.length(); j++)
                        if (al_read.charAt(j) == al_ref.charAt(j))
                            match++;
                    if (match > bestmatch)
                    {
                        bestmatch = match;
                        bestdirect = i;
                    }
                    File f = new File("allign_input.fas");
                    f.delete();
                    f = new File("allign_output.fas");
                    f.delete();
                    f = new File("allign_input.dnd");
                    f.delete();
                    nucl = r.getRevComp().getNucl();
                }
                if (bestdirect == 1)
                    r.nucl = r.getRevComp().getNucl();
            }
        }
        void PrintReadsSeparateFiles() throws IOException
        {
            int i = 1;
            for (Read r : reads)
            {
                FileWriter fw = new FileWriter(this.file_name + "_" + r.name + ".fas");
                fw.write(">" + r.name + "\n" + r.nucl + "\n");
                fw.close();
                i++;
            }
            
        }
        void PrintReadsSeparateFilesKeepName(String addr) throws FileNotFoundException, IOException
        // <editor-fold defaultstate="collapsed" desc=" DESCRIPTION ">
    {
        FileReader fr = new FileReader(addr);
        BufferedReader br = new BufferedReader(fr);
        int i = 1;
        String s = br.readLine();
        String nucl = new String();
        String name = s;
        s = br.readLine();
        while (s != null) {
            if (s.length() == 0) {
                s = br.readLine();
                continue;
            }
            if (s.charAt(0) == '>') {
                FileWriter fw = new FileWriter(addr + "_Read" + i + ".fas");
                fw.write(name + "\n" + nucl + "\n");      
                fw.close();
                nucl = "";
                name = s;
                i++;
            } else {
                nucl += s.toUpperCase();
            }
            s = br.readLine();
        }
        FileWriter fw = new FileWriter(addr + "_Read" + i + ".fas");
        fw.write(name + "\n" + nucl + "\n");
        fw.close();
        br.close();
        fr.close();
        
    }
// </editor-fold>
        void printMostFreqReads(String addr, int n) throws IOException
        {
            Collections.sort(reads, new ReadFreqComparator());
            FileWriter fw = new FileWriter(addr);
            for (int i = 1; i <= n; i++)
            {
                fw.write(">Read" + i + "_freq" + reads.get(i).frequency + "\n" + reads.get(i).getNucl() + "\n");
            }
            fw.close();
        }
        public Read findMostFreqRead()
                // <editor-fold defaultstate="collapsed" desc=" DESCRIPTION ">
        {
            Read maxr = new Read();
            double maxfreq = 0;
            for (Read r : this.reads)
                if (r.frequency > maxfreq)
                {
                    maxfreq = r.frequency;
                    maxr = r;
                }
            return maxr;
        }
        
        Haplotype findMostFreqHaplotype()
                // <editor-fold defaultstate="collapsed" desc=" DESCRIPTION ">
        {
            Haplotype maxh = this.haplotypes.get(0);
            double maxconc = 0;
            for (Haplotype h : this.haplotypes)
                if (h.concentration > maxconc)
                {
                    maxconc = h.concentration;
                    maxh = h;
                }
            return maxh;
        }
        
        // </editor-fold>
        void printpairdistReads(String addr) throws IOException
                // <editor-fold defaultstate="collapsed" desc=" DESCRIPTION ">
        {
            FileWriter pairdist = new FileWriter(addr);
            int[][] pd = new int[reads.size()][reads.size()];

            double gapop = 6;
            double gapext = 6;
            for (int i = 0; i < reads.size(); i++)
                for (int j = i+1; j < reads.size(); j++)
                {
                    DynamicOut.printStep("Calculating pair (" + i + "," + j + ") //" + reads.size());
                    Runtime run=Runtime.getRuntime();
                    Process p=null;
                    FileWriter fw_alin = new FileWriter("allign_input.fas");
                    fw_alin.write(">first" + "\n");
                    fw_alin.write(reads.get(i).getNucl() + "\n");
                    fw_alin.write(">second" + "\n");
                    fw_alin.write(reads.get(j).getNucl() + "\n");
                    fw_alin.close();
                    String param = "ClustalW2//clustalw2 -INFILE=" +"allign_input.fas" + " -OUTFILE=" + "allign_output.fas";
                    param+= " -OUTPUT=FASTA -DNAMATRIX=IUB -GAPOPEN=" + gapop +" -GAPEXT=" + gapext + " -TYPE=DNA -PWDNAMATRIX=IUB -PWGAPOPEN=" + gapop+ " -PWGAPEXT=" + gapext;

                    p = run.exec(param);
                    BufferedReader stdInput = new BufferedReader(new InputStreamReader(p.getInputStream()));
                    String s1;
                    while ((s1 = stdInput.readLine()) != null) {
        //              System.out.println(s1);
                    }
                    if (p.exitValue()!= 0)
                    {
                        System.out.println("Error in allgnment program");
                    }
                    DataSet allignment = new DataSet("allign_output.fas",'c');
                    int l = allignment.reads.get(0).getLength();
                    int mismatch = 0;
                    for (int k = 0; k < l; k++)
                        if (allignment.reads.get(0).getNucl().charAt(k) != allignment.reads.get(1).getNucl().charAt(k))
                            mismatch++;
                    pd[i][j] = mismatch;
                    pd[j][i] = mismatch;
                    File f = new File("allign_input.fas");
                    f.delete();
                    f = new File("allign_output.fas");
                    f.delete();
                    f = new File("allign_input.dnd");
                    f.delete();
                }
            DynamicOut.finishSteps();
            for (int i = 0; i < reads.size(); i++)
            {
                for (int j = 0; j < reads.size(); j++)
                    pairdist.write(pd[i][j] + " ");
                pairdist.write("\n");
            }
            pairdist.close();
        }
        // </editor-fold>
        public ArrayList<Read> FindMostFreqHaplotypes(int n) throws IOException
	{
		HaplComp hc = new HaplComp();
		Collections.sort(haplotypes, hc);
		
                ArrayList<Read> ar = new ArrayList<Read>();
		for (int  i = 0; i < n; i++)
		{
                    Haplotype h = haplotypes.get(i);
			ar.add(new Read(h.nucl,h.frequency));
		}
                return ar;
	}
        public ArrayList<Read> FindMostFreqReads(int n) throws IOException
	{
		ReadComp rc = new ReadComp();
		Collections.sort(reads, rc);
		
                ArrayList<Read> ar = new ArrayList<Read>();
		for (int  i = 0; i < n; i++)
		{
                    Read r = reads.get(i);
			ar.add(new Read(r.nucl,r.name,r.frequency));
		}
                return ar;
	}
        public boolean containRead(Read s)
        {
            for (Read r : reads)
                if (r.nucl.equalsIgnoreCase(s.nucl))
                    return true;
            
            return false;
        }
        public String containReadGetName(Read s)
        {
            for (Read r : reads)
                if (r.nucl.equalsIgnoreCase(s.nucl))
                    return r.name;
            
            return "";
        }
         public boolean containReadAlignment(Read s) throws IOException, CompoundNotFoundException
        {
            for (Read r : reads)
                if (r.calcEditDistAbsAlign(s, 15, 6) == 0)
                    return true;
            
            return false;
        }
         public void delGaps()
         {
             for (Read r : this.reads)
             {
//                 System.out.println(r.nucl);
                 String s = "";
                 for (int i = 0; i < r.nucl.length(); i++)
                     if (r.nucl.charAt(i) != '-')
                         s+=r.nucl.charAt(i);
                 r.nucl = s;
//                 System.out.println(r.nucl);
//                 System.out.println();
             }
         }
         public boolean equals(DataSet ds)
         {
             for (Read r : this.reads)
                 if (!ds.containRead(r))
                     return false;
             for (Read r : ds.reads)
                 if (!this.containRead(r))
                     return false;
             return true;
         }
         public int getTotalNReads()
         {
             int Nreads = 0;
             for (Read r : this.reads)
                 Nreads += r.frequency;
             return Nreads;
         }
         public int getFrequency(Read s)
         {
              for (Read r : reads)
                if (r.nucl.equalsIgnoreCase(s.nucl))
                    return r.frequency;
            
            return 0;
         }
          public int getFrequencyAlign(Read s) throws IOException
         {
             int freq = 0;
              for (Read r : reads)
                if (r.calcEditDistAbsAlignIgnoreIns(s) == 0)
                    freq += r.frequency;
             
            
            return freq;
         }
          public int getFrequencyAlign(Read s, int gapop, int gapext) throws IOException, CompoundNotFoundException
         {
             int freq = 0;
              for (Read r : reads)
              {
//                    System.out.println(s.name);
//                    System.out.println(r.name);
                    if (r.calcEditDistAbsAlign(s, gapop, gapext) == 0)
                        freq += r.frequency;
              }
             
            
            return freq;
         }
          public void delPrimers(ArrayList<Read> primers)
          {
              for (Read r : this.reads)
                  for (Read prim : primers)
                    if (r.nucl.endsWith(prim.nucl))
                    {
                        r.nucl = r.nucl.substring(0, r.getLength() - prim.getLength());
                        break;
                    }
          }
         public void clipToRefFillGaps(DataSet refs, int gapop, int gapext) throws IOException, CompoundNotFoundException
        {
                 int i = 1;
		 for (Read r : reads)
		 {
                    System.out.println("Clipping: read "+i + "/" + reads.size()); 
                    i++;
                     
                    int bestDist = Integer.MAX_VALUE;
                    Read bestRef = null;
                    for (Read ref : refs.reads)
                    {
                        int dDir = r.calcEditDistAbsAlignWithGaps(ref, gapop, gapext);
                        if (dDir < bestDist)
                        {
                            bestDist = dDir;
                            bestRef = ref;
                        }
                    }
                    r.clipToRefFillGaps(bestRef, gapop, gapext);
			 
		 }
        }
         public void fillGapsConsensusAligned()
         {
             int l = reads.get(0).getLength();
             char[] nucl = {'A','C','T','G'};
             
             String cons = "";
             for (int i=0; i < l; i++)
             {
                 int[] freqNuclPos = new int[nucl.length];
                 for (Read r : reads)
                 {
                     switch (r.getNucl(i))
                     {
                         case 'A':
                                 freqNuclPos[0] += r.frequency;
                                 break;
                         case 'C':
                                 freqNuclPos[1] += r.frequency;
                                 break;
                         case 'T':
                                 freqNuclPos[2] += r.frequency;
                                 break;
                         case 'G':
                                 freqNuclPos[3] += r.frequency;
                                 break;
                     }                                
                 }
                 int maxfreq = 0;
                 char maxnucl = 'X';
                 for (int j = 0; j < nucl.length; j++)
                     if (freqNuclPos[j] > maxfreq)
                     {
                         maxfreq = freqNuclPos[j];
                         maxnucl = nucl[j];
                     }
                 cons += maxnucl;
             }
             for (Read r : this.reads)
             {
                 String newnucl = "";
                 for (int i = 0; i < r.getLength(); i++)
                 {
                     if (r.getNucl(i) == '-')
                     {
                         if (cons.charAt(i) != 'X')
                            newnucl += cons.charAt(i);
                     }
                     else
                     {
                         newnucl += r.getNucl(i);
                     }
                 }
                 r.nucl = newnucl;
             }
                     
             
         }
         public void clipAligned(int stpos, int endpos)
         {
             HashMap<String,String> names = new HashMap();
             HashMap<String,Integer> seq = new HashMap();
             for (Read r : this.reads)
             {
                 String s = r.nucl.substring(stpos,endpos);
                 if (seq.containsKey(s))
                 {
                     seq.put(s, seq.get(s) + r.frequency);
                 }
                 else
                 {
                     seq.put(s, r.frequency);
                     names.put(s, r.name);
                 }
             }
             ArrayList<Read> newReads = new ArrayList();
             for (Map.Entry me : seq.entrySet())
             {
                 String nucl = (String) me.getKey();
                 int freq = (Integer) me.getValue();
                 String name = names.get(nucl);
                 newReads.add(new Read(nucl,name,freq));
             }
             this.reads = newReads;
         }
         public void sortByFreq()
         {
             Collections.sort(reads, new ReadFreqComparator());
         }
         public HashMap<String,DataSet> separateGenotypes()
         {
             HashMap<String,DataSet> hm = new HashMap();
             for (Read r : this.reads)
             {
                 if (hm.containsKey(r.genotype))
                 {
                     DataSet ds = hm.get(r.genotype);
                     ds.reads.add(r);
                     hm.put(r.genotype, ds);
                 }
                 else
                 {
                     DataSet ds = new DataSet();
                     ds.reads.add(r);
                     hm.put(r.genotype, ds);
                 }
             }
             return hm;
         }
         public void clipToRef(DataSet refs, int gapop, int gapext) throws IOException, CompoundNotFoundException
        {
                 int i = 1;
		 for (Read r : reads)
		 {
                    System.out.println("Clipping: read "+i + "/" + reads.size()); 
                    i++;
                     
                    int bestDist = Integer.MAX_VALUE;
                    Read bestRef = null;
                    for (Read ref : refs.reads)
                    {
                        int dDir = r.calcEditDistAbsAlignWithGaps(ref, gapop, gapext);
                        if (dDir < bestDist)
                        {
                            bestDist = dDir;
                            bestRef = ref;
                        }
                    }
                    r.clipToRef(bestRef, gapop, gapext);			 
		 }
                 
                HashMap<String,String> names = new HashMap();
                HashMap<String,Integer> seq = new HashMap();
                for (Read r : this.reads)
                {
                    String s = r.nucl;
                    if (seq.containsKey(s))
                    {
                        seq.put(s, seq.get(s) + r.frequency);
                    }
                    else
                    {
                        seq.put(s, r.frequency);
                        names.put(s, r.name);
                    }
             }
             ArrayList<Read> newReads = new ArrayList();
             for (Map.Entry me : seq.entrySet())
             {
                 String nucl = (String) me.getKey();
                 int freq = (Integer) me.getValue();
                 String name = names.get(nucl);
                 newReads.add(new Read(nucl,name,freq));
             }
             this.reads = newReads;
        }
         public void clipToRefParallel(DataSet refs, int gapop, int gapext) throws IOException, CompoundNotFoundException, InterruptedException, ExecutionException
        {
                List< Future > futuresList = new ArrayList< Future >();
//                 int nrOfProcessors = Runtime.getRuntime().availableProcessors();
                 int nrOfProcessors = this.avProc;
                 ExecutorService eservice = Executors.newFixedThreadPool(nrOfProcessors);
                 
		 for (Read r : reads)
		 {
                        ClipRefTask rt = new ClipRefTask(r,refs,gapop,gapext);
                        futuresList.add(eservice.submit(rt));
			 
		 }
                  Object taskResult;
                  int i = 0;
                  for(Future future:futuresList) 
                  {
                        i++;
                        taskResult = future.get();
                        System.out.println("Clipping: read "+i + "/" + reads.size());                            
                  }

                 
                HashMap<String,String> names = new HashMap();
                HashMap<String,Integer> seq = new HashMap();
                for (Read r : this.reads)
                {
                    String s = r.nucl;
                    if (seq.containsKey(s))
                    {
                        seq.put(s, seq.get(s) + r.frequency);
                    }
                    else
                    {
                        seq.put(s, r.frequency);
                        names.put(s, r.name);
                    }
             }
             ArrayList<Read> newReads = new ArrayList();
             for (Map.Entry me : seq.entrySet())
             {
                 String nucl = (String) me.getKey();
                 int freq = (Integer) me.getValue();
                 String name = names.get(nucl);
                 newReads.add(new Read(nucl,name,freq));
             }
             this.reads = newReads;
        }
        public  void delShortReads(int cutoff)
        {
            ArrayList<Read> toDel = new ArrayList();
            for (int i = 0; i < reads.size(); i++)
            {
                if (reads.get(i).getLength() < cutoff)
                {
                    toDel.add(reads.get(i));
                }
            }
            reads.removeAll(toDel);
        }
        public  void delNonStandardReadsLength(int len)
        {
            ArrayList<Read> toDel = new ArrayList();
            for (int i = 0; i < reads.size(); i++)
            {
                if (reads.get(i).getLength() < cutoff)
                {
                    toDel.add(reads.get(i));
                }
            }
            reads.removeAll(toDel);
        }
}
