/*
 * predExpLevel.cpp
 * This file is part of isoLasso library; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License as
 * published by the Free Software Foundation; version 2.1 of the License.
 *
 *  Created on: 2011-10-14
 *      Author: daniel
 */
#include "predExpLevel.h"
#include "mathlib.h"
#include "library.h"
#include "NewInstance.h"
#include "emfunction.h"
#include <algorithm>
#include <functional>
//#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_linalg.h>
#include <sstream>
#include <numeric>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <set>
#include <limits>
#include <iostream>
#include <cstdio>

using namespace std;

bool ReCalgetAllPossibleForm;
long exon;
vector<vector<long > > countmat;
vector<long> indeg, outdeg;
vector<vector<long > > conngraph;
vector<vector<long> > refinfo;
long *stack,*stackprev, *isvisited;


STAT predExpLevel(NewInstance &ins, double &lambda, vector<string> &varargin  ) {
  STAT stat;
  stat.negw = 0;
  stat.nw = 0;
  stat.ncandidate = 0;
  stat.tau = 0;
  stat.adjust = 1;

  vector<vector<long> > &pred=ins.Pred;
  vector<double> &predexplv=ins.PredExp;
  vector<long> &predsubinst=ins.PredSubInst;
  vector<char> &preddir=ins.PredDir;
  vector<string> & predid=ins.PredID;

  long minjuncreads=1, maxniso=10;
  double exonminfrac = 0.003, juncexpfrac = 0.002 , intronretentionfrac = 0.1 ,fracstep = 1.2;
  long useref=0;

  for (int i = 0; i < varargin.size(); ++i){
    if (varargin[i] == "--maxniso")
      maxniso = atol(varargin[i+1].c_str());
    if (varargin[i] == "--useref")
      useref = 1;
    if (varargin[i] == "--userefonly")
      useref = 2;
    if (varargin[i] == "--forceref")
      useref = 3;
  }

  if (ins.nExons == 0) return stat;
 
  vector<vector<long> > stap;
  vector<long> subinst;
  if(useref!=3){
    vector<string> var;
    string tmpST;
    if (testRunningTime) {
      checkpoint = clock();
      printf("Start loop to count instance at %ld\t", checkpoint);
      printf("Previous phase cost %ld\n", checkpoint - begintime);
      begintime = checkpoint;
    }
    long getcount = 0;
    ReCalgetAllPossibleForm = true;
    stack = new long[ins.rcount.size() * ins.rcount.size()];
    stackprev = new long[ins.rcount.size() * ins.rcount.size()];
    isvisited = new long[ins.rcount.size() * ins.rcount.size()];
    while (1){
      var.clear();
      stringstream ss (stringstream::in | stringstream::out);
      ss << " --minjunccount " << minjuncreads << " --maxniso " << maxniso  
         << " --countonly --exonexpfrac " << exonminfrac << " --juncexpfrac " 
         << juncexpfrac << " --intronretentionexpfrac " << intronretentionfrac 
        << " --useref " << useref;
      while (ss >> tmpST) var.push_back(tmpST);
  
      vector<vector<long> > currentniso;
      vector<long> tmpa, tmpb; 
      long isreachmax;
      ++getcount;
      /*printf("\nCall GetAllIsoform %ld:\t", getcount);
      if (testRunningTime) {
        checkpoint = clock();
        printf("Previous phase cost %ld\n", checkpoint - begintime);
        begintime = checkpoint;
      }*/
      getAllPossibleIsoform(ins, var, currentniso, tmpa, tmpb, isreachmax, subinst);
      if (isreachmax == false) break;
      if (exonminfrac > 1) {
        printf("Warning: exonminfrac > 1\n");
        break;
      }
      exonminfrac = exonminfrac * fracstep;
      intronretentionfrac = intronretentionfrac * fracstep;
      if (intronretentionfrac < exonminfrac)
        intronretentionfrac = exonminfrac;
      if (intronretentionfrac > 1)
        intronretentionfrac = 0.99;
      juncexpfrac = juncexpfrac * fracstep;
      if (juncexpfrac > 1) juncexpfrac = 0.99;
      //cout<<"exonminfrac:"<<exonminfrac<<",sz:"<<currentniso.size()<<endl;
    }
  
    vector<long> tmpa, tmpb; 
    long tmpc;
  
    stringstream ss2 (stringstream::in | stringstream::out);
    ss2 << " --minjunccount " << minjuncreads << " --maxniso " << maxniso  
       <<" --exonexpfrac " << exonminfrac << " --juncexpfrac " 
       << juncexpfrac << " --intronretentionexpfrac " << intronretentionfrac 
       << " --useref " << useref;
    var.clear();
    while (ss2 >> tmpST) var.push_back(tmpST);
  
    if (testRunningTime) {
      checkpoint = clock();
      printf("Begin to get all possible isofrom from %ld\t", checkpoint);
      printf("Previous phase cost %ld\n", checkpoint - begintime);
      begintime = checkpoint;
    }
    ReCalgetAllPossibleForm = true;
    //enumerate all possible isoforms
    getAllPossibleIsoform(ins, var, stap, tmpa, tmpb, tmpc, subinst);
    delete []stack;
    delete []stackprev;
    delete []isvisited;
  
    if (testRunningTime) {
      checkpoint = clock();
      printf("Finished get all possible isofrom at %ld\t", checkpoint);
      printf("Previous phase cost %ld\n", checkpoint - begintime);
      begintime = checkpoint;
    }
  }else{
    stap=ins.Refs;
    subinst=vector<long>(stap.size(),1);
  }

  vector<vector<long> > currenttap;
  if (stap.size() != 0) zeros(pred, stap.size(), stap[0].size());
  zeros(predexplv, stap.size());
  zeros(predsubinst, stap.size());

  long max_subinst;
  if (subinst.size() != 0) max_subinst = *max_element(subinst.begin(), subinst.end());
  else max_subinst = -1;
  long ntap = 0;
  if (testRunningTime) {
    checkpoint = clock();
    printf("Begin to predict all possible isofrom from %ld\t", checkpoint);
    printf("Previous phase cost %ld\n", checkpoint - begintime);
    begintime = checkpoint;
  }
  for (long i = 0; i <= max_subinst; ++i){
    currenttap.clear();
    for (long j = 0; j < subinst.size(); ++j)
      if (subinst[j] == i) {
        long isolen = 0;
        for (long k = 0; k < stap[j].size(); ++k){
          isolen += stap[j][k] * ins.exonLen[k];
        }
        if (isolen > ins.readLen|| useref==3) currenttap.push_back(stap[j]);
      }
    if (currenttap.size() == 0) continue;
    vector<vector<long> >npred;
    vector<double> npredexplv;
    STAT nstat;

    //predict subisoform here
    predsubisoform(currenttap, ins, lambda, varargin, npred, npredexplv, nstat);
    long currentp = npred.size();
    for (long j = 0; j < currentp; ++j){
      pred[ntap + j] = npred[j];
      predexplv[ntap + j] = npredexplv[j];
      predsubinst[ntap + j] = i;
    }
    ntap += currentp;
    stat.negw = stat.negw + nstat.negw;
    stat.nw = stat.nw += nstat.nw;
    if (nstat.tau != 0) stat.tau = nstat.tau;
    if (nstat.adjust != 0) stat.adjust = nstat.adjust;
  }
  if (testRunningTime) {
    checkpoint = clock();
    printf("End to predicate all possible isofrom at %ld\t", checkpoint);
    printf("Previous phase cost %ld\n", checkpoint - begintime);
    begintime = checkpoint;
  }

  pred.erase(pred.begin() + ntap, pred.end());
  predexplv.erase(predexplv.begin() + ntap, predexplv.end());
  predsubinst.erase(predsubinst.begin() + ntap, predsubinst.end());

  //preddir and predid
  if(useref!=3){
    //use '+' right now
    char cdir='.';
    if(ins.dir==1)cdir='+';
    if(ins.dir==-1)cdir='-';
    preddir=vector<char>(pred.size(),cdir);
    //ID
    for(int i=0;i<predsubinst.size();i++){
      stringstream ss;
      ss<<"Pred_"<<ins.id<<"_"<<predsubinst[i]<<"_"<<i;
      predid.push_back(ss.str());
    }
  }
  else{
    preddir=ins.RefDirs;
    predid=ins.RefIDs;
  }
  return stat;
}

