c_cpp/readscorr/sdleftTF.c

   1 #include <stdint.h> //uint64_t
   2 #include <stdlib.h> //calloc
   3 #include <math.h>
   4 #include "sdleft.h"
   5 //#include "sdleftTF.h"
   6 #include <stdio.h>
   7 /*
   8 Well, let's try "Template" in C with #define and pointer to function.
   9 Since the function body is the same, I prefer to call it "template" instead of "overload".
  10
  11 Here, we select typeof(Count) on SUM(Count)==dleftobj->ItemInsideAll
  12 */
  13
  14 #define TOKENPASTE(x, y) x ## y
  15 #define TOKENPASTE2(x, y) TOKENPASTE(x, y)
  16 #define CAT0(x) #x
  17 #define CAT(x) CAT0(x)
  18
  19 #ifdef USEUINT16
  20     #define THETYPE uint16_t
  21     #define DILTYPE uint32_t
  22     #define FLOTYPE double
  23 #elif defined USEUINT32
  24     #define THETYPE uint32_t
  25     #define DILTYPE uint64_t
  26     #define FLOTYPE double
  27 #elif defined USEUINT64
  28     #define THETYPE uint64_t
  29     #define DILTYPE uint128_t
  30     #define FLOTYPE float128
  31 #elif defined PUBLIC
  32     // NA
  33 #else
  34     #error Must define USEUINT{16,32,64} or PUBLIC before compilering !
  35 #endif
  36
  37 #define ADDSUFFIX(x) TOKENPASTE2(x ## _,THETYPE)
  38
  39 SDLeftStat_t *dleft_stat_uint16_t(SDLeftArray_t * const dleftobj, FILE *);
  40 SDLeftStat_t *dleft_stat_uint32_t(SDLeftArray_t * const dleftobj, FILE *);
  41 SDLeftStat_t *dleft_stat_uint64_t(SDLeftArray_t * const dleftobj, FILE *);
  42
  43 #ifndef PUBLIC  // USEUINT{16,32,64}
  44 SDLeftStat_t * ADDSUFFIX(dleft_stat) (SDLeftArray_t * const dleftobj, FILE *stream) {
  45 fprintf(stderr,"[!]Count with[%s]\n",CAT(THETYPE));
  46     SDLeftStat_t *pSDLeftStat = malloc(sizeof(SDLeftStat_t));
  47     THETYPE * const pCountHistArray = calloc(sizeof(THETYPE),1+dleftobj->maxCountSeen);
  48     size_t totalDLAsize = dleftobj->SubItemCount * dleftobj->itemByte * dleftobj->ArraySize;
  49     //size_t firstlevelDLAitemsize = SDLA_ITEMARRAY*dleftobj->itemByte;
  50     const unsigned char * const pDLA = dleftobj->pDLA;
  51     THETYPE Item_CountBits=0;  // set value in case SDLA_ITEMARRAY*dleftobj->itemByte == 0 (EP ?)
  52     uint128_t theItem;
  53     for (size_t i=0;i<totalDLAsize;i+=dleftobj->itemByte) {
  54         //const unsigned char * pChunk = pDLA + i;
  55         theItem = 0;
  56         for (uint_fast8_t j=0;j<dleftobj->itemByte;j++) {
  57             theItem |= ((uint128_t)*(pDLA + i + j)) << (j*8u);
  58         }
  59         Item_CountBits = theItem & dleftobj->Item_CountBitMask;
  60         ++pCountHistArray[Item_CountBits];
  61         //++HistSum;
  62     }
  63     //THETYPE HistSum=0;    // HistSum == dleftobj->ItemInsideAll
  64     DILTYPE HistSumSquared=0;
  65     double SStd;    // We need to return in a fixed type for printf
  66     for (size_t p=1;p<=dleftobj->maxCountSeen;p++) {
  67         //HistSum += pCountHistArray[p];
  68         HistSumSquared += pCountHistArray[p] * pCountHistArray[p];
  69     }
  70     //http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
  71     SStd = sqrtl( ( (long double)HistSumSquared-((long double)dleftobj->ItemInsideAll*(long double)dleftobj->ItemInsideAll/(long double)dleftobj->maxCountSeen) ) / (long double)(dleftobj->maxCountSeen -1) );
  72     pSDLeftStat->HistSStd = SStd;
  73     pSDLeftStat->HistMean = (double)dleftobj->ItemInsideAll / (double)dleftobj->maxCountSeen;
  74     pSDLeftStat->HistMaxCntVal = 1; //later
  75     pSDLeftStat->HistMaxHistVal = 1; //later
  76     fprintf(stream,"#Kmer_real_count: %ld\n#Kmer_count_hist: %ld\n#Kmer_depth_mean: %f\n#Kmer_depth_sStd: %f\n\n#Kmer_frequence\tHist_value\tKmer_count\tHist_ratio\n",
  77         dleftobj->ItemInsideAll,dleftobj->maxCountSeen,pSDLeftStat->HistMean,SStd);
  78     for (size_t p=1;p<=dleftobj->maxCountSeen;p++) {
  79         fprintf(stream,"%zu\t%lu\t%lu\t%g\n",p,(uint64_t)pCountHistArray[p],
  80             (uint64_t)pCountHistArray[p]*(uint64_t)p,(double)pCountHistArray[p]/(double)dleftobj->ItemInsideAll);
  81     }
  82     free(pCountHistArray);
  83     // deal(*pextree);
  84     return pSDLeftStat;
  85 }
  86 #else   // PUBLIC
  87 //G_SDLeftArray_IN *pf = dleft_stat_uint16_t;
  88 SDLeftStat_t * dleft_stat(SDLeftArray_t * const dleftobj, FILE *stream) {
  89     G_SDLeftArray_IN *pf;
  90     if (dleftobj->ItemInsideAll <= UINT16_MAX) {
  91         pf = dleft_stat_uint16_t;
  92     } else if (dleftobj->ItemInsideAll <= UINT32_MAX) {
  93         pf = dleft_stat_uint32_t;
  94     } else {    // uint64_t ItemInsideAll
  95         pf = dleft_stat_uint64_t;
  96     }
  97     return (*pf)(dleftobj,stream);
  98 }
  99 #endif
 100
 101 /*
 102 Not using now.
 103
 104 customobjects = $(addprefix $(OBJDIR), sdleftTFuint16.o sdleftTFuint32.o sdleftTFuint64.o sdleftTFpublic.o )
 105
 106 $(OBJDIR)sdleftTFuint16.o:
 107    $(CC) -std=gnu99 $(MAKEARG) -D USEUINT16 -c sdleftTF.c -o $@ > $(@:.o=.asm)
 108 $(OBJDIR)sdleftTFuint32.o:
 109    $(CC) -std=gnu99 $(MAKEARG) -D USEUINT32 -c sdleftTF.c -o $@ > $(@:.o=.asm)
 110 $(OBJDIR)sdleftTFuint64.o:
 111    $(CC) -std=gnu99 $(MAKEARG) -D USEUINT64 -c sdleftTF.c -o $@ > $(@:.o=.asm)
 112 $(OBJDIR)sdleftTFpublic.o:
 113    $(CC) -std=gnu99 $(MAKEARG) -D PUBLIC -c sdleftTF.c -o $@ > $(@:.o=.asm)
 114 */