/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
package org.cobi.kgg.ui.action;

import cern.colt.list.DoubleArrayList;
import cern.colt.list.IntArrayList;
import cern.colt.map.OpenIntIntHashMap;
import cern.colt.matrix.DoubleMatrix1D;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.EigenvalueDecomposition;
import cern.jet.stat.Gamma;
import cern.jet.stat.Probability;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.io.File;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.logging.Logger;
import javax.swing.JOptionPane;
import javax.swing.SwingWorker;
import net.sf.picard.liftover.LiftOver;
import org.cobi.kgg.business.GenotypeSetUtil;
import org.cobi.kgg.business.SetBasedTest;
import org.cobi.kgg.business.entity.Chromosome;
import org.cobi.kgg.business.entity.Constants;
import static org.cobi.kgg.business.entity.Constants.CHROM_NAMES;
import org.cobi.kgg.business.entity.CorrelationBasedByteLDSparseMatrix;
import org.cobi.kgg.business.entity.Gene;
import org.cobi.kgg.business.entity.Genome;
import org.cobi.kgg.business.entity.HaplotypeDataset;
import org.cobi.kgg.business.entity.PValueWeight;
import org.cobi.kgg.business.entity.PlinkDataset;
import org.cobi.kgg.business.entity.SNP;
import org.cobi.kgg.business.entity.SNPPosiComparator;
import org.cobi.kgg.business.entity.StatusGtySet;
import org.cobi.kgg.ui.ArrayListObjectArrayTableModel;
import org.cobi.kgg.ui.GlobalManager;
import org.cobi.util.download.stable.HttpClient4API;
import org.cobi.util.file.Zipper;
import org.cobi.util.math.MultipleTestingMethod;
import org.cobi.util.text.Util;
import org.ejml.data.DenseMatrix64F;
import org.netbeans.api.progress.ProgressHandle;
import org.netbeans.api.progress.ProgressHandleFactory;
import org.openide.ErrorManager;
import org.openide.awt.ActionID;
import org.openide.awt.ActionReference;
import org.openide.awt.ActionRegistration;
import org.openide.awt.StatusDisplayer;
import org.openide.util.Cancellable;
import org.openide.util.NbBundle.Messages;
import org.openide.util.RequestProcessor;
import org.rosuda.REngine.Rserve.RConnection;
import umontreal.iro.lecuyer.probdist.ChiSquareDist;

@ActionID(
        category = "Gene",
        id = "org.cobi.kgg.ui.action.ScanConditionalGeneBasedAssociation"
)
@ActionRegistration(
        displayName = "#CTL_ScanConditionalGeneBasedAssociation"
)
@ActionReference(path = "Actions/Gene", position = 550, separatorAfter = 600)
@Messages("CTL_ScanConditionalGeneBasedAssociation=ConditionalAssociation")
public final class ScanConditionalGeneBasedAssociation implements ActionListener, Constants {

    private final static Logger LOG = Logger.getLogger(ScanConditionalGeneBasedAssociation.class.getName());
    private final static RequestProcessor RP = new RequestProcessor("Scan conditional gene", 1, true);
    private RequestProcessor.Task buildTask = null;

    private Genome genome;
    private Map<String, int[]> geneGenomeIndexes;

    private IntArrayList chromIDs;
    private int snpPIndex;
    private List<Object[]> geneTable;

    @Override
    public void actionPerformed(ActionEvent e) {
        // TODO implement action body

    }
//Genome genome, Map<String, int[]> geneGenomeIndexes,  int snpPIndex

    public Genome getGenome() {
        return genome;
    }

    public void setGenome(Genome genome) {
        this.genome = genome;
    }

    public int getSnpPIndex() {
        return snpPIndex;
    }

    public void setSnpPIndex(int snpPIndex) {
        this.snpPIndex = snpPIndex;
    }

    public ScanConditionalGeneBasedAssociation(List<Object[]> geneTable) {
        this.geneTable = geneTable;
    }

    public Map<String, int[]> getGeneGenomeIndexes() {
        return geneGenomeIndexes;
    }

    public void setGeneGenomeIndexes(Map<String, int[]> geneGenomeIndexes) {
        this.geneGenomeIndexes = geneGenomeIndexes;
    }

    public void conditionScan(ArrayListObjectArrayTableModel listPValueTableModel) {
        //record the classification settings
        ScanConditionalGeneAssocSwingWorker worker = new ScanConditionalGeneAssocSwingWorker(listPValueTableModel);
        // buildTask = buildTask.create(buildingThread); //the task is not started yet
        buildTask = RP.create(worker); //the task is not started yet
        buildTask.schedule(0); //start the task
    }

    private boolean handleCancel() {
        if (null == buildTask) {
            return false;
        }
        return buildTask.cancel();
    }

    public void setChromIDs(IntArrayList chromIDs) {
        this.chromIDs = chromIDs;
    }

    class ScanConditionalGeneAssocSwingWorker extends SwingWorker<Void, String> {

        int runningThread = 0;
        boolean succeed = false;
        ProgressHandle ph = null;
        long time;
        ArrayListObjectArrayTableModel listPValueTableModel;