//int explvasmatrev(vector<vector<long> > &usecal, vector<double> &predexplv, vector<double>& yvalue);


void predsubisoform(vector<vector<long> >&tap, NewInstance &ins, double &lambda, vector<string> inputarg, vector<vector<long> > &pred, vector<double> &predexplv, STAT &stat){
  stat.negw = 0;
  stat.nw = 0;
  stat.ncandidate = 0;
  stat.tau = 0;
  stat.adjust = 0;
  bool verbose = false, outputallcandidate = false, useem = false;
  double steplen = 1.3;
  double smallestlambda = 1e-5;

  bool nofilter=false;

  for (int i = 0; i < inputarg.size(); ++i){
    if (inputarg[i] == "--verbose")
      verbose = true;
    if (inputarg[i] == "--steplen")
      steplen = atof(inputarg[i + 1].c_str());
    if (inputarg[i] == "--smallestlambda")
      smallestlambda = atof(inputarg[i + 1].c_str());
    if (inputarg[i] == "--outputall")
      outputallcandidate = true;
    if (inputarg[i] == "--useem")
      useem = true;
    if (inputarg[i] == "--no-filter")
      nofilter=true;
  }

  double lambdafornoambiguous = 0.0001, smallpositive = 1e-3, maxwfraction = 0.005;
  double minerrorsolution = 1.01;
  bool junccoverage = false, exoncoverage = false, addjunccvg = false, addexoncvg = false;
  bool usescalequadprog = false; // for quadprog, use scale version (true) or constraint norm version (false)
  long maxtapsize = 1;

  vector<double> exonexp(ins.rcount.size(), 0);
  for (long i = 0; i < ins.rcount.size(); ++i)
    exonexp[i] = ((double)ins.rcountColSum[i] + ins.rcountRowSum[i]) / 2;
  vector<long> exoncover;
  vector<vector<long> > oldtap;
  matrixSum(tap, 1, exoncover);
  oldtap.assign(tap.begin(), tap.end());
  tap.clear();
  tap.assign(oldtap.size(), vector<long>());
  vector<long> exl_exoncover;
  for (long j = 0; j < exoncover.size(); ++j)
    if (exoncover[j] > 0){
      for (long i = 0; i < oldtap.size(); ++i) tap[i].push_back(oldtap[i][j]);
      exl_exoncover.push_back(ins.exonLen[j]);
    }
  stat.ncandidate = tap.size();
  vector<double> weight, yvalue;
  getExonWeight(exl_exoncover, ins.readLen, weight);
  for (long j = 0; j < exoncover.size(); ++j)
    if (exoncover[j] > 0)
      yvalue.push_back(exonexp[j] / ins.exonLen[j]);
  for (long j = 0; j < weight.size(); ++j)
    yvalue[j] = yvalue[j] * weight[j];
  vector<vector<double> > diag_weight, xvalue;
  zeros(xvalue, weight.size(), tap.size());
  zeros(diag_weight, weight.size(), weight.size());
  for (long j = 0; j < weight.size(); ++j) diag_weight[j][j] = weight[j];
  for (long i = 0; i < diag_weight.size(); ++i)
    for (long j = 0; j < tap.size(); ++j){
      for (long r = 0; r < diag_weight[i].size(); ++r)
        xvalue[i][j] += diag_weight[i][r] * tap[j][r];
    }

  if (outputallcandidate) {
    //pred.assign(oldtap.begin(), oldtap.end());
    //use EM instead of mldivide
    emselect(ins, oldtap,  inputarg, pred, predexplv, stat);

    //mldivide(tap, yvalue, predexplv, true);
    //for (long k = 0; k < predexplv.size(); ++k)
    //  if (predexplv[k] < 0) ++stat.negw;
    //  else if (predexplv[k] > 0) ++stat.nw;
    //return;
  }

  //Use EM
  if (useem == true){
    emselect(ins, oldtap,  inputarg, pred, predexplv, stat);
    return;
  }

  //only one predictions
  if (tap.size() == 1){
    pred.assign(oldtap.begin(), oldtap.end());
    double exonexp_sum = 0, exl_sum = 0;
    for (long k = 0; k < exoncover.size(); ++k)
      if (exoncover[k] != 0){
        exonexp_sum += exonexp[k];
        exl_sum += ins.exonLen[k];
      }
    double rpkm=exonexp_sum / exl_sum;
    if(!nofilter && rpkm<smallpositive){ pred.clear();predexplv.clear(); }
    predexplv.push_back(rpkm);
    stat.negw = 0;
    stat.nw = 1;
    return;
  }

  long nlambda = ceil((double)(log((double)(lambda / smallestlambda)) / log ((double)steplen)));
  double newlambda;
  vector<double> lambdapool;
  zeros(lambdapool, nlambda);
  if (tap.size() > maxtapsize || usescalequadprog == false)
    newlambda = lambda;
  else newlambda = smallestlambda;
  for (long i = 0; i < nlambda; ++i){
    lambdapool[i] = newlambda;
    if (tap.size() > maxtapsize || usescalequadprog == false)
      newlambda = newlambda / steplen;
    else newlambda = newlambda * steplen;
  }

  vector<double> errpool,  ecvg, jcvg;
  vector<long> sizepool;
  zeros(errpool, nlambda);
  zeros(sizepool, nlambda);
  vector<vector<double> > wpool;
  zeros(ecvg, nlambda);
  zeros(jcvg, nlambda);
  double minerr = numeric_limits<double>::max();
  vector<string> var;

  stringstream ss (stringstream::in | stringstream::out);
  string tmpST;
  ss << "usejunccon " << junccoverage << " useexoncon " << exoncoverage 
    << " juncp 1e-2 usescale " << usescalequadprog;
  while (ss >> tmpST) var.push_back(tmpST);
  long exitflag;
  double fval;
  vector<double> w;
  long nl;
  if (testRunningTime) {
    checkpoint = clock();
    printf("Begin to start lasso with constraint from %ld\t", checkpoint);
    printf("Previous phase cost %ld\n", checkpoint - begintime);
    begintime = checkpoint;
  }

  int encountererror=0;
  for (nl = 0; nl < nlambda; ++nl){
    newlambda = lambdapool[nl];
    LassoWithConstraint(xvalue, yvalue, newlambda, ins, var, w, fval, exitflag);
  
    double err = 0, tmp;
    if (exitflag < 0) { //errors in quadratic program
      err = numeric_limits<double>::infinity();
      encountererror++;
      errpool[nl] = err;
      sizepool[nl] = 0;
      wpool.push_back(w);
      continue;
    }

    for (long i = 0; i < xvalue.size(); ++i){
      tmp = 0;
      for (long j = 0; j < xvalue[i].size(); ++j)
        tmp += xvalue[i][j] * w[j];
      tmp -= yvalue[i];
      tmp = tmp * tmp;
      err += tmp;
    }

    vector<long> wsel;
    if(!nofilter)
      selectw(tap, w, smallestlambda, maxwfraction, addjunccvg, addexoncvg, ins, wsel);
    else
      wsel=vector<long>(w.size(),1);
    double ijc, iec;
    vector<vector<long> > tap_wsel;
    for (long i = 0; i< wsel.size(); ++i)
      if (wsel[i] != 0) tap_wsel.push_back(tap[i]);
    checkCoverage(tap_wsel, ins, ijc, iec);
    ecvg[nl] = iec;
    jcvg[nl] = ijc;
    long t = 0;
    vector<long> ws;
    for (long i = 0; i < wsel.size(); ++i){
      t += wsel[i];
      if (wsel[i] > 0) ws.push_back(i);
    }

    if (verbose){
      printf("Lambda: %lf sqr err: %lf\n", newlambda, err);
      double wsum = 0,  wpsum = 0;
      long v = 0;
      for (long i = 0; i < w.size(); ++i){
        wsum += w[i];
        if (w[i] > smallpositive){
          ++v;
          wpsum += w[i];
        }
      }
      printf("prediction: %ld positive: %ld\n", t, v);
      printf("sum(w): %lf sum(w>0): %lf\n", wsum, wpsum);
      printf("junc cvg %lf, exonc cvg: %lf\n:", ijc, iec);
      printf("selection: ");
      for (long i = 0; i < ws.size(); ++i) printf("%ld ", ws[i]);
      printf("\n");
    }

    errpool[nl] = err;
    sizepool[nl] = t;
    wpool.push_back(w);

    if ((err != numeric_limits<double>::infinity()) && (err != numeric_limits<double>::signaling_NaN())){
      if (err < minerr) minerr = err;
      if (err > minerr * minerrorsolution && !verbose) break;
    }

  }
  
  if (encountererror==nlambda  ){
      printf("error happened when processing instance %s, switch to EM.\n", ins.id.c_str());
      emselect(ins, oldtap,  inputarg, pred, predexplv, stat);
      //stat.nw = -1;
      //stat.negw = 0;
      //pred.clear();
      //predexplv.clear();
      return;
  }


  if (testRunningTime) {
    checkpoint = clock();
    printf("End lasso with constraint from %ld\t", checkpoint);
    printf("Previous phase cost %ld\n", checkpoint - begintime);
    begintime = checkpoint;
  }
  //not support plotlambda, always assuming its false
  if(nl<errpool.size()-1){
    errpool.erase(errpool.begin() + nl + 1, errpool.end());
    sizepool.erase(sizepool.begin() + nl + 1, sizepool.end());
    wpool.erase(wpool.begin() + nl + 1, wpool.end());
  }
  if (verbose) printf("Possible isoforms %d", tap.size());

  if (minerr == numeric_limits<double>::infinity()){
    stat.nw = 0;
    stat.negw = 0;
    pred.clear();
    predexplv.clear();
    return;
  }
  vector<long> efferr, ssel;
  double msz = numeric_limits<double>::max();
  double merr = numeric_limits<double>::max();
  for (long i = 0; i < errpool.size(); ++i)
    if (errpool[i] != numeric_limits<double>::signaling_NaN() && errpool[i] <= minerrorsolution*minerr){
      efferr.push_back(1);
      if (sizepool[i] < msz) msz = sizepool[i];
    } else efferr.push_back(0);
  for (long i = 0; i <efferr.size(); ++i)
    if (efferr[i] != 0 && (fabs((double)sizepool[i] - (double)msz) < smalldelta)){
      ssel.push_back(1);
      if (errpool[i] < merr)
        merr = errpool[i];
    }
    else ssel.push_back(0);
  for (long i = 0; i < ssel.size(); ++i){
    if (ssel[i] != 0 && (fabs((double)errpool[i] - (double)merr) < smalldelta)){
      w.assign(wpool[i].begin(), wpool[i].end());
      break;
    }
  }
  pred.clear();
  predexplv.clear();
  for(int i=0;i<w.size();++i){
    if (fabs((double)w[i]) > smallpositive ||  nofilter){
      pred.push_back(oldtap[i]);
      predexplv.push_back(w[i]);
    }
  }
 
  //use (X^TX+lambda*I)^-1X^Ty to calculate expression level
  //Not used anymore, since it may create negative exp lvs
  if(false){
    //pred.clear(); vector<vector<long> >usecal;
    //for (long i = 0; i < w.size(); ++i)  if (fabs((double)w[i]) > smallpositive){
    //    usecal.push_back(tap[i]);   pred.push_back(oldtap[i]);
    //}
    //explvasmatrev(usecal, predexplv,  yvalue);
  }
}//end of function

