Skip to content

Commit 64bdc5b

Browse files
committed
Add the -n, --nbins option. Note var2 is interpolated, maybe could collect more experimental data points in future
1 parent a53419b commit 64bdc5b

File tree

1 file changed

+31
-4
lines changed

1 file changed

+31
-4
lines changed

plugins/afs.c

+31-4
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* The MIT License
22
3-
Copyright (c) 2024 Genome Research Ltd.
3+
Copyright (c) 2024-2025 Genome Research Ltd.
44
55
Author: Petr Danecek <[email protected]>
66
@@ -802,20 +802,47 @@ static void merge_add_batch(args_t *args, const char *fname)
802802
}
803803
static void init_var2(args_t *args, batch_t *batch)
804804
{
805+
805806
if ( !strcmp(args->recalc_type,"hc") )
806807
{
807-
if ( batch->nbins!=N_BINS ) error("todo: --nbins %d != %d\n",batch->nbins,N_BINS);
808808
batch->var2 = malloc(sizeof(*batch->var2)*batch->nbins);
809+
809810
int i;
810-
for (i=0; i<batch->nbins; i++) batch->var2[i] = var2[i];
811+
if ( batch->nbins==N_BINS )
812+
{
813+
for (i=0; i<batch->nbins; i++) batch->var2[i] = var2[i];
814+
}
815+
else
816+
{
817+
// intrapolate data points for nbins!=N_BINS
818+
819+
batch->var2[0] = var2[0];
820+
batch->var2[batch->nbins - 1] = var2[N_BINS - 1];
821+
822+
double dx = (double)1/(batch->nbins-1);
823+
double DX = (double)1/(N_BINS-1);
824+
int J=1;
825+
for (i=1; i<batch->nbins - 1; i++)
826+
{
827+
double x = i*dx;
828+
while ( J*DX < x ) J++;
829+
assert( J < N_BINS );
830+
double X = (J-1)*DX;
831+
assert( X <= x );
832+
if ( x==X )
833+
batch->var2[i] = var2[J-1];
834+
else
835+
batch->var2[i] = var2[J-1] + (var2[J]-var2[J-1])*(x-X)/DX;
836+
}
837+
}
811838
return;
812839
}
813840
if ( !strcmp(args->recalc_type,"data") )
814841
{
815842
// variance will be calculated from the data
816843
return;
817844
}
818-
error("todo: --recalc %s\n",args->recalc_type);
845+
error("todo: --recalc %s\n",args->recalc_type);
819846
}
820847
static int merge(args_t *args)
821848
{

0 commit comments

Comments
 (0)