        public ScanConditionalGeneAssocSwingWorker(ArrayListObjectArrayTableModel listPValueTableModel) {
            ph = ProgressHandleFactory.createHandle("task thats shows progress", new Cancellable() {
                @Override
                public boolean cancel() {
                    return handleCancel();
                }
            }); 
            this.listPValueTableModel=listPValueTableModel;
            time = System.nanoTime();
        }

        void condtionalTestOnChrom(int currChromIndex, List<Object[]> geneTableTmp) throws Exception {
            Set<Integer> positionsSet = new HashSet<Integer>();
            Chromosome currChrom = genome.readChromosomefromDisk(currChromIndex);

            Map<Double, Integer> orderMap = new HashMap<Double, Integer>();
            DoubleArrayList rankList = new DoubleArrayList();
            Set<Integer> groupID = new HashSet<Integer>();
            double order;
            String ordeValue;
            //pickup positions of  SNPs
            for (int i = 0; i < geneTableTmp.size(); i++) {
                ordeValue = geneTableTmp.get(i)[6].toString();
                if (ordeValue.equals("-")) {
                    continue;
                }
                order = Double.parseDouble(ordeValue);
                
                while (orderMap.containsKey(order)) {
                    order+=0.000001;
                }
                rankList.add(order);
                orderMap.put(order, i);

                String geneSymb = geneTableTmp.get(i)[1].toString();
                int[] indexes1 = geneGenomeIndexes.get(geneSymb);
                if (indexes1 == null) {
                    continue;
                }
                groupID.add(Integer.parseInt(geneTableTmp.get(i)[0].toString()));
                List<SNP> snps = currChrom.genes.get(indexes1[1]).snps;
                for (SNP snp : snps) {
                    positionsSet.add(snp.physicalPosition);
                }
            }

            File CHAIN_FILE = null;
            LiftOver liftOverLDGenome2pValueFile = null;
            Zipper ziper = new Zipper();
            String pValueFileGenomeVersion = genome.getFinalBuildGenomeVersion();

            Genome ldGenome;
            if (genome.getLdSourceCode() == -2) {
                ldGenome = GlobalManager.currentProject.getGenomeByName(genome.getSameLDGenome());
            } else {
                ldGenome = genome;
            }
            String info;
            String ldFileGenomeVersion = ldGenome.getLdFileGenomeVersion();

            try {
                if (ldFileGenomeVersion != null && !pValueFileGenomeVersion.equals(ldFileGenomeVersion)) {
                    CHAIN_FILE = new File(GlobalManager.RESOURCE_PATH + "liftOver/" + ldFileGenomeVersion + "ToH" + pValueFileGenomeVersion.substring(1) + ".over.chain.gz");
                    if (!CHAIN_FILE.exists()) {
                        if (!CHAIN_FILE.getParentFile().exists()) {
                            CHAIN_FILE.getParentFile().mkdirs();
                        }
                        String url = "http://hgdownload.cse.ucsc.edu/goldenPath/" + ldFileGenomeVersion + "/liftOver/" + CHAIN_FILE.getName();
                        //HttpClient4API.downloadAFile(url, CHAIN_FILE);
                        HttpClient4API.simpleRetriever(url, CHAIN_FILE.getCanonicalPath(), GlobalManager.proxyBean);
                        ziper.extractTarGz(CHAIN_FILE.getCanonicalPath(), CHAIN_FILE.getParent());
                    }
                    CHAIN_FILE = new File(GlobalManager.RESOURCE_PATH + "liftOver/" + ldFileGenomeVersion + "ToH" + pValueFileGenomeVersion.substring(1) + ".over.chain");
                    liftOverLDGenome2pValueFile = new LiftOver(CHAIN_FILE);
                }

                CorrelationBasedByteLDSparseMatrix ldRsMatrixes = null;
                info = "Reading LD information of Key SNPs";
                LOG.info(info);

                if (positionsSet == null || positionsSet.isEmpty()) {
                    return;
                }
                //ld source code
                //-2 others LD
                //0 genotype plink binary file
                //1 hapap ld
                //2 1kG haplomap
                //3 local LD calcualted by plink
                //4 1kG haplomap vcf format

                if (ldGenome.getLdSourceCode() == 0) {
                    if (ldGenome.getPlinkSet() == null || !ldGenome.getPlinkSet().avaibleFiles()) {
                        String infor = "No Plink binary file to account for LD for SNPs! ";
                        LOG.info(infor);
                        return;
                    }
                    ldRsMatrixes = calculateLocalLDRSquarebyGenotypesPositions(ldGenome.getPlinkSet(), positionsSet, CHROM_NAMES[currChromIndex], liftOverLDGenome2pValueFile);

                } else if (ldGenome.getLdSourceCode() == 4) {
                    if (ldGenome.getChromLDFiles()[currChromIndex] != null) {
                        File hapMapLDFile = new File(ldGenome.getChromLDFiles()[currChromIndex]);
                        if (!hapMapLDFile.exists()) {
                            return;
                        }
                    }
                    ldRsMatrixes = calculateLocalLDRSquarebyHaplotypeVCFByPositions(ldGenome.getChromLDFiles()[currChromIndex], CHROM_NAMES[currChromIndex], positionsSet, liftOverLDGenome2pValueFile);
                }

                boolean ingoreSNPNoLD = true;
                float geneWeight = 1;
                List<SNP> allSnpList = new ArrayList<SNP>();
                List<SNP> allLastSnpList = new ArrayList<SNP>();
                List<SNP> tmpSnpList = new ArrayList<SNP>();
                double[] tmpResult = new double[2];
                double[] accuResult = new double[2];
                Arrays.fill(accuResult, 0);
                int snpPVTypeIndex = 0;
                double df = 0;
                double Y = 0;
                boolean usingSquareChi = false;
                double p1 = 0;
                DoubleArrayList pValueArray = new DoubleArrayList();
                RConnection rcon = null;
                if (usingSquareChi) {
                    rcon = new RConnection();
                    rcon.eval("pack=\"survey\"; if (!require(pack,character.only = TRUE)) { install.packages(pack,dep=TRUE,repos='http://cran.us.r-project.org');if(!require(pack,character.only = TRUE)) stop(\"Package not found\")}");
                    rcon.eval("library(survey)");
                }

                double lastChisquare = 0;
                double currChisquare = 0;
                rankList.quickSort();

                double orgP;
                int rankSize = rankList.size();
                List<Integer> groupIDs = new ArrayList<Integer>(groupID);
                Collections.sort(groupIDs);
                int groupSize = groupIDs.size();
                String groupIDStr;
                int addedGeneNum = 0;
                double factor = 1.25,chi;
                for (int g = 0; g < groupSize; g++) {
                    allSnpList.clear();
                    allLastSnpList.clear();
                    groupIDStr = String.valueOf(groupIDs.get(g));
                    addedGeneNum = 0;
                    Arrays.fill(accuResult, 0);
                    for (int i = 0; i < rankSize; i++) {
                        int index = orderMap.get(rankList.getQuick(i));
                      
                        if (!geneTableTmp.get(index)[0].toString().equals(groupIDStr)) {
                            continue;
                        }
                        String geneSymb = geneTableTmp.get(index)[1].toString();
                         
                        int[] indexes1 = geneGenomeIndexes.get(geneSymb);
                        if (indexes1 == null) {
                            continue;
                        }
                        Gene gene = currChrom.genes.get(indexes1[1]);
                        addedGeneNum++;
                        allSnpList.addAll(gene.snps);
                        Collections.sort(allSnpList, new SNPPosiComparator());
                        if (usingSquareChi) {
                            currChisquare = snpSetPValuebyJohnnyChiSquare(allSnpList, pValueArray, ldRsMatrixes, snpPVTypeIndex, false, rcon);
                            if (currChisquare < 1E-0) {
                                currChisquare = MultipleTestingMethod.iterativeChisquareInverse(2, currChisquare);
                            } else {
                                currChisquare = ChiSquareDist.inverseF(2, 1 - currChisquare);
                            }
                            p1 = Probability.chiSquareComplemented(1, (currChisquare - lastChisquare) <= 0 ? 0 : currChisquare - lastChisquare);
                            lastChisquare = currChisquare;
                        } else {
                            boolean toPrint=false;                   
                            snpSetPValuebyMyChiSquareApproxEJML(allSnpList, ldRsMatrixes, snpPVTypeIndex, toPrint, tmpResult);
                            df = tmpResult[0] - accuResult[0];
                            Y = tmpResult[1] - accuResult[1];

                            if (Y < 0) {
                                Y = 0;
                            }
                            if (df < 1) {
                                df = 1;
                            }
                            // System.out.println(allSnpList.size() + "\t" + gene.getSymbol());
                            p1 = Gamma.incompleteGammaComplement(df / 2, Y / 2);
                            accuResult[0] = tmpResult[0];
                            accuResult[1] = tmpResult[1];
                        }
                        tmpSnpList.clear();

                        if (addedGeneNum > 1) {
                            tmpSnpList.addAll(gene.snps);
                            Collections.sort(tmpSnpList, new SNPPosiComparator());
                            orgP = snpSetPValuebyMyChiSquareApproxEJML(tmpSnpList, ldRsMatrixes, snpPVTypeIndex, false, tmpResult);
                             
                            if (p1*5>orgP) {
                                //if it is reasonable
                                allLastSnpList.addAll(tmpSnpList);                               
                                //because the conditional test has some inflation, we need an inflation factor to correct it
                                 chi = MultipleTestingMethod.zScore(p1 / 2);
                                 chi=chi*chi;
                                 chi = chi / factor;                                     
                                 p1= Gamma.incompleteGammaComplement(0.5, chi / 2);     
                                geneTableTmp.get(index)[8] = Util.formatPValue(p1);
                            } else {
                                geneTableTmp.get(index)[8] = "-";
                            }
                            allSnpList.clear();
                            allSnpList.addAll(allLastSnpList);
                        } else {
                            allLastSnpList.addAll(allSnpList);
                            geneTableTmp.get(index)[8] = Util.formatPValue(p1);
                        }
                    }
                }

                if (usingSquareChi) {
                    rcon.close();
                }
            } catch (Exception ex) {
                ex.printStackTrace();
            }
        }