/*
use (X^TX+lambda*I)^-1X^Ty to calculate expression level
*/
/*
int explvasmatrev(vector<vector<long> > &usecal, vector<double> &predexplv, vector<double>& yvalue){

    double lambdafornoambiguous = 0.0001;
    gsl_matrix *result = gsl_matrix_alloc(usecal.size(), usecal.size());
    gsl_matrix *inverse = gsl_matrix_alloc(usecal.size(), usecal.size());
    for (long i = 0; i < usecal.size(); ++i)
      for (long j = 0; j < usecal.size(); ++j){
        double tmp = 0;
        for (long k = 0; k < usecal[i].size(); ++k) tmp += usecal[i][k] * usecal[j][k];
        if (i == j) gsl_matrix_set(result, i, j, tmp + lambdafornoambiguous);
        else gsl_matrix_set(result, i, j, tmp);
      }
  
  
    gsl_permutation *p = gsl_permutation_alloc(usecal.size());
    int s;
    gsl_linalg_LU_decomp (result, p, &s);
    gsl_linalg_LU_invert(result, p, inverse);
    gsl_matrix_free(result);
    gsl_matrix *v = gsl_matrix_alloc(usecal.size(), usecal[0].size());  
    for (long i = 0; i < usecal.size(); ++i)
      for (long j = 0; j < usecal[0].size(); ++j){
        double tmp = 0;
        for (long k = 0; k < usecal.size(); ++k)
          tmp += gsl_matrix_get(inverse, i, k) * usecal[k][j];
        gsl_matrix_set(v, i, j, tmp);
      }
    gsl_matrix_free(inverse);
    gsl_permutation_free(p);
    predexplv.resize(usecal.size(), 0);
    for (long i = 0; i < usecal.size(); ++i){
      for (long k = 0; k < yvalue.size(); ++k) 
        predexplv[i] = predexplv[i] + gsl_matrix_get(v, i, k) * yvalue[k];
      //if (predexplv[i] < 0) ++stat.negw;
    }
  return 0;
}
*/

void selectw(vector<vector<long> >&tap, vector<double> &w, double smallpositive, double maxwfraction, bool requirejunccvg, bool requireexoncvg, NewInstance &ins, vector<long> &wsel){
  double max_w = *max_element(w.begin(), w.end());
  for (long i = 0; i < w.size(); ++i)
    if (w[i] > smallpositive && w[i] > max_w * maxwfraction)
      wsel.push_back(1);
    else wsel.push_back(0);
  vector<MATLAB_SORT_ENTRY<double> > newdata;
  matlabSortWithIndex(w, newdata);
  for (long i = newdata.size() - 1; i >= 0; --i){
    if (wsel[newdata[i].idx] == 1) continue;
    vector<vector<long> > tap_wsel;
    double jc, ec;
    for (long j = 0; j < wsel.size(); ++j)
      if (wsel[j] != 0) tap_wsel.push_back(tap[j]);
    checkCoverage(tap_wsel, ins, jc, ec);
    if ((!requirejunccvg || fabs(((double)jc - 1)) < smalldelta ) && (fabs(((double)ec - 1)) < smalldelta || !requireexoncvg)) return;
    wsel[i] = 1;
  }
}

