[med-svn] [eigensoft] 01/06: New upstream version 6.1.4+dfsg
Andreas Tille
tille at debian.org
Sun Oct 16 18:07:55 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository eigensoft.
commit aa59770db532929fa36b53ff25e3982776795e33
Author: Andreas Tille <tille at debian.org>
Date: Sun Oct 16 15:02:38 2016 +0200
New upstream version 6.1.4+dfsg
---
POPGEN/README | 2 +-
README | 17 +-
bin/ploteig | 209 +++++++
include/egsubs.h | 22 +-
include/strsubs.h | 213 +++----
include/vsubs.h | 507 ++++++---------
src/Makefile | 4 +-
src/egsubs.c | 182 +++---
src/nicksrc/strsubs.c | 1515 +++++++++++++++++++++++++--------------------
src/nicksrc/vsubs.c | 1642 ++++++++++++++++++++++++++++---------------------
src/pcaselection.c | 4 +-
11 files changed, 2366 insertions(+), 1951 deletions(-)
diff --git a/POPGEN/README b/POPGEN/README
index 64656a4..27be164 100644
--- a/POPGEN/README
+++ b/POPGEN/README
@@ -341,7 +341,7 @@ relthresh: threshhold for correlation coefficients. Only coefficients
------------------------------------------------------------------------------
Questions?
-See http://www.hsph.harvard.edu/faculty/alkes-price/files/eigensoftfaq.htm
+See https://www.hsph.harvard.edu/alkes-price/eigensoft-frequently-asked-questions/
or email Samuela Pollack, spollack at hsph.harvard.edu
SOFTWARE COPYRIGHT NOTICE AGREEMENT
diff --git a/README b/README
index 2ffa017..2fed8bb 100644
--- a/README
+++ b/README
@@ -1,12 +1,21 @@
-EIGENSOFT version 6.1.3, 7/29/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)
Price et al. 2006 Nat Genet (EIGENSTRAT stratification correction)
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
+-- 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
+-- Added updated ploteig script that disappeared from the repository
NEW features of EIGENSOFT version 6.1.2 include:
-- Updated license info to be GPL compliant required by linking the GSL
@@ -43,7 +52,7 @@ See POPGEN/README for documentation of population structure programs.
See EIGENSTRAT/README for documentation of EIGENSTRAT programs.
Questions?
-See http://www.hsph.harvard.edu/faculty/alkes-price/files/eigensoftfaq.htm
+See https://www.hsph.harvard.edu/alkes-price/eigensoft-frequently-asked-questions/
or email Samuela Pollack, spollack at hsph.harvard.edu
Executables and source code:
@@ -76,13 +85,15 @@ To remake the entire package:
"make clobber"
"make install"
-To remake EIG6.0 it is necessary to link to the OpenBLAS library. On orchestra,
+To remake EIG6.x it is necessary to link to the OpenBLAS library. On orchestra,
the path is /opt/openblas and should work automatically. On Broad institute machines,
the user should execute "use .openblas-0.2.8" and "use GCC-4.9" at the command
prompt before attempting to remake. All other users should install OpenBLAS and
set the variable OPENBLAS to the path at the make command line,
e.g. "make install OPENBLAS=/usr/local/openblas"
+The ploteig utility requires gnuplot to run.
+
----------------------------
Acknowledgements:
EIGENSOFT was written by Nick Patterson, Alkes Price, Samuela Pollack,
diff --git a/bin/ploteig b/bin/ploteig
new file mode 100755
index 0000000..d611b72
--- /dev/null
+++ b/bin/ploteig
@@ -0,0 +1,209 @@
+#!/usr/local/bin/perl -w
+
+### ploteig -i eigfile -p pops -c a:b [-t title] [-s stem] [-o outfile] [-x] [-k] [-y] [-z sep] -r colorstring -m xmul -n ymul
+use Getopt::Std ;
+use File::Basename ;
+
+## pops : separated -x = make postscript and pdf -z use another separator
+## -k keep intermediate files
+## NEW if pops is a file names are read one per line
+## scaling on x, y axes
+
+getopts('i:o:p:c:s:d:z:t:r:m:n:xkyq',\%opts) ;
+$postscmode = $opts{"x"} ;
+$oldkeystyle = $opts{"y"} ;
+$nopoptitle = $opts{"q"} ;
+$kflag = $opts{"k"} ;
+$keepflag = 1 if ($kflag) ;
+$keepflag = 1 unless ($postscmode) ;
+
+$zsep = ":" ;
+if (defined $opts{"z"}) {
+ $zsep = $opts{"z"} ;
+ $zsep = "\+" if ($zsep eq "+") ;
+}
+
+if (defined $opts{"r"}) {
+ $colorstr = $opts{"r"} ;
+ setcolor($colorstr) ;
+}
+$xmul = $opts{"m"} ;
+$xmul = 1 unless (defined $xmul) ;
+
+$ymul = $opts{"n"} ;
+$ymul = 1 unless (defined $ymul) ;
+
+$title = "" ;
+if (defined $opts{"t"}) {
+ $title = $opts{"t"} ;
+}
+if (defined $opts{"i"}) {
+ $infile = $opts{"i"} ;
+}
+else {
+ usage() ;
+ exit 0 ;
+}
+open (FF, $infile) || die "can't open $infile\n" ;
+ at L = (<FF>) ;
+chomp @L ;
+$nf = 0 ;
+foreach $line (@L) {
+ next if ($line =~ /\#/) ;
+ @Z = split " ", $line ;
+ $x = @Z ;
+ $nf = $x if ($nf < $x) ;
+}
+printf "## number of fields: %d\n", $nf ;
+$popcol = $nf-1 ;
+
+
+if (defined $opts{"p"}) {
+ $pops = $opts{"p"} ;
+}
+else {
+ die "p parameter compulsory\n" ;
+}
+
+$popsname = setpops ($pops) ;
+print "$popsname\n" ;
+
+$c1 = 1; $c2 =2 ;
+if (defined $opts{"c"}) {
+ $cols = $opts{"c"} ;
+ ($c1, $c2) = split ":", $cols ;
+ die "bad c param: $cols\n" unless (defined $cols) ;
+}
+
+$stem = "$infile.$c1:$c2" ;
+if (defined $opts{"s"}) {
+ $stem = $opts{"s"} ;
+}
+$gnfile = "$stem.$popsname.xtxt" ;
+
+if (defined $opts{"o"}) {
+ $gnfile = $opts{"o"} ;
+}
+
+
+ at T = () ; ## trash
+open (GG, ">$gnfile") || die "can't open $gnfile\n" ;
+print GG "## " unless ($postscmode) ;
+print GG "set terminal postscript color noenhanced\n" ;
+print GG "set title \"$title\" \n" ;
+print GG "set key outside\n" unless ($oldkeystyle) ;
+print GG "set xlabel \"eigenvector $c1\" \n" if ($xmul == 1) ;
+print GG "set xlabel \"eigenvector $c1 (* $xmul) \" \n" if ($xmul != 1) ;
+print GG "set ylabel \"eigenvector $c2\" \n" if ($ymul == 1) ;
+print GG "set ylabel \"eigenvector $c1 (* $ymul) \" \n" if ($ymul != 1) ;
+print GG "plot " ;
+$np = @P ;
+$lastpop = $P[$np-1] ;
+$d1 = $c1+1 ;
+$d2 = $c2+1 ;
+foreach $pop (@P) {
+ $dfile = "$stem:$pop" ;
+ push @T, $dfile ;
+ print GG " \"$dfile\" using (\$$d1)*$xmul:(\$$d2)*$ymul " ;
+ print GG "notitle " if (defined $nopoptitle) ;
+ print GG "title \"$pop\" " unless (defined $nopoptitle) ;
+ if (defined $COL{$pop}) {
+ $color = $COL{$pop} ;
+ print GG "lt rgb \"$color\" " ;
+ }
+ print GG ", \\\n" unless ($pop eq $lastpop) ;
+ open (YY, ">$dfile") || die "can't open $dfile\n" ;
+ foreach $line (@L) {
+ next if ($line =~ /\#/) ;
+ @Z = split " ", $line ;
+ next unless (defined $Z[$popcol]) ;
+ next unless ($Z[$popcol] eq $pop) ;
+ print YY "$line\n" ;
+ }
+ close YY ;
+}
+print GG "\n" ;
+print GG "## " if ($postscmode) ;
+print GG "pause 9999\n" ;
+close GG ;
+
+if ($postscmode) {
+$psfile = "$stem.ps" ;
+
+ if ($gnfile =~ /xtxt/) {
+ $psfile = $gnfile ;
+ $psfile =~ s/xtxt/ps/ ;
+ }
+system "gnuplot < $gnfile > $psfile" ;
+system "/home/np29/bin/fixgreen $psfile" ;
+system "ps2pdf $psfile " ;
+}
+unlink (@T) unless $keepflag ;
+
+sub setcolor {
+ my ($colorst) = @_ ;
+ local ($cp, $pop, $color, @CP, $line) ;
+ if (-r $colorst) {
+ open (C1, $colorst) || die "can't open $colorst\n" ;
+ foreach $line (<C1>) {
+ chomp $line ;
+ ($pop, $color) = split " ", $line ;
+ next if ($pop =~ /\#/) ;
+ next unless (defined $color) ;
+ print STDERR "setting color for $pop to $color\n" ;
+ $COL{$pop} = $color ;
+ }
+ close C1 ;
+ return ;
+ }
+
+ @CP = split " ", $colorst ;
+ foreach $cp (@CP) {
+ ($pop, $color) = split ":", $cp ;
+ $COL{$pop} = $color ;
+ }
+}
+
+sub usage {
+
+print "ploteig -i eigfile -p pops -c a:b [-t title] [-s stem] [-o outfile] [-x] [-k] -c colorstringh [-m xmul] [-n ymul]\n" ;
+print "-i eigfile input file first col indiv-id last col population\n" ;
+print "## as output by smartpca in outputvecs \n" ;
+print "-c a:b a, b columns to plot. 1:2 would be common and leading 2 eigenvectors\n" ;
+print "-p pops Populations to plot. : delimited. eg -p Bantu:San:French\n" ;
+print "## pops can also be a filename. List populations 1 per line\n" ;
+print "[-s stem] stem will start various output files\n" ;
+print "[-o ofile] ofile will be gnuplot control file. Should have xtxt suffix\n";
+print "[-x] make ps and pdf files\n" ;
+print "[-k] keep various intermediate files although -x set\n" ;
+print "## necessary if .xtxt file is to be hand edited\n" ;
+print "[-r colorstringpairs or colorstringfile]\n" ;
+print "[-y] put key at top right inside box (old mode)\n" ;
+print "[-t] title (legend)\n" ;
+
+print "The xtxt file is a gnuplot file and can be easily hand edited. Intermediate files
+needed if you want to make your own plot\n" ;
+
+}
+sub setpops {
+ my ($pops) = @_ ;
+ local (@a, $d, $b, $e) ;
+
+ if (-e $pops) {
+ open (FF1, $pops) || die "can't open $pops\n" ;
+ @P = () ;
+ foreach $line (<FF1>) {
+ ($a) = split " ", $line ;
+ next unless (defined $a) ;
+ next if ($a =~ /\#/) ;
+ push @P, $a ;
+ }
+ $out = join ":", @P ;
+ print "## pops: $out\n" ;
+ ($b, $d , $e) = fileparse($pops) ;
+ return $b ;
+ }
+ @P = split $zsep, $pops ;
+ return $pops ;
+
+}
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 <stdlib.h>
-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 <math.h>
#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/Makefile b/src/Makefile
index 319e8dc..18cda52 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,5 +1,5 @@
CFLAGS += -I../include
-LDLIBS += -lgsl -lopenblas -lgfortran -lrt
+LDLIBS += -lgfortran -lrt -Wl,-Bstatic -lgsl -lopenblas -Wl,-Bdynamic -fopenmp
ifeq ($(OPTIMIZE), 1)
CFLAGS += -O2
@@ -22,7 +22,7 @@ NLIB = $(ND)/libnick.a
# ----- phony targets
.PHONY: all clean clobber install
-EXE = baseprog convertf mergeit pca \
+EXE = baseprog convertf mergeit pca pcaselection \
$(ED)/pcatoy $(ED)/smartrel $(ED)/smarteigenstrat \
$(ED)/twstats $(ED)/eigenstrat $(ED)/eigenstratQTL $(ED)/smartpca
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 <errno.h>
#include <sys/types.h>
#include <sys/stat.h>
+#include <xsearch.h>
+
#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,36 +725,96 @@ 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 <stdio.h>
#include <string.h>
#include <unistd.h>
+#include <limits.h>
#include <math.h>
-#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 ;
+
+
+}
+
+
diff --git a/src/pcaselection.c b/src/pcaselection.c
index 1bc0a7a..a23a083 100644
--- a/src/pcaselection.c
+++ b/src/pcaselection.c
@@ -114,8 +114,8 @@ main (int argc, char **argv)
{
size_t i, j, k;
- double *vg = (double *) malloc (k * sizeof(double));
- double *v1 = (double *) malloc (k * sizeof(double));
+ double *vg = (double *) malloc (K * sizeof(double));
+ double *v1 = (double *) malloc (K * sizeof(double));
for (i = 0; i < numsnps; i++)
{
N = 0;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/eigensoft.git
More information about the debian-med-commit
mailing list