        @Override
        protected Void doInBackground() {
            try {

                long startTime = System.currentTimeMillis();

                String inforString = "Scan conditional gene-based association on the genome ...";

                StatusDisplayer.getDefault().setStatusText(inforString);

//pathwayBasedAssociation.getpValueSources();
                ph.start(); //we must start the PH before we swith to determinate
                ph.switchToIndeterminate();

                List<Object[]> geneTableTmp = new ArrayList<Object[]>();
                int size = geneTable.size();
                if (size < 1) {
                    return null;
                }
                int chrSize = chromIDs.size();

                for (int t = 0; t < chrSize; t++) {
                    int currChromIndex = chromIDs.get(t);
                    String chrC = CHROM_NAMES[currChromIndex];
                    geneTableTmp.clear();
                    for (int i = 0; i < size; i++) {
                        String chr = geneTable.get(i)[2].toString();
                        boolean sel = (Boolean) geneTable.get(i)[7];
                        if (!sel) {
                            continue;
                        }
                        if (chrC.equals(chr)) {
                            geneTableTmp.add(geneTable.get(i));
                        }
                    }
                    if (geneTableTmp.isEmpty()) {
                        continue;
                    }

                    condtionalTestOnChrom(currChromIndex, geneTableTmp);
                    System.gc();
                    publish("Chromsom " + CHROM_NAMES[currChromIndex] + " is processed!");
                }

                //GlobalManager.mainView.displayPathwayTree(pathwayTreeRoot);
                // InterfaceUtil.saveTreeNode2XMLFile(pathwayTreeRoot, storagePath);
                System.gc();

                succeed = true;
            } catch (Exception ex) {
                StatusDisplayer.getDefault().setStatusText("Scan conditional gene-based association task was CANCELLED!");
                java.util.logging.Logger.getLogger(BuildGenome.class.getName()).log(java.util.logging.Level.SEVERE, null, ex);
            }

            return null;
        }