void checkCoverage(vector<vector<long> > &pred, NewInstance &ins, double &isjunccov, double &isexoncov){
  long m = pred.size(), n = pred[0].size();
  vector<long> exoncov;
  vector<vector<long> > covcount;
  matrixSum(pred, 1, exoncov);
  isexoncov = 0;
  for (long i = 0; i < exoncov.size(); ++i) 
    if (exoncov[i] != 0) isexoncov += 1;
  isexoncov /= n;
  zeros(covcount, n, n);
  for (long i = 0; i < m; ++i){
    vector<long> isonz;
    for (long j = 0; j < pred[i].size(); ++j)
      if (pred[i][j] != 0) isonz.push_back(j);
    for (long j = 1; j < isonz.size(); ++j)
      ++covcount[isonz[j - 1]][isonz[j]];
  }
  long nallj = 0, ncov = 0;
  for (long i = 0; i < ins.rcount.size(); ++i)
    for (long j = i + 1; j < ins.rcount[i].size(); ++j)
      if (ins.rcount[i][j] != 0)
        ++nallj;
  for (long i = 0; i < covcount.size(); ++i)
    for (long j = 0; j < covcount[i].size(); ++j)
      if (covcount[i][j] != 0) ++ncov;
  isjunccov = (double) ncov / nallj;
}



