From f3559c3f82fba7f02038478942fda82e1bca8b25 Mon Sep 17 00:00:00 2001 From: Kevin Galinsky Date: Wed, 7 Sep 2016 11:00:02 -0400 Subject: [PATCH] merged nick's changes --- README | 9 +- include/egsubs.h | 22 +- include/strsubs.h | 213 ++---- include/vsubs.h | 507 +++++-------- src/egsubs.c | 182 ++--- src/nicksrc/strsubs.c | 1515 ++++++++++++++++++++----------------- src/nicksrc/vsubs.c | 1642 ++++++++++++++++++++++++----------------- 7 files changed, 2144 insertions(+), 1946 deletions(-) diff --git a/README b/README index 55ac1a9..2fed8bb 100644 --- a/README +++ b/README @@ -1,4 +1,4 @@ -EIGENSOFT version 6.1.4, 9/2/16 (for Linux only) +EIGENSOFT version 6.1.4, 9/7/16 (for Linux only) The EIGENSOFT package implements methods from the following 3 papers: Patterson et al. 2006 PLoS Genet (population structure) @@ -7,8 +7,11 @@ Galinsky et al. 2016 Am J Hum Genet (FastPCA and PC-based selection statistic) NEW features of EIGENSOFT version 6.1.4 include: -- pcaselection was omitted from 6.1.3 by accident --- statically linked GSL/openblas --- fixed memory allocation bug in pcaselection +-- Statically linked GSL/openblas +-- Fixed memory allocation bug in pcaselection +-- Some routines moved into nicklib +-- Error message on allocate failure now prints length as "%ld" + supporting long values. NEW features of EIGENSOFT version 6.1.3 include: -- Restored script file extensions to .perl instead of .pl diff --git a/include/egsubs.h b/include/egsubs.h index 8c58cdc..cf39c4a 100644 --- a/include/egsubs.h +++ b/include/egsubs.h @@ -1,14 +1,10 @@ -#include "admutils.h" +#include "admutils.h" -int -makeeglist (char **eglist, int maxnumeg, Indiv **indivmarkers, int numindivs); -int -mkeglist (Indiv **indm, int numindivs, char **eglist); -void -seteglist (Indiv **indm, int nindiv, char *eglistname); -void -seteglistv (Indiv **indm, int nindiv, char *eglistname, int val); -int -loadlist (char **list, char *listname); -int -loadlist_type (char **list, char *listname, int *ztypes, int off); + +int makeeglist (char **eglist, int maxnumeg, Indiv ** indivmarkers, + int numindivs); +int mkeglist (Indiv ** indm, int numindivs, char **eglist); +void seteglist (Indiv ** indm, int nindiv, char *eglistname); +void seteglistv (Indiv ** indm, int nindiv, char *eglistname, int val); +int loadlist (char **list, char *listname); +int loadlist_type (char **list, char *listname, int *ztypes, int off); diff --git a/include/strsubs.h b/include/strsubs.h index 697cda2..5c9f5a3 100644 --- a/include/strsubs.h +++ b/include/strsubs.h @@ -1,142 +1,87 @@ #include -int -splitup (char *strin, char *strpt[], int maxpt); -int -splitupx (char *strin, char **spt, int maxpt, char splitc); -int -splitupwxbuff (char *strin, char **spt, int maxpt, char *bigbuff, - int bigbufflen); -int -splitupxbuff (char *strin, char **spt, int maxpt, char splitc, char *bigbuff, - int bigbufflen); -int -oldsplitup (char *strin, char *strpt[], int maxpt); -void -freeup (char *strpt[], int numpt); -int -split1 (char *strin, char *strpt[], char splitc); -int -first_word (char *string, char *word, char *rest); -char * -fnwhite (char *ss); -char * -fwhite (char *ss); -char * -ftab (char *ss); -int -NPisnumber (char c); -int -isnumword (char *str); -void -fatalx (char *fmt, ...); -long -seednum (); -void -printbl (int n); -void -printnl (); -void -striptrail (char *sss, char c); -void -catx (char *sout, char **spt, int n); -void -catxx (char *sout, char **spt, int n); -void -catxc (char *sout, char **spt, int n, char c); -void -makedfn (char *dirname, char *fname, char *outname, int maxstr); -int -substring (char **ap, char *inx, char *outx); -int -numcols (char *name); -int -numlines (char *name); -void -openit (char *name, FILE **fff, char *type); -int -ftest (char *aname); -int -getxx (double **xx, int maxrow, int numcol, char *fname); -int -getss (char **ss, char *fname); -double -clocktime (); // cpu time in seconds -void -crevcomp (char *sout, char *sin); -int -indxstring (char **namelist, int len, char *strid); -int -indxstringr (char **namelist, int len, char *strid); -char * -strstrx (char *s1, char *s2); // case insensitive strstr -int -getxxnames (char ***pnames, double **xx, int maxrow, int numcol, char *fname); -int -getjjnames (char ***pnames, int **xx, int maxrow, int numcol, char *fname); -int -getxxnamesf (char ***pnames, double **xx, int maxrow, int numcol, FILE *fff); -int -getnameslohi (char ****pnames, int maxrow, int numcol, char *fname, int lo, - int hi); -int -getnames (char ****pnames, int maxrow, int numcol, char *fname); -char -num2iub (int num); -char -revchar (char c); -int -iub2num (char c); -char -num2base (int num); -int -base2num (char c); -char * -int_string (int a, int len, int base); -char * -binary_string (int a, int len); -int -string_binary (char *sx); -void -freestring (char **ss); -void -copystrings (char **sa, char **sb, int n); -void -printstringsw (char **ss, int n, int slen, int width); -void -printstrings (char **ss, int n); -int -ridfile (char *fname); -char -compbase (char x); -void -mkupper (char *sx); -void -mklower (char *sx); -int -iubdekode (char *a, char iub); -int -isiub (char iub); -int -isiub2 (char iub); -int -iubcbases (char *cbases, char iub); -int -ishet (char c); -int -char2int (char cc); -char -int2char (int x); -void -chomp (char *str); +int splitup (char *strin, char *strpt[],int maxpt) ; +int splitupx(char *strin, char **spt, int maxpt, char splitc) ; +int splitupwxbuff(char *strin, char **spt, int maxpt, char *bigbuff, int bigbufflen) ; +int splitupxbuff(char *strin, char **spt, int maxpt, char splitc, char *bigbuff, int bigbufflen) ; +int oldsplitup (char *strin, char *strpt[],int maxpt) ; +void freeup (char *strpt[],int numpt) ; +int split1 (char *strin, char *strpt[], char splitc); +int first_word(char *string, char *word, char *rest) ; +char *fnwhite (char *ss) ; +char *fwhite (char *ss) ; +char *ftab (char *ss) ; +int NPisnumber (char c) ; +int isnumword (char *str) ; +void fatalx( char *fmt, ...) ; +long seednum() ; +void printbl(int n) ; +void printnl() ; +void striptrail(char *sss, char c) ; +void catx(char *sout, char **spt, int n) ; +void catxx(char *sout, char **spt, int n) ; +void catxc(char *sout, char **spt, int n, char c) ; +void makedfn(char *dirname, char *fname, char *outname, int maxstr) ; +int substring (char **ap, char *inx, char *outx) ; +int numcols (char *name) ; +int numlines(char *name) ; +void openit(char *name, FILE **fff, char *type) ; +int ftest(char *aname) ; +int getxx(double **xx, int maxrow, int numcol, char *fname) ; +int getss(char **ss, char *fname) ; +int loadlist(char **list, char *listname) ; // with dup check +void printdups(char **list, int n) ; +int checkdup(char **list, int n) ; +double clocktime() ; // cpu time in seconds +void crevcomp(char *sout, char *sin) ; +int indxstring(char **namelist, int len, char *strid) ; +int indxstringr(char **namelist, int len, char *strid) ; +char *strstrx(char *s1, char *s2) ; // case insensitive strstr +int getxxnames(char ***pnames, double **xx, int maxrow, int numcol, char *fname); +int getjjnames(char ***pnames, int **xx, int maxrow, int numcol, char *fname); +int getxxnamesf(char ***pnames, double **xx, int maxrow, int numcol, FILE *fff) ; +int getnameslohi(char ****pnames, int maxrow, int numcol, char *fname, int lo, int hi) ; +int getnamesstripcolon(char ****pnames, int maxrow, int numcol, char *fname, int lo, int hi) ; +int getnames(char ****pnames, int maxrow, int numcol, char *fname) ; +char num2iub (int num) ; +char revchar(char c) ; +int iub2num(char c) ; +char num2base (int num) ; +int base2num(char c) ; +char *int_string(int a, int len, int base) ; +char *binary_string(int a, int len) ; +int string_binary(char *sx) ; +void freestring (char **ss) ; +void copystrings(char **sa, char **sb, int n) ; +void printstringsw(char **ss, int n, int slen, int width) ; +void printstrings(char **ss, int n) ; +int ridfile(char *fname) ; +char compbase(char x) ; +void mkupper(char *sx) ; +void mklower(char *sx) ; +int iubdekode(char *a, char iub) ; +int isiub(char iub) ; +int isiub2(char iub) ; +int iubcbases(char *cbases, char iub) ; +int ishet(char c) ; +int cttype(char cc) ; +int char2int(char cc) ; +char int2char(int x) ; +void chomp(char *str) ; + +int numcmatch(char *cc, int len, char c) ; +int numcnomatch(char *cc, int len, char c) ; +char *strnotchar(char *s, char c) ; +char *findupper(char *s) ; +char *fgetstrap(char *buff, int maxlen, FILE *fff, int *ret) ; +char readtonl(FILE *fff) ; +int filehash(char *name) ; + + -int -numcmatch (char *cc, int len, char c); -int -numcnomatch (char *cc, int len, char c); #define ZALLOC(item,n,type) if ((item = (type *)calloc((n),sizeof(type))) == NULL) \ - fatalx("Unable to allocate %d unit(s) for item \n",n) + fatalx("Unable to allocate %ld unit(s) for item \n", (long) n) #undef MAX #undef MIN diff --git a/include/vsubs.h b/include/vsubs.h index 993aa0f..29c4ec2 100644 --- a/include/vsubs.h +++ b/include/vsubs.h @@ -4,327 +4,194 @@ #include #include "strsubs.h" -void -vsp (double *a, double *b, double c, int n); -void -vst (double *a, double *b, double c, int n); -void -vvt (double *a, double *b, double *c, int n); -void -vvp (double *a, double *b, double *c, int n); -void -vvm (double *a, double *b, double *c, int n); -void -vvd (double *a, double *b, double *c, int n); -void -vsqrt (double *a, double *b, int n); -void -vinvert (double *a, double *b, int n); -void -vabs (double *a, double *b, int n); -void -vlog (double *a, double *b, int n); -void -vlog2 (double *a, double *b, int n); -void -vexp (double *a, double *b, int n); -void -vclear (double *a, double c, int n); -void -vzero (double *a, int n); -void -cpzero (char **a, int n); -void -ivvp (int *a, int *b, int *c, int n); -void -ivvm (int *a, int *b, int *c, int n); -void -ivsp (int *a, int *b, int c, int n); -void -ivst (int *a, int *b, int c, int n); -void -ivclear (int *a, int c, long n); -void -lvclear (long *a, long c, long n); -void -ivzero (int *a, int n); -void -cclear (unsigned char *a, unsigned char c, long n); - -double -clip (double x, double lo, double hi); -void -ivclip (int *a, int *b, int loval, int hival, int n); -void -vclip (double *a, double *b, double loval, double hival, int n); - -void -vmaxmin (double *a, int n, double *max, double *min); -void -vlmaxmin (double *a, int n, int *max, int *min); -void -ivmaxmin (int *a, int n, int *max, int *min); -int -minivec (int *a, int n); -int -maxivec (int *a, int n); -void -ivlmaxmin (int *a, int n, int *max, int *min); -void -getdiag (double *a, double *b, int n); -void -setdiag (double *a, double *diag, int n); -void -flipiarr (int *a, int *b, int n); -void -fliparr (double *a, double *b, int n); -int -ipow2 (int l); - -void -copyarr (double *a, double *b, int n); -void -revarr (double *a, double *b, int n); -void -reviarr (int *a, int *b, int n); -void -revuiarr (unsigned int *a, unsigned int *b, int n); -void -copyiarr (int *a, int *b, int n); -void -copyiparr (int **a, int **b, int n); - -void -dpermute (double *a, int *ind, int len); -void -ipermute (int *a, int *ind, int len); -void -dppermute (double **a, int *ind, int len); -void -ippermute (int **a, int *ind, int len); - -double -asum (double *a, int n); -double -asum2 (double *a, int n); -int -intsum (int *a, int n); -long -longsum (long *a, int n); -int -idot (int *a, int *b, int n); -int -iprod (int *a, int n); -double -aprod (double *a, int n); -double -vdot (double *a, double *b, int n); -double -corr (double *a, double *b, int n); -double -corrx (double *a, double *b, int n); -double -variance (double *a, int n); -double -trace (double *a, int n); -int -nnint (double a); -void -countcat (int *tags, int n, int *ncat, int nclass); -void -rowsum (double *a, double *rr, int n); -void -colsum (double *a, double *cc, int n); -void -rrsum (double *a, double *cc, int m, int n); -void -ccsum (double *a, double *cc, int m, int n); -void -printmatfile (double *a, int m, int n, FILE *fff); -void -printmatwfile (double *a, int m, int n, int w, FILE *fff); -void -printmat (double *a, int m, int n); -void -printmatw (double *a, int m, int n, int w); -void -printmatl (double *a, int m, int n); -void -printmatwl (double *a, int m, int n, int w); -void -printmatwf (double *a, int m, int n, int w, char *format); -void -int2c (char *cc, int *b, int n); -void -floatit (double *a, int *b, int n); -void -fixit (int *a, double *b, int n); -void -rndit (double *a, double *b, int n); -void -printimatw (int *a, int m, int n, int w); -void -printimatx (int *a, int m, int n); -void -printimat (int *a, int m, int n); -void -printimatl (int *a, int m, int n); -void -printimatlfile (int *a, int m, int n, FILE *fff); -void -printimatfile (int *a, int m, int n, FILE *fff); -void -printimatwfile (int *a, int m, int n, int w, FILE *fff); -void -printstring (char *ss, int width); -void -printstringbasepos (char *ss, int w, int basepos); -void -printstringf (char *ss, int width, FILE *fff); - -int -findfirst (int *a, int n, int val); -int -findfirstl (long *a, int n, long val); -int -findfirstu (unsigned int *a, int n, unsigned int val); -int -findlastu (unsigned int *a, int n, unsigned int val); - -int -findlast (int *a, int n, int val); -int -binsearch (int *a, int n, int val); -void -idperm (int *a, int n); -double -NPlog2 (double y); -double -log2fac (int n); -double -logfac (int n); -double -logbino (int n, int k); -double -loghprob (int n, int a, int m, int k); +void vsp(double *a, double *b, double c, int n); +void vst(double *a, double *b, double c, int n); +void vvt(double *a, double *b, double *c, int n); +void vvp(double *a, double *b, double *c, int n); +void vvm(double *a, double *b, double *c, int n); +void vvd(double *a, double *b, double *c, int n); +void vsqrt(double *a, double *b, int n) ; +void vinvert(double *a, double *b, int n) ; +void vabs(double *a, double *b, int n) ; +void vlog(double *a, double *b, int n) ; +void vlog2(double *a, double *b, int n) ; +void vexp(double *a, double *b, int n) ; +void vclear(double *a, double c, int n) ; +void vzero(double *a, int n) ; +void cpzero(char **a, int n) ; +void ivvp(int *a, int *b, int *c, int n); +void ivvm(int *a, int *b, int *c, int n); +void ivsp(int *a, int *b, int c, int n); +void ivst(int *a, int *b, int c, int n); +void ivclear(int *a, int c, long n) ; +void lvclear(long *a, long c, long n) ; +void ivzero(int *a, int n) ; +void lvzero(long *a, long n) ; +void cclear(unsigned char *a, unsigned char c, long n) ; +void charclear(char *a, unsigned char c, long n) ; + +double clip(double x, double lo, double hi) ; +void ivclip(int *a, int *b,int loval, int hival,int n) ; +void vclip(double *a, double *b,double loval, double hival,int n) ; + +void vmaxmin(double *a, int n, double *max, double *min) ; +void vlmaxmin(double *a, int n, int *max, int *min) ; +void ivmaxmin(int *a, int n, int *max, int *min) ; +int minivec(int *a, int n) ; +int maxivec(int *a, int n) ; +void ivlmaxmin(int *a, int n, int *max, int *min) ; +void getdiag(double *a, double *b, int n) ; +void setdiag(double *a, double *diag, int n) ; +void adddiag(double *a, double *diag, int n) ; +void flipiarr(int *a, int *b, int n) ; +void fliparr(double *a, double *b, int n) ; +int ipow2 (int l) ; + +void copyarr(double *a,double *b,int n) ; +void revarr(double *a, double *b,int n) ; +void reviarr(int *a,int *b,int n) ; +void revuiarr(unsigned int *a, unsigned int *b,int n) ; +void copyiarr(int *a,int *b,int n) ; +void copylarr(long *a, long *b, int n) ; +void copyiparr(int **a,int **b,int n) ; + +void dpermute(double *a, int *ind, int len) ; +void ipermute(int *a, int *ind, int len) ; +void dppermute(double **a, int *ind, int len) ; +void ippermute(int **a, int *ind, int len) ; + +double asum(double *a, int n) ; +double asum2(double *a, int n) ; +int intsum(int *a, int n) ; +long longsum(long *a, int n) ; +int idot(int *a, int *b, int n) ; +int iprod(int *a, int n) ; +double aprod(double *a, int n) ; +double vdot(double *a, double *b, int n) ; +double corr(double *a, double *b, int n) ; +double corrx(double *a, double *b, int n) ; +double variance(double *a, int n) ; +double trace(double *a, int n) ; +int nnint(double a) ; +void countcat(int *tags, int n,int *ncat,int nclass) ; +void rowsum(double *a, double *rr, int n) ; +void colsum(double *a, double *cc, int n) ; +void rrsum(double *a, double *cc, int m, int n) ; +void ccsum(double *a, double *cc, int m, int n) ; +void printmatfile(double *a, int m, int n, FILE *fff) ; +void printmatwfile(double *a, int m, int n, int w, FILE *fff) ; +void printmatx(double *a, int m, int n) ; +void printmat(double *a, int m, int n) ; +void printmatwx(double *a, int m, int n, int w) ; +void printmatw(double *a, int m, int n, int w) ; +void printmatl(double *a, int m, int n) ; +void printmatwl(double *a, int m, int n, int w) ; +void printmatwf(double *a, int m, int n, int w, char *format); +void int2c(char *cc, int *b, int n) ; +void floatit(double *a, int *b, int n) ; +void fixit(int *a, double *b, int n) ; +void rndit(double *a, double *b, int n) ; +void printimatw(int *a, int m, int n, int w) ; +void printimatx(int *a, int m, int n) ; +void printimat1(int *a, int m, int n) ; +void printimat1x(int *a, int m, int n) ; +void printimat(int *a, int m, int n) ; +void printimatl(int *a, int m, int n) ; +void printimatlfile(int *a, int m, int n, FILE *fff) ; +void printimatfile(int *a, int m, int n, FILE *fff) ; +void printimatwfile(int *a, int m, int n, int w, FILE *fff) ; +void printimat2D(int **a, int m, int n) ; +void printmat2D(double **a, int m, int n) ; +void printstring(char *ss, int width) ; +void printstringbasepos(char *ss, int w, int basepos) ; +void printstringf(char *ss, int width, FILE *fff) ; + +int findfirst(int *a, int n, int val) ; +int findfirstl(long *a, int n, long val) ; +int findfirstu(unsigned int *a, int n, unsigned int val) ; +int findlastu(unsigned int *a, int n, unsigned int val) ; + +int findlast(int *a, int n, int val) ; +int binsearch(int *a, int n, int val) ; +void idperm(int *a, int n) ; +double NPlog2(double y) ; +double log2fac(int n) ; +double logfac(int n) ; +double logbino(int n, int k) ; +double loghprob(int n, int a, int m, int k) ; /* hypergeometric probability */ -double -logmultinom (int *cc, int n); -double -addlog (double a, double b); -double -vldot (double *x, double *y, int n); -double -pow10 (double x); -double -vpow10 (double *a, double *b, int n); -double -vlog10 (double *a, double *b, int n); +double logmultinom(int *cc, int n) ; +double addlog(double a, double b) ; +double vldot(double *x, double *y, int n) ; +double pow10 (double x) ; +void vpow10 (double *a, double *b, int n) ; +void vlog10 (double *a, double *b, int n) ; /* matrix transpose */ -void -transpose (double *aout, double *ain, int m, int n); -void -addoutmul (double *out, double *a, double mul, int n); -void -addouter (double *out, double *a, int n); -void -subouter (double *out, double *a, int n); +void transpose(double *aout, double *ain, int m, int n) ; +void addoutmul(double *out, double *a, double mul, int n) ; +void addouter(double *out, double *a, int n) ; +void subouter(double *out, double *a, int n) ; +int mktriang(double *out, double *in, int n) ; +int mkfull(double *out, double *in, int n) ; /* storage allocation */ -int ** -initarray_2Dint (int numrows, int numcolumns, int initval); -long ** -initarray_2Dlong (int numrows, int numcolumns, int initval); -void -free2Dint (int ***xx, int numrows); -double ** -initarray_2Ddouble (int numrows, int numcolumns, double initval); -long double ** -initarray_2Dlongdouble (int numrows, int numcolumns, long double initval); -void -clear2D (double ***xx, int numrows, int numcols, double val); -void -iclear2D (int ***xx, int numrows, int numcols, int val); -void -free2D (double ***xx, int numrows); -void -free2Dlongdouble (long double ***xx, int numrows); -void -free_darray (double **xx); -void -free_iarray (int **xx); - -double -bal1 (double *a, int n); -void -vcompl (double *a, double *b, int n); -void -setidmat (double *a, int n); - -int -stripit (double *a, double *b, int *x, int len); -int -istripit (int *a, int *b, int *x, int len); -int -cstripit (char **a, char **b, int *x, int len); - -void -mapit (int *a, int *b, int n, int inval, int outval); -int -ifall (int n, int k); // falling factorial = n (n-1) (n-2) ... (n-k+1) -double -hlife (double val); -void * -topheap (); - -void -swap (double *pa, double *pb); -void -iswap (int *pa, int *pb); -void -cswap (char *c1, char *c2); - -void -copyarr2D (double **a, double **b, int nrows, int ncols); // a input b output -void -copyiarr2D (int **a, int **b, int nrows, int ncols); // a input b output -void -plus2Dint (int **a, int **b, int **c, int nrows, int ncols); -void -minus2Dint (int **a, int **b, int **c, int nrows, int ncols); - -void -plus2D (double **a, double **b, double **c, int nrows, int ncols); -void -minus2D (double **a, double **b, double **c, int nrows, int ncols); -void -sum2D (double *a, double **b, int nrows, int ncols); -int -total2D (double **a, int nrows, int ncols); -int -total2Dint (int **a, int nrows, int ncols); - -int -kodeitb (int *xx, int len, int base); -int -dekodeitb (int *xx, int kode, int len, int base); -int -kodeitbb (int *xx, int len, int *baselist); -int -dekodeitbb (int *xx, int kode, int len, int *baselist); - -int -isprime (long num); -long -nextprime (long num); -int -irevcomp (int xx, int stringlen); -long -lrevcomp (long xx, int stringlen); -void -ismatch (int *a, int *b, int n, int val); -int -pmult (double *a, double *b, double *c, int na, int nb); -void -pdiff (double *a, double *b, int deg); - +int **initarray_2Dint(int numrows, int numcolumns, int initval); +long **initarray_2Dlong(int numrows, int numcolumns, int initval); +void free2Dint(int ***xx, int numrows) ; +void free2Dlong(long ***xx, int numrows) ; +double **initarray_2Ddouble(int numrows, int numcolumns, double initval); +long double **initarray_2Dlongdouble(int numrows, int numcolumns, long double initval); +void clear2D(double ***xx, int numrows, int numcols, double val) ; +void iclear2D(int ***xx, int numrows, int numcols, int val) ; +void lclear2D(long ***xx, int numrows, int numcols, long val) ; +void free2D(double ***xx, int numrows) ; +void free2Dlongdouble(long double ***xx, int numrows) ; +void free_darray (double **xx) ; +void free_iarray (int **xx) ; + +double bal1 (double *a, int n) ; +void vcompl(double *a, double *b, int n) ; +void setidmat(double *a, int n) ; + +int stripit(double *a, double *b, int *x, int len) ; +int istripit(int *a, int *b, int *x, int len) ; +int cstripit(char **a, char **b, int *x, int len) ; + +void mapit(int *a, int *b, int n, int inval, int outval) ; +int ifall(int n, int k) ; // falling factorial = n (n-1) (n-2) ... (n-k+1) +double hlife(double val) ; +void *topheap () ; + +void swap (double *pa, double *pb) ; +void iswap (int *pa, int *pb) ; +void cswap(char *c1, char *c2) ; + + +void floatit2D(double **a, int **b, int nrows, int ncols) ; +void copyarr2D(double **a, double **b, int nrows, int ncols) ; // a input b output +void copyiarr2D(int **a, int **b, int nrows, int ncols) ; // a input b output +void plus2Dint(int **a, int **b, int **c, int nrows, int ncols) ; +void minus2Dint(int **a, int **b, int **c, int nrows, int ncols) ; + +void plus2D(double **a, double **b, double **c, int nrows, int ncols) ; +void minus2D(double **a, double **b, double **c, int nrows, int ncols) ; +void sum2D(double *a, double **b, int nrows, int ncols) ; +double total2D(double **a, int nrows, int ncols) ; +int total2Dint(int **a, int nrows, int ncols) ; + +int kodeitb(int *xx, int len, int base) ; +int dekodeitb(int *xx, int kode, int len, int base) ; +long lkodeitbb(int *xx, int len, int *baselist) ; +int ldekodeitbb(int *xx, long kode, int len, int *baselist) ; +int kodeitbb(int *xx, int len, int *baselist) ; +int dekodeitbb(int *xx, int kode, int len, int *baselist) ; + +int isprime(long num) ; +long nextprime(long num) ; +int irevcomp (int xx, int stringlen) ; +long lrevcomp (long xx, int stringlen) ; +void ismatch(int *a, int *b, int n, int val) ; +int pmult(double *a, double *b, double *c, int na, int nb) ; +void pdiff(double *a, double *b, int deg) ; +void vswap(double *a, double *b, int n) ; +void setlong(long *pplen, long a, long b) ; diff --git a/src/egsubs.c b/src/egsubs.c index 32a800e..a283f8c 100644 --- a/src/egsubs.c +++ b/src/egsubs.c @@ -1,84 +1,51 @@ -#include "mcio.h" -#include "egsubs.h" +#include "mcio.h" +#include "egsubs.h" + int -makeeglist (char **eglist, int maxnumeg, Indiv **indivmarkers, int numindivs) +makeeglist (char **eglist, int maxnumeg, Indiv ** indivmarkers, int numindivs) // old routine mkeglist { Indiv *indx; int i, k, numeg = 0; - for (i = 0; i < numindivs; i++) - { - indx = indivmarkers[i]; - if (indx->ignore) - continue; - k = indxindex (eglist, numeg, indx->egroup); - if (k < 0) - { - if (numeg >= maxnumeg) - { - printf ( - "number of populations too large. Increase maxpops if you wish\n"); - fatalx ( - "(makeeglist) You really want to analyse more than %d populations?\n", - maxnumeg); - } - eglist[numeg] = strdup (indx->egroup); - ++numeg; - } + for (i = 0; i < numindivs; i++) { + indx = indivmarkers[i]; + if (indx->ignore) + continue; + k = indxindex (eglist, numeg, indx->egroup); + if (k < 0) { + if (numeg >= maxnumeg) { + printf + ("number of populations too large. Increase maxpops if you wish\n"); + fatalx + ("(makeeglist) You really want to analyse more than %d populations?\n", + maxnumeg); + } + eglist[numeg] = strdup (indx->egroup); + ++numeg; } + } return numeg; } + int -mkeglist (Indiv **indm, int numindivs, char **eglist) +mkeglist (Indiv ** indm, int numindivs, char **eglist) { Indiv *indx; int i, k, numeg = 0; - for (i = 0; i < numindivs; i++) - { - indx = indm[i]; - if (indx->ignore) - continue; - k = indxindex (eglist, numeg, indx->egroup); - if (k < 0) - { - eglist[numeg] = strdup (indx->egroup); - ++numeg; - } + for (i = 0; i < numindivs; i++) { + indx = indm[i]; + if (indx->ignore) + continue; + k = indxindex (eglist, numeg, indx->egroup); + if (k < 0) { + eglist[numeg] = strdup (indx->egroup); + ++numeg; } + } return numeg; } -int -loadlist (char **list, char *listname) -// listname is just a list of names ... -{ - FILE *lfile; - char line[MAXSTR]; - char *spt[MAXFF]; - char *sx; - int nsplit, i, n = 0; - - if (listname == NULL) - return 0; - openit (listname, &lfile, "r"); - while (fgets (line, MAXSTR, lfile) != NULL) - { - nsplit = splitup (line, spt, MAXFF); - if (nsplit == 0) - continue; - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - list[n] = strdup (sx); - ++n; - freeup (spt, nsplit); - } - return n; -} int loadlist_type (char **list, char *listname, int *ztypes, int off) @@ -94,30 +61,29 @@ loadlist_type (char **list, char *listname, int *ztypes, int off) if (listname == NULL) return 0; openit (listname, &lfile, "r"); - while (fgets (line, MAXSTR, lfile) != NULL) - { - nsplit = splitup (line, spt, MAXFF); - if (nsplit == 0) - continue; - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (nsplit < 2) - fatalx ("bad listname: %s\n", sx); - list[n] = strdup (sx); - tt = atoi (spt[1]); - ztypes[n] = tt + off; - ++n; + while (fgets (line, MAXSTR, lfile) != NULL) { + nsplit = splitup (line, spt, MAXFF); + if (nsplit == 0) + continue; + sx = spt[0]; + if (sx[0] == '#') { freeup (spt, nsplit); + continue; } + if (nsplit < 2) + fatalx ("bad listname: %s\n", sx); + list[n] = strdup (sx); + tt = atoi (spt[1]); + ztypes[n] = tt + off; + ++n; + freeup (spt, nsplit); + } return n; } + void -seteglist (Indiv **indm, int nindiv, char *eglistname) +seteglist (Indiv ** indm, int nindiv, char *eglistname) { FILE *egfile; char line[MAXSTR]; @@ -129,22 +95,21 @@ seteglist (Indiv **indm, int nindiv, char *eglistname) if (eglistname == NULL) return; openit (eglistname, &egfile, "r"); - while (fgets (line, MAXSTR, egfile) != NULL) - { - nsplit = splitup (line, spt, MAXFF); - if (nsplit == 0) - continue; - sx = spt[0]; - if (sx[0] == '#') - continue; - setstatus (indm, nindiv, sx); - freeup (spt, nsplit); - } + while (fgets (line, MAXSTR, egfile) != NULL) { + nsplit = splitup (line, spt, MAXFF); + if (nsplit == 0) + continue; + sx = spt[0]; + if (sx[0] == '#') + continue; + setstatus (indm, nindiv, sx); + freeup (spt, nsplit); + } fclose (egfile); } void -seteglistv (Indiv **indm, int nindiv, char *eglistname, int val) +seteglistv (Indiv ** indm, int nindiv, char *eglistname, int val) { FILE *egfile; char line[MAXSTR]; @@ -153,23 +118,20 @@ seteglistv (Indiv **indm, int nindiv, char *eglistname, int val) Indiv *indx; int nsplit, i; - if (eglistname == NULL) - { - setstatusv (indm, nindiv, NULL, val); - } + if (eglistname == NULL) { + setstatusv (indm, nindiv, NULL, val); + } openit (eglistname, &egfile, "r"); - while (fgets (line, MAXSTR, egfile) != NULL) - { - nsplit = splitup (line, spt, MAXFF); - if (nsplit == 0) - continue; - sx = spt[0]; - if (sx[0] == '#') - continue; - setstatusv (indm, nindiv, sx, val); - freeup (spt, nsplit); - } + while (fgets (line, MAXSTR, egfile) != NULL) { + nsplit = splitup (line, spt, MAXFF); + if (nsplit == 0) + continue; + sx = spt[0]; + if (sx[0] == '#') + continue; + setstatusv (indm, nindiv, sx, val); + freeup (spt, nsplit); + } fclose (egfile); } - diff --git a/src/nicksrc/strsubs.c b/src/nicksrc/strsubs.c index 58ab14f..b96fac7 100644 --- a/src/nicksrc/strsubs.c +++ b/src/nicksrc/strsubs.c @@ -8,20 +8,25 @@ #include #include #include +#include + #define MAXSTR 10000 #define MAXFF 50 -#include "strsubs.h" -#include "vsubs.h" +#include "strsubs.h" +#include "vsubs.h" +#include "getpars.h" extern int errno; + int -oldsplitup (char *strin, char**spt, int maxpt) +oldsplitup (char *strin, char **spt, int maxpt) + /** retained in case there are compatibility problems - */ +*/ { char *s1, *s2, *sx; char *str; @@ -30,43 +35,40 @@ oldsplitup (char *strin, char**spt, int maxpt) len = strlen (strin); if (len == 0) return 0; - ZALLOC(str, 2*len, char); + ZALLOC (str, 2 * len, char); num = 0; sx = strin; - for (i = 0; i < maxpt; i++) - { - s1 = fnwhite (sx); - if (s1 == NULL) - { - break; - } - s2 = fwhite (s1); - if (s2 == NULL) - { - s2 = s1 + strlen (s1); - } - s2--; /* now points at last character of next word */ - len = s2 - s1 + 1; - strncpy (str, s1, len); - str[len] = '\0'; - spt[num] = strdup (str); - ++num; - sx = s2 + 1; + for (i = 0; i < maxpt; i++) { + s1 = fnwhite (sx); + if (s1 == NULL) { + break; } + s2 = fwhite (s1); + if (s2 == NULL) { + s2 = s1 + strlen (s1); + } + s2--; /* now points at last character of next word */ + len = s2 - s1 + 1; + strncpy (str, s1, len); + str[len] = '\0'; + spt[num] = strdup (str); + ++num; + sx = s2 + 1; + } freestring (&str); return num; } void freeup (char *strpt[], int numpt) + /** free up array of strings */ { int i; - for (i = numpt - 1; i >= 0; i--) - { - if (strpt[i] != NULL) - freestring (&strpt[i]); - } + for (i = numpt - 1; i >= 0; i--) { + if (strpt[i] != NULL) + freestring (&strpt[i]); + } } int @@ -74,104 +76,99 @@ first_word (char *string, char *xword, char *xrest) /* first_word(string, *word, *rest) - Break the string into the first word and the rest. Both word and - rest begin with non-white space, unless rest is null. - Return: - 0 means string is all white - 1 means word is non-white, but rest is white - 2 means word and rest are non-white + Break the string into the first word and the rest. Both word and +rest begin with non-white space, unless rest is null. + Return: + 0 means string is all white + 1 means word is non-white, but rest is white + 2 means word and rest are non-white - If string and rest coincide, string will be overwritten + If string and rest coincide, string will be overwritten - */ +*/ { char *spt, x; char *ss = NULL, *sx; int l1, l2; ss = strdup (string); - if (ss == NULL) - { - printf ("strdup fails\n"); - printf ("%s\n", string); - fatalx ("first_word... strdup fails\n"); - } + if (ss == NULL) { + printf ("strdup fails\n"); + printf ("%s\n", string); + fatalx ("first_word... strdup fails\n"); + } fflush (stdout); spt = ss; xword[0] = xrest[0] = '\0'; - if ((spt = fnwhite (ss)) == NULL) - { - free (ss); - return 0; - } + if ((spt = fnwhite (ss)) == NULL) { + free (ss); + return 0; + } sx = fwhite (spt); - if (sx == NULL) - { - strcpy (xword, spt); - free (ss); - return 1; - } + if (sx == NULL) { + strcpy (xword, spt); + free (ss); + return 1; + } l1 = sx - spt; l2 = strlen (sx) - 1; *sx = '\0'; strcpy (xword, spt); - if (l2 <= 0) - { - free (ss); - return 1; - } + if (l2 <= 0) { + free (ss); + return 1; + } sx = fnwhite (sx + 1); - if (sx == NULL) - { - free (ss); - return 1; - } + if (sx == NULL) { + free (ss); + return 1; + } strcpy (xrest, sx); free (ss); return 2; } char * -fnwhite (char * ss) +fnwhite (char *ss) + /* return first non white space */ { char *x; if (ss == NULL) fatalx ("fnwhite: logic bug\n"); - for (x = ss; *x != '\0'; ++x) - { - if (!isspace(*x)) - return x; - } + for (x = ss; *x != '\0'; ++x) { + if (!isspace (*x)) + return x; + } return NULL; } char * ftab (char *ss) + /* return first tab */ { char *x; int n; - for (x = ss; *x != '\0'; ++x) - { - if (*x == CTAB) - return x; - } + for (x = ss; *x != '\0'; ++x) { + if (*x == CTAB) + return x; + } return NULL; } char * -fwhite (char * ss) +fwhite (char *ss) + /* return first white space */ { char *x; int n; - for (x = ss; *x != '\0'; ++x) - { - if (isspace(*x)) - return x; - } + for (x = ss; *x != '\0'; ++x) { + if (isspace (*x)) + return x; + } return NULL; } @@ -182,22 +179,24 @@ fatalx (char *fmt, ...) { va_list args; - va_start(args, fmt); + va_start (args, fmt); vsprintf (Estr, fmt, args); - va_end(args); + va_end (args); fflush (stdout); fprintf (stderr, "fatalx:\n%s", Estr); fflush (stderr); abort (); } + int NPisnumber (char c) + /** returns 1 if - + or digit - */ +*/ { - if (isdigit(c)) + if (isdigit (c)) return 1; if (c == '+') return 1; @@ -206,6 +205,7 @@ NPisnumber (char c) return 0; } + int isnumword (char *str) { @@ -215,21 +215,19 @@ isnumword (char *str) len = strlen (str); numpt = 0; - for (i = 0; i < len; i++) - { - c = str[i]; - - if ((c == '.') && (numpt == 0)) - { - ++numpt; - continue; - } - - if (!NPisnumber (c)) - return NO; - if (!isdigit(c) && (i > 0)) - return NO; + for (i = 0; i < len; i++) { + c = str[i]; + + if ((c == '.') && (numpt == 0)) { + ++numpt; + continue; } + + if (!NPisnumber (c)) + return NO; + if (!isdigit (c) && (i > 0)) + return NO; + } return YES; } @@ -246,9 +244,11 @@ seednum () c = d ^ ((a + b) << 15); + return c; } + int splitupwxbuff (char *strin, char **spt, int maxpt, char *bigbuff, int bigbufflen) @@ -263,31 +263,27 @@ splitupwxbuff (char *strin, char **spt, int maxpt, char *bigbuff, fatalx ("(splitupwxbuff) overflow\n%s", strin); strcpy (bigbuff, strin); num = 0; - for (k = 0; k < len; ++k) - { - if (!isspace(bigbuff[k])) - { - empty = NO; - klo = k; - sx = bigbuff + k; - break; - } + for (k = 0; k < len; ++k) { + if (!isspace (bigbuff[k])) { + empty = NO; + klo = k; + sx = bigbuff + k; + break; } + } if (empty) return 0; - for (k = klo; k < len; ++k) - { - if (isspace(strin[k])) - { - bigbuff[k] = CNULL; - if (num >= maxpt) - break; - spt[num] = sx; - if (strlen (sx) > 0) - ++num; - sx = bigbuff + k + 1; - } + for (k = klo; k < len; ++k) { + if (isspace (strin[k])) { + bigbuff[k] = CNULL; + if (num >= maxpt) + break; + spt[num] = sx; + if (strlen (sx) > 0) + ++num; + sx = bigbuff + k + 1; } + } if (num >= maxpt) return num; spt[num] = sx; @@ -295,6 +291,7 @@ splitupwxbuff (char *strin, char **spt, int maxpt, char *bigbuff, ++num; return num; } + int splitupxbuff (char *strin, char **spt, int maxpt, char splitc, char *bigbuff, int bigbufflen) @@ -308,30 +305,26 @@ splitupxbuff (char *strin, char **spt, int maxpt, char splitc, char *bigbuff, fatalx ("(splitupxbuff) overflow \n%s\n", strin); strcpy (bigbuff, strin); num = 0; - for (k = 0; k < len; ++k) - { - if (strin[k] != splitc) - { - empty = NO; - klo = k; - sx = bigbuff + k; - break; - } + for (k = 0; k < len; ++k) { + if (strin[k] != splitc) { + empty = NO; + klo = k; + sx = bigbuff + k; + break; } + } if (empty) return 0; - for (k = klo; k < len; ++k) - { - if (strin[k] == splitc) - { - bigbuff[k] = CNULL; - if (num >= maxpt) - fatalx ("overflow\n"); - spt[num] = sx; - sx = bigbuff + k + 1; - ++num; - } + for (k = klo; k < len; ++k) { + if (strin[k] == splitc) { + bigbuff[k] = CNULL; + if (num >= maxpt) + fatalx ("overflow\n"); + spt[num] = sx; + sx = bigbuff + k + 1; + ++num; } + } if (num >= maxpt) fatalx ("overflow\n"); spt[num] = sx; @@ -348,17 +341,17 @@ splitup (char *strin, char **spt, int maxpt) if (strin == NULL) return 0; len = strlen (strin); - ZALLOC(bigb, len+1, char); - ZALLOC(qpt, maxpt, char *); + ZALLOC (bigb, len + 1, char); + ZALLOC (qpt, maxpt + 10, char *); num = splitupwxbuff (strin, qpt, maxpt, bigb, len + 1); - for (k = 0; k < num; ++k) - { - spt[k] = strdup (qpt[k]); - } + for (k = 0; k < num; ++k) { + spt[k] = strdup (qpt[k]); + } free (bigb); free (qpt); return num; } + int splitupx (char *strin, char **spt, int maxpt, char splitc) { @@ -368,47 +361,45 @@ splitupx (char *strin, char **spt, int maxpt, char splitc) if (strin == NULL) return 0; len = strlen (strin); - ZALLOC(bigb, len+1, char); - ZALLOC(qpt, maxpt, char *); + ZALLOC (bigb, len + 1, char); + ZALLOC (qpt, maxpt + 10, char *); num = splitupxbuff (strin, qpt, maxpt, splitc, bigb, len + 1); - for (k = 0; k < num; ++k) - { - spt[k] = strdup (qpt[k]); - } - free (bigb); + for (k = 0; k < num; ++k) { + spt[k] = strdup (qpt[k]); + } free (qpt); + free (bigb); return num; } int split1 (char *strin, char *strpt[], char splitc) + /* - take a string and break it into 2 substrings separated by splitc ; - numpt is number of words returned (1 or 2) - */ +take a string and break it into 2 substrings separated by splitc ; +numpt is number of words returned (1 or 2) +*/ { char rest[MAXSTR], str[MAXSTR], ww[MAXSTR]; int len, i, l; strncpy (str, strin, MAXSTR); len = strlen (strin); - for (i = 0; i < len; i++) - { - if (str[i] == splitc) - { - l = i; - strncpy (ww, str, l); - ww[l] = '\0'; - strpt[0] = strdup (ww); - l = len - (i + 1); - if (l <= 0) - return 1; - strncpy (rest, str + i + 1, l); - rest[l] = '\0'; - strpt[1] = strdup (rest); - return 2; - } + for (i = 0; i < len; i++) { + if (str[i] == splitc) { + l = i; + strncpy (ww, str, l); + ww[l] = '\0'; + strpt[0] = strdup (ww); + l = len - (i + 1); + if (l <= 0) + return 1; + strncpy (rest, str + i + 1, l); + rest[l] = '\0'; + strpt[1] = strdup (rest); + return 2; } + } strpt[0] = strdup (strin); strpt[1] = NULL; return 1; @@ -418,10 +409,9 @@ void printbl (int n) { int i; - for (i = 0; i < n; i++) - { - printf (" "); - } + for (i = 0; i < n; i++) { + printf (" "); + } } void @@ -432,19 +422,19 @@ printnl () void striptrail (char *sss, char c) + /** strip out trailing characters c will usually be ' ' - */ +*/ { int len, i; len = strlen (sss); - for (i = len - 1; i >= 0; --i) - { - if (sss[i] != c) - return; - sss[i] = '\0'; - } + for (i = len - 1; i >= 0; --i) { + if (sss[i] != c) + return; + sss[i] = '\0'; + } } void @@ -453,34 +443,36 @@ catx (char *sxout, char **spt, int n) int i; sxout[0] = CNULL; - for (i = 0; i < n; i++) - { - strcat (sxout, spt[i]); - } + for (i = 0; i < n; i++) { + strcat (sxout, spt[i]); + } + } void catxx (char *sxout, char **spt, int n) + /** like catx but with space between items - */ +*/ { int i; sxout[0] = CNULL; - for (i = 0; i < n; i++) - { - strcat (sxout, spt[i]); - if (i < (n - 1)) - strcat (sxout, " "); - } + for (i = 0; i < n; i++) { + strcat (sxout, spt[i]); + if (i < (n - 1)) + strcat (sxout, " "); + } } + void catxc (char *sxout, char **spt, int n, char c) + /** like catx but with char c between items - */ +*/ { int i; char cc[2]; @@ -490,34 +482,34 @@ catxc (char *sxout, char **spt, int n, char c) cc[0] = c; cc[1] = CNULL; - for (i = 0; i < n; i++) - { - strcat (sxout, spt[i]); - if (i < (n - 1)) - strcat (sxout, cc); - } + for (i = 0; i < n; i++) { + strcat (sxout, spt[i]); + if (i < (n - 1)) + strcat (sxout, cc); + } } void makedfn (char *dirname, char *fname, char *outname, int maxstr) + /** makes full path name. - If fname starts with '/' or dirname = NULL we - so nothing. - outname MUST be allocated of length at least maxstr - */ + If fname starts with '/' or dirname = NULL we + so nothing. + outname MUST be allocated of length at least maxstr +*/ { char *ss; int len; - if ((dirname == NULL) || (fname[0] == '/')) - { - /* if fname starts with / we assume absolute pathname */ - len = strlen (fname); - if (len >= maxstr) - fatalx ("(makedfn) maxstr too short\n"); - strcpy (outname, fname); - return; - } + if ((dirname == NULL) || (fname[0] == '/')) { + +/* if fname starts with / we assume absolute pathname */ + len = strlen (fname); + if (len >= maxstr) + fatalx ("(makedfn) maxstr too short\n"); + strcpy (outname, fname); + return; + } len = strlen (dirname) + strlen (fname) + 1; if (len >= maxstr) fatalx ("(makedfn) maxstr too short\n"); @@ -532,13 +524,14 @@ makedfn (char *dirname, char *fname, char *outname, int maxstr) int substringx (char **ap, char *inx, char *outx, int niter) + /** *ap is original string all occurrences of inx are substituted with outx can loop so be careful !! NB. ap must be on heap. Fixed allocation not supported - */ +*/ { char *a, *pt; char *str; @@ -550,11 +543,10 @@ substringx (char **ap, char *inx, char *outx, int niter) a = *ap; len = strlen (a) + strlen (inx) + strlen (outx) + 1; pt = strstr (a, inx); - if (pt == NULL) - { - return 0; - } - ZALLOC(str, len, char); + if (pt == NULL) { + return 0; + } + ZALLOC (str, len, char); off = pt - a; strncpy (str, a, off); strcpy (str + off, outx); @@ -570,6 +562,7 @@ substringx (char **ap, char *inx, char *outx, int niter) int substring (char **ap, char *inx, char *outx) + /** *ap is original string all occurrences of inx are substituted with outx @@ -577,11 +570,13 @@ substring (char **ap, char *inx, char *outx) NB. ap must be on heap. Fixed allocation not supported - */ +*/ { return (substringx (ap, inx, outx, 0)); } + + int numcols (char *name) // number of cols @@ -595,21 +590,20 @@ numcols (char *name) if (name == NULL) fatalx ("(numlines) no name"); openit (name, &fff, "r"); - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, MAXFF); - if (nsplit == 0) - continue; - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, MAXFF); + if (nsplit == 0) + continue; + sx = spt[0]; + if (sx[0] == '#') { freeup (spt, nsplit); - fclose (fff); - return nsplit; + continue; } + freeup (spt, nsplit); + fclose (fff); + return nsplit; + } + return -1 ; } int @@ -626,20 +620,18 @@ numlines (char *name) if (name == NULL) fatalx ("(numlines) no name"); openit (name, &fff, "r"); - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, MAXFF); - if (nsplit == 0) - continue; - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - ++num; + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, MAXFF); + if (nsplit == 0) + continue; + sx = spt[0]; + if (sx[0] == '#') { freeup (spt, nsplit); + continue; } + ++num; + freeup (spt, nsplit); + } fclose (fff); return num; } @@ -659,20 +651,19 @@ ftest (char *sss) } void -openit (char *name, FILE **fff, char *type) +openit (char *name, FILE ** fff, char *type) { char *ss; if (name == NULL) fatalx ("\n(openit) null name\n"); *fff = fopen (name, type); - if (*fff == NULL) - { - ss = strerror (errno); - printf ("bad open %s\n", name); + if (*fff == NULL) { + ss = strerror (errno); + printf ("bad open %s\n", name); // system("lsof | fgrep np29") ; - fatalx ("can't open file %s of type %s\n error info: %s\n", name, type, - ss); - } + fatalx ("can't open file %s of type %s\n error info: %s\n", name, type, + ss); + } } int @@ -688,43 +679,37 @@ getxx (double **xx, int maxrow, int numcol, char *fname) if (fname == NULL) fff = stdin; - else - { - openit (fname, &fff, "r"); + else { + openit (fname, &fff, "r"); + } + maxff = MAX (MAXFF, numcol); + + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { + freeup (spt, nsplit); + continue; } - maxff = MAX(MAXFF, numcol); - - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (nsplit < numcol) - { - ++nbad; - if (nbad < 10) - printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, - line); - continue; - } - if (num >= maxrow) - fatalx ("too much data\n"); - for (i = 0; i < numcol; i++) - { - xx[i][num] = atof (spt[i]); - } + sx = spt[0]; + if (sx[0] == '#') { freeup (spt, nsplit); - ++num; + continue; + } + if (nsplit < numcol) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; + } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + xx[i][num] = atof (spt[i]); } + freeup (spt, nsplit); + ++num; + } if (fname != NULL) fclose (fff); return num; @@ -740,35 +725,95 @@ clocktime () y = xtime / (double) CLOCKS_PER_SEC; return y; } + int indxstring (char **namelist, int len, char *strid) // look for string in list. Was called indxindex { int k; - for (k = 0; k < len; k++) - { - if (namelist[k] == NULL) - continue; - if (strcmp (namelist[k], strid) == 0) - return k; - } + for (k = 0; k < len; k++) { + if (namelist[k] == NULL) + continue; + if (strcmp (namelist[k], strid) == 0) + return k; + } return -1; } + int indxstringr (char **namelist, int len, char *strid) // look for string in list. Searches array in reverse ; { int k; - for (k = len - 1; k >= 0; k--) - { - if (namelist[k] == NULL) - continue; - if (strcmp (namelist[k], strid) == 0) - return k; - } + for (k = len - 1; k >= 0; k--) { + if (namelist[k] == NULL) + continue; + if (strcmp (namelist[k], strid) == 0) + return k; + } return -1; } +int +getnamesstripcolon (char ****pnames, int maxrow, int numcol, char *fname, + int lo, int hi) +{ + +// count is base 1 + char line[MAXSTR]; + char *spt[MAXFF]; + char *sx; + int nsplit, i, j, num = 0, maxff, numcolp, lcount = 0; + FILE *fff; + int nbad = 0; + char ***names; + + names = *pnames; + if (fname == NULL) + fff = stdin; + else { + openit (fname, &fff, "r"); + } + numcolp = numcol + 1; + maxff = MAX (MAXFF, numcolp); + + while (fgets (line, MAXSTR, fff) != NULL) { + subcolon (line); + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { + freeup (spt, nsplit); + continue; + } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; + } + if (nsplit < numcol) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; + } + ++lcount; + if ((lcount < lo) || (lcount > hi)) { + freeup (spt, nsplit); + continue; + } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + names[i][num] = strdup (spt[i]); + } + freeup (spt, nsplit); + ++num; + } + if (fname != NULL) + fclose (fff); + return num; +} + int getnameslohi (char ****pnames, int maxrow, int numcol, char *fname, int lo, int hi) @@ -786,50 +831,43 @@ getnameslohi (char ****pnames, int maxrow, int numcol, char *fname, int lo, names = *pnames; if (fname == NULL) fff = stdin; - else - { - openit (fname, &fff, "r"); - } + else { + openit (fname, &fff, "r"); + } numcolp = numcol + 1; - maxff = MAX(MAXFF, numcolp); - - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (nsplit < numcol) - { - ++nbad; - if (nbad < 10) - printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, - line); - continue; - } - ++lcount; - if ((lcount < lo) || (lcount > hi)) - { - freeup (spt, nsplit); - continue; - } - if (num >= maxrow) - fatalx ("too much data\n"); - for (i = 0; i < numcol; i++) - { - names[i][num] = strdup (spt[i]); - } + maxff = MAX (MAXFF, numcolp); + + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { freeup (spt, nsplit); - ++num; + continue; + } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; + } + if (nsplit < numcol) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; } + ++lcount; + if ((lcount < lo) || (lcount > hi)) { + freeup (spt, nsplit); + continue; + } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + names[i][num] = strdup (spt[i]); + } + freeup (spt, nsplit); + ++num; + } if (fname != NULL) fclose (fff); return num; @@ -850,44 +888,38 @@ getnames (char ****pnames, int maxrow, int numcol, char *fname) names = *pnames; if (fname == NULL) fff = stdin; - else - { - openit (fname, &fff, "r"); - } + else { + openit (fname, &fff, "r"); + } numcolp = numcol + 1; - maxff = MAX(MAXFF, numcolp); - - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (nsplit < numcol) - { - ++nbad; - if (nbad < 10) - printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, - line); - continue; - } - if (num >= maxrow) - fatalx ("too much data\n"); - for (i = 0; i < numcol; i++) - { - names[i][num] = strdup (spt[i]); - } + maxff = MAX (MAXFF, numcolp); + + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { freeup (spt, nsplit); - ++num; + continue; + } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; + } + if (nsplit < numcol) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + names[i][num] = strdup (spt[i]); + } + freeup (spt, nsplit); + ++num; + } if (fname != NULL) fclose (fff); return num; @@ -909,56 +941,51 @@ getxxnames (char ***pnames, double **xx, int maxrow, int numcol, char *fname) names = *pnames; if (fname == NULL) fff = stdin; - else - { - openit (fname, &fff, "r"); - } + else { + openit (fname, &fff, "r"); + } numcolp = numcol + 1; - maxff = MAX(MAXFF, numcolp); - - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (names != NULL) - names[num] = strdup (sx); - if (nsplit < numcolp) - { - ++nbad; - if (nbad < 10) - printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, - line); - continue; - } - if (num >= maxrow) - fatalx ("too much data\n"); - for (i = 0; i < numcol; i++) - { - xx[i][num] = atof (spt[i + 1]); - } + maxff = MAX (MAXFF, numcolp); + + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { freeup (spt, nsplit); - ++num; + continue; + } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; + } + if (names != NULL) + names[num] = strdup (sx); + if (nsplit < numcolp) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + xx[i][num] = atof (spt[i + 1]); + } + freeup (spt, nsplit); + ++num; + } if (fname != NULL) fclose (fff); return num; } int -getxxnamesf (char ***pnames, double **xx, int maxrow, int numcol, FILE *fff) +getxxnamesf (char ***pnames, double **xx, int maxrow, int numcol, FILE * fff) + /** - like getxxnames but file already open - */ +like getxxnames but file already open +*/ { #define MAXFF 50 @@ -974,49 +1001,46 @@ getxxnamesf (char ***pnames, double **xx, int maxrow, int numcol, FILE *fff) names = *pnames; numcolp = numcol + 1; - maxff = MAX(MAXFF, numcolp); - - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (names != NULL) - names[num] = strdup (sx); - if (nsplit < numcolp) - { - ++nbad; - if (nbad < 10) - printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, - line); - continue; - } - if (num >= maxrow) - fatalx ("too much data\n"); - for (i = 0; i < numcol; i++) - { - xx[i][num] = atof (spt[i + 1]); - } + maxff = MAX (MAXFF, numcolp); + + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { freeup (spt, nsplit); - ++num; + continue; + } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; + } + if (names != NULL) + names[num] = strdup (sx); + if (nsplit < numcolp) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + xx[i][num] = atof (spt[i + 1]); + } + freeup (spt, nsplit); + ++num; + } return num; } + int getss (char **ss, char *fname) + /** get list of names - */ +*/ { char line[MAXSTR]; @@ -1026,46 +1050,91 @@ getss (char **ss, char *fname) int nsplit, i, j, num = 0, maxff; FILE *fff; + if (fname == NULL) fff = stdin; - else - { - openit (fname, &fff, "r"); - } + else { + openit (fname, &fff, "r"); + } maxff = MAXFF; - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - if (nsplit < 1) - { - continue; - } - ss[num] = strdup (spt[0]); + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { freeup (spt, nsplit); - ++num; + continue; + } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; } + if (nsplit < 1) { + continue; + } + ss[num] = strdup (spt[0]); + freeup (spt, nsplit); + ++num; + } if (fname != NULL) fclose (fff); return num; } + +int +checkdup (char **list, int n) +{ + int a, b, t; + for (a = 0; a < n; ++a) { + for (b = a + 1; b < n; ++b) { + t = strcmp (list[a], list[b]); + if (t == 0) + return YES; // dup + } + } + return NO; +} + +void +printdups (char **list, int n) +{ + int a, b, t; + for (a = 0; a < n; ++a) { + for (b = a + 1; b < n; ++b) { + t = strcmp (list[a], list[b]); + if (t == 0) + printf ("dup: %s\n", list[a]); + } + } +} + + +int +loadlist (char **list, char *listname) +// listname is just a list of names ... +// dup check made +{ + int n, a, b, t; + n = getss (list, listname); + for (a = 0; a < n; ++a) { + for (b = a + 1; b < n; ++b) { + t = strcmp (list[a], list[b]); + if (t == 0) { + printf ("(loadlist) duplicate in list: %s\n", list[a]); + fflush (stdout); + fatalx ("(loadlist) duplicate in list: %s\n", list[a]); + } + } + } + return n; +} + char revchar (char c) { char cc; - cc = toupper(c); + cc = toupper (c); if (cc == 'A') return 'T'; if (cc == 'C') @@ -1086,24 +1155,22 @@ crevcomp (char *sout, char *sin) int i, j, t; len = strlen (sin); - ZALLOC(sss, len+1, char); + ZALLOC (sss, len + 1, char); sss[len] = CNULL; - for (i = 0; i < len; ++i) - { - j = len - i - 1; - c = sin[i]; - t = base2num (c); - if (t < 0) - { - sss[j] = c; - continue; - } - cout = num2base (3 - t); - if (islower(c)) - cout = tolower(cout); - sss[j] = cout; + for (i = 0; i < len; ++i) { + j = len - i - 1; + c = sin[i]; + t = base2num (c); + if (t < 0) { + sss[j] = c; + continue; } + cout = num2base (3 - t); + if (islower (c)) + cout = tolower (cout); + sss[j] = cout; + } strcpy (sout, sss); free (sss); } @@ -1116,15 +1183,15 @@ int_string (int a, int len, int base) char *binary = "01"; ss[len] = CNULL; - for (i = 0; i < len; i++) - { - k = t % base; - ss[len - i - 1] = '0' + k; - t = t / base; - } + for (i = 0; i < len; i++) { + k = t % base; + ss[len - i - 1] = '0' + k; + t = t / base; + } return ss; // fragile } + char * binary_string (int a, int len) { @@ -1133,12 +1200,11 @@ binary_string (int a, int len) char *binary = "01"; ss[len] = CNULL; - for (i = 0; i < len; i++) - { - k = t % 2; - ss[len - i - 1] = binary[k]; - t = t / 2; - } + for (i = 0; i < len; i++) { + k = t % 2; + ss[len - i - 1] = binary[k]; + t = t / 2; + } return ss; // fragile } @@ -1188,31 +1254,30 @@ num2base (int num) return bases[num]; } + int base2num (char c) - { char cc; - cc = toupper(c); + cc = toupper (c); - switch (cc) - { - case 'A': - return 0; - break; - case 'C': - return 1; - break; - case 'G': - return 2; - break; - case 'T': - return 3; - break; - default: - return -1; - } + switch (cc) { + case 'A': + return 0; + break; + case 'C': + return 1; + break; + case 'G': + return 2; + break; + case 'T': + return 3; + break; + default: + return -1; + } } int @@ -1222,25 +1287,26 @@ string_binary (char *sx) char c; len = strlen (sx); - ZALLOC(aa, len, int); + ZALLOC (aa, len, int); - for (i = 0; i < len; i++) - { + for (i = 0; i < len; i++) { - c = sx[i]; - if (c == '0') - continue; - if (c != '1') - fatalx ("bad string: %s\n", sx); - aa[i] = 1; - } + c = sx[i]; + if (c == '0') + continue; + if (c != '1') + fatalx ("bad string: %s\n", sx); + aa[i] = 1; + } t = kodeitb (aa, len, 2); free (aa); return t; } + void freestring (char **ss) + /* note extra indirection */ { if (*ss == NULL) @@ -1253,10 +1319,9 @@ void copystrings (char **sa, char **sb, int n) { int i; - for (i = 0; i < n; ++i) - { - sb[i] = strdup (sa[i]); - } + for (i = 0; i < n; ++i) { + sb[i] = strdup (sa[i]); + } } void @@ -1265,39 +1330,35 @@ printstringsw (char **ss, int n, int slen, int width) int k, kmod; char fmt[10], s1[5]; - sprintf (s1, "%ds", slen); + sprintf (s1, "%ds ", slen); strcpy (fmt, "%"); strcat (fmt, s1); - for (k = 0; k < n; ++k) - { - if (ss[k] != NULL) - printf (fmt, ss[k]); - else - printf (fmt, "NULL"); - kmod = (k + 1) % width; - if ((kmod == 0) && (k < (n - 1))) - { - printnl (); - } + for (k = 0; k < n; ++k) { + if (ss[k] != NULL) + printf (fmt, ss[k]); + else + printf (fmt, "NULL"); + kmod = (k + 1) % width; + if ((kmod == 0) && (k < (n - 1))) { + printnl (); } + } printnl (); } void printstrings (char **ss, int n) - { int k; - for (k = 0; k < n; ++k) - { - if (ss[k] != NULL) - printf ("%s", ss[k]); - else - printf ("%s", "NULL"); - printnl (); - } + for (k = 0; k < n; ++k) { + if (ss[k] != NULL) + printf ("%s", ss[k]); + else + printf ("%s", "NULL"); + printnl (); + } } int @@ -1327,16 +1388,16 @@ compbase (char x) return x; } + void mkupper (char *sx) { int len, k; len = strlen (sx); - for (k = 0; k < len; ++k) - { - sx[k] = toupper(sx[k]); - } + for (k = 0; k < len; ++k) { + sx[k] = toupper (sx[k]); + } } void @@ -1345,12 +1406,12 @@ mklower (char *sx) int len, k; len = strlen (sx); - for (k = 0; k < len; ++k) - { - sx[k] = tolower(sx[k]); - } + for (k = 0; k < len; ++k) { + sx[k] = tolower (sx[k]); + } } + char * strstrx (char *s1, char *s2) // like strstr but case insensitive @@ -1361,14 +1422,14 @@ strstrx (char *s1, char *s2) ss1 = strdup (s1); ss2 = strdup (s2); + mkupper (ss1); mkupper (ss2); spt = strstr (ss1, ss2); - if (spt != NULL) - { - spt = s1 + (spt - ss1); - } + if (spt != NULL) { + spt = s1 + (spt - ss1); + } freestring (&ss1); freestring (&ss2); @@ -1377,6 +1438,7 @@ strstrx (char *s1, char *s2) } + int getjjnames (char ***pnames, int **jj, int maxrow, int numcol, char *fname) { @@ -1392,49 +1454,44 @@ getjjnames (char ***pnames, int **jj, int maxrow, int numcol, char *fname) names = *pnames; if (fname == NULL) fff = stdin; - else - { - openit (fname, &fff, "r"); - } + else { + openit (fname, &fff, "r"); + } numcolp = numcol + 1; - maxff = MAX(MAXFF, numcolp); - - while (fgets (line, MAXSTR, fff) != NULL) - { - nsplit = splitup (line, spt, maxff); - if (nsplit == 0) - { - freeup (spt, nsplit); - continue; - } - sx = spt[0]; - if (sx[0] == '#') - { - freeup (spt, nsplit); - continue; - } - names[num] = strdup (sx); - if (nsplit < numcolp) - { - ++nbad; - if (nbad < 10) - printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, - line); - continue; - } - if (num >= maxrow) - fatalx ("too much data\n"); - for (i = 0; i < numcol; i++) - { - jj[i][num] = atoi (spt[i + 1]); - } + maxff = MAX (MAXFF, numcolp); + + while (fgets (line, MAXSTR, fff) != NULL) { + nsplit = splitup (line, spt, maxff); + if (nsplit == 0) { freeup (spt, nsplit); - ++num; + continue; } + sx = spt[0]; + if (sx[0] == '#') { + freeup (spt, nsplit); + continue; + } + names[num] = strdup (sx); + if (nsplit < numcolp) { + ++nbad; + if (nbad < 10) + printf ("+++ bad line: nsplit: %d numcol: %d\n%s\n", nsplit, numcol, + line); + continue; + } + if (num >= maxrow) + fatalx ("too much data\n"); + for (i = 0; i < numcol; i++) { + jj[i][num] = atoi (spt[i + 1]); + } + freeup (spt, nsplit); + ++num; + } if (fname != NULL) fclose (fff); return num; } + int isiub (char iub) { @@ -1471,65 +1528,65 @@ iubdekode (char *aa, char iub) char a[5]; - switch (iub) - { - - case 'A': - strcpy (a, "A"); - break; - case 'C': - strcpy (a, "C"); - break; - case 'G': - strcpy (a, "G"); - break; - case 'T': - strcpy (a, "T"); - break; - case 'M': - strcpy (a, "AC"); - break; - case 'R': - strcpy (a, "AG"); - break; - case 'W': - strcpy (a, "AT"); - break; - case 'S': - strcpy (a, "CG"); - break; - case 'Y': - strcpy (a, "CT"); - break; - case 'K': - strcpy (a, "GT"); - break; - case 'V': - strcpy (a, "ACG"); - break; - case 'H': - strcpy (a, "ACT"); - break; - case 'D': - strcpy (a, "AGT"); - break; - case 'B': - strcpy (a, "CGT"); - break; - case 'X': - strcpy (a, "ACGT"); - break; - case 'N': - strcpy (a, "ACGT"); - break; - - default: - a[0] = CNULL; - } + switch (iub) { + + case 'A': + strcpy (a, "A"); + break; + case 'C': + strcpy (a, "C"); + break; + case 'G': + strcpy (a, "G"); + break; + case 'T': + strcpy (a, "T"); + break; + case 'M': + strcpy (a, "AC"); + break; + case 'R': + strcpy (a, "AG"); + break; + case 'W': + strcpy (a, "AT"); + break; + case 'S': + strcpy (a, "CG"); + break; + case 'Y': + strcpy (a, "CT"); + break; + case 'K': + strcpy (a, "GT"); + break; + case 'V': + strcpy (a, "ACG"); + break; + case 'H': + strcpy (a, "ACT"); + break; + case 'D': + strcpy (a, "AGT"); + break; + case 'B': + strcpy (a, "CGT"); + break; + case 'X': + strcpy (a, "ACGT"); + break; + case 'N': + strcpy (a, "ACGT"); + break; + + default: + a[0] = CNULL; + } if (aa != NULL) strcpy (aa, a); return strlen (a); } + int iubcbases (char *cbases, char iub) // crack iub into 2 bases (which may agree) @@ -1540,6 +1597,7 @@ iubcbases (char *cbases, char iub) int nuu; nuu = iubdekode (uu, iub); + if (nuu < 1) return -1; if (nuu > 2) @@ -1566,9 +1624,26 @@ ishet (char c) return NO; } +int +cttype (char cx) +// useful for CPG detection +{ + char cc; + cc = toupper (cx); + + if (cc == 'A') + return 2; + if (cc == 'C') + return 1; + if (cc == 'G') + return 2; + if (cc == 'T') + return 1; + return -1; +} + char * lastff (char *sss) - { char *sx; sx = strrchr (sss, '/'); @@ -1586,12 +1661,15 @@ char2int (char cc) return x; } + char int2char (int x) { char c; c = (char) ('0' + x); + + return c ; } void @@ -1612,11 +1690,10 @@ numcmatch (char *cc, int len, char c) { int k, t = 0; - for (k = 0; k < len; ++k) - { - if (cc[k] == c) - ++t; - } + for (k = 0; k < len; ++k) { + if (cc[k] == c) + ++t; + } return t; @@ -1627,13 +1704,107 @@ numcnomatch (char *cc, int len, char c) { int k, t = 0; - for (k = 0; k < len; ++k) - { - if (cc[k] != c) - ++t; - } + for (k = 0; k < len; ++k) { + if (cc[k] != c) + ++t; + } return t; } +char * +findupper (char *s) +// CNULL not tested +{ + char x; + int len, k; + + len = strlen (s); + + for (k = 0; k < len; ++k) { + if (isupper (s[k])) + return s + k; + } + + return NULL; + +} + +char * +strnotchar (char *s, char c) +// CNULL not tested +{ + char x; + int len, k; + + len = strlen (s); + + for (k = 0; k < len; ++k) { + if (s[k] != c) + return s + k; + } + + return NULL; + +} + +char +readtonl (FILE * fff) +{ + char c; + + for (;;) { + c = fgetc (fff); + if (c == EOF) + break; + if (c == CNL) + break; + } + return c; +} + +char * +fgetstrap (char *buff, int maxlen, FILE * fff, int *ret) +// fgets with long lines trapped +{ + int len; + char c; + + if (fgets (buff, maxlen, fff) == NULL) + return NULL; + + *ret = 1; + len = strlen (buff); + if (buff[len - 1] == CNL) + return buff; + + *ret = 0; + + c = readtonl (fff); + buff[0] = c; + buff[1] = CNULL; + return buff; + +} +int filehash(char *name) +{ +#define MAXKL 256 + FILE *fff; + char line[MAXKL]; + int num = 0; + int hash = 0, thash ; + + num = 0; + if (name == NULL) + fatalx ("(filehash) no name"); + openit (name, &fff, "r"); + while (fgets (line, MAXKL, fff) != NULL) { + thash = stringhash(line) ; + hash += xhash1(thash ^ num) ; + ++num ; + } + fclose (fff); + return abs(hash) ; +} + diff --git a/src/nicksrc/vsubs.c b/src/nicksrc/vsubs.c index 453d366..52930e8 100644 --- a/src/nicksrc/vsubs.c +++ b/src/nicksrc/vsubs.c @@ -1,15 +1,16 @@ #include #include #include +#include #include -#include "strsubs.h" -#include "vsubs.h" +#include "strsubs.h" +#include "vsubs.h" /** tiny routines BLAS? a small library to do simple arithmetic on 1D vectors with no skips - */ +*/ void vsp (double *a, double *b, double c, int n) { @@ -17,6 +18,7 @@ vsp (double *a, double *b, double c, int n) for (i = 0; i < n; i++) a[i] = b[i] + c; } + void vst (double *a, double *b, double c, int n) { @@ -24,6 +26,7 @@ vst (double *a, double *b, double c, int n) for (i = 0; i < n; i++) a[i] = b[i] * c; } + void vvt (double *a, double *b, double *c, int n) { @@ -31,6 +34,7 @@ vvt (double *a, double *b, double *c, int n) for (i = 0; i < n; i++) a[i] = b[i] * c[i]; } + void vvp (double *a, double *b, double *c, int n) { @@ -38,6 +42,7 @@ vvp (double *a, double *b, double *c, int n) for (i = 0; i < n; i++) a[i] = b[i] + c[i]; } + void vvm (double *a, double *b, double *c, int n) { @@ -45,85 +50,84 @@ vvm (double *a, double *b, double *c, int n) for (i = 0; i < n; i++) a[i] = b[i] - c[i]; } + void vvd (double *a, double *b, double *c, int n) { int i; - for (i = 0; i < n; i++) - { - if (c[i] == 0.0) - fatalx ("(vvd): zero value in denominator\n"); - a[i] = b[i] / c[i]; - } + for (i = 0; i < n; i++) { + if (c[i] == 0.0) + fatalx ("(vvd): zero value in denominator\n"); + a[i] = b[i] / c[i]; + } } + void vsqrt (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - if (b[i] < 0.0) - fatalx ("(vsqrt): negative value %g\n", b[i]); - if (b[i] == 0.0) - { - a[i] = 0.0; - continue; - } - a[i] = sqrt (b[i]); + for (i = 0; i < n; i++) { + if (b[i] < 0.0) + fatalx ("(vsqrt): negative value %g\n", b[i]); + if (b[i] == 0.0) { + a[i] = 0.0; + continue; } + a[i] = sqrt (b[i]); + } } + void vinvert (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - if (b[i] == 0.0) - fatalx ("(vinvert): zero value\n"); - a[i] = 1.0 / b[i]; - } + for (i = 0; i < n; i++) { + if (b[i] == 0.0) + fatalx ("(vinvert): zero value\n"); + a[i] = 1.0 / b[i]; + } } + void vabs (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - a[i] = fabs (b[i]); - } + for (i = 0; i < n; i++) { + a[i] = fabs (b[i]); + } } + void vlog (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - if (b[i] <= 0.0) - fatalx ("(vlog): negative or zero value %g\n", b[i]); - a[i] = log (b[i]); - } + for (i = 0; i < n; i++) { + if (b[i] <= 0.0) + fatalx ("(vlog): negative or zero value %g\n", b[i]); + a[i] = log (b[i]); + } } + void vlog2 (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - if (b[i] <= 0.0) - fatalx ("(vlog2): negative or zero value %g\n", b[i]); - a[i] = NPlog2 (b[i]); - } + for (i = 0; i < n; i++) { + if (b[i] <= 0.0) + fatalx ("(vlog2): negative or zero value %g\n", b[i]); + a[i] = NPlog2 (b[i]); + } } void vexp (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - a[i] = exp (b[i]); - } + for (i = 0; i < n; i++) { + a[i] = exp (b[i]); + } } + void vclear (double *a, double c, int n) { @@ -131,32 +135,44 @@ vclear (double *a, double c, int n) for (i = 0; i < n; i++) a[i] = c; } + void vzero (double *a, int n) { vclear (a, 0.0, n); } + void cpzero (char **a, int n) { int i; - for (i = 0; i < n; ++i) - { - a[i] = NULL; - } + for (i = 0; i < n; ++i) { + a[i] = NULL; + } } + void cclear (unsigned char *a, unsigned char c, long n) + /** be careful nothing done about NULL at end - */ +*/ { long i; - for (i = 0; i < n; i++) - { - a[i] = c; - } + for (i = 0; i < n; i++) { + a[i] = c; + } +} + +void +charclear (char *a, unsigned char c, long n) +// fussy compiler warnigns about unsigned char conversions avoided +{ + + cclear( (unsigned char *) a, c, n) ; + } + void ivvp (int *a, int *b, int *c, int n) { @@ -164,6 +180,7 @@ ivvp (int *a, int *b, int *c, int n) for (i = 0; i < n; i++) a[i] = b[i] + c[i]; } + void ivvm (int *a, int *b, int *c, int n) { @@ -171,6 +188,7 @@ ivvm (int *a, int *b, int *c, int n) for (i = 0; i < n; i++) a[i] = b[i] - c[i]; } + void ivsp (int *a, int *b, int c, int n) { @@ -178,6 +196,7 @@ ivsp (int *a, int *b, int c, int n) for (i = 0; i < n; i++) a[i] = b[i] + c; } + void ivst (int *a, int *b, int c, int n) { @@ -185,6 +204,7 @@ ivst (int *a, int *b, int c, int n) for (i = 0; i < n; i++) a[i] = b[i] * c; } + void ivclear (int *a, int c, long n) { @@ -192,6 +212,7 @@ ivclear (int *a, int c, long n) for (i = 0; i < n; i++) a[i] = c; } + void lvclear (long *a, long c, long n) { @@ -205,8 +226,16 @@ ivzero (int *a, int n) { ivclear (a, 0, n); } + +void +lvzero (long *a, long n) +{ + lvclear (a, 0, n); +} + double clip (double x, double lo, double hi) + /* clip off values to range [lo,hi] */ { if (x < lo) @@ -219,30 +248,31 @@ clip (double x, double lo, double hi) void ivclip (int *a, int *b, int loval, int hival, int n) { - /* clip off values to range [loval,hival] */ + +/* clip off values to range [loval,hival] */ int i; int t; - for (i = 0; i < n; i++) - { - t = MAX(b[i], loval); - a[i] = MIN(t, hival); - } + for (i = 0; i < n; i++) { + t = MAX (b[i], loval); + a[i] = MIN (t, hival); + } } void vclip (double *a, double *b, double loval, double hival, int n) { - /* clip off values to range [loval,hival] */ + +/* clip off values to range [loval,hival] */ int i; double t; - for (i = 0; i < n; i++) - { - t = MAX(b[i], loval); - a[i] = MIN(t, hival); - } + for (i = 0; i < n; i++) { + t = MAX (b[i], loval); + a[i] = MIN (t, hival); + } } + void vmaxmin (double *a, int n, double *max, double *min) { @@ -251,21 +281,22 @@ vmaxmin (double *a, int n, double *max, double *min) double tmax, tmin; tmax = tmin = a[0]; - for (i = 1; i < n; i++) - { - tmax = MAX(tmax, a[i]); - tmin = MIN(tmin, a[i]); - } + for (i = 1; i < n; i++) { + tmax = MAX (tmax, a[i]); + tmin = MIN (tmin, a[i]); + } if (max != NULL) *max = tmax; if (min != NULL) *min = tmin; } + void vlmaxmin (double *a, int n, int *pmax, int *pmin) + /** return location - */ +*/ { int i; @@ -274,24 +305,22 @@ vlmaxmin (double *a, int n, int *pmax, int *pmin) tmax = tmin = a[0]; lmax = lmin = 0; - for (i = 1; i < n; i++) - { - if (a[i] > tmax) - { - tmax = a[i]; - lmax = i; - } - if (a[i] < tmin) - { - tmin = a[i]; - lmin = i; - } + for (i = 1; i < n; i++) { + if (a[i] > tmax) { + tmax = a[i]; + lmax = i; + } + if (a[i] < tmin) { + tmin = a[i]; + lmin = i; } + } if (pmax != NULL) *pmax = lmax; if (pmin != NULL) *pmin = lmin; } + void ivmaxmin (int *a, int n, int *max, int *min) { @@ -300,16 +329,16 @@ ivmaxmin (int *a, int n, int *max, int *min) int tmax, tmin; tmax = tmin = a[0]; - for (i = 1; i < n; i++) - { - tmax = MAX(tmax, a[i]); - tmin = MIN(tmin, a[i]); - } + for (i = 1; i < n; i++) { + tmax = MAX (tmax, a[i]); + tmin = MIN (tmin, a[i]); + } if (max != NULL) *max = tmax; if (min != NULL) *min = tmin; } + int minivec (int *a, int n) { @@ -332,9 +361,10 @@ maxivec (int *a, int n) void ivlmaxmin (int *a, int n, int *pmax, int *pmin) + /** return location - */ +*/ { int i; @@ -343,24 +373,22 @@ ivlmaxmin (int *a, int n, int *pmax, int *pmin) tmax = tmin = a[0]; lmax = lmin = 0; - for (i = 1; i < n; i++) - { - if (a[i] > tmax) - { - tmax = a[i]; - lmax = i; - } - if (a[i] < tmin) - { - tmin = a[i]; - lmin = i; - } + for (i = 1; i < n; i++) { + if (a[i] > tmax) { + tmax = a[i]; + lmax = i; } + if (a[i] < tmin) { + tmin = a[i]; + lmin = i; + } + } if (pmax != NULL) *pmax = lmax; if (pmin != NULL) *pmin = lmin; } + double vdot (double *a, double *b, int n) { @@ -372,13 +400,14 @@ vdot (double *a, double *b, int n) return ans; } + double corr (double *a, double *b, int n) { double v12, v11, v22, y1, y2, y; double *aa, *bb; - ZALLOC(aa, n, double); - ZALLOC(bb, n, double); + ZALLOC (aa, n, double); + ZALLOC (bb, n, double); y1 = asum (a, n) / (double) n; y2 = asum (b, n) / (double) n; @@ -393,11 +422,13 @@ corr (double *a, double *b, int n) if (y == 0.0) fatalx ("(corr) constant vector\n"); + free (aa); free (bb); return (v12 / sqrt (y)); } + double corrx (double *a, double *b, int n) // like corr but constant vec returns 0 @@ -405,8 +436,8 @@ corrx (double *a, double *b, int n) double v12, v11, v22, y1, y2, y; double *aa, *bb; - ZALLOC(aa, n, double); - ZALLOC(bb, n, double); + ZALLOC (aa, n, double); + ZALLOC (bb, n, double); y1 = asum (a, n) / (double) n; y2 = asum (b, n) / (double) n; @@ -427,6 +458,7 @@ corrx (double *a, double *b, int n) } + double variance (double *a, int n) { @@ -434,7 +466,7 @@ variance (double *a, int n) double *aa; double y1, y2; - ZALLOC(aa, n, double); + ZALLOC (aa, n, double); y1 = asum (a, n) / (double) n; vsp (aa, a, -y1, n); @@ -447,39 +479,53 @@ variance (double *a, int n) void getdiag (double *a, double *b, int n) + /* extract diagonal */ { int i, k; - for (i = 0; i < n; i++) - { - k = i * n + i; - a[i] = b[k]; - } + for (i = 0; i < n; i++) { + k = i * n + i; + a[i] = b[k]; + } } void setdiag (double *a, double *diag, int n) + /* set diagonal matrix */ { int i, k; vzero (a, n * n); - for (i = 0; i < n; i++) - { - k = i * n + i; - a[k] = diag[i]; - } + for (i = 0; i < n; i++) { + k = i * n + i; + a[k] = diag[i]; + } } +void +adddiag (double *a, double *diag, int n) + +/* add diagonal matrix */ +{ + int i, k; + + for (i = 0; i < n; i++) { + k = i * n + i; + a[k] += diag[i]; + } +} + + + void copyarr (double *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - b[i] = a[i]; - } + for (i = 0; i < n; i++) { + b[i] = a[i]; + } } void @@ -487,28 +533,26 @@ revarr (double *b, double *a, int n) { int i; double *x; - ZALLOC(x, n, double); - for (i = 0; i < n; i++) - { - x[n - i - 1] = a[i]; - } + ZALLOC (x, n, double); + for (i = 0; i < n; i++) { + x[n - i - 1] = a[i]; + } copyarr (x, b, n); free (x); } + void revuiarr (unsigned int *b, unsigned int *a, int n) { int i; unsigned int *x; - ZALLOC(x, n, unsigned int); - for (i = 0; i < n; i++) - { - x[n - i - 1] = a[i]; - } - for (i = 0; i < n; i++) - { - b[i] = x[i]; - } + ZALLOC (x, n, unsigned int); + for (i = 0; i < n; i++) { + x[n - i - 1] = a[i]; + } + for (i = 0; i < n; i++) { + b[i] = x[i]; + } free (x); } @@ -517,33 +561,42 @@ reviarr (int *b, int *a, int n) { int i; int *x; - ZALLOC(x, n, int); - for (i = 0; i < n; i++) - { - x[n - i - 1] = a[i]; - } + ZALLOC (x, n, int); + for (i = 0; i < n; i++) { + x[n - i - 1] = a[i]; + } copyiarr (x, b, n); free (x); } + void copyiarr (int *a, int *b, int n) { int i; - for (i = 0; i < n; i++) - { - b[i] = a[i]; - } + for (i = 0; i < n; i++) { + b[i] = a[i]; + } } + +void +copylarr (long *a, long *b, int n) +{ + int i; + for (i = 0; i < n; i++) { + b[i] = a[i]; + } +} + void copyiparr (int **a, int **b, int n) { int i; - for (i = 0; i < n; i++) - { - b[i] = a[i]; - } + for (i = 0; i < n; i++) { + b[i] = a[i]; + } } + void dpermute (double *a, int *ind, int len) { @@ -551,21 +604,20 @@ dpermute (double *a, int *ind, int len) int i, k; double *rrr; - ZALLOC(rrr, len, double); + ZALLOC (rrr, len, double); - for (i = 0; i < len; i++) - { - rrr[i] = a[i]; - } + for (i = 0; i < len; i++) { + rrr[i] = a[i]; + } - for (i = 0; i < len; i++) - { - k = ind[i]; - a[i] = rrr[k]; - } + for (i = 0; i < len; i++) { + k = ind[i]; + a[i] = rrr[k]; + } free (rrr); } + void ipermute (int *a, int *ind, int len) { @@ -573,18 +625,18 @@ ipermute (int *a, int *ind, int len) int i, k; int *rrr; - ZALLOC(rrr, len, int); + ZALLOC (rrr, len, int); copyiarr (a, rrr, len); - for (i = 0; i < len; i++) - { - k = ind[i]; - a[i] = rrr[k]; - } + for (i = 0; i < len; i++) { + k = ind[i]; + a[i] = rrr[k]; + } free (rrr); } + void dppermute (double **a, int *ind, int len) { @@ -592,21 +644,20 @@ dppermute (double **a, int *ind, int len) int i, k; double **rrr; - ZALLOC(rrr, len, double *); + ZALLOC (rrr, len, double *); - for (i = 0; i < len; i++) - { - rrr[i] = a[i]; - } + for (i = 0; i < len; i++) { + rrr[i] = a[i]; + } - for (i = 0; i < len; i++) - { - k = ind[i]; - a[i] = rrr[k]; - } + for (i = 0; i < len; i++) { + k = ind[i]; + a[i] = rrr[k]; + } free (rrr); } + void ippermute (int **a, int *ind, int len) { @@ -614,18 +665,16 @@ ippermute (int **a, int *ind, int len) int i, k; int **rrr; - ZALLOC(rrr, len, int *); + ZALLOC (rrr, len, int *); - for (i = 0; i < len; i++) - { - rrr[i] = a[i]; - } + for (i = 0; i < len; i++) { + rrr[i] = a[i]; + } - for (i = 0; i < len; i++) - { - k = ind[i]; - a[i] = rrr[k]; - } + for (i = 0; i < len; i++) { + k = ind[i]; + a[i] = rrr[k]; + } free (rrr); } @@ -640,6 +689,7 @@ asum (double *a, int n) return ans; } + int intsum (int *a, int n) { @@ -650,6 +700,7 @@ intsum (int *a, int n) return ans; } + long longsum (long *a, int n) { @@ -660,6 +711,7 @@ longsum (long *a, int n) return ans; } + int idot (int *a, int *b, int n) { @@ -674,6 +726,7 @@ idot (int *a, int *b, int n) int iprod (int *a, int n) + /* overflow not checked */ { int i; @@ -684,8 +737,10 @@ iprod (int *a, int n) return ans; } + double aprod (double *a, int n) + /* overflow not checked */ { int i; @@ -695,6 +750,7 @@ aprod (double *a, int n) return ans; } + double asum2 (double *a, int n) { @@ -710,8 +766,8 @@ double trace (double *a, int n) { double *diags, t; - ZALLOC(diags,n,double); - getdiag (diags, a, n); /* extract diagonal */ + ZALLOC (diags, n, double); + getdiag (diags, a, n); /* extract diagonal */ t = asum (diags, n); free (diags); return t; @@ -720,202 +776,239 @@ trace (double *a, int n) int nnint (double x) { - long int - lrint (double x); + long int lrint (double x); // double round(double x) ; return (int) lrint (x); } + void countcat (int *tags, int n, int *ncat, int nclass) -/* simple frequency count of integer array */ +/* simple frequency count of integer array */ { int i, k; ivzero (ncat, nclass); - for (i = 0; i < n; i++) - { - k = tags[i]; - if ((k < 0) || (k >= nclass)) - fatalx ("(countcat) bounds error\n"); - ++ncat[k]; - } + for (i = 0; i < n; i++) { + k = tags[i]; + if ((k < 0) || (k >= nclass)) + fatalx ("(countcat) bounds error\n"); + ++ncat[k]; + } } + void rowsum (double *a, double *rr, int n) { int i, j; vclear (rr, 0.0, n); - for (i = 0; i < n; i++) - { - for (j = 0; j < n; j++) - { - rr[j] += a[i + j * n]; - } + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + rr[j] += a[i + j * n]; } + } } + void colsum (double *a, double *cc, int n) { int i, j; vclear (cc, 0.0, n); - for (i = 0; i < n; i++) - { - for (j = 0; j < n; j++) - { - cc[i] += a[i + j * n]; - } + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + cc[i] += a[i + j * n]; } + } } + void rrsum (double *a, double *rr, int m, int n) { int i, j; vclear (rr, 0.0, n); - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - rr[j] += a[i + j * m]; - } + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + rr[j] += a[i + j * m]; } + } } + void ccsum (double *a, double *cc, int m, int n) { int i, j; vclear (cc, 0.0, m); - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - cc[i] += a[i + j * m]; - } + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + cc[i] += a[i + j * m]; } + } } + void -printmatfile (double *a, int m, int n, FILE *fff) +printmatfile (double *a, int m, int n, FILE * fff) + /** print a matrix n wide m rows - */ +*/ { printmatwfile (a, m, n, 5, fff); } + void -printmatwfile (double *a, int m, int n, int w, FILE *fff) +printmatwfile (double *a, int m, int n, int w, FILE * fff) + /** print a matrix n wide m rows w to a row - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - fprintf (fff, "%9.3f ", a[i * n + j]); - jmod = (j + 1) % w; - if ((jmod == 0) && (j < (n - 1))) - { - fprintf (fff, " ...\n"); - } - } - fprintf (fff, "\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + fprintf (fff, "%9.3f ", a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + fprintf (fff, " ...\n"); + } } + fprintf (fff, "\n"); + } +} + +void +printmatx (double *a, int m, int n) + +/** + print a matrix n wide m rows no final nl +*/ +{ + printmatwx (a, m, n, 5); } + void printmat (double *a, int m, int n) + /** print a matrix n wide m rows - */ +*/ { printmatw (a, m, n, 5); } + void -printmatw (double *a, int m, int n, int w) +printmatwx (double *a, int m, int n, int w) + /** print a matrix n wide m rows w to a row - */ + no final nl +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf ("%9.3f ", a[i * n + j]); - jmod = (j + 1) % w; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%9.3f ", a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } + } + if (i < (m - 1)) printf ("\n"); + } +} + +void +printmatw (double *a, int m, int n, int w) + +/** + print a matrix n wide m rows w to a row +*/ +{ + int i, j, jmod; + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%9.3f ", a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + printf ("\n"); + } } + void printmatl (double *a, int m, int n) + /** print a matrix n wide m rows - */ +*/ { printmatwl (a, m, n, 5); } + void printmatwl (double *a, int m, int n, int w) + /** print a matrix n wide m rows w to a row 15.9f format - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf ("%15.9f ", a[i * n + j]); - jmod = (j + 1) % w; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } - printf ("\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%15.9f ", a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + printf ("\n"); + } } + void printmatwf (double *a, int m, int n, int w, char *format) + /** print a matrix n wide m rows w to a row with format no spacing introduced here. User must supply - */ +*/ { int i, j, jmod; - if (format == NULL) - { - printmatw (a, m, n, w); - return; - } - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf (format, a[i * n + j]); - jmod = (j + 1) % w; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } - printf ("\n"); + if (format == NULL) { + printmatw (a, m, n, w); + return; + } + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf (format, a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + printf ("\n"); + } +} + +void +printmat2D (double **a, int m, int n) +{ + int k; + for (k = 0; k < m; ++k) { + printf ("%3d: ", k); + printmat (a[k], 1, n); + } } void int2c (char *cc, int *b, int n) { int i; - for (i = 0; i < n; i++) - { - cc[i] = (char) b[i]; - } + for (i = 0; i < n; i++) { + cc[i] = (char) b[i]; + } cc[n] = '\0'; } @@ -923,170 +1016,215 @@ void floatit (double *a, int *b, int n) { int i; - for (i = 0; i < n; i++) - { - a[i] = (double) b[i]; - } + for (i = 0; i < n; i++) { + a[i] = (double) b[i]; + } } + void -printimatwfile (int *a, int m, int n, int w, FILE *fff) +printimatwfile (int *a, int m, int n, int w, FILE * fff) + /** print a matrix n wide m rows w to a row - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - fprintf (fff, "%5d ", a[i * n + j]); - jmod = (j + 1) % w; - if ((jmod == 0) && (j < (n - 1))) - { - fprintf (fff, " ...\n"); - } - } - fprintf (fff, "\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + fprintf (fff, "%5d ", a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + fprintf (fff, " ...\n"); + } } + fprintf (fff, "\n"); + } } + void printimatw (int *a, int m, int n, int w) + /** print a matrix n wide m rows w to a row - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf ("%5d ", a[i * n + j]); - jmod = (j + 1) % w; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } - printf ("\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%5d ", a[i * n + j]); + jmod = (j + 1) % w; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + printf ("\n"); + } } + void printimatx (int *a, int m, int n) + /** print a matrix n wide m rows no final new line - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf ("%5d ", a[i * n + j]); - jmod = (j + 1) % 10; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%5d ", a[i * n + j]); + jmod = (j + 1) % 10; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + } } + void -printimatfile (int *a, int m, int n, FILE *fff) +printimatfile (int *a, int m, int n, FILE * fff) + /** print a matrix n wide m rows - */ +*/ { printimatwfile (a, m, n, 10, fff); } +void +printimat2D (int **a, int m, int n) +{ + int k; + + for (k = 0; k < m; ++k) { + printimat (a[k], 1, n); + } +} + +void +printimat1 (int *a, int m, int n) + +/** + print a matrix n wide m rows, %1d format +*/ +{ + int i, j, jmod; + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%1d", a[i * n + j]); + jmod = (j + 1) % 50; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } + } + printf ("\n"); + } +} + void printimat (int *a, int m, int n) + /** print a matrix n wide m rows - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf ("%5d ", a[i * n + j]); - jmod = (j + 1) % 10; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } - printf ("\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%5d ", a[i * n + j]); + jmod = (j + 1) % 10; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + printf ("\n"); + } } + void -printimatlfile (int *a, int m, int n, FILE *fff) +printimatlfile (int *a, int m, int n, FILE * fff) + /** print a matrix n wide m rows %10d format - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - fprintf (fff, "%10d ", a[i * n + j]); - jmod = (j + 1) % 10; - if ((jmod == 0) && (j < (n - 1))) - { - fprintf (fff, " ...\n"); - } - } - fprintf (fff, "\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + fprintf (fff, "%10d ", a[i * n + j]); + jmod = (j + 1) % 10; + if ((jmod == 0) && (j < (n - 1))) { + fprintf (fff, " ...\n"); + } } + fprintf (fff, "\n"); + } } void printimatl (int *a, int m, int n) + /** print a matrix n wide m rows %10d format - */ +*/ { int i, j, jmod; - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - { - printf ("%10d ", a[i * n + j]); - jmod = (j + 1) % 10; - if ((jmod == 0) && (j < (n - 1))) - { - printf (" ...\n"); - } - } - printf ("\n"); + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%10d ", a[i * n + j]); + jmod = (j + 1) % 10; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } + } + printf ("\n"); + } +} + +void +printimatlx (int *a, int m, int n) + +/** + print a matrix n wide m rows %10d format + no final newline +*/ +{ + int i, j, jmod; + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) { + printf ("%10d ", a[i * n + j]); + jmod = (j + 1) % 10; + if ((jmod == 0) && (j < (n - 1))) { + printf (" ...\n"); + } } + if (i < (m - 1)) + printf ("\n"); + } } void -printstringf (char *ss, int w, FILE *fff) +printstringf (char *ss, int w, FILE * fff) { char *sss; char *sx; - ZALLOC(sss, w+1, char); - cclear (sss, CNULL, w + 1); + ZALLOC (sss, w + 1, char); + cclear ((unsigned char *) sss, CNULL, w + 1); sx = ss; - for (;;) - { - strncpy (sss, sx, w); - if (strlen (sss) <= 0) - break; - sx += w; - fprintf (fff, "%s\n", sss); - } + for (;;) { + strncpy (sss, sx, w); + if (strlen (sss) <= 0) + break; + sx += w; + fprintf (fff, "%s\n", sss); + } free (sss); } + void printstringbasepos (char *ss, int w, int basepos) { @@ -1094,81 +1232,82 @@ printstringbasepos (char *ss, int w, int basepos) char *sx; int pos = basepos; - ZALLOC(sss, w+1, char); - cclear (sss, CNULL, w + 1); + ZALLOC (sss, w + 1, char); + cclear ((unsigned char *) sss, CNULL, w + 1); sx = ss; - for (;;) - { - strncpy (sss, sx, w); - if (strlen (sss) <= 0) - break; - printf ("%12d ", pos); - printf ("%s\n", sss); - sx += w; - pos += w; - } + for (;;) { + strncpy (sss, sx, w); + if (strlen (sss) <= 0) + break; + printf ("%12d ", pos); + printf ("%s\n", sss); + sx += w; + pos += w; + } free (sss); } + + void printstring (char *ss, int w) { printstringf (ss, w, stdout); } + void rndit (double *a, double *b, int n) { int i; - for (i = 0; i < n; ++i) - { - a[i] = nearbyint (b[i]); - } + for (i = 0; i < n; ++i) { + a[i] = nearbyint (b[i]); + } } + void fixit (int *a, double *b, int n) { int i; - for (i = 0; i < n; i++) - { - a[i] = nnint (b[i]); - } + for (i = 0; i < n; i++) { + a[i] = nnint (b[i]); + } } + int findfirst (int *a, int n, int val) { int i; - for (i = 0; i < n; i++) - { - if (a[i] == val) - return i; - } + for (i = 0; i < n; i++) { + if (a[i] == val) + return i; + } return -1; } + int findfirstl (long *a, int n, long val) { int i; - for (i = 0; i < n; i++) - { - if (a[i] == val) - return i; - } + for (i = 0; i < n; i++) { + if (a[i] == val) + return i; + } return -1; } + int findfirstu (unsigned int *a, int n, unsigned int val) { int i; - for (i = 0; i < n; i++) - { - if (a[i] == val) - return i; - } + for (i = 0; i < n; i++) { + if (a[i] == val) + return i; + } return -1; } @@ -1176,11 +1315,10 @@ int findlastu (unsigned int *a, int n, unsigned int val) { int i; - for (i = n - 1; i >= 0; i--) - { - if (a[i] == val) - return i; - } + for (i = n - 1; i >= 0; i--) { + if (a[i] == val) + return i; + } return -1; } @@ -1188,13 +1326,13 @@ int findlast (int *a, int n, int val) { int i; - for (i = n - 1; i >= 0; i--) - { - if (a[i] == val) - return i; - } + for (i = n - 1; i >= 0; i--) { + if (a[i] == val) + return i; + } return -1; } + int binsearch (int *a, int n, int val) // binary search. a sorted in ascending order @@ -1228,6 +1366,7 @@ idperm (int *a, int n) for (i = 0; i < n; i++) a[i] = i; } + double NPlog2 (double y) { @@ -1238,9 +1377,10 @@ NPlog2 (double y) double logfac (int n) + /** log (factorial n)) - */ +*/ { double y, x; x = (double) (n + 1); @@ -1250,6 +1390,7 @@ logfac (int n) double logbino (int n, int k) + /* log n choose k */ { double top, bot; @@ -1264,10 +1405,11 @@ double loghprob (int n, int a, int m, int k) // http://www.math.uah.edu/stat/urn/Hypergeometric.xhtml { - /** - n balls a black. Pick m without replacement - return log prob (k black) - */ + +/** + n balls a black. Pick m without replacement + return log prob (k black) +*/ double ytop, ybot; @@ -1286,11 +1428,13 @@ loghprob (int n, int a, int m, int k) } + double log2fac (int n) + /** log base2 (factorial n)) - */ +*/ { double y, x; x = (double) (n + 1); @@ -1302,18 +1446,18 @@ double addlog (double a, double b) { /* given a = log(A) - b = log(B) - returns log(A+B) - with precautions for overflow etc + b = log(B) + returns log(A+B) + with precautions for overflow etc */ double x, y, z; - x = MIN(a, b); - y = MAX(a, b); + x = MIN (a, b); + y = MAX (a, b); - /** - answer is log(1 + A/B) + log (B) - */ +/** + answer is log(1 + A/B) + log (B) +*/ z = x - y; if (z < -50.0) return y; @@ -1323,26 +1467,28 @@ addlog (double a, double b) } + + double vldot (double *x, double *y, int n) + /** x. log(y) - */ +*/ { double *z, ans; double tiny = 1.0e-19; int i; - ZALLOC(z, n, double); + ZALLOC (z, n, double); vsp (z, y, 1.0e-20, n); vlog (z, z, n); ans = 0.0; - for (i = 0; i < n; i++) - { - if (x[i] > tiny) - ans += x[i] * z[i]; - } + for (i = 0; i < n; i++) { + if (x[i] > tiny) + ans += x[i] * z[i]; + } free (z); return ans; } @@ -1359,7 +1505,8 @@ pow10 (double x) return exp (x * log (10.0)); } -double + +void vpow10 (double *a, double *b, int n) { int i; @@ -1367,7 +1514,7 @@ vpow10 (double *a, double *b, int n) a[i] = exp (b[i] * log (10.0)); } -double +void vlog10 (double *a, double *b, int n) { int i; @@ -1377,67 +1524,68 @@ vlog10 (double *a, double *b, int n) void transpose (double *aout, double *ain, int m, int n) + /** aout and ain must be identical or not overlap does matrix transpose input m vectors of length n (m x n) output n vectors of length m - */ +*/ { double *ttt; int i, j, k1, k2; - if (aout == ain) - { - ZALLOC(ttt, m*n, double); - } + if (aout == ain) { + ZALLOC (ttt, m * n, double); + } else ttt = aout; for (i = 0; i < m; i++) - for (j = 0; j < n; j++) - { - k1 = i * n + j; - k2 = j * m + i; - ttt[k2] = ain[k1]; - } - if (aout == ain) - { - copyarr (ttt, aout, m * n); - free (ttt); + for (j = 0; j < n; j++) { + k1 = i * n + j; + k2 = j * m + i; + ttt[k2] = ain[k1]; } + if (aout == ain) { + copyarr (ttt, aout, m * n); + free (ttt); + } } + int ** initarray_2Dint (int numrows, int numcolumns, int initval) { int i, j; int **array; - ZALLOC(array, numrows, int *); - for (i = 0; i < numrows; i++) - { - ZALLOC(array[i], numcolumns, int); - if (initval != 0) - ivclear (array[i], initval, numcolumns); - } + + ZALLOC (array, numrows, int *); + for (i = 0; i < numrows; i++) { + ZALLOC (array[i], numcolumns, int); + if (initval != 0) + ivclear (array[i], initval, numcolumns); + } return array; } + long ** initarray_2Dlong (int numrows, int numcolumns, int initval) { int i, j; long **array; - ZALLOC(array, numrows, long *); - for (i = 0; i < numrows; i++) - { - ZALLOC(array[i], numcolumns, long); - if (initval != 0) - lvclear (array[i], initval, numcolumns); - } + + ZALLOC (array, numrows, long *); + for (i = 0; i < numrows; i++) { + ZALLOC (array[i], numcolumns, long); + if (initval != 0) + lvclear (array[i], initval, numcolumns); + } return array; } + void free2Dint (int ***xx, int numrows) { @@ -1445,10 +1593,23 @@ free2Dint (int ***xx, int numrows) int i; array = *xx; - for (i = numrows - 1; i >= 0; i--) - { - free (array[i]); - } + for (i = numrows - 1; i >= 0; i--) { + free (array[i]); + } + free (array); + *xx = NULL; +} + +void +free2Dlong (long ***xx, int numrows) +{ + long **array; + int i; + array = *xx; + + for (i = numrows - 1; i >= 0; i--) { + free (array[i]); + } free (array); *xx = NULL; } @@ -1467,19 +1628,20 @@ free_iarray (int **xx) *xx = NULL; } + double ** initarray_2Ddouble (int numrows, int numcolumns, double initval) { int i, j; double **array; - ZALLOC(array, numrows, double *); - for (i = 0; i < numrows; i++) - { - ZALLOC(array[i], numcolumns, double); - if (initval != 0.0) - vclear (array[i], initval, numcolumns); - } + + ZALLOC (array, numrows, double *); + for (i = 0; i < numrows; i++) { + ZALLOC (array[i], numcolumns, double); + if (initval != 0.0) + vclear (array[i], initval, numcolumns); + } return array; } @@ -1489,22 +1651,21 @@ initarray_2Dlongdouble (int numrows, int numcolumns, long double initval) int i, j; long double **array, *bb; - ZALLOC(array, numrows, long double *); - for (i = 0; i < numrows; i++) - { - ZALLOC(array[i], numcolumns, long double); - if (initval != 0.0) - { - bb = array[i]; - for (j = 0; j < numcolumns; ++j) - { - bb[j] = initval; - } - } + + ZALLOC (array, numrows, long double *); + for (i = 0; i < numrows; i++) { + ZALLOC (array[i], numcolumns, long double); + if (initval != 0.0) { + bb = array[i]; + for (j = 0; j < numcolumns; ++j) { + bb[j] = initval; + } } + } return array; } + void clear2D (double ***xx, int numrows, int numcols, double val) { @@ -1512,10 +1673,9 @@ clear2D (double ***xx, int numrows, int numcols, double val) int i; array = *xx; - for (i = numrows - 1; i >= 0; i--) - { - vclear (array[i], val, numcols); - } + for (i = numrows - 1; i >= 0; i--) { + vclear (array[i], val, numcols); + } } @@ -1527,10 +1687,23 @@ iclear2D (int ***xx, int numrows, int numcols, int val) array = *xx; - for (i = numrows - 1; i >= 0; i--) - { - ivclear (array[i], val, numcols); - } + for (i = numrows - 1; i >= 0; i--) { + ivclear (array[i], val, numcols); + } + +} + +void +lclear2D (long ***xx, int numrows, int numcols, long val) +{ + long **array; + int i; + + array = *xx; + + for (i = numrows - 1; i >= 0; i--) { + lvclear (array[i], val, numcols); + } } @@ -1541,10 +1714,9 @@ free2D (double ***xx, int numrows) int i; array = *xx; - for (i = numrows - 1; i >= 0; i--) - { - free (array[i]); - } + for (i = numrows - 1; i >= 0; i--) { + free (array[i]); + } free (array); *xx = NULL; } @@ -1557,10 +1729,9 @@ free2Dlongdouble (long double ***xx, int numrows) array = *xx; - for (i = numrows - 1; i >= 0; i--) - { - free (array[i]); - } + for (i = numrows - 1; i >= 0; i--) { + free (array[i]); + } free (array); *xx = NULL; } @@ -1569,21 +1740,23 @@ void addoutmul (double *mat, double *v, double mul, int n) { int a, b; - for (a = 0; a < n; ++a) - { - for (b = 0; b < n; ++b) - { - mat[a * n + b] += v[a] * v[b] * mul; - } + for (a = 0; a < n; ++a) { + for (b = 0; b < n; ++b) { + mat[a * n + b] += v[a] * v[b] * mul; } + } } + + + void addouter (double *out, double *a, int n) + /* add outerprod(a) to out trival to recode to make ~ 2 * faster - */ +*/ { addoutmul (out, a, 1.0, n); @@ -1592,10 +1765,11 @@ addouter (double *out, double *a, int n) void subouter (double *out, double *a, int n) + /* subtract outerprod(a) to out trival to recode to make ~ 2 * faster - */ +*/ { addoutmul (out, a, -1.0, n); @@ -1617,6 +1791,7 @@ bal1 (double *a, int n) double logmultinom (int *cc, int n) + /* log multinomial */ { int t, k, i; @@ -1628,31 +1803,31 @@ logmultinom (int *cc, int n) if (t == 0) return 0.0; ytot = 0; - for (i = 0; i < n - 1; i++) - { - k = cc[i]; - y = logbino (t, k); - ytot += y; - t -= k; - } + for (i = 0; i < n - 1; i++) { + k = cc[i]; + y = logbino (t, k); + ytot += y; + t -= k; + } return ytot; } + void flipiarr (int *a, int *b, int n) // reverse array { int *x, k; - ZALLOC(x, n, int); + ZALLOC (x, n, int); - for (k = 0; k < n; ++k) - { - x[n - 1 - k] = b[k]; - } + for (k = 0; k < n; ++k) { + x[n - 1 - k] = b[k]; + } copyiarr (x, a, n); free (x); + } void @@ -1661,24 +1836,24 @@ fliparr (double *a, double *b, int n) double *x; int k; - ZALLOC(x, n, double); + ZALLOC (x, n, double); - for (k = 0; k < n; ++k) - { - x[n - 1 - k] = b[k]; - } + for (k = 0; k < n; ++k) { + x[n - 1 - k] = b[k]; + } copyarr (x, a, n); free (x); } + void vcompl (double *a, double *b, int n) // a <- 1 - b { double *x; - ZALLOC(x, n, double); + ZALLOC (x, n, double); vvm (x, x, b, n); vsp (x, x, 1.0, n); @@ -1687,18 +1862,19 @@ vcompl (double *a, double *b, int n) free (x); } + void setidmat (double *a, int n) // a <- identity matrix { int i; vzero (a, n * n); - for (i = 0; i < n; i++) - { - a[i * n + i] = 1.0; - } + for (i = 0; i < n; i++) { + a[i * n + i] = 1.0; + } } + int stripit (double *a, double *b, int *x, int len) // copy b to a leave out elems where x < 0 @@ -1706,14 +1882,12 @@ stripit (double *a, double *b, int *x, int len) int k, n; n = 0; - for (k = 0; k < len; ++k) - { - if (x[k] >= 0) - { - a[n] = b[k]; - ++n; - } + for (k = 0; k < len; ++k) { + if (x[k] >= 0) { + a[n] = b[k]; + ++n; } + } return n; } @@ -1724,17 +1898,16 @@ istripit (int *a, int *b, int *x, int len) int k, n; n = 0; - for (k = 0; k < len; ++k) - { - if (x[k] >= 0) - { - a[n] = b[k]; - ++n; - } + for (k = 0; k < len; ++k) { + if (x[k] >= 0) { + a[n] = b[k]; + ++n; } + } return n; } + int cstripit (char **a, char **b, int *x, int len) // copy b to a leave out elems where x < 0 @@ -1742,14 +1915,12 @@ cstripit (char **a, char **b, int *x, int len) int k, n; n = 0; - for (k = 0; k < len; ++k) - { - if (x[k] >= 0) - { - a[n] = b[k]; - ++n; - } + for (k = 0; k < len; ++k) { + if (x[k] >= 0) { + a[n] = b[k]; + ++n; } + } return n; } @@ -1759,11 +1930,10 @@ mapit (int *a, int *b, int n, int inval, int outval) int k; copyiarr (b, a, n); - for (k = 0; k < n; ++k) - { - if (a[k] == inval) - a[k] = outval; - } + for (k = 0; k < n; ++k) { + if (a[k] == inval) + a[k] = outval; + } } int @@ -1773,11 +1943,10 @@ ifall (int n, int k) int prod = 1, t = n, j; - for (j = 0; j < k; ++j) - { - prod *= t; - --t; - } + for (j = 0; j < k; ++j) { + prod *= t; + --t; + } return prod; } @@ -1788,6 +1957,7 @@ hlife (double val) return -log (2.0) / log (val); } + void * topheap () // find top of heap (address). Useful for finding memory leaks @@ -1835,20 +2005,22 @@ cswap (char *c1, char *c2) *c1 = *c2; *c2 = cc; + } + int kodeitb (int *xx, int len, int base) { int t = 0, i; - for (i = 0; i < len; ++i) - { - t *= base; - t += xx[i]; - } + for (i = 0; i < len; ++i) { + t *= base; + t += xx[i]; + } return t; } + int dekodeitb (int *xx, int kode, int len, int base) { @@ -1856,25 +2028,34 @@ dekodeitb (int *xx, int kode, int len, int base) int i, t; t = kode; - for (i = len - 1; i >= 0; --i) - { - xx[i] = t % base; - t /= base; - } - return intsum (xx, len); // weight + for (i = len - 1; i >= 0; --i) { + xx[i] = t % base; + t /= base; + } + return intsum (xx, len); // weight } +void +floatit2D (double **a, int **b, int nrows, int ncols) +{ + + int x; + + for (x = 0; x < nrows; ++x) { + floatit (a[x], b[x], ncols); + } +} + void copyarr2D (double **a, double **b, int nrows, int ncols) { int x; - for (x = 0; x < nrows; ++x) - { - copyarr (a[x], b[x], ncols); - } + for (x = 0; x < nrows; ++x) { + copyarr (a[x], b[x], ncols); + } } void @@ -1883,21 +2064,20 @@ copyiarr2D (int **a, int **b, int nrows, int ncols) int x; - for (x = 0; x < nrows; ++x) - { - copyiarr (a[x], b[x], ncols); - } + for (x = 0; x < nrows; ++x) { + copyiarr (a[x], b[x], ncols); + } } + void plus2Dint (int **a, int **b, int **c, int nrows, int ncols) { int x; - for (x = 0; x < nrows; ++x) - { - ivvp (a[x], b[x], c[x], ncols); - } + for (x = 0; x < nrows; ++x) { + ivvp (a[x], b[x], c[x], ncols); + } } void @@ -1905,10 +2085,9 @@ minus2Dint (int **a, int **b, int **c, int nrows, int ncols) { int x; - for (x = 0; x < nrows; ++x) - { - ivvm (a[x], b[x], c[x], ncols); - } + for (x = 0; x < nrows; ++x) { + ivvm (a[x], b[x], c[x], ncols); + } } void @@ -1916,10 +2095,9 @@ plus2D (double **a, double **b, double **c, int nrows, int ncols) { int x; - for (x = 0; x < nrows; ++x) - { - vvp (a[x], b[x], c[x], ncols); - } + for (x = 0; x < nrows; ++x) { + vvp (a[x], b[x], c[x], ncols); + } } void @@ -1927,10 +2105,9 @@ minus2D (double **a, double **b, double **c, int nrows, int ncols) { int x; - for (x = 0; x < nrows; ++x) - { - vvm (a[x], b[x], c[x], ncols); - } + for (x = 0; x < nrows; ++x) { + vvm (a[x], b[x], c[x], ncols); + } } void @@ -1939,21 +2116,20 @@ sum2D (double *a, double **b, int nrows, int ncols) int x; vzero (a, ncols); - for (x = 0; x < nrows; ++x) - { - vvp (a, a, b[x], ncols); - } + for (x = 0; x < nrows; ++x) { + vvp (a, a, b[x], ncols); + } } -int + +double total2D (double **a, int nrows, int ncols) { int x; double sum = 0; - for (x = 0; x < nrows; ++x) - { - sum += asum (a[x], ncols); - } + for (x = 0; x < nrows; ++x) { + sum += asum (a[x], ncols); + } return sum; @@ -1964,31 +2140,64 @@ total2Dint (int **a, int nrows, int ncols) { int x, sum = 0; - for (x = 0; x < nrows; ++x) - { - sum += intsum (a[x], ncols); - } + for (x = 0; x < nrows; ++x) { + sum += intsum (a[x], ncols); + } return sum; } + /** mixed modulus coding (see .../popgen/kimfitdir - */ +*/ +long +lkodeitbb (int *xx, int len, int *baselist) +{ + int i, base; + long t = 0; + + for (i = 0; i < len; ++i) { + base = baselist[i]; + t *= base; + t += xx[i]; + if (t < 0) + fatalx ("(lkodeitbb) overflow\n"); + } + return t; +} + +int +ldekodeitbb (int *xx, long kode, int len, int *baselist) +{ +// return weight + + int i, base; + long t; + + t = kode; + for (i = len - 1; i >= 0; --i) { + base = baselist[i]; + xx[i] = t % base; + t /= base; + } + return intsum (xx, len); + +} + int kodeitbb (int *xx, int len, int *baselist) { int t = 0, i, base; - for (i = 0; i < len; ++i) - { - base = baselist[i]; - t *= base; - t += xx[i]; - if (t < 0) - fatalx ("(kodeitbb) overflow\n"); - } + for (i = 0; i < len; ++i) { + base = baselist[i]; + t *= base; + t += xx[i]; + if (t < 0) + fatalx ("(kodeitbb) overflow\n"); + } return t; } @@ -2000,15 +2209,15 @@ dekodeitbb (int *xx, int kode, int len, int *baselist) int i, t, base; t = kode; - for (i = len - 1; i >= 0; --i) - { - base = baselist[i]; - xx[i] = t % base; - t /= base; - } + for (i = len - 1; i >= 0; --i) { + base = baselist[i]; + xx[i] = t % base; + t /= base; + } return intsum (xx, len); } + long nextprime (long num) // return nextprime >= num @@ -2016,12 +2225,11 @@ nextprime (long num) long x; int t; - for (x = num;; ++x) - { - t = isprime (x); - if (t == YES) - return x; - } + for (x = num;; ++x) { + t = isprime (x); + if (t == YES) + return x; + } } int @@ -2036,16 +2244,16 @@ isprime (long num) return YES; top = nnint (sqrt (num)); - for (x = 2; x <= top; ++x) - { - t = num % x; - if (t == 0) - return NO; - } + for (x = 2; x <= top; ++x) { + t = num % x; + if (t == 0) + return NO; + } return YES; } + int irevcomp (int xx, int stringlen) // consists of stringlen "mininibbles" (2 bits) @@ -2055,19 +2263,18 @@ irevcomp (int xx, int stringlen) if (stringlen > 16) fatalx ("stringlen > 16\n"); xxx = xx; - for (k = 0; k < stringlen; ++k) - { - aa[k] = (xxx & 3) ^ 3; - xxx = xxx >> 2; - } + for (k = 0; k < stringlen; ++k) { + aa[k] = (xxx & 3) ^ 3; + xxx = xxx >> 2; + } xxx = 0; - for (k = 0; k < stringlen; ++k) - { - t = aa[k]; - xxx = (xxx << 2) | t; - } + for (k = 0; k < stringlen; ++k) { + t = aa[k]; + xxx = (xxx << 2) | t; + } return xxx; } + long lrevcomp (long xx, int stringlen) // consists of stringlen "mininibbles" (2 bits) @@ -2079,34 +2286,33 @@ lrevcomp (long xx, int stringlen) if (stringlen > 32) fatalx ("stringlen > 32\n"); xxx = xx; - for (k = 0; k < stringlen; ++k) - { - aa[k] = (xxx & 3) ^ 3; - xxx = xxx >> 2; - } + for (k = 0; k < stringlen; ++k) { + aa[k] = (xxx & 3) ^ 3; + xxx = xxx >> 2; + } xxx = 0; - for (k = 0; k < stringlen; ++k) - { - t = aa[k]; - xxx = (xxx << 2) | t; - } + for (k = 0; k < stringlen; ++k) { + t = aa[k]; + xxx = (xxx << 2) | t; + } return xxx; } + void ismatch (int *a, int *b, int n, int val) { int i; - for (i = 0; i < n; i++) - { - if (b[i] == val) - a[i] = YES; - else - a[i] = NO; + for (i = 0; i < n; i++) { + if (b[i] == val) + a[i] = YES; + else + a[i] = NO; - } + } } + int pmult (double *a, double *b, double *c, int nb, int nc) // polynomial multiplication @@ -2114,15 +2320,13 @@ pmult (double *a, double *b, double *c, int nb, int nc) double *ww; int i, j; - ZALLOC(ww, nb+nc+1, double); + ZALLOC (ww, nb + nc + 1, double); - for (i = 0; i <= nb; ++i) - { - for (j = 0; j <= nc; ++j) - { - ww[i + j] += b[i] * c[j]; - } + for (i = 0; i <= nb; ++i) { + for (j = 0; j <= nc; ++j) { + ww[i + j] += b[i] * c[j]; } + } copyarr (ww, a, nb + nc + 1); free (ww); @@ -2138,15 +2342,65 @@ pdiff (double *a, double *b, int deg) double *ww, y; int k; - ZALLOC(ww, deg+1, double); - for (k = 1; k <= deg; ++k) - { - y = (double) k; - ww[k - 1] = y * b[k]; - } + ZALLOC (ww, deg + 1, double); + for (k = 1; k <= deg; ++k) { + y = (double) k; + ww[k - 1] = y * b[k]; + } copyarr (ww, a, deg + 1); free (ww); +} + +int +mktriang (double *out, double *in, int n) +{ + int a, b, x = 0; + for (a = 0; a < n; ++a) { + for (b = a; b < n; ++b) { + out[x] = in[a * n + b]; + ++x; + } + } + return x; +} + +int +mkfull (double *out, double *in, int n) +// inverse to mktriang +{ + int a, b, x = 0; + for (a = 0; a < n; ++a) { + for (b = a; b < n; ++b) { + out[a * n + b] = out[b * n + a] = in[x];; + ++x; + } + } + return x; +} +void vswap(double *a, double *b, int n) +{ + double *w ; + ZALLOC(w, n, double) ; + + copyarr(a, w, n) ; + copyarr(b, a, n) ; + copyarr(w, b, n) ; + + free(w) ; } +void setlong(long *pplen, long a, long b) +// *pplen is a*b with check for overflow +{ + long long int xx ; + + xx = a*b ; + if (xx > LONG_MAX) fatalx("overflow: Are you on a 32 bit machine?\n") ; + *pplen = xx ; + + +} + +