        @Override
        protected void process(List<String> chunks) {
            // TODO Auto-generated method stub  
            for (String message : chunks) {
                LOG.info(message);
                StatusDisplayer.getDefault().setStatusText(message);
            }

        }

        @Override
        protected void done() {
            try {
                String message;
                if (!succeed) {
                    message = ("Scan conditional gene-based association failed!");
                    LOG.info(message);
                    return;
                }

                message = ("Conditional gene-based association scan has been finished!");
                LOG.info(message);
                StatusDisplayer.getDefault().setStatusText(message);
                ph.finish();

                String prjName = GlobalManager.currentProject.getName();
                String prjWorkingPath = GlobalManager.currentProject.getWorkingPath();

                time = System.nanoTime() - time;
                time = time / 1000000000;
                long min = time / 60;
                long sec = time % 60;
                String info = ("Elapsed time: " + min + " min. " + sec + " sec.");
                listPValueTableModel.fireTableDataChanged();
                JOptionPane.showMessageDialog(null, "The conditional gene-based test has been finished!", "Message", JOptionPane.INFORMATION_MESSAGE);
            } catch (Exception e) {
                ErrorManager.getDefault().notify(e);
            }
        }
    }

    // this function attempts to calculate SNP LD only within genes because it is too slow to do others
    public CorrelationBasedByteLDSparseMatrix calculateLocalLDRSquarebyGenotypesPositions(PlinkDataset plinkSet, Set<Integer> positionsSet,
            String chromName, LiftOver liftOver) throws Exception {
        if (positionsSet == null || positionsSet.isEmpty()) {
            return null;
        }
        if (!plinkSet.avaibleFiles()) {
            String infor = "No Plink binary file on chromosome " + chromName;

            LOG.info(infor);
            return null;
        }
        Map<String, StatusGtySet> indivGtyMap = new HashMap<String, StatusGtySet>();
        List<SNP> mappedSNP = plinkSet.readSNPsinPlinkBinaryMapFileByPositions(positionsSet, chromName, liftOver);
        IntArrayList snpPositionList = new IntArrayList();
        snpPositionList.addAllOf(positionsSet);
        snpPositionList.quickSort();
        int ldSNPNum = snpPositionList.size();
        OpenIntIntHashMap allIndexMap = new OpenIntIntHashMap(ldSNPNum);
        for (int i = 0; i < ldSNPNum; i++) {
            allIndexMap.put(snpPositionList.getQuick(i), i);
        }
        snpPositionList.clear();

        CorrelationBasedByteLDSparseMatrix ldRsMatrix = new CorrelationBasedByteLDSparseMatrix(allIndexMap);

        String infor = "Reading local genotypes on chromosome " + chromName + "...";
        plinkSet.readPlinkBinaryGenotypeinPedigreeFile(mappedSNP, indivGtyMap);
        LOG.info(infor);
        List<StatusGtySet> chromGtys = new ArrayList<StatusGtySet>(indivGtyMap.values());
        infor = "Calculating local pair-wise LD of SNP within genes on chromosome " + chromName + "...";

        LOG.info(infor);
        if (!GenotypeSetUtil.calculateGenotypeCorrelationSquareFastBit(mappedSNP, chromGtys, ldRsMatrix, genome.getMinEffectiveR2())) {
            infor = "Using slow function due to mssing genotypes!!!";

            LOG.info(infor);
            GenotypeSetUtil.calculateGenotypeCorrelationSquareFast(mappedSNP, chromGtys, ldRsMatrix, genome.getMinEffectiveR2());
        }

        return ldRsMatrix;
    }