void getAllPossibleIsoform(NewInstance &ins, 
  vector<string> varargin, 
  vector<vector<long> >&isf, //where results stored
  vector<long> &isambigious, 
  vector<long> &fullysupport, 
  long &isreachmax, 
  vector<long> &subinst //the sub-instance
){
  //1st optional parameters
  long minjunccount = 1; // the minimum number of reads to be considered junctions
  // 2nd optional parameters
  bool enumerateadjacentjunction = false; // try to connect adjacent junctions even if they're not covered by junction reads
  long maxniso = 5000;
  double exonexpfrac = 0.0001; // used to decide whether exons should be covered, or whether exons be encountered
  double juncexpfrac = 0.0001; // used to decide whether one exon with low exp lv but high junction exp lv should be counted as start
  double intronretentionexpfrac = 0.03; // for possible intron retention exons, should expect higher value of exonexpfrac
  bool countonly = false; // only return counts?
  subinst.push_back(1);
  isambigious.push_back(0);
  fullysupport.push_back(1);
  long subid = 0;
  long useref = 0; // enumerate the exon junction with the help of ref seq?

  // if all the expressed segs of one isoform  have exp lv greater than this
  //  value, then add it to the list
  long minrefexonexp=1000000; 

  for (int i = 0; i < varargin.size(); ++i){
    if (varargin[i] == "--minjunccount")
      minjunccount = atol(varargin[i + 1].c_str());
    //TODO check what will happen for boolean
    if (varargin[i] == "--enumerateadjacentjunction")
      enumerateadjacentjunction = getBoolenValue(varargin[i + 1]);
    if (varargin[i] == "--maxniso")
      maxniso = atol(varargin[i + 1].c_str());
    if (varargin[i] == "--exonexpfrac")
      exonexpfrac = atof(varargin[i + 1].c_str());
    if (varargin[i] == "--juncexpfrac")
      juncexpfrac = atof(varargin[i + 1].c_str());
    if (varargin[i] == "--countonly")
      countonly = true;
    if (varargin[i] == "--intronretentionexpfrac")
      intronretentionexpfrac = atof(varargin[i + 1].c_str());
    if (varargin[i] == "--useref")
      useref = atol(varargin[i + 1].c_str());
  }

  if (useref > 2){
    subinst.clear();
    for (int i = 0; i < ins.Refs.size(); ++i)
      subinst.push_back(0);
    if (countonly != 0 ){
      isf.assign(ins.Refs.begin(), ins.Refs.end());
      while (1) {
        long sum = 0;
        for(int i = 0; i < subinst.size(); ++i)
          if (subinst[i] != 0) ++sum;
        if (sum <= 0) break;
        subid += 1;
        vector<long> sel;
        for (int i = 0; i < subinst.size(); ++i)
          if (subinst[i] == 0)
            sel.push_back(i);
        vector<long> subbin;
        subbin.assign(isf[sel[0]].begin(), isf[sel[0]].end());
        subinst[sel[0]] = subid;
        while (1) {
          bool isupdated = false;
          for (int i = 1; i < sel.size(); ++i)
            if (subinst[sel[i]] == 0){
              long sumtmp;      
              for (int j = 0; j <subbin.size(); ++j)
                sumtmp += isf[sel[i]][j] * subbin[j];
              if (sumtmp >0) {
                isupdated = true;
                for (int j = 0; j < isf[sel[i]].size(); ++j)
                  if (isf[sel[i]][j] > 0)
                    subbin[j] = 1;
                subinst[sel[i]] = subid;
              }
          }
            if (isupdated == false) break;
        }
      }
    }
    else {
      vector<long> tmp;
      tmp.push_back(ins.Refs.size());
      isf.push_back(tmp);
    }
    zeros(isambigious, subinst.size());
    fullysupport.assign(subinst.begin(), subinst.end());
    isreachmax = 0;
  }
  // Connect adjacent junctions
  if (ReCalgetAllPossibleForm == true){
    ReCalgetAllPossibleForm = false;
    exon = ins.nExons;
    countmat.assign(ins.rcount.begin(), ins.rcount.end());
    if (enumerateadjacentjunction) {
      for (long i = 0; i < countmat.size() - 1; ++i){
        long j = i + 1;
        if (countmat[i][i] > 0 && countmat[i][j] < minjunccount && (ins.exonBoundary[j][0] - ins.exonBoundary[i][1] == 1))
          countmat[i][j] = minjunccount;
      }
    }
    //calculate the in degree and out degree
    zeros(indeg, exon);
    zeros(outdeg, exon);
    for (long i = 0; i < countmat.size(); ++i)
      for (long j = i + 1; j < countmat.size(); ++j)
        if (countmat[i][j] >= minjunccount) {
          outdeg[i] = outdeg[i] + 1;
          indeg[j] = indeg[j] + 1;
        }
    
    // get the connectivity graph
    zeros(conngraph, exon, exon);
    for (long i = 0; i < countmat.size(); ++i)
      for (long j = i + 1; j < countmat.size(); ++j)
        if (countmat[i][j] >= minjunccount) {
          conngraph[i][j] = 1;
          for (long k = 0; k <= i; ++k)
            if (conngraph[k][i] == 1)
              conngraph[k][j] = 1;
        }
  
    // build a matrix of existing annotation junction map
    zeros(refinfo, exon, exon);
    if (useref != 0) {
      for (long i = 0; i < ins.Refs.size(); ++i){
        vector<long> allj;
        for (long j = 0; j < ins.Refs[i].size(); ++j)
          if (ins.Refs[i][j] != 0) allj.push_back(j);
        if (allj.size() == 1)
          refinfo[allj[0]][allj[0]] = 1;
        else
          for (long j = 0; j < allj.size(); ++j)
            refinfo[allj[j]][allj[j+1]] = refinfo[allj[j]][allj[j+1]] + 1;
      }
    }
  }
  long nsegs, segid, lastsegid; 
  vector<long> covermap;
  zeros(covermap, exon);
  isreachmax = 0;
  vector<vector<long> > alliso;
  zeros(alliso, maxniso, exon);
  zeros(subinst, maxniso);

  long niso = 0, nsubinst = -1;
  vector<vector<long> >  ctypes(exon, vector<long>());
  long currentctypes = -1;

  //exon expression level and junction expression level
  vector<double> explv, exonweight,exonexplv,juncexplv, rightexplv, leftexplv;
  getExonExpLv(ins, explv);
  getExonWeight(ins.exonLen, ins.readLen, exonweight);
  double max_exonexplv = 0, max_juncexplv = 0;
  for (long i = 0; i < explv.size(); ++i){
    exonexplv.push_back(explv[i] * exonweight[i]);
    if (exonexplv[i] > max_exonexplv) max_exonexplv = exonexplv[i];
  }
  getExonJunctionExpLv(ins, juncexplv, rightexplv, leftexplv);
  if (juncexplv.size() > 0) max_juncexplv = juncexplv[0];
  else max_juncexplv = 0;
  for (long i= 0; i < juncexplv.size(); ++i)
    if (juncexplv[i] > max_juncexplv) max_juncexplv = juncexplv[i];

  vector<long> isintronretentionleft, isintronretentionright;
  zeros(isintronretentionleft, exon);
  zeros(isintronretentionright, exon);

  for (long i = 1; i < exon; ++i)
    if (ins.exonBoundary[i][0] == ins.exonBoundary[i-1][1] + 1)
      isintronretentionleft[i] = 1;
  for (long i = 0; i < exon - 1; ++i)
    if (ins.exonBoundary[i][1] == ins.exonBoundary[i+1][0] - 1)
      isintronretentionright[i] = 1;
  vector<long> isintronretention;
  for (long i = 0; i < exon; ++i)
    if (isintronretentionleft[i] + isintronretentionright[i] > 0)
      isintronretention.push_back(1);
    else isintronretention.push_back(0);

  vector<long> selectedexons; 
  for (long i = 0; i < isintronretention.size(); ++i) {
    long value = ((isintronretention[i] == 0  && exonexplv[i] >= max_exonexplv * exonexpfrac) || (isintronretention[i] == 1 && exonexplv[i] >= max_exonexplv * intronretentionexpfrac) || (juncexplv[i] >= max_juncexplv * juncexpfrac));
    selectedexons.push_back(value);
  }

  selectedexons[0] = 1;
  for (long i = 0; i < ins.exonStat.size(); ++i)
    if (ins.exonStat[i][3] > 0.7)
      selectedexons[i] = 0;
  
  if (useref != 0 ){
    for (long i = 0; i < ins.Refs.size(); ++i){
      getlastoffirstexon(ins.Refs[i], ins.exonBoundary, nsegs, segid, lastsegid);
      selectedexons[segid] = 1;
    }
  }
  long stackcount = 0; 
  long sumcounts = 0;
  while (1){
    long sum = 0;
    for (long i = 0; i < selectedexons.size(); ++i)
      if (selectedexons[i] != 0 && covermap[i] != 1)
        ++sum;
    if (sum <= 0) break;
    vector<vector<long> > siso, imap;
    long reachmax, ind;
    for (long i = 0; i < covermap.size(); ++i){
      long t = ((covermap[i] != 1) && selectedexons[i]);
      if (t != 0) {
        ind = i;
        break;
      }
    }
    ++stackcount;
    //printf("\nstack count %ld", stackcount);
    stackiso(countmat, ind, minjunccount, maxniso, exonexplv, exonexpfrac, countonly, isintronretentionleft, intronretentionexpfrac, refinfo, siso, reachmax, imap);
    if (countonly){
      sumcounts += siso[0][0];
      if (sumcounts >= maxniso || reachmax){
        isreachmax = true;
        isf.clear();
        isf.push_back(vector<long>(1, sumcounts));
        return;
      }
    }
    else {
      long sum_covermap = 0;
      vector<long>cvgmap;
      for (long i = 0; i < imap[0].size(); ++i){
        if (imap[0][i] > 0){
            cvgmap.push_back(1);
            sum_covermap += covermap[i];
          }
        else cvgmap.push_back(0);
      }
      long currentnsub;
      if (sum_covermap == 0){
        ++nsubinst;
        currentnsub = nsubinst;
        currentctypes = currentctypes + 1;
        ctypes[currentctypes].assign(cvgmap.begin(), cvgmap.end());
      }
      else {
        for (long i = 0; i <= currentctypes; ++i){
          long testsum = 0;
          for (long j = 0; j < cvgmap.size(); ++j){
            if (cvgmap[j] != 0) testsum += ctypes[i][j];
          }
          if (testsum > 0) currentnsub = i;
        }
      }
      long remainisoc = min(long(alliso.size()) - niso, long(siso.size()));
      if (remainisoc > 0){
        for (long i = 0; i < remainisoc; ++i){
          alliso[niso + i].assign(siso[i].begin(), siso[i].end());
          subinst[niso + i] = currentnsub;
        }
        niso += remainisoc;
      }
      if (reachmax == 1 || niso>=alliso.size()){
        isreachmax = true;
        break;
      }
    }
    for (long i = 0; i < covermap.size(); ++i)
      covermap[i] = covermap[i] ||  imap[0][i];
  }

  if (countonly){
    isf.clear();
    isf.push_back(vector<long>(1, sumcounts));
    if (useref > 1){
      isf[0][0] = ins.Refs.size();
    }
    return;
  }

  alliso.erase(alliso.begin() + niso, alliso.end());
  subinst.erase(subinst.begin() + niso, subinst.end());
  vector<long> tempSubinst;
  isf.clear();
  for (long i = 0; i < alliso.size(); ++i){
    long rowSum = 0;
    for (long j  = 0; j < alliso[i].size(); ++j) rowSum += alliso[i][j];
    if (rowSum > 0) {
      isf.push_back(alliso[i]);
      tempSubinst.push_back(subinst[i]);
    }
  }
  subinst.assign(tempSubinst.begin(), tempSubinst.end());
  
  vector<long> orvalid, isplsupported;
  if (testRunningTime) {
    checkpoint = clock();
    printf("Begin to filter impossible isofrom from %ld\t", checkpoint);
    printf("Previous phase cost %ld\n", checkpoint - begintime);
    begintime = checkpoint;
  }
  if (isf.size() != 0) filterImpossilbeCandidates(isf, ins, isf, orvalid, isplsupported);
  long index = 0;
  fullysupport[0] = alliso.size();
  for (long i = 0; i < subinst.size(); ++i)
    if (orvalid[i] != 0) subinst[index++] = subinst[i];
  subinst.erase(subinst.begin() + index, subinst.end());
  zeros(isambigious, isf.size());

  for (long i = 0; i < isf.size(); ++i){
    long isindg1 = 0;
    for (long j = 0; j < exon; ++j)
      if (isf[i][j] == 1){
        if (indeg[j] > 1)
          isindg1 = 1;
        if (isindg1 == 1 && outdeg[j] > 1){
          isambigious[i] = 1;
          break;
        }
      }
  }

  if (useref){
    vector<long> invalid;
    zeros(invalid, ins.Refs.size());
    for (long i = 0; i < ins.Refs.size(); ++i){
      vector<long> thisref;
      thisref.assign(ins.Refs[i].begin(), ins.Refs[i].end());
      getlastoffirstexon(thisref, ins.exonBoundary, nsegs, segid, lastsegid);
      if (nsegs > 2 && segid < lastsegid){
        for (long j = 0; j < segid - 1; ++j) thisref[j] = 0;
        for (long j = lastsegid; j < thisref.size(); ++j) thisref[j] = 0;
      }
      for (long j = 0; j < isf.size(); ++j){
        if (checkcompatibleofisoforms(thisref, isf[j], ins.exonBoundary)){
          invalid[i] = 1;
          break;
        }
      }
      long ts = 0, ts2 = 0;
      for (long j = 0; j < thisref.size(); ++j){
        if (thisref[j] > 0 && exonexplv[j] > minrefexonexp) ++ts;
        ts2 += thisref[j];
      }
      if (ts < ts2)
        invalid[i] = 1;
    }
    bool f = false;
    for (long i = 0; i < invalid.size(); ++i)
      if (invalid[i] == 0){
        f = true;
        break;
      }
    
    if (f){
      vector<vector<long> > toadd;
      for (long i = 0; i < invalid.size(); ++i)
        if (invalid[i] == 0)
          toadd.push_back(ins.Refs[i]);
      vector<long> subinsttoadd;
      zeros(subinsttoadd, toadd.size());
      for (long i = 0; i < toadd.size(); ++i) {
        long maxelement = *max_element(subinst.begin(), subinst.end());
        for (long s = 1; s < maxelement; ++s){
          vector<vector<long> > tcan;
          for (long j = 0; j < subinst.size(); ++j)
            if (subinst[j] == s)
              tcan.push_back(isf[j]);
          if (tcan.size() == 0) continue;
          vector<long> sign,tsign;
          matrixSum(tcan, 1, tsign);
          for (long j = 0; j < tsign.size(); ++j)
            if (tsign[j] > 0) sign.push_back(1);
            else sign.push_back(0);
          long ttsign = 0;
          for (long j = 0; j < toadd[i].size(); ++i)
            ttsign += toadd[i][j] * sign[j];
          if (ttsign > 0){
            subinsttoadd[i] = s;
            break;
          }
        }

        if (subinsttoadd[i] == 0){
          for (long j = 0; j <= i - 1; ++j){
            long toaddtt = 0;
            for (long k = 0; k < toadd[j].size(); ++k)
              toaddtt += toadd[j][k] * toadd[i][k];
            if (toaddtt > 0){
              subinsttoadd[i] = subinsttoadd[j];
              break;
            }
          }
        }

        if (subinsttoadd[i] == 0){
          long ms = *max_element(subinst.begin(), subinst.end());
          subinsttoadd[i] = max(long (1), (*max_element(subinsttoadd.begin(), subinsttoadd.end())) + 1);

        }
      }
      if (useref > 1){
        isf.assign(toadd.begin(),toadd.end());
        subinst.assign(subinsttoadd.begin(), subinsttoadd.end());
      }
      else {
        isf.insert(isf.end(), toadd.begin(), toadd.end());
        subinst.insert(subinst.end(), subinsttoadd.begin(), subinsttoadd.end());
      }
    }

  }

}

