[perf] Reuse memory for simulated tag arrays

bcf_utils.cpp

@@ -1,5 +1,148 @@
 #include "bcf_utils.h"
 
+#include "io.h"
+
+sim_rec::sim_rec(bcf_hdr_t* in_hdr){
+
+	this->nSamples=bcf_hdr_nsamples(in_hdr);
+
+	this->set_hdr(in_hdr);
+
+	if(args->in_mps_depths!=NULL){
+		this->mps_depths=read_depthsFile(args->in_mps_depths, this->nSamples);
+
+		fprintf(stderr, "\n");
+		for (int sample_i=0; sample_i<this->nSamples; sample_i++) {
+			fprintf(stderr, "Individual %d mean per-site depth is set to %f\n", sample_i,this->mps_depths[sample_i]);
+		}
+	}
+
+
+	this->gt_arr=NULL; // allocated in-place
+
+	this->dp_arr= bcf_tag_alloc<int32_t>(DP, 0);
+
+	this->gl_vals = (double*) malloc(this->nSamples*SIM_NGTS*sizeof(double));
+	ASSERT(NULL!=this->gl_vals);
+	for (int i=0;i<nSamples*SIM_NGTS; ++i){
+		this->gl_vals[i]=-0.0;
+	}
+
+
+	this->gl_arr= bcf_tag_alloc<float>(GL,-0.0,nSamples*SIM_NGTS);
+
+	if(1==args->addGP){
+		this->gp_arr= bcf_tag_alloc<float>(GL,-0.0,nSamples*SIM_NGTS);
+	}
+
+	if(1==args->addPL){
+		this->pl_arr = bcf_tag_alloc<int32_t>(PL,0, nSamples*SIM_NGTS);
+	}
+
+	if(1==args->addQS){
+		this->qs_arr= bcf_tag_alloc<float>(QS,0.0);
+	}
+
+	if(1==args->addI16){
+		this->i16_arr= bcf_tag_alloc<float>(I16,0.0);
+	}
+
+
+}
+
+sim_rec::~sim_rec(){
+
+	bcf_hdr_destroy(this->hdr);
+
+	free(this->gt_arr);
+	this->gt_arr=NULL;
+
+	//TODO delme
+	free(this->gl_vals);
+	this->gl_vals=NULL;
+
+	free(this->gl_arr);
+	this->gl_arr=NULL;
+
+	if(1==args->addGP){
+		free(this->gp_arr);
+		this->gp_arr=NULL;
+	}
+
+	if(1==args->addPL){
+		free(this->pl_arr);
+		this->pl_arr=NULL;
+	}
+
+	if(1==args->addQS){
+		free(this->qs_arr);
+		this->qs_arr=NULL;
+	}
+
+	if(1==args->addI16){
+		free(this->i16_arr);
+		this->i16_arr=NULL;
+	}
+
+
+	free(this->dp_arr);
+	this->dp_arr=NULL;
+
+}
+
+
+void sim_rec::set_hdr(bcf_hdr_t* in_hdr){
+
+	this->hdr = bcf_hdr_dup(in_hdr);
+
+	char *DATE_TAG=NULL;
+
+	ASSERT(asprintf(&DATE_TAG, "##fileDate=%s", args->datetime)>0);
+	ASSERT(0==bcf_hdr_append(this->hdr, DATE_TAG));
+	free(DATE_TAG);
+	DATE_TAG=NULL;
+
+	char *SOURCE_TAG=NULL;
+	ASSERT(asprintf(&SOURCE_TAG, "##source=%s", args->command)>0);
+	ASSERT(0==bcf_hdr_append(this->hdr, SOURCE_TAG));
+	free(SOURCE_TAG);
+	SOURCE_TAG=NULL;
+
+	char *SOURCE_VERSION_TAG;
+	ASSERT(asprintf(&SOURCE_VERSION_TAG, "##source=vcfgl version: %s",VCFGL_VERSION)>0);
+	ASSERT(0==bcf_hdr_append(this->hdr, SOURCE_VERSION_TAG));
+	free(SOURCE_VERSION_TAG);
+	SOURCE_VERSION_TAG=NULL;
+
+
+	ASSERT(0==bcf_hdr_append(this->hdr,bcf_tags[GL].hdr));
+
+
+	if(1==args->addGP){
+		ASSERT(0==bcf_hdr_append(this->hdr,bcf_tags[GP].hdr));
+	}
+
+	if(1==args->addPL){
+		ASSERT(0==bcf_hdr_append(this->hdr,bcf_tags[PL].hdr));
+	}
+
+	if(-999!=args->mps_depth){
+		bcf_tag_set_size(DP, nSamples);
+		ASSERT(0==bcf_hdr_append(this->hdr,bcf_tags[DP].hdr));
+
+		if(1==args->addI16){
+			ASSERT(0==bcf_hdr_append(this->hdr,bcf_tags[I16].hdr));
+		}
+
+		if(1==args->addQS){
+			ASSERT(0==bcf_hdr_append(this->hdr,bcf_tags[QS].hdr));
+		}
+	}
+
+
+
+}
+
 void bcf_tag_set_size(enum bcf_tag t, const int size){
 	int bcf_tag_size = bcf_tags[t].n;
 

bcf_utils.h

@@ -1,12 +1,66 @@
 #ifndef __BCF_UTILS__
 #define __BCF_UTILS__
 
+
+
+#include <htslib/kstring.h> // kstring_t
+
+#include <htslib/vcf.h>
+#include <htslib/vcfutils.h>
+#include <limits>
+
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <ctype.h>
+#include <sys/stat.h>
+#include <stdio.h>
+
 #include "shared.h"
+#include "version.h"
 
 
 /* ========================================================================== */
 /* /BEGIN/ BCF UTILS ======================================================== */
 
+
+typedef struct sim_rec{
+
+	bcf_hdr_t* hdr=NULL;
+
+	bcf1_t* rec=NULL;
+	bcf1_t* blank_rec=NULL;
+
+	kstring_t* alleles_str=NULL;
+
+	int nSamples=0;
+	int nSites=0;
+
+	int site_i=-1;
+
+	//TODO delme
+	double* gl_vals=NULL;
+
+	double* mps_depths=NULL;
+
+	int32_t* dp_arr=NULL;
+	int32_t* gt_arr=NULL;
+
+	float* gl_arr=NULL;
+	float* gp_arr=NULL;
+	int32_t* pl_arr=NULL;
+
+	float* qs_arr=NULL;
+	float* i16_arr=NULL;
+
+	sim_rec(bcf_hdr_t* in_hdr);
+	~sim_rec();
+
+	void set_hdr(bcf_hdr_t* in_hdr);
+
+
+}sim_rec;
+
 /* -> BCF TAGS ---------------------------------------------------------------*/
 
 // modified from source: bcftools/tag2tag
@@ -22,8 +76,7 @@ typedef struct
 	int type;
 	const char *str=NULL;
 	const char *hdr=NULL;
-}
-bcf_tag_t;
+}bcf_tag_t;
 
 enum bcf_tag{DP, GT, GL, GP, PL, QS, I16};
 extern bcf_tag_t bcf_tags[7];