    // this function attempts to calculate SNP LD only within genes because it is too slow to do others
    public CorrelationBasedByteLDSparseMatrix calculateLocalLDRSquarebyHaplotypeVCFByPositions(
            String vcfFilePath, String chromName, Set<Integer> positionsSet, LiftOver liftOver) throws Exception {

        if (vcfFilePath == null) {
            return null;
        }
        File vcfFile = new File(vcfFilePath);

        HaplotypeDataset haplotypeDataset = new HaplotypeDataset();
        if (!vcfFile.exists()) {
            String infor = "The VCF file for chromosome " + chromName + " does not exist!";

            LOG.info(infor);
            return null;
        }

        String infor = "Reading haplotpyes on chromosome " + chromName + "...";

        LOG.info(infor);

        List<StatusGtySet> chromGtys = new ArrayList< StatusGtySet>();
        List<SNP> mappedSNP = haplotypeDataset.readSNPMapHapVCFFileByPositions(vcfFile, chromName, positionsSet, liftOver, chromGtys);

        Set<Integer> addedPositionSet = new HashSet<Integer>();
        for (SNP snp : mappedSNP) {
            addedPositionSet.add(snp.physicalPosition);
        }
        //count SNPs have not genotype information
        StringBuilder sb = new StringBuilder();
        int unmappedSNPinGeneNum = 0;
        int totalSNPinGeneNum = positionsSet.size();
        IntArrayList snpPositionList = new IntArrayList();

        sb.append("Position\n");
        for (Integer pos1 : positionsSet) {
            if (!addedPositionSet.contains(pos1)) {
                sb.append(pos1);
                sb.append('\n');
                unmappedSNPinGeneNum++;
            } else {
                snpPositionList.add(pos1);
            }
        }
        /*
         if (unmappedSNPinGeneNum > 0) {
         String prjName = GlobalManager.currentProject.getName();
         String prjWorkingPath = GlobalManager.currentProject.getWorkingPath();
         File outPath = new File(prjWorkingPath + File.separator + prjName + File.separator + geneScan.getName() + ".NoLDSNPs." + chromName + ".txt");
         BufferedWriter bw = new BufferedWriter(new FileWriter(outPath));
         bw.write(sb.toString());
         bw.close();

         sb.delete(0, sb.length());
         infor = "Warning!!! " + unmappedSNPinGeneNum + " variants within genes out of " + totalSNPinGeneNum + " on chromosome " + chromName + " have NO haplotype information and will be assumed to be independent of others! Detailed information of these variants is saved in " + outPath.getCanonicalPath() + ".";
            
         LOG.info(infor);
         }
         */
        infor = "Calculating local pair-wise LD of SNP within genes on chromosome " + chromName + " ...";

        LOG.info(infor);
        snpPositionList.quickSort();
        int ldSNPNum = snpPositionList.size();
        OpenIntIntHashMap allIndexMap = new OpenIntIntHashMap(ldSNPNum);
        for (int i = 0; i < ldSNPNum; i++) {
            allIndexMap.put(snpPositionList.getQuick(i), i);
        }
        snpPositionList.clear();
        CorrelationBasedByteLDSparseMatrix ldRsMatrix = new CorrelationBasedByteLDSparseMatrix(allIndexMap);
        //GenotypeSetUtil.calculateGenotypeCorrelationSquare(mappedSNP, chromGtys, ldRsMatrix);

        if (!GenotypeSetUtil.calculateLDRSequareByHaplotypeFast(mappedSNP, chromGtys, ldRsMatrix, genome.getMinEffectiveR2())) {
            infor = "Using slow function due to mssing genotypes!!!";

            LOG.info(infor);
            GenotypeSetUtil.calculateLDRSequareByHaplotype(mappedSNP, chromGtys, ldRsMatrix, genome.getMinEffectiveR2());
        }
        mappedSNP.clear();
        chromGtys.clear();
        addedPositionSet.clear();

        return ldRsMatrix;
    }

    public String LDPruning(List<SNP> mainSnpMap, CorrelationBasedByteLDSparseMatrix ldCorr, double maxCorr, boolean ingoreNOGty) throws Exception {
        List<SNP> tmpSNPMap = new ArrayList<SNP>();
        tmpSNPMap.addAll(mainSnpMap);
        mainSnpMap.clear();
        int listSize = tmpSNPMap.size();
        Set<Integer> highlyCorrIndexes = new HashSet<Integer>();
        int windowSize = 50;
        int stepLen = 5;
        double r, c;
        int[] counts = new int[2];
        for (int s = 0; s < listSize; s += stepLen) {
            for (int i = s; (i - s <= windowSize) && (i < listSize); i++) {
                SNP snp1 = tmpSNPMap.get(i);
                if (ingoreNOGty) {
                    if (snp1.genotypeOrder < 0) {
                        highlyCorrIndexes.add(i);
                        continue;
                    }
                }
                if (highlyCorrIndexes.contains(i)) {
                    continue;
                }
                for (int j = i + 1; (j - i <= windowSize) && (j < listSize); j++) {
                    if (highlyCorrIndexes.contains(j)) {
                        continue;
                    }
                    SNP snp2 = tmpSNPMap.get(j);
                    if (ingoreNOGty) {
                        if (snp2.genotypeOrder < 0) {
                            highlyCorrIndexes.add(j);
                            continue;
                        }
                    }

                    r = ldCorr.getLDAt(snp1.physicalPosition, snp2.physicalPosition);
                    /*
                     r = Math.sqrt(ldRMatrix.getQuick(t, j));
                     //for R
                     c = (0.6065 * r - 1.033) * r + 1.7351;
                     if (c > 2) {
                     c = 2;
                     }
                     */

                    //R2 
                    //y = -35.741x6 + 111.16x5 - 128.42x4 + 66.906x3 - 14.641x2 + 0.6075x + 0.8596
                    //c = (((((-35.741 * r + 111.16) * r - 128.42) * r + 66.906) * r - 14.641) * r + 0.6075) * r + 0.8596;
                    //y = 0.2725x2 - 0.3759x + 0.8508
                    //c = (0.2725 * r - 0.3759) * r + 0.8508;
                    // y = 0.2814x2 - 0.4308x + 0.86
                    //c = (0.2814 * r - 0.4308) * r + 0.86;
                    //y = -0.155x + 0.8172
                    //c = -0.155 * r + 0.8172;
                   // c = 0.9;

                   // r = Math.pow(r, c);

                    if (r >= maxCorr) {
                        highlyCorrIndexes.add(j);
                    }
                }
            }
        }

        counts[0] = listSize - highlyCorrIndexes.size();
        counts[1] = listSize;
        String info = (listSize - (highlyCorrIndexes.size()) + " SNPs (out of " + listSize + ") passed LD pruning (r2>=" + maxCorr + ").");

        for (int s = 0; s < listSize; s++) {
            SNP snp1 = tmpSNPMap.get(s);
            if (!highlyCorrIndexes.contains(s)) {
                mainSnpMap.add(snp1);
            }
        }
        return info;
    }