bool checkcompatibleofisoforms(vector<long> &ref, vector<long> &candidate, vector<vector<long> > &segboundary){
  vector<vector<long> > refbound, canbound;
  getbound(ref, segboundary, refbound);
  getbound(candidate, segboundary, canbound);

  if (refbound.size() != canbound.size()) return false;
  for (long i = 0; i < refbound.size() - 1; ++i)
    if (refbound[i][1] != canbound[i][1]) return false;
  for (long i = 1; i < refbound.size(); ++i)
    if (refbound[i][0] != canbound[i][0]) return false;
  return true;
}

void getbound(vector<long> &iso, vector<vector<long> > &segboundary, vector<vector<long> > &bound){
  vector<vector<long> > canbound;
  for (long i = 0; i < segboundary.size(); ++i)
    if (iso[i] == 1)
      canbound.push_back(segboundary[i]);
  vector<long> valid;
  ones(valid, canbound.size());
  for (long i = 0; i < valid.size() - 1; ++i)
    if (canbound[i][1] + 1 == canbound[i + 1][0]){
      canbound[i + 1][0] = canbound[i][0];
      valid[i] = 0;
    }
    else bound.push_back(canbound[i]);
  if (valid[valid.size() - 1] == 1)
    bound.push_back(canbound[valid.size() - 1]);
}