@@ -32,20 +85,43 @@ extern bcf_tag_t bcf_tags[7];
 template <typename T> T* bcf_tag_alloc(enum bcf_tag t){
 	const int bcf_tag_size = bcf_tags[t].n;
 	ASSERT(bcf_tag_size>0);
-	T* obj = (T*) malloc(bcf_tag_size * sizeof(T));
-	ASSERT(NULL!=obj);
-	return obj;
+	T* arr = (T*) malloc(bcf_tag_size * sizeof(T));
+	ASSERT(NULL!=arr);
+	return arr;
 }
 
 template <typename T> T* bcf_tag_alloc(enum bcf_tag t, T init_val){
 	const int bcf_tag_size = bcf_tags[t].n;
 	ASSERT(bcf_tag_size>0);
-	T* obj = (T*) malloc(bcf_tag_size * sizeof(T));
+	T* arr = (T*) malloc(bcf_tag_size * sizeof(T));
 	for(int i=0;i<bcf_tag_size;++i){
-		obj[i]=init_val;
+		arr[i]=init_val;
+	}
+	ASSERT(NULL!=arr);
+	return arr;
+}
+
+template <typename T> void bcf_tag_reset(T* arr, enum bcf_tag t, T init_val, const int size){
+
+	ASSERT(size>0);
+	bcf_tags[t].n=size;
+
+	ASSERT(NULL!=arr);
+	for(int i=0;i<size;++i){
+		arr[i]=init_val;
+	}
+	return;
+}
+
+template <typename T> void bcf_tag_reset(T* arr, enum bcf_tag t, T init_val){
+
+	const int size=bcf_tags[t].n;
+
+	ASSERT(NULL!=arr);
+	for(int i=0;i<size;++i){
+		arr[i]=init_val;
 	}
-	ASSERT(NULL!=obj);
-	return obj;
+	return;
 }
 
 template <typename T> T* bcf_tag_alloc(enum bcf_tag t, T init_val, const int size){
@@ -57,12 +133,12 @@ template <typename T> T* bcf_tag_alloc(enum bcf_tag t, T init_val, const int siz
 	}
 
 	ASSERT(size>0);
-	T* obj = (T*) malloc(size * sizeof(T));
+	T* arr = (T*) malloc(size * sizeof(T));
 	for(int i=0;i<size;++i){
-		obj[i]=init_val;
+		arr[i]=init_val;
 	}
-	ASSERT(NULL!=obj);
-	return obj;
+	ASSERT(NULL!=arr);
+	return arr;
 }
 
 

estimator.cpp

@@ -1,9 +1,5 @@
 #include "estimator.h"
 
-#include <stdio.h>
-#include <math.h>
-
-
 
 void gl_log10(int base, double errate, double *like){
 

estimator.h

@@ -3,6 +3,10 @@
 
 #include "lut.h"
 
+#include <stdio.h>
+#include <math.h>
+
+
 
 
 void gl_log10(int base, double errate, double *like);