Skip to content

Commit

Permalink
Synchronize with upstream haplo.stats, fix some redundant checks (#196
Browse files Browse the repository at this point in the history
)

* fix null checks

* sync with upstream `haplo.stats`

* fix type when setting of `printf` widths

* fix missing argument in `gthwe`

* `error` apparently missing from MacOS C library
  • Loading branch information
alexlancaster authored Feb 25, 2024
1 parent 718f675 commit 09aba6b
Show file tree
Hide file tree
Showing 5 changed files with 46 additions and 50 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/codeql.yml
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ name: "CodeQL"

on:
push:
branches: [ "main" ]
branches: [ "main", ]
pull_request:
branches: [ "main" ]
schedule:
Expand Down
2 changes: 1 addition & 1 deletion src/emhaplofreq/emhaplofreq.c
Original file line number Diff line number Diff line change
Expand Up @@ -973,7 +973,7 @@ int main_proc(
/* remove the trailing haplotype separator character in haplotype */
/* FIXME: the need for a trailing character should be removed when creating haplotype */
#ifdef XML_OUTPUT
xmlfprintf(fp_out, "<haplotype name=\"%.*s\"><frequency>%.5f</frequency><numCopies>%.1f</numCopies></haplotype>\n", (strlen(haplo[i]) - 1), haplo[i], freq_zero[i], freq_zero[i]*2.0*n_recs);
xmlfprintf(fp_out, "<haplotype name=\"%.*s\"><frequency>%.5f</frequency><numCopies>%.1f</numCopies></haplotype>\n", ((int) strlen(haplo[i]) - 1), haplo[i], freq_zero[i], freq_zero[i]*2.0*n_recs);
#else
xmlfprintf(fp_out, "%3d %12.5f %8.1f %.*s\n", j, freq_zero[i], freq_zero[i]*2.0*n_recs, (strlen(haplo[i]) - 1), haplo[i]);
#endif
Expand Down
2 changes: 1 addition & 1 deletion src/gthwe/hwe.c
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ int main(int argc, char *argv[])
size = sample.size;

/* pass the parsed variables to do the main processing */
run_data(genotypes, allele_array, no_allele, total, step, group, size, title, outfile, 1);
run_data(genotypes, allele_array, no_allele, total, step, group, size, title, outfile, 1, 0);

free(genotypes);
free(allele_array);
Expand Down
84 changes: 41 additions & 43 deletions src/haplo-stats/haplo_em_pin.c
Original file line number Diff line number Diff line change
@@ -1,11 +1,5 @@
/* $Author: sinnwell $ */
/* $Date: 2013/01/14 19:10:42 $ */
/* $Header: /projects/genetics/cvs/cvsroot/haplo.stats/src/haplo_em_pin.c,v 1.18 2013/01/14 19:10:42 sinnwell Exp $ */
/* $Locker: $ */
/*
* $Log:
*/
/*
/* $Date: 2013/01/14 19:10:42 $
*License:
*
*Copyright 2003 Mayo Foundation for Medical Education and Research.
Expand All @@ -27,7 +21,7 @@
*
*Daniel J. Schaid, Ph.D.
*Division of Biostatistics
*Harwick Building û Room 775
*Harwick Building Room 775
*Mayo Clinic
*200 First St., SW
*Rochester, MN 55905
Expand All @@ -47,10 +41,8 @@
#include <nmath.h>
#include "haplo_em_pin.h"


/* Progressive insertion of loci into haplotypes with EM algorithm */


/*************** Global vars ******************************************************/

static int n_loci, *loci_used; /* used for qsort functions */
Expand Down Expand Up @@ -107,7 +99,6 @@ void haplo_em_pin(

HAP **hap_list; /* List of all haplotypes = array of pointers to hap structs */
HAPUNIQUE **u_hap_list; /* List of unique haplotypes */
HAP *h1, *h2;

/* convert from S vecs to C structures */

Expand Down Expand Up @@ -227,7 +218,7 @@ void haplo_em_pin(


if(*verbose){
REprintf("\nhap_list after insert batch %d & set_post, before code haplo\n\n",n_batch);
REprintf("\nhap_list after insert batch %i & set_post, before code haplo\n\n",n_batch);
write_hap_list(hap_list, n_hap);
}

Expand All @@ -243,7 +234,7 @@ void haplo_em_pin(
qsort(hap_list, n_hap, sizeof(HAP *), cmp_subId_hapPairId);

if(*verbose){
REprintf("\nhap_list after code haplo, before EM: %d\n\n",n_batch);
REprintf("\nhap_list: %4i after code haplo, before EM\n",n_batch);
write_hap_list(hap_list, n_hap);
}

Expand Down Expand Up @@ -287,31 +278,32 @@ void haplo_em_pin(

} /* end of EM loop */

if( (*converge)==0){
/* FIXME: had to remove the PROBLEM macro to get to compile on MacOS */
REprintf("failed to converge for batch %i after %i iterations", n_batch, n_iter);
// PROBLEM "failed to converge for batch %i after %i iterations", n_batch, n_iter
// RECOVER(NULL_ENTRY);
}
if( (*converge)==0){
/* Replaced with call to REprintf for 4.1.x
PROBLEM "failed to converge for batch %i after %i iterations", n_batch, n_iter
RECOVER(NULL_ENTRY);
*/
REprintf("failed to converge for batch %i after %i iterations", n_batch, n_iter);
}


divideKeep(hap_list, n_hap, &len_hap_list);
n_hap = len_hap_list;
divideKeep(hap_list, n_hap, &len_hap_list);
n_hap = len_hap_list;


if(*verbose){
REprintf("\nhap_list after EM and after divideKeep: %d \n\n",n_trim);
write_hap_list(hap_list, n_hap);
}
if(*verbose){
/*REprintf("\nhap_list after EM and after divideKeep \n\n",n_trim);*/
write_hap_list(hap_list, n_hap);
}


/* update priors, in case haplos were trimmed during posteior calculations */

hap_prior(n_hap, hap_list, prior, n_u_hap, *min_prior);

if(*verbose){
if( (*converge)==1){
REprintf("\n\nConverged after batch insertion, lnlike = %8.5f, n_iter = %i\n\n", lnlike, n_iter);
/* update priors, in case haplos were trimmed during posteior calculations */
hap_prior(n_hap, hap_list, prior, n_u_hap, *min_prior);
if(*verbose){
if( (*converge)==1){
REprintf("\n\nConverged after batch insertion, lnlike = %8.5f, n_iter = %i\n\n", lnlike, n_iter);
}
}

Expand Down Expand Up @@ -385,7 +377,7 @@ void haplo_em_pin(
}
Free(geno);
geno = NULL;

}

/***********************************************************************************/
Expand Down Expand Up @@ -1140,7 +1132,7 @@ static int ranAS183_seed(int iseed1, int iseed2, int iseed3)

/***********************************************************************************/

static double ranAS183()
static double ranAS183(void)
{
double u;

Expand All @@ -1156,8 +1148,12 @@ static void errmsg(char *string){

/* Function to emulate "stop" of S+ - see page 134, S Programing, by
Venables and Ripley */
REprintf("%s");
// PROBLEM "%s", string RECOVER(NULL_ENTRY);

/* PROBLEM "%s", string RECOVER(NULL_ENTRY);
Replace with call to Rf_error for R 4.1.x
*/
//Rf_error(string, "%s");
REprintf("%s", string);
}


Expand Down Expand Up @@ -1199,7 +1195,7 @@ static void add_more_memory(HAP ***hap_list, double **prior,int *max_haps){


/* check that max_haps will not exceed max limit for an int on a 32-bit processor */


if(*max_haps == INT_MAX)
{
Expand All @@ -1216,16 +1212,18 @@ static void add_more_memory(HAP ***hap_list, double **prior,int *max_haps){
}


*prior = (double *) Realloc(*prior, *max_haps, double);
if(prior==NULL){
double *new_prior = (double *) Realloc(*prior, *max_haps, double);
if(new_prior==NULL){
errmsg("could not realloc mem for prior");
}

*hap_list = (HAP **) Realloc(*hap_list, *max_haps, HAP* );
if(hap_list==NULL){
*prior = new_prior;

HAP **new_hap_list = (HAP **) Realloc(*hap_list, *max_haps, HAP* );
if(new_hap_list==NULL){
errmsg("could not realloc mem for hap_list");
}

*hap_list = new_hap_list;

return;
}

Expand Down
6 changes: 2 additions & 4 deletions src/haplo-stats/haplo_em_pin.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
*
*Daniel J. Schaid, Ph.D.
*Division of Biostatistics
*Harwick Building û Room 775
*Harwick Building û Room 775
*Mayo Clinic
*200 First St., SW
*Rochester, MN 55905
Expand Down Expand Up @@ -109,12 +109,10 @@ static void set_posterior(int n_hap, HAP **hap_list, int *random_start);

static int ranAS183_seed(int iseed1, int iseed2, int iseed3);

static double ranAS183();
static double ranAS183(void);

static void errmsg(char *string);

static void compact(HAP **hap_list, int n, int *nReturn);

static HAPUNIQUE* copy_hap_unique(HAP *old, double *prior);

static void unique_haps(int n_hap, HAP **hap_list, HAPUNIQUE **u_hap_list, double *prior);
Expand Down

0 comments on commit 09aba6b

Please sign in to comment.