void filterImpossilbeCandidates(vector<vector<long> >&tap, NewInstance &oneins, vector<vector<long> > &filtered, vector<long> &orvalid, vector<long> &validf){
  bool checkpairend = false;
  double supportalpha = 0.01;
  long ncan = tap.size();
  long nExons = tap[0].size();
  if (nExons != oneins.nExons)
    printf("Error: incompatible number of segments!\n");

  if (nExons == 1){
    ones(orvalid, ncan);
    ones(validf, ncan);
    return;
  }

  vector<long> valid;
  zeros(valid, ncan);

  vector<vector<long> > juncread;
  vector<long> junccnt;
  for (long i = 0; i < oneins.SGTypes.size(); ++i){
    long tmpSum = 0;
    for (long j = 0; j < oneins.SGTypes[i].size(); ++j)
      tmpSum += oneins.SGTypes[i][j];
    if (tmpSum > 1) {
      juncread.push_back(oneins.SGTypes[i]);
      junccnt.push_back(oneins.SGCounts[i]);
    }
  }
  vector<long> candidate;
  for (long i = 0; i < ncan; ++i){
    candidate.assign(tap[i].begin(), tap[i].end());
    vector<long> coveredseg;
    zeros(coveredseg, nExons);
    if (accumulate(tap[i].begin(), tap[i].end(), 0) == 1){
      continue;
    }
    
    for (long j = 0; j < juncread.size(); ++j){
      if (checkcompatible(candidate, juncread[j]) == false)
        continue;
      if (junccnt[j] <= 0)
        continue;
      for (long k = 0; k < coveredseg.size(); ++k)
        coveredseg[k] += juncread[j][k];
    }

    long tmpSum = 0;
    for (long j = 0; j < coveredseg.size(); ++j)
      if (candidate[j] - coveredseg[j] == 1) ++tmpSum;
    if (tmpSum > 0){
      for (long k = 0; k < coveredseg.size(); ++k)
        if (coveredseg[k] > 0) candidate[k] = 1;
        else candidate[k] = 0;
      long checkSum = 0;
      for (long l1 = 0; l1 < tap.size(); ++l1){
        long rowSum = 0;
        for (long l2 = 0; l2 < tap[l1].size(); ++l2)
          rowSum = tap[l1][l2] - candidate[l2];
        if (rowSum == nExons) {checkSum = 1; break;}
      }
      if (checkSum > 0)
        valid[i] = 0;
      else {
        valid[i] = 1;
        tap[i].assign(candidate.begin(), candidate.end());
      }
    }
    else valid[i] = 1;
  }

  vector<long> segcount, tmpVectorSum, seloneexon, coveredsingleexon;
  matrixSum(tap, 2, segcount);
  matrixSum(oneins.SGTypes, 2, tmpVectorSum);
  for (long i = 0; i < oneins.SGCounts.size(); ++i)
    if (tmpVectorSum[i] == 1 && oneins.SGCounts[i] > 0)
      seloneexon.push_back(1);
    else seloneexon.push_back(0);
  vector<vector<long> > tmpMatrix;
  for (long i = 0; i < seloneexon.size(); ++i)
    if (seloneexon[i] != 0)
      tmpMatrix.push_back(oneins.SGTypes[i]);
  if (tmpMatrix.size() > 0){
    matrixSum(tmpMatrix,1, tmpVectorSum);
    for (long i = 0; i < tmpVectorSum.size(); ++i)
      if (tmpVectorSum[i] > 0)
        coveredsingleexon.push_back(1);
      else coveredsingleexon.push_back(0);
  }
  else coveredsingleexon.assign(oneins.SGTypes[0].size(), 0);

  vector<long> validone;
  zeros(validone, tap.size());
  for (long i = 0 ; i < validone.size(); ++i){
    long tmpSum = 0;
    for (long j = 0; j < tap[i].size(); ++j)
      if (tap[i][j] > 0) tmpSum += coveredsingleexon[j];
    if (tmpSum == 1 && segcount[i] == 1)
      validone[i] = 1;
  }

  if (oneins.peReadDis == 0){
    oneins.peReadDis = 200;
    oneins.peReadDisSTD = 20;
  }
  vector<double> alpharange; 
  alpharange.push_back(oneins.peReadDis - 3 * oneins.peReadDisSTD);
  alpharange.push_back(oneins.peReadDis + 3 * oneins.peReadDisSTD);

  ones(validf, ncan);
  if (checkpairend && oneins.PETypes.size() > 0){
    for (long i = 0; i < ncan; ++i){
      candidate.assign(tap[i].begin(), tap[i].end());
      if (valid[i] == 0 ) continue;
      vector<long> sgcompatible;
      zeros(sgcompatible, oneins.SGTypes.size());
      for (long j = 0; j < oneins.SGTypes.size(); ++j)
        if (checkcompatible(candidate,oneins.SGTypes[j]))
          sgcompatible[j] = 1;
        else sgcompatible[j] = 0;
      // paired-end patterns
      vector<long> pecompatible, pepossible;
      for (long j = 0; j < oneins.PETypes.size(); ++j){
        pecompatible.push_back(sgcompatible[oneins.PETypes[j][0]] || sgcompatible[oneins.PETypes[j][1]]);
        pepossible.push_back(pecompatible[j] && (oneins.PETypes[j][0] != oneins.PETypes[j][1]));
      }
      vector<long> coverseg;
      zeros(coverseg, nExons);
      for (long j = 0; j < pepossible.size(); ++j){
        if (pepossible[j] == 0) continue;
        long sumFexon= 0, sumLexon = 0;
        vector<long> nfexon, nlexon;
        for (long k = 0; k <= oneins.SGTypes[oneins.PETypes[j][0]].size(); ++k){
          sumFexon += oneins.SGTypes[oneins.PETypes[j][0]][k];
          sumLexon += oneins.SGTypes[oneins.PETypes[j][1]][k];
          if (oneins.SGTypes[oneins.PETypes[j][0]][k] != 0) nfexon.push_back(oneins.SGTypes[oneins.PETypes[j][0]][k]);
          if (oneins.SGTypes[oneins.PETypes[j][1]][k] != 0) nlexon.push_back(oneins.SGTypes[oneins.PETypes[j][1]][k]);
        }

        if (sumFexon == 0 || sumLexon == 0) continue;
        long nlast = nfexon[nfexon.size() - 1];
        long nfirst = nlexon[0];
        vector<long> thispetheoryrange;
        for (long k = 0; k < alpharange.size(); ++k)
          thispetheoryrange.push_back(alpharange[k]);
        
        if (nfirst > nlast) {
          vector<long> middlemark;
          zeros(middlemark, nExons);
          for (long k = nlast + 1; k <= nfirst - 1; ++k)
            middlemark[k] = candidate[k];

          long gap = oneins.exonBoundary[nfirst][0] - oneins.exonBoundary[nlast][1] - 1;
          for (long k = 0;k < nExons; ++k)
            gap -= oneins.exonLen[k] * middlemark[k];
          for (long k = 0;k < thispetheoryrange.size(); ++k)
            thispetheoryrange[k] += gap;
        }
        long npossibledis;
        for (long k = 0; k < oneins.PETypesDistance[j].size(); ++k)
          if (oneins.PETypesDistance[j][k] >= thispetheoryrange[0] && oneins.PETypesDistance[j][k] <= thispetheoryrange[1])
            npossibledis += 1;
        if (npossibledis > 0){
          for (long k = 0; k < coverseg.size(); ++k)
            coverseg[k] = coverseg[k] + oneins.SGTypes[oneins.PETypes[j][0]][k] + oneins.SGTypes[oneins.PETypes[j][1]][k];
        }
      }
      long sumC = 0;
      for (long k = 0; k < candidate.size(); ++k)
        if (candidate[k] - coverseg[k] == 1)
          ++sumC;
      if (sumC > 0)
        validf[i] = 0;
      else validf[i] =1;
    }
  }
  orvalid.clear();
  for (long i = 0; i < valid.size(); ++i)
    orvalid.push_back((((valid[i] == 1) && validf[i] && (segcount[i] > 0)) || ((segcount[i] == 1) && (validone[i] > 0))));
  vector<vector<long> > tmp;
  tmp.assign(tap.begin(), tap.end());
  filtered.clear();
  for (long i = 0; i < orvalid.size(); ++i)
    if (orvalid[i] != 0) filtered.push_back(tmp[i]);
}

bool checkcompatible(vector<long> &candidate, vector<long> &juncpat){
  long nExons = candidate.size();
  long tmpSum1 = 0, tmpSum2 = 0;
  for (long i =0; i < candidate.size(); ++i){
    if (candidate[i] -juncpat[i] == -1) ++tmpSum1;
    tmpSum2 += juncpat[i];
  }
  if (tmpSum1 > 0 || tmpSum2 == 0) return false;
  vector<long> juncnz, juncmask;
  long tmpSum = 0;
  for (long i = 0; i < juncpat.size(); ++i)
    if (juncpat[i] != 0) juncnz.push_back(i);
  zeros(juncmask, nExons);
  for (long i = juncnz[0]; i <= juncnz[juncnz.size() - 1]; ++i)
    juncmask[i] = 1;
  for (long i = 0; i < juncpat.size(); ++i)
    if (juncpat[i] - juncmask[i] == -1)
      tmpSum += candidate[i];
  if (tmpSum > 0) return false;
  else return true;
}