    //use a faster matrix invers function DenseMatrix64F
    private double snpSetPValuebyMyChiSquareApproxEJML(List<SNP> snpList, CorrelationBasedByteLDSparseMatrix ldCorr, int snpPVTypeIndex,
            boolean printMatirx, double[] results) throws Exception {
        boolean ignoreNoLDSNP = true;
        int snpNum = snpList.size();
        Set<Integer> highlyCorrIndexes = new HashSet<Integer>();
        double maxCorr = 0.98;
        double r, c;
        LDPruning(snpList, ldCorr, maxCorr, ignoreNoLDSNP);

        double p1 = Double.NaN;
        //CALB2

        List<PValueWeight> pvalueWeightList = new ArrayList<PValueWeight>();
        snpNum = snpList.size();
        //here I think the only uesfulness of the pvalueWeightList is to filter out the null p1-value SNPs
        for (int k = 0; k < snpNum; k++) {
            SNP snp = snpList.get(k);
            if (ignoreNoLDSNP) {
                if (snp.genotypeOrder < 0) {
                    continue;
                }
            }
            double[] pValues = snp.getpValues();
            if (pValues == null) {
                continue;
            }

            if (!Double.isNaN(pValues[snpPVTypeIndex])) //
            {
                PValueWeight pv = new PValueWeight();
                pv.physicalPos = snp.physicalPosition;
                pv.pValue = pValues[snpPVTypeIndex];
                pv.chiSquare = pv.pValue / 2;
                pv.chiSquare = MultipleTestingMethod.zScore(pv.chiSquare);
                pv.chiSquare = pv.chiSquare * pv.chiSquare;

                pv.weight = 1;
                pvalueWeightList.add(pv);
            }
        }

        snpNum = pvalueWeightList.size();

        if (snpNum == 0) {
            results[0] = 1;
            results[1] = pvalueWeightList.get(0).chiSquare;
            return Double.NaN;
        } else if (snpNum == 1) {
            results[0] = 1;
            results[1] = pvalueWeightList.get(0).chiSquare;
            return pvalueWeightList.get(0).pValue;
        }
//min (Ax - b)'(Ax-b) = min x'A'Ax + b'Ax
//s.t. x > 0 
//Can therefore consider the packages quadprog and nnls on CRAN.

        DoubleMatrix2D ldRMatrix = new DenseDoubleMatrix2D(snpNum, snpNum);

        for (int i = 0; i < snpNum; i++) {
            PValueWeight pv = pvalueWeightList.get(i);
            ldRMatrix.setQuick(i, i, pv.weight * pv.weight);
            for (int j = i + 1; j < snpNum; j++) {
                //Math.sqrt(pValueArray[t].var * pValueArray[j].var) 
                r = 0;
                ldRMatrix.setQuick(i, j, pv.weight * pvalueWeightList.get(j).weight * ldCorr.getLDAt(pv.physicalPos, pvalueWeightList.get(j).physicalPos));
                if (ldRMatrix.getQuick(i, j) > 0) {
                    /*
                     r = Math.sqrt(ldRMatrix.getQuick(t, j));
                     //for R
                     c = (0.6065 * r - 1.033) * r + 1.7351;
                     if (c > 2) {
                     c = 2;
                     }
                     */

                    r = (ldRMatrix.getQuick(i, j));
                    //R2
                     //R2
                    //y = -35.741x6 + 111.16x5 - 128.42x4 + 66.906x3 - 14.641x2 + 0.6075x + 0.8596
                    //c = (((((-35.741 * r + 111.16) * r - 128.42) * r + 66.906) * r - 14.641) * r + 0.6075) * r + 0.8596;
                    //y = 0.2725x2 - 0.3759x + 0.8508
                    //c = (0.2725 * r - 0.3759) * r + 0.8508;
                    // y = 0.2814x2 - 0.4308x + 0.86
                    //c = (0.2814 * r - 0.4308) * r + 0.86;
                    //y = -0.155x + 0.8172
                    //c = -0.155 * r + 0.8172;
                   // c = 0.9;

                   // r = Math.pow(r, c);
                    ldRMatrix.setQuick(i, j, r);
                } else {
                    ldRMatrix.setQuick(i, j, 0);
                }
                ldRMatrix.setQuick(j, i, r);
            }
        }

        //remove redundant SNPs according to LD
        int originalSampleSize = snpNum;
        maxCorr = 0.98;
        for (int i = 0; i < originalSampleSize; i++) {
            for (int j = i + 1; j < originalSampleSize; j++) {
                if (Math.abs(ldRMatrix.getQuick(i, j)) >= maxCorr) {
                    if (!highlyCorrIndexes.contains(j) && !highlyCorrIndexes.contains(i)) {
                        highlyCorrIndexes.add(j);
                        // System.out.println(t + " <-> " + j);
                    }
                }
            }
        }

        if (highlyCorrIndexes.size() > 0) {
            // System.out.println("Removed columns and rows: " + highlyCorrIndexes.toString());
            snpNum = originalSampleSize - highlyCorrIndexes.size();
            List<PValueWeight> tmpPvalueWeightList = new ArrayList<PValueWeight>(pvalueWeightList);
            pvalueWeightList.clear();
            DoubleMatrix2D tmpCorMat = new DenseDoubleMatrix2D(snpNum, snpNum);
            int incRow = 0;
            int incCol = 0;
            for (int i = 0; i < originalSampleSize; i++) {
                if (highlyCorrIndexes.contains(i)) {
                    continue;
                }
                pvalueWeightList.add(tmpPvalueWeightList.get(i));
                incCol = 0;
                for (int j = 0; j < originalSampleSize; j++) {
                    if (highlyCorrIndexes.contains(j)) {
                        continue;
                    }
                    tmpCorMat.setQuick(incRow, incCol, ldRMatrix.getQuick(i, j));
                    incCol++;
                }
                incRow++;
            }
            ldRMatrix = tmpCorMat;
            tmpPvalueWeightList.clear();
            // System.out.println(corrMat.toString());
        }
        if (snpNum == 1) {
            results[0] = 1;
            results[1] = pvalueWeightList.get(0).chiSquare;
            return pvalueWeightList.get(0).pValue;
        }
//System.out.println(gene.+"\t"+snpNum);
         

        /*
         do {
         blockCluster = new MarixDensity().getCluster(ldRMatrix, threshold, maxBlockSize, scale, pairThreshold);
         threshold += 0.05;
         pairThreshold += 0.05;
         //  snpNum < 500 ||
         } while (blockCluster.size() > 3);
         */
        double df = 0;
        double Y = 0;
        double dft = 0;
        double Yt = 0;

    
        SetBasedTest sbt = new SetBasedTest();
        double[] result = new double[2];
        if(printMatirx){
            System.out.println(ldRMatrix.toString());
        }
     
            DenseMatrix64F A = new DenseMatrix64F(snpNum, snpNum);
         
            DenseMatrix64F b1 = new DenseMatrix64F(snpNum, 1);
            for (int i = 0; i < snpNum; i++) {
                A.set(i, i, 1);
                for (int j = i + 1; j < snpNum; j++) {
                    A.set(i, j, ldRMatrix.getQuick(i, j));
                    A.set(j, i, ldRMatrix.getQuick(i, j));
                }
            }

            for (int k = 0; k < snpNum; k++) {
                b1.set(k, 0, pvalueWeightList.get(k).chiSquare);
            }

           
            DenseMatrix64F b2 = new DenseMatrix64F(snpNum, 1);
            for (int i = 0; i < snpNum; i++) {
                b2.set(i, 0, 1);
            }

            sbt.mySudoSVDSolverOverlappedWindow(A, b1, b2, result);
            dft += result[0];
            Yt += result[1];
        
        //   System.out.println(ldRMatrix.toString());

        p1 = Gamma.incompleteGammaComplement(dft / 2, Yt / 2);

        results[0] = dft;
        results[1] = Yt;
        return p1;
    }

