c_cpp/etc/tiny/duptc.c

   1 #include <stdlib.h> //EXIT_FAILURE
   2 #include <stdio.h>
   3 #include <float.h>
   4 #include <math.h>
   5
   6 //#define SQRT2 1.4142135623730950488016887242097
   7 //const long double Sqrt2=sqrtl(2);
   8
   9 double InsMean, InsSTD, ReadsDep;
  10
  11 /* http://en.wikipedia.org/wiki/Normal_distribution#Numerical_approximations_for_the_normal_CDF
  12 如果用Ncdf(x+1)-Ncdf(x)，则当std<1时，最终结果会减半。
  13 还原回Ncdf(x+0.5)-Ncdf(x-0.5)，则如同预期地趋近SE的结果。
  14  */
  15 static inline long double NormalCDFdiff(double x, double m, double s) {
  16         //return 0.5 + 0.5 * erfl( (x - m) / (s * Sqrt2) );
  17         //return 0.5 * ( erfl( (x+1 - m) / (s * Sqrt2) ) - erfl( (x - m) / (s * Sqrt2) ) );
  18         return 0.5 * ( erfcl( (m-x-0.5) / (s * sqrtl(2)) ) - erfcl( (m-x+0.5) / (s * sqrtl(2)) ) );
  19 }
  20
  21 int main (int argc, char *argv[]) {
  22         if (argc!=4) {
  23                 fprintf(stderr,"Usage: %s <Reads_Depth> <InsertSize_Mean> <InsertSize_STD>\n",argv[0]);
  24                 exit(EXIT_FAILURE);
  25         }
  26         ReadsDep=atof(argv[1]);
  27         InsMean=atof(argv[2]);
  28         InsSTD=atof(argv[3]);
  29         if (ReadsDep<0 || InsMean<0 || InsSTD<=0) {
  30                 fputs("[x]Input Error !\n",stderr);
  31                 exit(EXIT_FAILURE);
  32         }
  33         double SEdup = 1-exp(-ReadsDep);
  34         printf("ReadsDep=(%.9g) InsMean=(%.12g) InsSTD=(%.12g)\nSE theory dup. rate: %.20e\n",
  35                 ReadsDep,InsMean,InsSTD,SEdup);
  36         long double sumdup=0;
  37         long double thisdup,Ci;
  38         long long int ins=1;
  39         do {
  40                 Ci = (long double)ReadsDep * NormalCDFdiff(ins,InsMean,InsSTD);
  41                 thisdup=Ci*(1-exp(-Ci));
  42                 sumdup += thisdup;
  43                 //printf("%d %.16Le %.16Le\n",ins,thisdup,sumdup);
  44                 ++ins;
  45         } while (ins<=2*InsMean+6*InsSTD || thisdup>LDBL_MIN);
  46         sumdup /= (long double)ReadsDep;
  47         printf("PE theory dup. rate: %.20Le (%.6Lf %% of SE)\n",sumdup,100*sumdup/SEdup);
  48         exit(EXIT_SUCCESS);
  49 }
  50