void stackiso(vector<vector<long > > &rcount, long &index, long &minjunccount, long &nmaxiso, vector<double> &exonexp, double &minexpfrac, bool &countsonly, vector<long> &isintronretention, double &intronretentionfrac, vector<vector<long> > &refinfo, vector<vector<long> >&alliso, long &reachmax, vector<vector<long> >&imap){
  long nexon = rcount.size();
  memset(stack, 0, sizeof(long) * rcount.size());
  memset(stackprev, 0, sizeof(long) * rcount.size());
  memset(isvisited, 0, sizeof(long) * rcount.size());
  long stackpt = 0;

  stack[stackpt] = index;
  stackprev[stackpt] = -1;
  long rightnowsize = nmaxiso;
  long rightnowexpand = 2;

  if (countsonly == false)
    zeros(alliso, rightnowsize, nexon);
  else
    zeros(alliso, 1, nexon);

  long nisocount = 0;
  zeros(imap, 1, nexon);
  reachmax = 0;
  long loopcount = 0;
  while (stackpt >= 0){
    ++loopcount;
    //if (loopcount % 10 == 0) printf("\nloop count %ld", loopcount);
    long current = stack[stackpt];
    double thisexonexp = exonexp[current];
    imap[0][current] = 1;

    if (isvisited[stackpt] == 0) {
      isvisited[stackpt] = 1;
      bool haschildren = false;
      long oldstackpt = stackpt;
      vector<long> cdt;
      for (long i = 0; i < rcount[current].size(); ++i)
        if (i <= current) cdt.push_back(0);
        else cdt.push_back(rcount[current][i]);

      vector<MATLAB_SORT_ENTRY<long> > tmp;
      matlabSortWithIndex(cdt,tmp);
      vector<long> seled;
      for (long i = 0; i < tmp.size(); ++i){
        if (tmp[i].value >= minjunccount)
          seled.push_back(tmp[i].idx);
      }
      long currentmjc = minjunccount;
      while (seled.size() == 0 && currentmjc > 1){
        --currentmjc;
        if (currentmjc == 0) printf("zero here!");
        for (long i = 0; i < tmp.size(); ++i){
          if (tmp[i].value >= currentmjc)
            seled.push_back(tmp[i].idx);
        }
      }
      double maxcexp = 0, maxjexp = 0;
      if (seled.size() > 0){
        for (long i = 0 ; i < seled.size(); ++i){
          if (exonexp[seled[i]] > maxcexp) maxcexp = exonexp[seled[i]];
          if (rcount[current][seled[i]] > maxjexp) maxjexp = rcount[current][seled[i]];
        }
        vector<long> selbin;
        long sumselbin = 0;
        for (long i = 0 ; i < seled.size(); ++i){
          if ((exonexp[seled[i]] >= maxcexp * minexpfrac) && (rcount[current][seled[i]] >= maxjexp*minexpfrac)){
            selbin.push_back(1);
            sumselbin += 1;
          }
          else selbin.push_back(0);
        }
        if (sumselbin == 0){
          for (long i = 0 ; i < seled.size(); ++i)
            if (exonexp[seled[i]] >= maxcexp * minexpfrac)
              selbin[i] = 1;
            else selbin[i] = 0;
        }
        long nextid = current + 1;
        if ((find(seled.begin(), seled.end(), nextid) != seled.end()) && isintronretention[nextid] != 0){
          long rightmostrd = 0;
          for (long i = current + 1; i < rcount[current].size(); ++i)
            rightmostrd += rcount[current][i];
          long leftrd = rcount[current][nextid];
          if (exonexp[nextid] < maxcexp*intronretentionfrac){
            for (long i = 0; i < selbin.size(); ++i)
              selbin[i] = selbin[i] && (!(seled[i] == nextid));
          }
        }
        vector<long> possiblejunc, unionj;
        possiblejunc.assign(refinfo[current].begin(), refinfo[current].end());
        possiblejunc[current] = 0;
        for (long i = 0; i < possiblejunc.size(); ++i)
          if (possiblejunc[i] && cdt[i])
            unionj.push_back(i);
        set<long> selseled(unionj.begin(), unionj.end());
        for (long i = 0; i < selbin.size(); ++i)
          if (selbin[i] != 0)
            selseled.insert(seled[i]);
        if (selseled.size() >0){
          haschildren = true;
          long offset = 1;
          for (set<long>::iterator iter = selseled.begin(); iter != selseled.end(); ++iter){
            stack[stackpt + offset] = *iter;
            stackprev[stackpt + offset] = oldstackpt;
            isvisited[stackpt + offset] = 0;
            ++offset;
          }
          stackpt = stackpt + selseled.size();
        }
      }
      if (!haschildren){
        if (!countsonly && nisocount >= alliso.size()){
          long oldrightnowsize = rightnowsize;
          rightnowsize = min(rightnowsize* rightnowexpand, nmaxiso);
          if (rightnowsize == oldrightnowsize){
            reachmax = true;
            return;
          }
          printf("Expanding Matrix size to %ld\n", rightnowsize);
          vector<vector<long> >alliso2;
          alliso2.assign(alliso.begin(), alliso.end());
          zeros(alliso, rightnowsize, nexon);
          for (long i = 0; i < alliso2.size(); ++i)
            alliso[i].assign(alliso2[i].begin(), alliso[2].end());
        }
        if (nisocount >= nmaxiso){
          reachmax = true;
          break;
        }
        else {
          ++nisocount;
          if (!countsonly){
            long stcurpt = stackpt;
            long stcur = current;
            while (1){
              alliso[nisocount - 1][stcur] = 1;
              stcurpt = stackprev[stcurpt];
              if (stcurpt == -1) break;
              stcur = stack[stcurpt];
            }
          }
        }
      }

    }
    else {
      --stackpt;
    }

  }
  //printf("\nloop count %ld", loopcount);
  if (!countsonly){
    alliso.erase(alliso.begin() + nisocount, alliso.end());
  }
  else {
    alliso.clear();
    vector<long> tmp3;
    tmp3.push_back(nisocount);
    alliso.push_back(tmp3);
  }
}

void getlastoffirstexon(vector<long> &candidate,vector<vector<long> > &bound, long &nsegs, long &segid, long &lastsegid) {
  vector<long> none;
  for (long i = 0; i < candidate.size(); ++i)
    if (candidate[i] !=0) none.push_back(i);
  if (none.size() == 0){
    nsegs = 0; 
    segid = 0; 
    lastsegid = 0;
  }
  if (none.size() == 1){
    nsegs = 1;
    segid = none[0];
    lastsegid = none[0];
  }
  nsegs = 1;
  segid = none[0];
  for (long i = 0; i < none.size() - 1; ++i) {
    if (bound[none[i]][1] == bound[none[i + 1]][0] - 1){
      if (nsegs == 1) segid = none[i + 1];
    }
    else nsegs = nsegs + 1;
  }

  if (nsegs == 1){
    lastsegid = none[0];
    return;
  }

  lastsegid = none[none.size() - 1];
  for (long i = none.size() - 2; i >=0; --i){
    if (bound[none[i]][1] == bound[none[i + 1]][0] - 1)
      lastsegid = none[i + 1];
    else {
      lastsegid = none[i + 1];
      break;
    }
  }

}

void getExonExpLv(NewInstance &ins, vector<double> &explv){
  for (long i = 0; i < ins.rcount.size(); ++i){
    explv.push_back(((double)ins.rcountRowSum[i] + ins.rcountColSum[i] ) / 2);
    explv[i] /= ins.exonLen[i];
  }
}

void getExonWeight(vector<long> &exonlength, long &readlength, vector<double> &exonweight){
  double maxreadlenfraction = 2, sigmoidcoef = 5;
  double a = maxreadlenfraction * readlength;
  for (long i = 0; i < exonlength.size(); ++i){
    if (exonlength[i] >= a) exonweight.push_back(1);
    else exonweight.push_back(0);
    if (exonlength[i] < a) exonweight[i] += double(exonlength[i]) / a;
  }
}

void getExonJunctionExpLv(NewInstance &ins, vector<double> &explv, vector<double> &rightexplv, vector<double> &leftexplv) {
  double nlen = ins.readLen * 2;
  zeros(rightexplv, ins.rcount[0].size());
  zeros(leftexplv, ins.rcount.size());
  for (long i = 0; i < ins.rcount.size(); ++i)
    for (long j = i + 1; j < ins.rcount[i].size(); ++j){
      rightexplv[i] += ins.rcount[i][j];
      leftexplv[j] += ins.rcount[i][j];
  }
  explv.clear();
  for (long i = 0; i < ins.rcount.size(); ++i){
    explv.push_back((rightexplv[i] + leftexplv[i]) / 2.0 / nlen);
    rightexplv[i] /= nlen;
    leftexplv[i] /= nlen;
  }

}