    //use a faster matrix invers function DenseMatrix64F
    private double snpSetPValuebyJohnnyChiSquare(List<SNP> snpList, DoubleArrayList pValueArray, CorrelationBasedByteLDSparseMatrix ldCorr, int snpPVTypeIndex,
            boolean needWeight, RConnection rcon) throws Exception {

        int snpNum = snpList.size();
        // System.out.println(gene.getSymbol() + " " + snpNum);
        //because it is very time consumming, when SNP number is over 200. I try to spl

        double p1 = Double.NaN;

        for (int k = 0; k < snpNum; k++) {
            SNP snp = snpList.get(k);
            double[] pValues = snp.getpValues();
            if (pValues == null) {
                continue;
            }
            if (!Double.isNaN(pValues[snpPVTypeIndex])) //
            {
                pValueArray.add(pValues[snpPVTypeIndex]);
            }
        }

        int[] keySNPPosition = new int[1];
        keySNPPosition[0] = -1;
        List<PValueWeight> pvalueWeightList = new ArrayList<PValueWeight>();

        boolean ignoreNoLDSNP = true;

        //here I think the only uesfulness of the pvalueWeightList is to filter out the null p1-value SNPs
        for (int k = 0; k < snpNum; k++) {
            SNP snp = snpList.get(k);
            if (ignoreNoLDSNP) {
                if (snp.genotypeOrder < 0) {
                    continue;
                }
            }
            double[] pValues = snp.getpValues();
            if (pValues == null) {
                continue;
            }

            if (!Double.isNaN(pValues[snpPVTypeIndex])) //
            {
                PValueWeight pv = new PValueWeight();
                pv.physicalPos = snp.physicalPosition;
                pv.pValue = pValues[snpPVTypeIndex];
                pv.chiSquare = pv.pValue / 2;
                pv.chiSquare = MultipleTestingMethod.zScore(pv.chiSquare);
                pv.chiSquare = pv.chiSquare * pv.chiSquare;

                pv.weight = 1;
                pvalueWeightList.add(pv);
            }
        }

        snpNum = pvalueWeightList.size();

        if (snpNum == 0) {
            keySNPPosition[0] = -1;
            return Double.NaN;
        } else if (snpNum == 1) {
            keySNPPosition[0] = pvalueWeightList.get(0).physicalPos;
            return pvalueWeightList.get(0).pValue;
        }

        DoubleMatrix2D ldRMatrix = new DenseDoubleMatrix2D(snpNum, snpNum);

        for (int i = 0; i < snpNum; i++) {
            PValueWeight pv = pvalueWeightList.get(i);
            ldRMatrix.setQuick(i, i, pv.weight * pv.weight);
            for (int j = i + 1; j < snpNum; j++) {
                //Math.sqrt(pValueArray[t].var * pValueArray[j].var) 
                ldRMatrix.setQuick(i, j, pv.weight * pvalueWeightList.get(j).weight * ldCorr.getLDAt(pv.physicalPos, pvalueWeightList.get(j).physicalPos));
                if (ldRMatrix.getQuick(i, j) > 0) {
                    ldRMatrix.setQuick(i, j, Math.sqrt(ldRMatrix.getQuick(i, j)));
                } else {
                    ldRMatrix.setQuick(i, j, 0);
                }
                ldRMatrix.setQuick(j, i, ldRMatrix.getQuick(i, j));
            }
        }
//remove redundant SNPs according to LD
        int originalSampleSize = snpNum;
        Set<Integer> highlyCorrIndexes = new HashSet<Integer>();
        double maxCorr = 0.99;
        for (int i = 0; i < originalSampleSize; i++) {
            for (int j = i + 1; j < originalSampleSize; j++) {
                if (Math.abs(ldRMatrix.getQuick(i, j)) >= maxCorr) {
                    if (!highlyCorrIndexes.contains(j) && !highlyCorrIndexes.contains(i)) {
                        highlyCorrIndexes.add(j);
                        // System.out.println(t + " <-> " + j);
                    }
                }
            }
        }

        if (highlyCorrIndexes.size() > 0) {
            // System.out.println("Removed columns and rows: " + highlyCorrIndexes.toString());
            snpNum = originalSampleSize - highlyCorrIndexes.size();
            List<PValueWeight> tmpPvalueWeightList = new ArrayList<PValueWeight>(pvalueWeightList);
            pvalueWeightList.clear();
            DoubleMatrix2D tmpCorMat = new DenseDoubleMatrix2D(snpNum, snpNum);
            int incRow = 0;
            int incCol = 0;
            for (int i = 0; i < originalSampleSize; i++) {
                if (highlyCorrIndexes.contains(i)) {
                    continue;
                }
                pvalueWeightList.add(tmpPvalueWeightList.get(i));
                incCol = 0;
                for (int j = 0; j < originalSampleSize; j++) {
                    if (highlyCorrIndexes.contains(j)) {
                        continue;
                    }
                    tmpCorMat.setQuick(incRow, incCol, ldRMatrix.getQuick(i, j));
                    incCol++;
                }
                incRow++;
            }
            ldRMatrix = tmpCorMat;
            tmpPvalueWeightList.clear();
            // System.out.println(corrMat.toString());
        }

        if (snpNum == 1) {
            keySNPPosition[0] = pvalueWeightList.get(0).physicalPos;
            return pvalueWeightList.get(0).pValue;
        }

        double df = 0;
        double Y = 0;
        for (int i = 0; i < snpNum; i++) {
            PValueWeight pv = pvalueWeightList.get(i);
            Y += (pv.weight * pv.chiSquare);
        }
        EigenvalueDecomposition ed = new EigenvalueDecomposition(ldRMatrix);
        DoubleMatrix1D eVR = ed.getRealEigenvalues();
        double[] weights = eVR.toArray();
        boolean hasNeg = false;
        for (int i = 0; i < weights.length; i++) {
            if (weights[i] < 0) {
                hasNeg = true;
                //  weights[t] = 0;
            }
        }
        if (hasNeg) {
            int sss = 0;
            // System.out.println(ed.getD().toString());
        }
        rcon.assign("weights", weights);

        double p = rcon.eval("pchisqsum(" + Y + ",df=rep(1," + snpNum + "),method=\"int\",lower=FALSE,a=weights)").asDouble();//This step will cost a lot of time. 
        if (p < 1E-8) {
            int sss = 0;
            System.out.println(eVR.toString());
            System.out.println(ed.getImagEigenvalues());
        }

        return p;
    }

}
