Saddle approx speedup

python/PascalX/genescorer.py

@@ -616,6 +616,15 @@ def plot_Manhattan(self,ScoringResult=None,region=None,sigLine=0,logsigThreshold
             plt.axhline(y=-np.log10(sigLine), color='r', linestyle='dotted')
 
 
+
+    def clean(self):
+        """
+        Removes scores obtained from previous runs
+        
+        """
+        self._SCORES = {}
+        self._SKIPPED = {}
+
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 
@@ -816,9 +825,9 @@ def _calcAndFilterEV(self,C):
             return None
 
         F = L > 0
+        L = L[F][::-1]
 
-        if len(F) > 0:
-            L = L[F][::-1]
+        if len(L) > 0:
             N_L = []
 
             # Leading EV
@@ -943,10 +952,10 @@ def _scoremain(self,gene,unloadRef=False,method='saddle',mode='auto',reqacc=1e-1
                     pbar.update(1)
 
 
-        with lock:
+        with lock:  
+            pbar.set_postfix_str("done")
             pbar.close()
-            #pbar.set_description(label+"(done          )")
-
+
         return RESULT,FAIL,TOTALFAIL
 
     def score(self,gene,parallel=1,unloadRef=False,method='saddle',mode='auto',reqacc=1e-100,intlimit=1000000,nobar=False,autorescore=False,keep_idx=None):

python/PascalX/wchissum.py

@@ -336,7 +336,9 @@ def onemin_cdf_satterthwaite(X, lb, mode='auto'):
     toc = time.time()
 
     return [max(rf,res), 0, round(toc-tic,5)]
-
+
+
+
 def onemin_cdf_pearson(X, lb, mode='auto'):
     tic = time.time()
 
@@ -359,9 +361,9 @@ def onemin_cdf_pearson(X, lb, mode='auto'):
         rf = 1e-15
 
     toc = time.time()
-
+     
     return [max(rf,res), 0, round(toc-tic,5)]
-
+    
 def onemin_cdf_saddle(X, lb, mode='auto'):
     tic = time.time()
 

python/wchissum.cpp

@@ -1246,6 +1246,7 @@ double oneminwchissum_m1nc0_saddle_500d(double* lambda, int N, double X) {
     }   
 }
 
+/*
 extern "C"
 double oneminwchissum_m1nc0_saddle_auto(double* lambda, int N, double X) {
     double res = oneminwchissum_m1nc0_saddle(lambda,N,X);
@@ -1268,6 +1269,82 @@ double oneminwchissum_m1nc0_saddle_auto(double* lambda, int N, double X) {
     
     return res;
 }
+*/
+
+
+extern "C"
+double oneminwchissum_m1nc0_saddle_auto(double* lambda, int N, double X) {
+    double sum = lambda[0];
+
+    // find maxb
+    double ma = 1./lambda[0];
+    for(int i=1;i<N;i++) {
+        double tmp = 1./lambda[i];
+        sum += lambda[i];
+
+        if(tmp < ma) {
+            ma = tmp;
+        }
+    }
+    ma *= 0.5;
+
+    // Do not use in unstable regime
+    if (abs(sum-X)/X < 1e-5) {
+        return -1;
+    }
+
+    const int digits = std::numeric_limits<double>::digits; 
+    int get_digits = static_cast<int>(digits * 0.6);
+
+    // Solve for zeta 
+    const boost::uintmax_t maxit = 10000;
+    boost::uintmax_t it = maxit;
+    try {
+        double zeta = boost::math::tools::newton_raphson_iterate(
+        [lambda,N,X](double x) {
+            return std::make_pair(K1(lambda,N,x)-X,K2(lambda,N,x));
+        },
+        -0.5*N/X,-0.5*N/X-1,ma,
+        get_digits, it
+        );
+
+        // Calc parameters    
+        double v = zeta*sqrt(K2(lambda,N,zeta));
+        double w = sign(zeta)*sqrt(2*(zeta*X - K(lambda,N,zeta)));
+        double z = (w + log(v/w)/w)/sqrt(2.);
+
+        // Calc cdf using error function
+        double Z = z;
+
+        double res = ( 0.5*(1 - erf( Z )) );
+
+        if (res < 1e-15) {
+            float128 Z(z);
+            res = ( 0.5*(1 - erf( Z )) ).convert_to<double>();
+
+            if (res < 1e-32) {
+                number<cpp_bin_float<100>> Z(z);
+                res = ( 0.5*(1 - erf( Z )) ).convert_to<double>();
+
+                if (res < 1e-98) {
+                    number<cpp_bin_float<200>> Z(z);
+                    res = ( 0.5*(1 - erf( Z )) ).convert_to<double>();
+
+                    if (res < 1e-195) {
+                        number<cpp_bin_float<300>> Z(z);
+                        res = ( 0.5*(1 - erf( Z )) ).convert_to<double>();
+                    }
+                }
+            }
+        }
+
+        // Calc and return solution
+        return res;
+
+    } catch(...) {
+        return -1;
+    }   
+}