# Date: Thu, 11 Jan 90 11:14:59 PST # From: fessler@isl.Stanford.EDU (Jeffrey A. Fessler) # To: ehg@research.att.com # # This is a bundled file containg the VSPLINE library. # (Non-parametric fixed-interval smoothing with vector-splines.) # To unbundle, remove everything above this line, then sh this file mkdir Doc Test matrix smooth Data numrecC echo License 1>&2 cat >License <<'End of License' STANFORD UNIVERSITY SOFTWARE LICENSE AGREEMENT BY ACCESSING THE SOFTWARE VSPLINE, YOU ("LICENSEE") ARE AGREEING TO BECOME BOUND BY THE TERMS OF THIS AGREEMENT. 1. STANFORD grants to LICENSEE a nonexclusive and nontransferable license to use the software VSPLINE furnished hereunder, upon the terms and conditions set out below. 2. LICENSEE acknowledges that VSPLINE is a research tool still in the development stage, that it is being supplied "as is," without any accompanying services or improvements from STANFORD, and that this license is entered into in order to encourage scientific collaboration aimed at further development and application of the Program. 3. STANFORD MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. By way of example, but not limitation, STANFORD MAKES NO REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS. STANFORD shall not be held liable for any liability nor for any direct, indirect or consequential damages with respect to any claim by LICENSEE or any third party on account of or arising from this Agreement or use of VPSLINE. 4. "LICENSEE agrees not to sell or sublicense VSPLINE to another location or to any other person without prior written permission from STANFORD." 5. Title to copyright to VSPLINE and to any associated documentation shall at all times remain with STANFORD, and LICENSEE agrees to preserve same. LICENSEE agrees to maintain the copyright notice and the license agreement in any copies LICENSEE makes of VSPLINE. 6. This Agreement shall be construed, interpreted and applied in accordance laws of the State of California. 7. Nothing in this Agreement shall be construed as conferring rights to use advertising, publicity or otherwise any trademark or the name of "STANFORD." End of License echo README 1>&2 cat >README <<'End of README' INSTALLATION INSTRUCTIONS FOR VSPLINE LIBRARY To install this code, # Edit Vlist and List to give the proper path names. # # Change defs.h if necessary for your compiler. # e.g., on a sun 4, add "#include ". # # Make sure that you have installed the # Numerical Recipes in C library, and have compiled # it with compatible floating point options : # Change the FLOATFLAGS definition in List to your # favorite option (or comment it out). # # Type make. # # To use, simply include libvspline.a and libnumrecC.a # when linking. # # For documentation, there is a Latex file in Doc/ # # The numrecC library is only needed for finding the # minimum cross-validation score. If you want to # specify the smoothing parameter yourself, then type # "make libvsmooth.a" to get a reduced library. # # This software can also be used with MATLAB mex files # by adding the flag -DMATLAB to CFLAGS in matrix/Makefile. End of README echo List 1>&2 cat >List <<'End of List' # master makefile definitions for vspline library #BIN=$(HOME)/bin VSPLINE=$(HOME)/recon/global/vspline DEFDIR=$(VSPLINE) BIN=$(VSPLINE)/Data MYNR=$(VSPLINE)/numrecC MATRIX=$(VSPLINE)/matrix SMOOTH=$(VSPLINE)/smooth NUMRECDIR=/usr/local/lib PREC=Double FLOATFLAGS=-fswitch -fsingle End of List echo Vlist 1>&2 cat >Vlist <<'End of Vlist' include /user/fessler/recon/global/vspline/numrecC/List include /user/fessler/recon/global/vspline/matrix/List include /user/fessler/recon/global/vspline/smooth/List End of Vlist echo Makefile 1>&2 cat >Makefile <<'End of Makefile' # Master Makefile for VSPLINE library all: \ makeall \ libvspline.a include List DIRS= \ matrix smooth DIR2= $(DIRS) Test # numrecC makeall: for i in $(DIRS) ; do echo $$i ; cd $(VSPLINE)/$$i ; make ; done LL: for i in $(DIR2) ; do echo $$i ; cd $(VSPLINE)/$$i ; make LL ; done .PRECIOUS: .print .print: for i in $(DIR2) ; do (echo .PRECIOUS: .print ; echo $(VSPLINE)/$$i/.print: $(VSPLINE)/$$i/*.[ch]; echo ' enscript -f Courier7 -2rG `/user/fessler/bin/sh/linesplit $$? | sort` ; touch $$@') | make -f - ; done .check: for i in $(DIR2) ; do (echo $(VSPLINE)/$$i/.print: $(VSPLINE)/$$i/*.[ch]; echo ' THESE $$?') | make -n -f - ; done include Vlist #$(MYNRO) libvspline.a: $(IODEPO) $(MATRIXO) $(SMOOTHEXO) ar ruv $@ $? ; ranlib $@ libvsmooth.a: $(IODEPO) $(MATRIXO) $(SCOREO) $(SMOOTH)/example1.o ar ruv $@ $? ; ranlib $@ SOURCE= \ matrix/*.c \ smooth/*.c # numrecC/*.c tags: $(SOURCE) ctags $(SOURCE) clean: rm */*.o End of Makefile echo defs.h 1>&2 cat >defs.h <<'End of defs.h' /***************************** defs.h ************************************ * Definitions that may need to be changed * to make the code more transportable. *************************************************************************/ #define Void char /* for (void *) declarations */ extern char *calloc(); extern char *sprintf(); extern double atof(); extern char *alloca(); #define TMP_ALLOC(n,s) alloca((int) ((n)*(s))) #define TMP_FREE(p) { } #define BZERO(p, n) bzero((char *) (p), (int) ((n) * sizeof(*p))); #define BCOPY(s, d, n) bcopy((char *) (s), (char *) (d), (int) ((n) * sizeof(*s))); extern char *mycalloc(), *alloc_check(); extern void myfree(); #define ALLOC_CHECK(n,s) alloc_check(n,s) #define ALLOC(n,s) mycalloc(n,s) #define FREE(p) myfree((char *) (p)) #define ABORT jumpship #define EXITBUG(msg) { char str[512]; \ sprintf(str, "ERROR: %s in %s line %d\n", msg, __FILE__, __LINE__); \ ABORT(str); } #define FOPEN(fp,n,t) {if ( !(fp = fopen(n, t)) ) EXITBUG("fopen") } #define FREAD(p,s,n,fp) {if ((n) != fread((char *) (p), s, n, fp)) EXITBUG("fread") } #define FWRITE(p,s,n,fp) {if ((n) != fwrite((char *) (p), s, n, fp)) EXITBUG("fwrite") } #define FFLUSH(fp) {if (fflush(fp)) EXITBUG("fflush") } #define FCLOSE(fp) {if (fclose(fp)) EXITBUG("fclose") } End of defs.h echo Doc/README 1>&2 cat >Doc/README <<'End of Doc/README' Don't remove root.bb* or you'll lose the bibliography! End of Doc/README echo Doc/defs_math.tex 1>&2 cat >Doc/defs_math.tex <<'End of Doc/defs_math.tex' % special definitions of math symbols % DONT FORGET : CANNOT USE NUMERALS IN COMMAND NAMES! \newcommand{\BA}{\begin{array}} \newcommand{\EA}{\end{array}} \newcommand{\BE}{\begin{equation}} \newcommand{\EE}{\end{equation}} \newcommand{\BM}[1]{ \left[ \begin{array}{#1}} \newcommand{\EM}{ \end{array} \right] } \newcommand{\cA}{ {\cal A}} \newcommand{\cB}{ {\cal B}} \newcommand{\cC}{ {\cal C}} \newcommand{\cD}{ {\cal D}} \newcommand{\cE}{ {\cal E}} \newcommand{\cF}{ {\cal F}} \newcommand{\cG}{ {\cal G}} \newcommand{\cH}{ {\cal H}} \newcommand{\cI}{ {\cal I}} \newcommand{\cJ}{ {\cal J}} \newcommand{\cK}{ {\cal K}} \newcommand{\cL}{ {\cal L}} \newcommand{\cM}{ {\cal M}} \newcommand{\cN}{ {\cal N}} \newcommand{\cO}{ {\cal O}} \newcommand{\cP}{ {\cal P}} \newcommand{\cQ}{ {\cal Q}} \newcommand{\cR}{ {\cal R}} \newcommand{\cS}{ {\cal S}} \newcommand{\cT}{ {\cal T}} \newcommand{\cU}{ {\cal U}} \newcommand{\cV}{ {\cal V}} \newcommand{\cW}{ {\cal W}} \newcommand{\cX}{ {\cal X}} \newcommand{\cY}{ {\cal Y}} \newcommand{\cZ}{ {\cal Z}} \newcommand{\bcY}{{\boldmath ${\cal Y}$}} \newcommand{\boldzero}{{\bf 0}} \newcommand{\bolda}{{\bf a}} \newcommand{\boldb}{{\bf b}} \newcommand{\boldc}{{\bf c}} \newcommand{\boldd}{{\bf d}} \newcommand{\bolde}{{\bf e}} \newcommand{\boldf}{{\bf f}} \newcommand{\boldg}{{\bf g}} \newcommand{\boldh}{{\bf h}} \newcommand{\boldi}{{\bf i}} \newcommand{\boldj}{{\bf j}} \newcommand{\boldk}{{\bf k}} \newcommand{\boldl}{{\bf l}} \newcommand{\boldm}{{\bf m}} \newcommand{\boldn}{{\bf n}} \newcommand{\boldo}{{\bf o}} \newcommand{\boldp}{{\bf p}} \newcommand{\boldq}{{\bf q}} \newcommand{\boldr}{{\bf r}} \newcommand{\bolds}{{\bf s}} \newcommand{\boldt}{{\bf t}} \newcommand{\boldu}{{\bf u}} \newcommand{\boldv}{{\bf v}} \newcommand{\boldw}{{\bf w}} \newcommand{\boldx}{{\bf x}} \newcommand{\boldy}{{\bf y}} \newcommand{\boldz}{{\bf z}} \newcommand{\boldA}{{\bf A}} \newcommand{\boldB}{{\bf B}} \newcommand{\boldC}{{\bf C}} \newcommand{\boldD}{{\bf D}} \newcommand{\boldE}{{\bf E}} \newcommand{\boldF}{{\bf F}} \newcommand{\boldG}{{\bf G}} \newcommand{\boldH}{{\bf H}} \newcommand{\boldI}{{\bf I}} \newcommand{\boldJ}{{\bf H}} \newcommand{\boldK}{{\bf K}} \newcommand{\boldL}{{\bf L}} \newcommand{\boldM}{{\bf M}} \newcommand{\boldN}{{\bf N}} \newcommand{\boldO}{{\bf O}} \newcommand{\boldP}{{\bf P}} \newcommand{\boldQ}{{\bf Q}} \newcommand{\boldR}{{\bf R}} \newcommand{\boldS}{{\bf S}} \newcommand{\boldT}{{\bf T}} \newcommand{\boldU}{{\bf U}} \newcommand{\boldV}{{\bf V}} \newcommand{\boldW}{{\bf W}} \newcommand{\boldX}{{\bf X}} \newcommand{\boldY}{{\bf Y}} \newcommand{\boldZ}{{\bf Z}} % lower case greek are "symbols" so \bf wont work \newcommand{\balph} {\mbox{\boldmath $\alpha$}} \newcommand{\balpha} {\mbox{\boldmath $\alpha$}} \newcommand{\bbeta} {\mbox{\boldmath $\beta$}} \newcommand{\bchi} {\mbox{\boldmath $\chi$}} \newcommand{\bdelt} {\mbox{\boldmath $\delta$}} \newcommand{\bdelta} {\mbox{\boldmath $\delta$}} \newcommand{\beps} {\mbox{\boldmath $\epsilon$}} \newcommand{\bveps} {\mbox{\boldmath $\varepsilon$}} \newcommand{\bgamm} {\mbox{\boldmath $\gamma$}} \newcommand{\bgamma} {\mbox{\boldmath $\gamma$}} \newcommand{\blamb} {\mbox{\boldmath $\lambda$}} \newcommand{\blambda} {\mbox{\boldmath $\lambda$}} \newcommand{\bmu} {\mbox{\boldmath $\mu$}} \newcommand{\bpsi} {\mbox{\boldmath $\psi$}} \newcommand{\btau} {\mbox{\boldmath $\tau$}} \newcommand{\bthet} {\mbox{\boldmath $\theta$}} \newcommand{\btheta} {\mbox{\boldmath $\theta$}} \newcommand{\brho} {\mbox{\boldmath $\rho$}} \newcommand{\bzeta} {\mbox{\boldmath $\zeta$}} \newcommand{\bGamm} {{\bf \Gamma}} \newcommand{\bLamb} {{\bf \Lambda}} \newcommand{\bLambda} {{\bf \Lambda}} \newcommand{\bPi} {{\bf \Pi}} \newcommand{\bSig} {{\bf \Sigma}} \newcommand{\bSigma} {{\bf \Sigma}} \newcommand{\eps}{\epsilon} \newcommand{\veps}{\varepsilon} \newcommand{\LI}{\mbox{\boldmath \langle}} \newcommand{\RI}{\mbox{\boldmath \rangle}} \newcommand{\defequ}{ \ \stackrel{\bigtriangleup}{=} \ } \newcommand{\iid}{\stackrel{iid}{\sim}} \newcommand{\ul}[1]{\underline{ #1 }} \newcommand{\degr}[1]{#1^{\circ}} \newcommand{\hathat}[1]{ \hat{ \hat{ \mbox{ $ #1 $ } } } } \newcommand{\myvect}[2]{\left[ \begin{array}{c} #1 \\ \vdots \\ #2 \end{array} \right]} End of Doc/defs_math.tex echo Doc/local_defs.tex 1>&2 cat >Doc/local_defs.tex <<'End of Doc/local_defs.tex' \newcommand{\Ad}[1]{ \boldA_{(#1)}(\balpha) } \newcommand{\ghata}{ \hat{g}_{\alpha} } \newcommand{\bghata}{ \hat{\boldg}_{\balpha} } \newcommand{\bghatan}{ \hat{\boldg}_{\balpha,-n} } \newcommand{\hatcalph}{ \hat{\boldc}_{\balpha} } \newcommand{\calph}{ \boldc_{\balpha} } \newcommand{\RE}{\eta} \newcommand{\CV}{ {\rm CV} } \newcommand{\RUF}{ {\rm R} } \newcommand{\diag}{ {\rm diag} } \newcommand{\trace}{ {\rm trace} } \newcommand{\meas}{\boldy} \newcommand{\Meas}{\boldy} \newcommand{\mcoviid}{\ul{\bSigma}} \newcommand{\mcov}{\bSigma} End of Doc/local_defs.tex echo Doc/root.aux 1>&2 cat >Doc/root.aux <<'End of Doc/root.aux' \relax \citation{netlib:87} \citation{Fessler:89assp} \@writefile{toc}{\string\contentsline\space {section}{\string\numberline\space {\uppercase {i}}Introduction}{1}} \@writefile{toc}{\string\contentsline\space {section}{\string\numberline\space {\uppercase {ii}}Spline smoothing of vector measurements}{1}} \newlabel{sec:vector}{{\uppercase {ii}}{1}} \newlabel{eq:multi-model}{{1}{1}} \citation{Fessler:89assp} \citation{Fessler:89assp} \citation{HutchHoog:85} \citation{Fessler:89assp} \citation{numrecc:88} \newlabel{eq:multi-min}{{2}{2}} \@writefile{toc}{\string\contentsline\space {section}{\string\numberline\space {\uppercase {iii}}Choosing the smoothing parameters - Cross Validation}{2}} \newlabel{sec:auto}{{\uppercase {iii}}{2}} \newlabel{eq:cv-score}{{3}{2}} \newlabel{eq:cv-score2}{{4}{2}} \@writefile{toc}{\string\contentsline\space {section}{\string\numberline\space {\uppercase {iv}}Software Interface}{3}} \@writefile{toc}{\string\contentsline\space {subsection}{\string\numberline\space {A}Conventions}{3}} \@writefile{toc}{\string\contentsline\space {subsection}{\string\numberline\space {B}Subroutine {\string\ptt\space vspline}}{3}} \@writefile{toc}{\string\contentsline\space {subsection}{\string\numberline\space {C}Subroutine {\string\ptt\space bestcv}}{3}} \@writefile{toc}{\string\contentsline\space {subsection}{\string\numberline\space {D}Subroutine {\string\ptt\space cubic\unhbox \voidb@x \kern .06em \vbox {\hrule width.3em}interp}}{3}} \citation{numrecc:88} \bibstyle{ieeetr} \bibdata{preamble,abbrev,kalman,misc,mybib,spline} \bibcite{netlib:87}{1} \bibcite{Fessler:89assp}{2} \bibcite{numrecc:88}{3} \bibcite{HutchHoog:85}{4} \@writefile{toc}{\string\contentsline\space {section}{\string\numberline\space {\uppercase {v}}Installation}{4}} End of Doc/root.aux echo Doc/root.bbl 1>&2 cat >Doc/root.bbl <<'End of Doc/root.bbl' \begin{thebibliography}{1} \bibitem{netlib:87} J.~J.~Dongarra and E.~Grosse, ``Distibution of mathematical software via electronic mail,'' {\it Comm. ACM}, vol.~30, pp.~403--407, Oct. 1987. \bibitem{Fessler:89assp} J.~Fessler, ``Non-parametric fixed-interval smoothing with vector splines,'' {\it IEEE Transactions on Acoustics Speech and Signal Processing}. \newblock vol.~39:4, pp.~852--859, Apr. 1991. \bibitem{numrecc:88} W.~H.~Press, B.~P.~Flannery, S.~A.~Teukolsky, and W.~T.~Vetterling, {\it Numerical Recipes in {C}}. \newblock Cambridge Univ. Press, 1988. \bibitem{HutchHoog:85} M.~Hutchinson and F.~deHoog, ``Smoothing noisy data with spline functions,'' {\it Numerische Mathematik}, vol.~47, pp.~99--106, 1985. \end{thebibliography} End of Doc/root.bbl echo Doc/root.blg 1>&2 cat >Doc/root.blg <<'End of Doc/root.blg' This is BibTeX, Version 0.98i for Berkeley UNIX The top-level auxiliary file: root.aux The style file: ieeetr.bst Database file #1: preamble.bib Database file #2: abbrev.bib Database file #3: kalman.bib Database file #4: misc.bib Database file #5: mybib.bib Database file #6: spline.bib : Warning: the pages shouldn't be empty in Fessler:89assp : Warning: the year shouldn't be empty in Fessler:89assp End of Doc/root.blg echo Doc/root.log 1>&2 cat >Doc/root.log <<'End of Doc/root.log' This is TeX, C Version 2.991 (preloaded format=lplain 89.10.13) 28 NOV 1989 16:05 **&lplain root (root.tex LaTeX Version 2.09 <24 May 1989> (/usr/local/lib/tex82/macros/article.sty Document Style `article' <16 Mar 88>. (/usr/local/lib/tex82/macros/art11.sty) \c@part=\count78 \c@section=\count79 \c@subsection=\count80 \c@subsubsection=\count81 \c@paragraph=\count82 \c@subparagraph=\count83 \c@figure=\count84 \c@table=\count85 ) (size.tex) (root.aux) (defs_math.tex) (local_defs.tex) (s0.tex) (s1.tex) (s2.tex [1 ]) (s3.tex) (s4.tex [2] Overfull \hbox (13.15578pt too wide) in paragraph at lines 33--52 \elvtt (double *) NULL\elvrm . The mea-sure-ments are stored in \elvtt yn \elvr m in the or-der $\elvsy f\elvmi y[]; [] ; y[]; [] ; y[]; [] ; y[]\elvsy g$\elvr m . \hbox(8.2125+3.25056)x469.75499, glue set - 1.0 .\elvtt ( .\elvtt d .\elvtt o .\elvtt u .\elvtt b .etc. ) (s5.tex [3]) (root.bbl) [4] (root.aux) Here is how much of TeX's memory you used: 359 strings out of 4441 2500 string characters out of 62996 33646 words of memory out of 262141 2318 multiletter control sequences out of 9500 43015 words of font info for 157 fonts, out of 72000 for 255 14 hyphenation exceptions out of 607 12i,12n,15p,222b,266s stack positions out of 300i,40n,60p,3000b,4000s Output written on root.dvi (4 pages, 14464 bytes). End of Doc/root.log echo Doc/root.tex 1>&2 cat >Doc/root.tex <<'End of Doc/root.tex' % documentation for vspline code %\batchmode \documentstyle[11pt]{article} \input{size.tex} \renewcommand{\baselinestretch}{1.1} % double spacing \renewcommand{\thesection}{\Roman{section}} \renewcommand{\thesubsection}{\Alph{subsection}} \begin{document} \input{defs_math.tex} \input{local_defs.tex} \renewcommand{\baselinestretch}{1.0} % double spacing \input{s0.tex} \renewcommand{\baselinestretch}{1.1} % double spacing \input{s1.tex} \input{s2.tex} \input{s3.tex} \input{s4.tex} \input{s5.tex} \bibliographystyle{ieeetr} \bibliography{preamble,abbrev,kalman,misc,mybib,spline} \end{document} End of Doc/root.tex echo Doc/s0.tex 1>&2 cat >Doc/s0.tex <<'End of Doc/s0.tex' \title{VSPLINE: a subroutine library for non-parametric fixed-interval smoothing with vector splines} \author{Jeffrey A. Fessler% \thanks{ This work was supported in part by National Institute of Health contract NO1-HV-38045 and grant R01-HL-39045, National Science Foundation contract ECS-8213959, and GE Medical Systems Group contract 22-84. }\\ Information Systems Laboratory \\ Department of Electrical Engineering \\ Stanford University $\cdot$ Stanford, CA 94305 \\ E-mail: fessler@isl.stanford.edu \\ (415) 723-1904 %\\ \vspace*{0.5in} \\ VSPLINE Software and Documentation \\ Copyright \copyright 1989 \\ The Board of Trustees of the Leland Stanford Junior University.\\ All Rights Reserved. } \date{\today} \maketitle End of Doc/s0.tex echo Doc/s1.tex 1>&2 cat >Doc/s1.tex <<'End of Doc/s1.tex' \section{Introduction} This document describes a software library called VSPLINE that is available via electronic-mail from NETLIB~\cite{netlib:87}. The VSPLINE library contains a set of subroutines that computes a non-parametric estimate of a smooth vector-valued function from noisy measurements. The algorithms are derived and described in detail in~\cite{Fessler:89assp}; this brief outlines the main ideas and explains the source code interface. End of Doc/s1.tex echo Doc/s2.tex 1>&2 cat >Doc/s2.tex <<'End of Doc/s2.tex' \section{Spline smoothing of vector measurements} \label{sec:vector} Consider the {\em vector} measurement model: \BE \boldy_n = \boldg(t_n) + \bveps_n, \ \ n=1, \ldots, N, \label{eq:multi-model} \EE \[ \boldg(t_n), \bveps_n, \boldy_n \in \Re^M, \ \ \bveps_n \sim N(0, \mcov_n), \ \ E \{ \bveps_n \bveps_m' \} = 0, \ n \neq m, \] with known symmetric and positive definite error covariances $\mcov_n$. The goal of smoothing is to estimate the smooth vector-valued function $\boldg$, given the measurements $\{\boldy_n\}_{n=1}^N$. The non-parametric approach to this problem prescribes a compromise between the conflicting goals of fit to the data and smoothness of the estimated functions. An estimate is obtained by using the following criterion, which is the natural generalization of a similar criterion used for scalar measurements: \BE \bghata = \arg \min_{\boldg} \sum_{n=1}^N (\boldy_n - \boldg(t_n))' \mcov_n^{-1} (\boldy_n - \boldg(t_n)) + R_{\balpha}(\boldg), \label{eq:multi-min} \EE where \[ R_{\balpha}(\boldg) = \sum_{m=1}^M \alpha_m \int (\ddot{g}_m(t))^2 \ dt . \] The smoothing parameter $\balpha$ controls the tradeoff between residual error and smoothness. Although $\bghata(t)$ is a continuous function, it is uniquely determined by its values at the knots $t_1, \ldots, t_N$, denoted $\hat{\boldy}_n = \bghata(t_n)$. Subroutine {\tt vspline}, described below, computes $\{ \hat{\boldy}_n \}_{n=1}^N$ given $\balpha$, %$\alpha_1, \ldots, \alpha_M$, $\{ \boldy_n \}_{n=1}^N $, $\{ t_n \}_{n=1}^N $, and $\{ \mcov_n^{-1} \}_{n=1}^N $. The computational requirements~\cite{Fessler:89assp} are $O(M^3\cdot N)$. End of Doc/s2.tex echo Doc/s3.tex 1>&2 cat >Doc/s3.tex <<'End of Doc/s3.tex' \section{Choosing the smoothing parameters - Cross Validation} \label{sec:auto} The smoothing parameter $\balpha$ can be chosen automatically~\cite{Fessler:89assp} by minimizing the cross validation (CV) score: \[ \balpha_{\CV} \defequ \arg \min_{\balpha} \CV(\balpha), \] \BE \CV(\balpha) \defequ \frac{1}{N} \sum_{n=1}^N (\meas_n - \hat{\boldg}_{\balpha, -n}(t_n))' \mcov_n^{-1} (\meas_n - \hat{\boldg}_{\balpha, -n}(t_n)). \label{eq:cv-score} \EE $\hat{\boldg}_{\balpha, -n}$ is the solution to the smoothing problem (\ref{eq:multi-min}) with N-1 data points, posed while excluding the pair $(t_n, \meas_n)$. Each data pair is dropped in turn, the smoothed curve $\bghatan$ is estimated, and the predicted value $\bghatan(t_n)$ is compared with the unused measurement $\boldy_n$. If the CV score is small, then we have chosen the smoothing parameter that makes the estimated curve a good self predictor. Although equation~(\ref{eq:cv-score}) illustrates the idea behind cross validation, it is computationally inefficient. %requires O($M^3N^2$) operations. We have shown that~(\ref{eq:cv-score}) can be rewritten \BE \CV(\balpha) = \frac{1}{N} \sum_{n=1}^N \| \mcov_n^{-\frac{1}{2}} (\boldI_M - \Ad{nn})^{-1} (\meas_n - \bghata(t_n)) \|^2, \label{eq:cv-score2} \EE where $\Ad{nn}$ is the $n$'th $M \times M$ block diagonal submatrix of the influence matrix. By using the Hutchinson and deHoog algorithm~\cite{HutchHoog:85}, (\ref{eq:cv-score2}) is computed in only O($M^3N$) operations. Subroutine {\tt bestcv}, described below, computes $\balpha_{\CV}$ and $\hat{\boldg}_{\balpha_{\CV}}$. The VSPLINE library also contains routines for computing Generalized Cross Validation and Unbiased Risk scores~\cite{Fessler:89assp}. End of Doc/s3.tex echo Doc/s4.tex 1>&2 cat >Doc/s4.tex <<'End of Doc/s4.tex' \section{Software Interface} \subsection{Conventions} We have adopted the Numerical Recipes in C~\cite{numrecc:88} convention of using unit offset arrays where convenient (almost everywhere). If an array is declared {\tt double vect[100];} then {\tt (\&vect[0]-1)} or {\tt (vect-1)} should be passed to the subroutines below. \subsection{Subroutine {\tt vspline}} The syntax for {\tt vspline} is: \begin{verbatim} int vspline(ys, alpha, tn, yn, inv_sigma, iid, M, N) double *ys, *alpha, *tn, *yn, *inv_sigma; int iid, M, N; \end{verbatim} %I repeat that the arrays are all unit offset, %i.e., {\tt alpha[1]} is the first element of {\tt alpha}, %and {\tt alpha[M]} is the last. {\tt M} is the dimension (length) the measurement vectors $\boldy_n$. {\tt N} is the number of samples. {\tt alpha} is the $M$ smoothing parameters. {\tt tn} is the $N$ sample points. If the samples are uniform, then set {\tt tn} to {\tt (double *) NULL}. The measurements are stored in {\tt yn} in the order $\{ y_{1,1},\ldots,y_{1,M}, \ldots, y_{N,1},\ldots,y_{N,M} \}$. If {\tt iid} is 1, the covariances are assumed to be identical, and {\tt inv\_sigma} is the $M^2$ elements of $\mcov^{-1}$, stored by columns or rows (it is symmetric). However, if {\tt iid} is 0, the covariances can vary with $n$, and {\tt inv\_sigma} is a $N \cdot M^2$ array of the concatenation of the elements of $\bSigma_1^{-1}, \ldots, \bSigma_N^{-1}$ in the natural order. On return, {\tt ys} is the $N\cdot M$ array of $\hat{\boldy}_1,\ldots,\hat{\boldy}_N$, stored in the same order as {\tt yn}. {\tt vspline} returns 1 if all goes well, 0 if there is an error. The possible errors are: running out of memory, {\tt tn} not strictly increasing, or singular covariance matrices. \subsection{Subroutine {\tt bestcv}} The syntax for {\tt bestcv} is identical to that of {\tt vspline}. The difference is that {\tt alpha} is used as an initial estimate for $\balpha_{\CV}$. On return, {\tt alpha} is set to $\balpha_{\CV}$, and {\tt ys} is the smoothed estimate at that value of $\balpha$. \subsection{Subroutine {\tt cubic\_interp}} Each component function of $\bghata(t)$ is a scalar cubic spline. To compute values of $\hat{g}_{\balpha,m}(t)$ for values of $t$ other than the knots, call {\tt cubic\_interp}. %Currently, %{\tt cubic\_interp} %is implemented only for scalar functions. %The syntax is: \begin{verbatim} int cubic_interp(yk, tn, ys, tk, N, K) double *yk, *tn, *ys, *tk; int N, K; \end{verbatim} {\tt tn} is the $N$ sample points, and {\tt ys} is the $N$ smoothed estimates for one of the component functions, e.g. $\{\hat{y}_{n,m}\}_{n=1}^N$. \ {\tt tk} is $K$ new points at which $\hat{g}_{\balpha,m}(t)$ is to be evaluated. On return, {\tt yk} is the $K$ interpolated samples. {\tt cubic\_interp} returns 1 if successful, 0 otherwise. {\tt cubic\_interp} is not very efficient, and may not be exact for values of $t$ outside of $[t_1,t_N]$. End of Doc/s4.tex echo Doc/s5.tex 1>&2 cat >Doc/s5.tex <<'End of Doc/s5.tex' \section{Installation} VSPLINE is distributed as a Bourne Shell ``bundle''. Create a directory called {\tt vspline}, move the bundled file to that directory, and type {\tt sh} {\em filename} at a UNIX prompt. The file {\tt README} contains important installation instructions. When installed, the C source code resides in several subdirectories. {\tt matrix} contains very general vector and matrix utility subroutines. Simply by changing the definitions in {\tt matrix/matrix.h}, one could use a different library. {\tt smooth} contains utilities for banded matrices, and the smoothing subroutines themselves. Finally, {\tt Test} contains programs we used to test the library. The test programs were applied to the data files in {\tt Data}; to verify your installation, follow the instructions in {\tt Data/README}. One {\tt makefile} in the parent directory controls compilation of all three subdirectories. This makefile creates a library called {\tt libvspline.a}. The {\tt powell} subroutine of the Numerical Recipes in C library~\cite{numrecc:88} is used to minimize the CV score. Therefore, you must load in both {\tt libvspline.a} and (e.g.) {\tt libnumrecC.a}. If you do not have this library, type {\tt make libvsmooth.a}. This latter library is a subset of the VSPLINE library that contains the basic smoothing routines (sufficient to use the {\tt vspline} routine), but does not contain the routines for minimizing the CV score. {\tt vspline} resides in {\tt smooth/example1.c}, {\tt bestcv} in {\tt smooth/example2.c}. A subroutine called {\tt vect\_smooth} in {\tt smooth/example1.c} gives another example of how the VSPLINE library can used. {\tt cubic\_interp} resides in {\tt smooth/interp.c}. The \LaTeX source for this document is in {\tt smooth/Doc}. End of Doc/s5.tex echo Doc/size.tex 1>&2 cat >Doc/size.tex <<'End of Doc/size.tex' \setlength{\headsep}{0.0in} \setlength{\headheight}{0.0in} \setlength{\textwidth}{6.5in} \setlength{\textheight}{8.5in} \setlength{\topmargin}{0.0in} \setlength{\oddsidemargin}{0.0in} \setlength{\evensidemargin}{0.0in} End of Doc/size.tex echo Test/data.h 1>&2 cat >Test/data.h <<'End of Test/data.h' /***************************** data.h ************************************ * DELCARATIONS FOR DATA THAT IS READ FROM A FILE * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ typedef struct { ALPHA_SET alpha; VECT ti; CMAT yi; COVAR *covar; int M, N; int flag_iid, flag_equ; } FILEDATA; End of Test/data.h echo Test/read_data.c 1>&2 cat >Test/read_data.c <<'End of Test/read_data.c' /*********************** read_data.c ************************************* * FOR TESTING: READ Data TO BE SMOOTHED FROM stdin * FILE FORMAT: * N M then alpha1 ... alphaM then * if (!flag_equ): t1, ..., tN * for n=1 to N y1n y2n ... yMn * for n=1 to N isig11n .. isigMMn * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" #include "data.h" void read_data(); extern int inv_to_covar(); extern COVAR *alloc_covar(); /* * INPUT: * fd->file_iid, file_equ */ void read_data(fd) FILEDATA *fd; { int MM, NN, nn, Nc; CMAT isig; if (2 != scanf("%d %d", &NN, &MM)) ABORT("Bad data file.\n"); if (NN < 3 || MM < 1) ABORT("Bad data file.\n"); fd->N = NN; fd->M = MM; Nc = fd->flag_iid ? 1 : NN; if (!fd->flag_equ) fd->ti = ALLOC_VECTOR(NN); else fd->ti = (VECT) 0; if (!(fd->alpha = ALLOC_ALPHA(MM))) ABORT("Mem 0"); if (!(fd->yi = ALLOC_VECTOR(MM*NN))) ABORT("Mem 1"); if (!(isig = ALLOC_VECTOR(MM*MM))) ABORT("Mem 2"); if (!(fd->covar = alloc_covar(MM, Nc))) ABORT("Mem 3"); SCAN_ALPHA(fd->alpha, MM); if (!fd->flag_equ) SCAN_VECTOR(fd->ti, NN); SCAN_VECTOR(fd->yi, NN*MM); for (nn = 1; nn <= Nc; nn++) { SCAN_VECTOR(isig, MM*MM); if (!inv_to_covar(&(fd->covar)[nn], isig, MM)) ABORT("Singular covariance matrix\n"); } FREE_VECTOR(isig); } End of Test/read_data.c echo Test/test_ci.c 1>&2 cat >Test/test_ci.c <<'End of Test/test_ci.c' /*********************** test_ci.c *************************************** * Test cubic_interp(). * Read data, interpolate, print out. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #define USAGE { fprintf(stderr, "Usage: %s Ni Nx\n", argv[0]); exit(-1);} #include #include "defs.h" #include "smooth.h" static void read_data(); extern int cubic_interp(); main(argc, argv) int argc; char *argv[]; { VECT ti, yi, xs, ys; int Ni, Nx; if (argc != 3) USAGE; Ni = atoi(argv[1]); Nx = atoi(argv[2]); TMP_ALLOC_VECT(ti, Ni); TMP_ALLOC_VECT(yi, Ni); TMP_ALLOC_VECT(xs, Nx); TMP_ALLOC_VECT(ys, Nx); SCAN_VECTOR(ti, Ni); SCAN_VECTOR(yi, Ni); SCAN_VECTOR(xs, Nx); cubic_interp(ys, ti, yi, xs, Ni, Nx); PRINT_VECTOR(ys, Nx); } End of Test/test_ci.c echo Test/test_ss.c 1>&2 cat >Test/test_ss.c <<'End of Test/test_ss.c' /*********************** test_ss.c *************************************** * TEST best_indiv() * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #define USAGE { fprintf(stderr, "Usage: %s score_id\n", argv[0]); exit(-1);} #include #include "smooth.h" #include "data.h" extern void read_data(); extern double best_indiv(); #define TOL 0.01 main(argc, argv) int argc; char *argv[]; { FILEDATA fd; VECT trace; CMAT ys; int score_id; if (argc < 2) USAGE; score_id = atoi(argv[1]); fd.flag_iid = 1; fd.flag_equ = 0; fd.ti = (VECT) 0; read_data(&fd); if (!(ys = ALLOC_VECTOR(fd.N*fd.M))) ABORT("Mem 98"); if (!(trace = ALLOC_VECTOR(fd.M))) ABORT("Mem 99"); if (!best_indiv(ys, trace, fd.alpha, fd.ti, fd.yi, score_id, TOL, fd.M, fd.N)) ABORT("Problem in best_ind"); PRINT_ALPHA(fd.alpha, fd.M); PRINT_VECTOR(trace, fd.M); } End of Test/test_ss.c echo Test/test_vs.c 1>&2 cat >Test/test_vs.c <<'End of Test/test_vs.c' /*********************** test_vs.c *************************************** * Read data from file (see read_data.c), smooth it, compute scores. * USES atof() * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #define USAGE { fprintf(stderr, "Usage: %s iid equi [score tol]\n", argv[0]); exit(-1);} #include #include #include "smooth.h" #include "data.h" #define TOL 0.01 extern void read_data(); extern int vect_smooth(); extern double best_score(); main(argc, argv) int argc; char *argv[]; { FILEDATA fd; double ur, cv, gcv, penin, penout, rough; CMAT ys; int flag_best; if (argc < 3) USAGE; flag_best = argc > 3; fd.flag_iid = atoi(argv[1]); fd.flag_equ = atoi(argv[2]); fd.ti = (VECT) 0; read_data(&fd); PRINT_ALPHA(fd.alpha, fd.M); if (!fd.flag_equ) PRINT_VECTOR(fd.ti, fd.N); PRINT_VECTOR(fd.yi, fd.N*fd.M); PRINT_BAND(fd.covar[1].inv, fd.M, fd.M-1); if (!(ys = ALLOC_VECTOR(fd.N*fd.M))) ABORT("Mem 99"); if (flag_best) { int score_id = atoi(argv[3]); double tol = (argc > 4) ? atof(argv[4]) : TOL; double bs = best_score(ys, fd.alpha, fd.ti, fd.yi, fd.covar, fd.flag_iid, fd.M, fd.N, score_id, tol); message0("best score ", bs); PRINT_ALPHA(fd.alpha, fd.M); PRINT_VECTOR(ys, fd.M*fd.N); } else { if ( !vect_smooth(ys, &ur, &cv, &gcv, &penin, &penout, &rough, fd.alpha, fd.ti, fd.yi, fd.covar, fd.flag_iid, fd.M, fd.N) ) ABORT("E 7"); PRINT_VECTOR(ys, fd.N*fd.M); message("ur score", ur); message("cv score", cv); message("gcv score", gcv); message("roughness", rough); message("penalty yi (in)", penin); message("penalty ys (out)", penout); } } End of Test/test_vs.c echo matrix/alloc.c 1>&2 cat >matrix/alloc.c <<'End of matrix/alloc.c' /***************************** alloc.c *********************************** * SUBROUTINES FOR ALLOCATING AND FREEING MATRICES AND VECTORS * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "defs.h" double *alloc_vector_d(), **alloc_matrix_d(); float *alloc_vector_f(), **alloc_matrix_f(), ***alloc_matrix3_f(); /* DOUBLE PRECISION SUBROUTINES */ /* * ALLOCATE SPACE FOR A UNIT OFFSET VECTOR: V[1..N] */ double *alloc_vector_d(N) int N; { double *p; if ( !(p = (double *) ALLOC(N, sizeof(*p))) ) return 0; return p-1; } /* * ALLOCATE SPACE FOR A [1..rows][1..cols] MATRIX * SUCH THAT MATRIX ELEMENTS ARE CONSECTIVE */ double **alloc_matrix_d(rows, cols) int rows, cols; { double **p; if ( !(p = (double **) ALLOC(rows, sizeof(*p))) ) return 0; p--; if ( !(p[1] = alloc_vector_d(rows*cols)) ) return 0; for (; rows > 1; rows--) p[rows] = (rows-1)*cols + p[1]; return p; } /* SINGLE PRECISION SUBROUTINES */ /* * ALLOCATE SPACE FOR A UNIT OFFSET VECTOR: V[1..N] */ float *alloc_vector_f(N) int N; { float *p; if ( !(p = (float *) ALLOC(N, sizeof(*p))) ) return 0; return p-1; } /* * ALLOCATE SPACE FOR A [1..rows][1..cols] MATRIX * SUCH THAT MATRIX ELEMENTS ARE CONSECTIVE */ float **alloc_matrix_f(rows, cols) int rows, cols; { float **p; if ( !(p = (float **) ALLOC(rows, sizeof(*p))) ) return 0; p--; if ( !(p[1] = alloc_vector_f(rows*cols)) ) return 0; for (; rows > 1; rows--) p[rows] = (rows-1)*cols + p[1]; return p; } /* * ALLOCATE SPACE FOR A [1..rows][1..cols][1..last] MATRIX * SUCH THAT MATRIX ELEMENTS ARE CONSECTIVE */ float ***alloc_matrix3_f(rows, cols, last) int rows, cols, last; { int r; float ***p; if ( !(p = (float ***) ALLOC(rows, sizeof(*p))) ) return 0; p--; if ( !(p[1] = (float **) ALLOC(rows*cols, sizeof(**p))) ) return 0; p[1]--; for (r=rows; r > 1; r--) p[r] = (r-1)*cols + p[1]; if ( !(p[1][1] = (float *) ALLOC(rows*cols*last, sizeof(***p))) ) return 0; p[1][1]--; for (r=rows*cols; r > 1; r--) p[1][r] = (r-1) * last + p[1][1]; return p; } End of matrix/alloc.c echo matrix/debug.c 1>&2 cat >matrix/debug.c <<'End of matrix/debug.c' /***************************** debug.c *********************************** * MATRIX AND VECTOR DEBUGGING AND IOAND IO SUBROUTINES * print_xxx ascii to stdout * scan_xxx ascii from stdin * USES: sprintf(), fwrite(), fopen(), fclose(), fflush(), etc * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "mexdep.h" message(str, val) char *str; double val; { printf("%s %g\n", str, val); #ifndef MATLAB FFLUSH(stdout); #endif } message0(str, val) char *str; double val; { printf("%s%g ", str, val); } /* DOUBLE PRECISION SUBROUTINES */ /* * vect[1..N] */ print_vector_d(vect, N) double *vect; int N; { int n; for (n=1; n <= N; n++) { printf("%g ", *++vect); if (!(n % 10)) printf("\n"); } if (N % 10) printf("\n"); } /* * vect[1..N] from stdin */ scan_vector_d(vect, N) double *vect; int N; { int n; for (n=1; n <= N; n++) if (1 != scanf("%lg", &vect[n])) EXITBUG("scanf") } /* SINGLE PRECISION SUBROUTINES */ /* * vect[1..N] */ print_vector_f(vect, N) float *vect; int N; { int n; for (n=1; n <= N; n++) { printf("%g ", *++vect); if (!(n % 10)) printf("\n"); } if (N % 10) printf("\n"); } /* * vect[1..N] from stdin */ scan_vector_f(vect, N) float *vect; int N; { int n; for (n=1; n <= N; n++) if (1 != scanf("%g", &vect[n])) EXITBUG("scanf") } End of matrix/debug.c echo matrix/ident.c 1>&2 cat >matrix/ident.c <<'End of matrix/ident.c' /***************************** ident.c *************************** * Return identity matrix. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *****************************************************************/ #include "defs.h" #include "matrix.h" float **identity_f(); double **identity_d(); /* * IDENTITY MATRIX * OUTPUT: * ident[1..M][1..M] M*M identity matrix */ float **identity_f(M) register int M; { float **ident; if (!(ident = ALLOC_MATRIX_F(M, M))) EXITBUG("mem alloc"); for (; M; M--) ident[M][M] = 1.; return ident; } /* * IDENTITY MATRIX * OUTPUT: * ident[1..M][1..M] M*M identity matrix */ double **identity_d(M) register int M; { double **ident; if (!(ident = ALLOC_MATRIX_D(M, M))) EXITBUG("mem alloc"); for (; M; M--) ident[M][M] = 1.; return ident; } End of matrix/ident.c echo matrix/io_d.c 1>&2 cat >matrix/io_d.c <<'End of matrix/io_d.c' /***************************** matrix_io_d.c ***************************** * Double Precision IO SUBROUTINES (mex independent) * write_xxx binary to file * fprint_xxx ascii to file.ext (or file if !ext) * fscan_xxx ascii from file * read_xxx binary from file * USES: sprintf(), fwrite(), fopen(), fclose(), fflush(), etc * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include #include "defs.h" #define TYPE double /* * mat[1..N][1..M] */ print_matrix_d(mat, N, M) TYPE **mat; int N, M; { int n; for (n=1; n <= N; n++) print_vector_d(mat[n], M); } /* * vect[1..N] */ write_vector_d(name, vect, N) char *name; TYPE *vect; int N; { FILE *fp; FOPEN(fp, name, "w"); FWRITE(&vect[1], sizeof(*vect), N, fp); FCLOSE(fp); } /* * mat[1..N][1..M] -> file "name"(.ext) */ fprint_matrix_d(name, ext, mat, N, M) char *name; int *ext; TYPE **mat; int N, M; { int n, m; FILE *fp; char fname[256]; if (ext) sprintf(fname, "%s.%02d", name, *ext); else sprintf(fname, "%s", name); FOPEN(fp, fname, "w"); for (n=1; n <= N; n++) { for (m=1; m <= M; m++) if (EOF == fprintf(fp, "%g ", mat[n][m])) EXITBUG("fprintf") if (EOF == fprintf(fp, "\n")) EXITBUG("fprintf") } FCLOSE(fp); } /* * vect[1..N] from file "name" */ fscan_vector_d(vect, name, N) TYPE *vect; char *name; int N; { int n; FILE *fp; FOPEN(fp, name, "r"); for (n=1; n <= N; n++) if (1 != fscanf(fp, "%lg", &vect[n])) EXITBUG("fscanf") FCLOSE(fp); } /* * vect[1..N] from file "name" (binary) */ read_vector_d(vect, name, N) TYPE *vect; char *name; int N; { FILE *fp; FOPEN(fp, name, "r"); FREAD(&vect[1], sizeof(*vect), N, fp); FCLOSE(fp); } End of matrix/io_d.c echo matrix/io_f.c 1>&2 cat >matrix/io_f.c <<'End of matrix/io_f.c' /***************************** matrix_io_f.c ***************************** * Single Precision IO SUBROUTINES (mex independent) * write_xxx binary to file * fprint_xxx ascii to file.ext (or file if !ext) * fscan_xxx ascii from file * read_xxx binary from file * USES: sprintf(), fwrite(), fopen(), fclose(), fflush(), etc * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include #include "defs.h" #define TYPE float /* * mat[1..N][1..M] */ print_matrix_f(mat, N, M) TYPE **mat; int N, M; { int n; for (n=1; n <= N; n++) print_vector_f(mat[n], M); } /* * vect[1..N] */ write_vector_f(name, vect, N) char *name; TYPE *vect; int N; { FILE *fp; FOPEN(fp, name, "w"); FWRITE(&vect[1], sizeof(*vect), N, fp); FCLOSE(fp); } /* * mat[1..N][1..M] -> file "name"(.ext) */ fprint_matrix_f(name, ext, mat, N, M) char *name; int *ext; TYPE **mat; int N, M; { int n, m; FILE *fp; char fname[256]; if (ext) sprintf(fname, "%s.%02d", name, *ext); else sprintf(fname, "%s", name); FOPEN(fp, fname, "w"); for (n=1; n <= N; n++) { for (m=1; m <= M; m++) if (EOF == fprintf(fp, "%g ", mat[n][m])) EXITBUG("fprintf") if (EOF == fprintf(fp, "\n")) EXITBUG("fprintf") } FCLOSE(fp); } /* * vect[1..N] from file "name" */ fscan_vector_f(vect, name, N) TYPE *vect; char *name; int N; { int n; FILE *fp; FOPEN(fp, name, "r"); for (n=1; n <= N; n++) if (1 != fscanf(fp, "%g", &vect[n])) EXITBUG("fscanf") FCLOSE(fp); } /* * vect[1..N] from file "name" (binary) */ read_vector_f(vect, name, N) TYPE *vect; char *name; int N; { FILE *fp; FOPEN(fp, name, "r"); FREAD(&vect[1], sizeof(*vect), N, fp); FCLOSE(fp); } End of matrix/io_f.c echo matrix/matrix.h 1>&2 cat >matrix/matrix.h <<'End of matrix/matrix.h' /***************************** matrix.h ********************************** * General purpose vector and matrix utilites. * I use the "Numerical Recipes in C" convention for allocating * vectors and matrices, so * A[a..b][c..d] is a matrix with b-a+1 rows and d-c+1 columns. * A[1..nxm]: an nxm "column" matrix stored 1,1 2,1 ... n,1 ... n,m. * However, not only are they allocated as a doubly subscripted array, * i.e., mat[1..N][1..M], but the data is stored sequentially * as well. So &mat[1][1] is a pointer to NxM elements stored * in the order: * mat[1][1], mat[1][2], ..., mat[1][M], mat[2][1], ... mat[N][M]. * Thus we can easily operate on the individual rows of the matrix, * or on the entire matrix. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ extern double *alloc_vector_d(), **alloc_matrix_d(); extern float *alloc_vector_f(), **alloc_matrix_f(), ***alloc_matrix3_f(); extern void scale_vector_d(), inc_scaled_vector_d(); extern void scale_vector_f(), inc_scaled_vector_f(); extern void vector_add_d(), vector_sub_d(), vector_inc_d(), vector_dec_d(); extern void vector_add_f(), vector_sub_f(), vector_inc_f(), vector_dec_f(); extern double norm_vector_d(), dot_vector_d(); extern double norm_vector_f(), dot_vector_f(); extern float **identity_f(); extern double **identity_d(); /* Precision Independent */ #define COPY_VECTOR(src,dst,n) BCOPY(&src[1],&dst[1],n) #define COPY_MATRIX(src,dst,r,c) COPY_VECTOR(src[1],dst[1],(r)*(c)) #define COPY_MATRIX3(src,dst,r,c,l) COPY_VECTOR(src[1][1],dst[1][1],(r)*(c)*(l)) #define ZERO_VECTOR(v,n) BZERO(&v[1],n) #define ZERO_MATRIX(m,r,c) ZERO_VECTOR(m[1],(r)*(c)) #define ZERO_MATRIX3(m,r,c,l) ZERO_VECTOR(m[1][1],(r)*(c)*(l)) #define FREE_VECTOR(v) { FREE(&v[1]); } #define FREE_MATRIX(m) { FREE(&m[1][1]); FREE(&m[1]); } #define FREE_MATRIX3(m) { FREE(&m[1][1][1]); FREE(&m[1][1]); FREE(&m[1]); } /* Double Precision */ #define ALLOC_VECTOR_D alloc_vector_d #define ALLOC_MATRIX_D alloc_matrix_d #define ALLOC_MATRIX3_D alloc_matrix3_d #define SCALE_VECTOR_D scale_vector_d #define NORM_VECTOR_D norm_vector_d #define DOT_VECTOR_D dot_vector_d #define INC_SCALED_VECTOR_D inc_scaled_vector_d #define ADD_VECTOR_D vector_add_d #define INC_VECTOR_D vector_inc_d #define SUB_VECTOR_D vector_sub_d #define DEC_VECTOR_D vector_dec_d #define SCAN_VECTOR_D scan_vector_d #define READ_VECTOR_D read_vector_d #define FSCAN_VECTOR_D fscan_vector_d #define PRINT_VECTOR_D print_vector_d #define PRINT_MATRIX_D print_matrix_d #define WRITE_VECTOR_D write_vector_d #define FPRINT_MATRIX_D fprint_matrix_d #define SCALE_MATRIX_D(m,s,r,c) SCALE_VECTOR_D(m[1],s,(r)*(c)) #define SCAN_MATRIX_D(m,r,c) SCAN_VECTOR_D(m[1],(r)*(c)) #define READ_MATRIX_D(m,f,r,c) READ_VECTOR_D(m[1],f,(r)*(c)) #define FSCAN_MATRIX_D(m,f,r,c) FSCAN_VECTOR_D(m[1],f,(r)*(c)) #define WRITE_MATRIX_D(name,m,r,c) WRITE_VECTOR_D(name,m[1],(r)*(c)) /* Single Precision */ #define ALLOC_VECTOR_F alloc_vector_f #define ALLOC_MATRIX_F alloc_matrix_f #define ALLOC_MATRIX3_F alloc_matrix3_f #define SCALE_VECTOR_F scale_vector_f #define NORM_VECTOR_F norm_vector_f #define DOT_VECTOR_F dot_vector_f #define INC_SCALED_VECTOR_F inc_scaled_vector_f #define ADD_VECTOR_F vector_add_f #define INC_VECTOR_F vector_inc_f #define SUB_VECTOR_F vector_sub_f #define DEC_VECTOR_F vector_dec_f #define SCAN_VECTOR_F scan_vector_f #define READ_VECTOR_F read_vector_f #define FSCAN_VECTOR_F fscan_vector_f #define PRINT_VECTOR_F print_vector_f #define PRINT_MATRIX_F print_matrix_f #define WRITE_VECTOR_F write_vector_f #define FPRINT_MATRIX_F fprint_matrix_f #define SCALE_MATRIX_F(m,s,r,c) SCALE_VECTOR_F(m[1],s,(r)*(c)) #define SCAN_MATRIX_F(m,r,c) SCAN_VECTOR_F(m[1],(r)*(c)) #define READ_MATRIX_F(m,f,r,c) READ_VECTOR_F(m[1],f,(r)*(c)) #define FSCAN_MATRIX_F(m,f,r,c) FSCAN_VECTOR_F(m[1],f,(r)*(c)) #define WRITE_MATRIX_F(name,m,r,c) WRITE_VECTOR_F(name,m[1],(r)*(c)) /* MACROS FOR ALLOCATING TEMPORARY WORK SPACE */ /* Double Precision */ #define TMP_ALLOC_VECTOR_D(v,n) {v = (double *) TMP_ALLOC((n),sizeof(*v)); \ if (!v) ABORT("alloca failed"); \ BZERO(v, (n)); \ v--; \ } #define TMP_FREE_VECTOR_D(v) TMP_FREE(&v[1]) #define TMP_ALLOC_MATRIX_D(m,r,c) { int i; \ m = (double **) TMP_ALLOC((r),sizeof(*m)); \ if (!m) ABORT("alloca failed"); \ m--; \ TMP_ALLOC_VECTOR_D(m[1], (r)*(c)) \ for (i=(r); i > 1; i--) \ m[i] = (i-1)*(c) + m[1]; \ } #define TMP_FREE_MATRIX_D(m) { TMP_FREE(&m[1][1]); TMP_FREE(&m[1]); } /* Single Precision */ #define TMP_ALLOC_VECTOR_F(v,n) {v = (float *) TMP_ALLOC((n),sizeof(*v)); \ if (!v) ABORT("alloca failed"); \ BZERO(v, (n)); \ v--; \ } #define TMP_FREE_VECTOR_F(v) TMP_FREE(&v[1]) #define TMP_ALLOC_MATRIX_F(m,r,c) { int i; \ m = (float **) TMP_ALLOC((r),sizeof(*m)); \ if (!m) ABORT("alloca failed"); \ m--; \ TMP_ALLOC_VECTOR_F(m[1], (r)*(c)) \ for (i=(r); i > 1; i--) \ m[i] = (i-1)*(c) + m[1]; \ } #define TMP_FREE_MATRIX_F(m) { TMP_FREE(&m[1][1]); TMP_FREE(&m[1]); } #define TMP_ALLOC_MATRIX3_F(m,r,c,l) { int i; \ m = (float ***) TMP_ALLOC((r),sizeof(*m)); \ if (!m) EXITBUG("alloca failed"); \ m--; \ m[1] = (float **) TMP_ALLOC((r)*(c),sizeof(**m)); \ if (!m[1]) EXITBUG("alloca failed"); \ m[1]--; \ TMP_ALLOC_VECTOR_F(m[1][1], (r)*(c)*(l)); \ for (i=(r); i > 1; i--) \ m[i] = (i-1)*(c) + m[1]; \ for (i=(r)*(c); i > 1; i--) \ m[1][i] = (i-1)*(l) + m[1][1]; \ } #define TMP_FREE_MATRIX3_F(m) { TMP_FREE(&m[1][1][1]); TMP_FREE(&m[1][1]); TMP_FREE(&m[1]); } End of matrix/matrix.h echo matrix/mexdep.c 1>&2 cat >matrix/mexdep.c <<'End of matrix/mexdep.c' /***************************** mexdep.c ********************************** * SUBROUTINES WHICH DEPEND ON WHETHER OR NOT MEX FILES ARE USED. * USES fprintf(), calloc(), free() * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "mexdep.h" char *mycalloc(), *alloc_check(); void myfree(); jumpship(); jumpship(x) char *x; { #ifdef MATLAB mex_error(x); #else if (EOF == fprintf(stderr, "%s\n", x)) exit(-1); /* pretty silly, huh? */ exit(-1); #endif } /* * A GENERIC INTERFACE SO THAT OTHER SUBROUTINES * NEED NOT WORRY ABOUT mex VS reg. */ char *mycalloc(n1, n2) int n1, n2; { return calloc((unsigned) n1, (unsigned) n2); } /* * ALLOCATION WITH CHECKING. */ char *alloc_check(n1, n2) int n1, n2; { char *p; p = calloc((unsigned) n1, (unsigned) n2); if (!p) ABORT("Can't allocate memory"); return p; } void myfree(p) char *p; { free(p); } End of matrix/mexdep.c echo matrix/mexdep.h 1>&2 cat >matrix/mexdep.h <<'End of matrix/mexdep.h' /***************************** mexdep.h ********************************** * Header files for code that uses MEX dependent subroutines. * For compiling Matlab MEX files, must use the MEX library * rather than stdio. Quite a pain, frankly. * Affected functions: malloc(), calloc(), free(), printf() * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "defs.h" #ifdef MATLAB #include "cmex.h" #else #include #endif #ifdef MATLAB #ifdef DoublePrec /* This is the right way... #define MEX_TO_VECTOR(v, mat) BCOPY(mat->pr-1, &v[1], mat->n*mat->m) #define VECTOR_TO_MEX(mat, v) BCOPY(&v[1], mat->pr-1, mat->n*mat->m) */ #define MEX_TO_VECTOR(v, mat) {v = (mat->pr - 1);} #define VECTOR_TO_MEX(mat, v) #else #define MEX_TO_VECTOR(v, mat) VECTOR_D_TO_F(mat->pr-1, &v[1], mat->n*mat->m) #define VECTOR_TO_MEX(mat, v) VECTOR_F_TO_D(&v[1], mat->pr-1, mat->n*mat->m) #endif DoublePrec #endif MATLAB End of matrix/mexdep.h echo matrix/sub_mat.c 1>&2 cat >matrix/sub_mat.c <<'End of matrix/sub_mat.c' /***************************** sub_mat.c ********************************* * UTILITY MATRIX SUBROUTINES * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #ifdef NOTUSED /* * MULTIPLY TWO MATRICES (OR MATRIX AND VECTOR) * INPUT: * in1[1..r1xc1] * in2[1..c1xc2] * OUTPUT: * out[1..r1xc2] out = in1 * in2 * FLOPS: r1 c2 (2 c1 - 1) */ void cmat_mul(out, in1, r1, c1, c2, in2) CMAT out, in1, in2; int r1, c1, c2; { register int i,j,k; register DATA *t; for (i = r1; i >= 1; i--) { for (j = c2; j >= 1; j--) { t = &out[(j-1)*r1 + i]; *t += in1[(c1-1)*r1 + i] * in2[(j-1)*c1 + c1]; for (k = c1-1; k >= 1; k--) *t += in1[(k-1)*r1 + i] * in2[(j-1)*c1 + k]; } } } /* * TRACE OF A MATRIX * INPUT: * in[1..NxN] * RETURNS: * trace(in) * FLOPS: N-1 */ DATA cmat_trace(in, N) CMAT in; int N; { register int n; register DATA tr = in[N*N]; for (n = N-1; n >= 1; n--) tr += in[(n-1)*N + n]; return tr; } /* * TAKE THE "TRANSPOSE" OF A SET OF SIGNALS * INPUT: * in[0..NM-1] * OUTPUT: * out[0..MN-1] */ void cmat_transpose(out, in, N, M) VECT out, in; int N, M; { register int n, m; for (n=N-1; n >= 0; n--) for (m=M-1; m >= 0; m--) out[n*M+m] = in[m*N+n]; } #endif End of matrix/sub_mat.c echo matrix/vector_d.c 1>&2 cat >matrix/vector_d.c <<'End of matrix/vector_d.c' /***************************** vector_d.c ******************************** * Double Precision UTILITY VECTOR SUBROUTINES * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ void vector_add_d(), vector_sub_d(), vector_inc_d(), vector_dec_d(); void scale_vector_d(), vector_f_to_d(), inc_scaled_vector_d(); double norm_vector_d(), dot_vector_d(); typedef double *VECTOR; /* * ADD TWO VECTORS * INPUT: * in1[1..N] * in2[1..N] * OUTPUT: * out[1..N] out = in1 - in2 * FLOPS: N */ void vector_add_d(out, in1, in2, N) register VECTOR out, in1, in2; register int N; { for (; N; N--) *(++out) = *(++in1) + *(++in2); } /* * SUBTRACT TWO VECTORS * INPUT: * in1[1..N] * in2[1..N] * OUTPUT: * out[1..N] out = in1 - in2 * FLOPS: N */ void vector_sub_d(out, in1, in2, N) register VECTOR out, in1, in2; register int N; { for (; N; N--) *(++out) = *(++in1) - *(++in2); } /* * ACCUMULATE A VECTOR * INPUT: * inX[1..N] * OUTPUT: * in1[1..N] in1 = in1 + in2 * FLOPS: N */ void vector_inc_d(in1, in2, N) register VECTOR in1, in2; register int N; { for (; N; N--) *(++in1) += *(++in2); } /* * DECREMENT A VECTOR * INPUT: * inX[1..N] * OUTPUT: * in1[1..N] in1 = in1 - in2 * FLOPS: N */ void vector_dec_d(in1, in2, N) register VECTOR in1, in2; register int N; { for (; N; N--) *(++in1) -= *(++in2); } /* * SCALE v[1..N] BY s */ void scale_vector_d(v, s, N) register VECTOR v; register double s; register int N; { for ( ; N; N--) *++v *= s; } /* * Convert v[1..n] from single to double */ void vector_f_to_d(v, vf, n) register VECTOR v; register float *vf; register int n; { for (; n; n--) *++v = *++vf; } /* * v1[1..n] += scale * v2[1..n] */ void inc_scaled_vector_d(v1, s, v2, n) register VECTOR v1, v2; register double s; register int n; { for (; n; n--) *++v1 += s * *++v2; } /* * square of Euclidian norm of v[1..n] */ double norm_vector_d(v, n) register VECTOR v; register int n; { register double d = 0.; v++; for (; n; n--, v++) d += *v * *v; return d; } /* * dot product of v1[1..n] and v2[1..n] */ double dot_vector_d(v1, v2, n) register VECTOR v1, v2; register int n; { register double d = 0.; for (; n; n--) d += *++v1 * *++v2; return d; } End of matrix/vector_d.c echo matrix/vector_f.c 1>&2 cat >matrix/vector_f.c <<'End of matrix/vector_f.c' /***************************** vector_f.c ******************************** * Single Precision UTILITY VECTOR SUBROUTINES * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ void vector_add_f(), vector_sub_f(), vector_inc_f(), vector_dec_f(); void scale_vector_f(), vector_d_to_f(), inc_scaled_vector_f(); double norm_vector_f(), dot_vector_f(); typedef float *VECTOR; /* * ADD TWO VECTORS * INPUT: * in1[1..N] * in2[1..N] * OUTPUT: * out[1..N] out = in1 - in2 * FLOPS: N */ void vector_add_f(out, in1, in2, N) register VECTOR out, in1, in2; register int N; { for (; N; N--) *(++out) = *(++in1) + *(++in2); } /* * SUBTRACT TWO VECTORS * INPUT: * in1[1..N] * in2[1..N] * OUTPUT: * out[1..N] out = in1 - in2 * FLOPS: N */ void vector_sub_f(out, in1, in2, N) register VECTOR out, in1, in2; register int N; { for (; N; N--) *(++out) = *(++in1) - *(++in2); } /* * ACCUMULATE A VECTOR * INPUT: * inX[1..N] * OUTPUT: * in1[1..N] in1 = in1 + in2 * FLOPS: N */ void vector_inc_f(in1, in2, N) register VECTOR in1, in2; register int N; { for (; N; N--) *(++in1) += *(++in2); } /* * DECREMENT A VECTOR * INPUT: * inX[1..N] * OUTPUT: * in1[1..N] in1 = in1 - in2 * FLOPS: N */ void vector_dec_f(in1, in2, N) register VECTOR in1, in2; register int N; { for (; N; N--) *(++in1) -= *(++in2); } /* * SCALE v[1..N] BY s */ void scale_vector_f(v, s, N) register VECTOR v; register double s; register int N; { for ( ; N; N--) *++v *= s; } /* * Convert v[1..n] from double to single */ void vector_d_to_f(v, vd, n) register VECTOR v; register double *vd; register int n; { for (; n; n--) *++v = *++vd; } /* * v1[1..n] += scale * v2[1..n] */ void inc_scaled_vector_f(v1, s, v2, n) register VECTOR v1, v2; register double s; register int n; { for (; n; n--) *++v1 += s * *++v2; } /* * square of Euclidian norm of v[1..n] */ double norm_vector_f(v, n) register VECTOR v; register int n; { register double d = 0.; v++; for (; n; n--, v++) d += *v * *v; return d; } /* * dot product of v1[1..n] and v2[1..n] */ double dot_vector_f(v1, v2, n) register VECTOR v1, v2; register int n; { register double d = 0.; for (; n; n--) d += *++v1 * *++v2; return d; } End of matrix/vector_f.c echo smooth/band_util.c 1>&2 cat >smooth/band_util.c <<'End of smooth/band_util.c' /****************************** band_util.c ******************************* * UTILITY SUBROUTINES FOR BAND MATRICES * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" double band_trace(), trace_prod(); void mult_band_vect(), mult_full_vect(), mult_full_band(); void sym_band_to_mat(), mat_to_sym_band(); /****************************** band_trace() ****************************** * TRACE OF A SYMMETRIC BAND MATRIX * INPUT: * A: SYMBAND(J,M) * RETURNS: * trace(A) *************************************************************************/ double band_trace(A, J, M) register SYMBAND A; register int J, M; { register double trace = 0.; J++; for (; M; M--) trace += A[M][J]; return trace; } /****************************** trace_prod() ****************************** * TRACE OF PRODUCT OF TWO SYMMETRIC BAND MATRICES * INPUT: * A: SYMBAND(M-1,M) * B: SYMBAND(M-1,M) * RETURNS: * trace( A * B ) * FLOPS: **************************************************************************/ double trace_prod(A, B, M) register SYMBAND A, B; register int M; { register double trace = 0.; register int m, i; for (m=M; m >= 1; m--) { for (i=M; i >= m; i--) trace += A[i][m-i+M] * B[i][m-i+M]; for ( ; i >= 1; i--) trace += A[m][i-m+M] * B[m][i-m+M]; } return trace; } /****************************** mult_band_vect() ************************** * MULTIPLY A BAND MATRIX TIMES A VECTOR * INPUT: * A: SYMBAND(J,N) * x[1..N] * OUTPUT: * y[1..N] y = A x * FLOPS: N(2J+1)2J - ? **************************************************************************/ void mult_band_vect(y, A, x, J, N) register VECT y, x; register SYMBAND A; register int J, N; { register int n, k; /* ASSUMES J+J <= N */ if (J+J > N) EXITBUG("Band matrix too fat for this algorithm"); for (n=J; n >= 1; n--) { y[n] = A[n][J+1] * x[n]; for (k=J; k >= 1; k--) y[n] += A[n+k][-k+J+1] * x[n+k]; /* A(n,n+k) A(n+k,n) */ for (k=n-1; k >= 1; k--) y[n] += A[n][-k+ J+1] * x[n-k]; /* A(n, n-k) */ } for (n=J+1; n <= N; n++) { int k0 = (n+J > N) ? N-n : J; y[n] = A[n][J+1] * x[n]; for (k=k0; k >= 1; k--) y[n] += A[n+k][-k+J+1] * x[n+k]; /* A(n,n+k) A(n+k,n) */ for (k=J; k >= 1; k--) y[n] += A[n][-k+ J+1] * x[n-k]; /* A(n, n-k) */ } } /****************************** mult_full_vect() ************************** * MULTIPLY A FULL BAND MATRIX TIMES A VECTOR * INPUT: * A: SYMBAND(M-1,M) * z[1..M] * OUTPUT: * v[1..M] v = A z * FLOPS: M(2M-1) **************************************************************************/ void mult_full_vect(v, A, z, M) register VECT v, z; register SYMBAND A; register int M; { register int m, k; for (m=M; m >= 1; m--) { v[m] = A[m][M] * z[m]; for (k=m-1; k >= 1; k--) v[m] += A[m][k-m + M] * z[k]; for (k=m+1; k <= M; k++) v[m] += A[k][m-k + M] * z[k]; } } /****************************** mult_full_band() ************************** * MULTIPLY TWO FULL M*M MATRICES STORED IN BAND FORMAT * INPUT: * A: SYMBAND(M-1, M) * B: SYMBAND(M-1, M) * OUTPUT: * C: SYMBAND(M-1, M) C = A * B *************************************************************************/ void mult_full_band(C, A, B, M) register SYMBAND C, A, B; register int M; { register int i, j, k; register DATA *sum; for (i=M; i >= 1; i--) { for (j=i; j >= 1; j--) { *(sum = &C[i][j-i + M]) = 0.; for (k=1; k <= j; k++) *sum += A[i][k-i + M] * B[j][k-j + M]; for (k=j+1; k <= i; k++) *sum += A[i][k-i + M] * B[k][j-k + M]; for (k=i+1; k <= M; k++) *sum += A[k][i-k + M] * B[k][j-k + M]; } } } /************************************************************************* * SUBROUTINES FOR CONVERTING BETWEEN BAND AND COLUMN-WISE FORMATS. *************************************************************************/ /****************************** sym_band_to_mat() ************************ * CONVERT A SYMMETRIC BAND MATRIX B TO CMAT FORMAT MATRIX A * INPUT: * B: SYMBAND(J, N) * OUTPUT: * A: CMAT[1..NxN] *************************************************************************/ void sym_band_to_mat(A, B, J, N) CMAT A; SYMBAND B; int J, N; { register int i, j; for (i=1; i <= J; i++) for (j=1; j <= i; j++) A[(i-1)*N + j] = A[(j-1)*N + i] = B[i][j-i+J+1]; for (i=J+1; i <= N; i++) for (j=i-J; j <= i; j++) A[(i-1)*N + j] = A[(j-1)*N + i] = B[i][j-i+J+1]; } /****************************** mat_to_sym_band() ************************ * CONVERT A SYMMETRIC COLUMN-WISE MATRIX A TO BAND FORMAT MATRIX B * INPUT: * A: CMAT[1..NxN] * OUTPUT: * B: SYMBAND(J, N) *************************************************************************/ void mat_to_sym_band(B, A, J, N) CMAT A; SYMBAND B; int J, N; { register int i, j; for (j=J+1; j >= 1; j--) for (i=J+2-j; i <= N; i++) B[i][j] = A[(i-1)*N + i+j-(J+1)]; } End of smooth/band_util.c echo smooth/best_indiv.c 1>&2 cat >smooth/best_indiv.c <<'End of smooth/best_indiv.c' /*********************** best_indiv.c ************************************ * SMOOTH CURVES OPTIMALLY - INDIVIDUALLY * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" double best_indiv(); static COVAR *fake_covar(); static double best_single(); static NRCTYPE score_func(); extern int vs_setup_alloc(), vs_setup_ti(), vs_setup_yi(); extern int vs_setup_alpha(), vs_setup_cov(), vs_smooth(); extern int cv_setup_alloc(), cv_setup_covar(), cross_val(); extern void vs_cleanup(), cv_cleanup(); extern void mnbrak(); extern NRCTYPE brent(); extern COVAR *alloc_covar(); extern void covar_free(); static int score_id; /* * CURRENTLY ONLY IMPLEMENTED FOR iid CASE * INPUT: * alpha[1..M] initial values * ti[1..N] sample points * yi[1..NxM] NOISY MEASUREMENTS TO BE SMOOTHED * score_id: WHICH SCORE * tol: CONVERGENCE TOLERANCE * OUTPUT: * ys[1..NxM] best estimates * trace[1..M] trace of influence matrices * EFFECT: * alpha: best alpha values * RETURNS: sum of scores */ double best_indiv(ys, trace, alpha, ti, yi, l_score_id, tol, MM, NN) VECT ys, trace, ti, yi; ALPHA_SET alpha; int l_score_id; double tol; int MM, NN; { int m; COVAR *covar = fake_covar(); double score_sum = 0.; /* SETUP: INDEPENDENT OF ALPHA */ if (!vs_setup_alloc(1, NN)) EXITBUG("bi 1"); if (!vs_setup_ti(ti, NN)) EXITBUG("bi 2"); if (!vs_setup_yi(yi, 1, NN)) EXITBUG("bi 3"); if (!vs_setup_cov(covar, 1)) EXITBUG("bi 4"); if (!cv_setup_alloc()) EXITBUG("bi 5"); if (!cv_setup_covar()) EXITBUG("bi 6"); score_id = l_score_id; /* LOOP THROUGH INDIVIDUAL VECTOR COMPONENTS */ for (m=1; m <= MM; m++) { /* alpha+m-1 HERE SINCE best_score WANTS alpha[1..1] */ /* best_score(ys+(m-1)*NN, ti, yi+(m-1)*NN, covar, alpha+m-1, 1, 1, NN, score_id, tol); */ score_sum += best_single(ys+(m-1)*NN, &alpha[m], yi+(m-1)*NN, NN, tol); if (!trace_i_a(&trace[m], alpha+m-1)) EXITBUG("bi 7"); } covar_free(covar, 1); vs_cleanup(); cv_cleanup(); return score_sum; } /* * INPUT: * score_id: which score to do * alpha: initial value * tol: termination tolerance * OUTPUT: * alpha: optimal alpha * ys: best estimates * RETURNS: * best score */ static double best_single(ys, alpha, yi, NN, tol) ALPHA *alpha; VECT ys, yi; int NN; double tol; { double score; NRCTYPE ax, bx, cx, fa, fb, fc; /* "SAVE" THE PASSED VARIABLES */ if (!vs_setup_yi(yi, 1, NN) ) EXITBUG("bs 1"); /* THE ROUTINES BELOW ALL DEPEND ON ALPHA */ ax = *alpha-1; bx = *alpha; mnbrak(&ax, &bx, &cx, &fa, &fb, &fc, score_func); score = brent(ax, bx, cx, score_func, tol, alpha); if ( !vs_setup_alpha(alpha-1) ) EXITBUG("bs 2"); if ( !vs_smooth(ys, 1) ) EXITBUG("bs 3"); return score; } /* * RETURN PARTICULAR SCORE */ static NRCTYPE score_func(pass_alpha) ALPHA pass_alpha; { double ur, cv, gcv; NRCTYPE score = 0.; ALPHA alpha = pass_alpha; ALPHA_SET pseudo_alpha = &alpha-1; /* tricky... */ if ( !vs_setup_alpha(pseudo_alpha) ) EXITBUG("sf 1"); if ( !vs_smooth((VECT) 0, 0) ) EXITBUG("sf 2"); if ( !cross_val(&ur, &cv, &gcv, score_id) ) EXITBUG("sf 3"); if (score_id == JOBUR) score = ur; else if (score_id == JOBCV) score = cv; else if (score_id == JOBGCV) score = gcv; else EXITBUG("Invalid score id"); return score; } /* * SETUP FAKE 1x1 COVARIANCE MATRIX * RETURNS: * covar[1] */ static COVAR *fake_covar() { int n; COVAR *covar; /* FIRST SETUP A FAKE COVARIANCE MATRIX */ if (! (covar = alloc_covar(1, 1)) ) EXITBUG("fc 1"); for (n=1; n >= 1; n--) { DATA one = 1.0; VECT isig = &one - 1; if (!inv_to_covar(&covar[n], isig, 1)) EXITBUG("fc 2"); } return covar; } End of smooth/best_indiv.c echo smooth/best_mse.c 1>&2 cat >smooth/best_mse.c <<'End of smooth/best_mse.c' /*********************** best_mse.c ************************************** * VECTOR SPLINE SMOOTHING WITH SMALLEST MSE * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" double best_mse(); static NRCTYPE func_mse(); extern int vs_setup_alloc(), vs_setup_ti(), vs_setup_yi(); extern int vs_setup_cov(), vs_setup_alpha(), vs_smooth(); extern void vs_cleanup(); extern void powell(); static int MM, NN; static VECT ys, yr; /* * INPUT: * yr[1..NxM] true values * alpha[1..M] initial value * tol: termination tolerance * OUTPUT: * ys: best estimates * EFFECT: * alpha: optimal alpha * RETURNS: * best mse */ double best_mse(l_ys, l_yr, alpha, ti, yi, covar, flag_iid, M, N, tol) CMAT l_ys, l_yr; ALPHA_SET alpha; VECT ti; CMAT yi; COVAR *covar; int flag_iid, M, N; double tol; { NRCTYPE mse; DIRSET dirset; int iter; /* "SAVE" THE PASSED VARIABLES */ ys = l_ys; yr = l_yr; MM = M; NN = N; /* SETUP: INDEPENDENT OF ALPHA */ if (!vs_setup_alloc(MM, NN)) EXITBUG("1."); if (!vs_setup_ti(ti, NN)) EXITBUG("2."); if (!vs_setup_yi(yi, MM, NN)) EXITBUG("3."); if (!vs_setup_cov(covar, flag_iid)) EXITBUG("4."); /* THE ROUTINES BELOW ALL DEPEND ON ALPHA */ dirset = SETUP_DIRSET(MM); powell(alpha, dirset, MM, (NRCTYPE) tol, &iter, &mse, func_mse); vs_cleanup(); return (double) mse; } /* * RETURN PARTICULAR MSE */ static NRCTYPE func_mse(alpha) ALPHA_SET alpha; { int n; VECT pys = ys+1, pyr = yr+1; double mse = 0.; if ( !vs_setup_alpha(alpha) ) EXITBUG("8."); if ( !vs_smooth(ys, 1) ) EXITBUG("9."); for (n=NN*MM; n >= 1; n--, pys++, pyr++) mse += (*pys - *pyr) * (*pys - *pyr); return (NRCTYPE) mse / NN; } End of smooth/best_mse.c echo smooth/best_score.c 1>&2 cat >smooth/best_score.c <<'End of smooth/best_score.c' /*********************** best_score.c ************************************ * SMOOTH VECTOR MEASUREMENTS WITH AUTOMATIC SCORE CRITERIA * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" double best_score(); static NRCTYPE score_func(); extern int vs_setup_alloc(), vs_setup_ti(), vs_setup_yi(); extern int vs_setup_alpha(), vs_setup_cov(), vs_smooth(); extern int cv_setup_alloc(), cv_setup_covar(), cross_val(); extern void vs_cleanup(), cv_cleanup(); extern void powell(); static int MM, NN, score_id; /* * INPUT: * score_id: which score to do * alpha: initial value * tol: termination tolerance * OUTPUT: * alpha: optimal alpha * ys: best estimates * RETURNS: * best score */ double best_score(ys, alpha, ti, yi, covar, flag_iid, M, N, l_score_id, tol) CMAT ys; ALPHA_SET alpha; VECT ti; CMAT yi; COVAR *covar; int flag_iid, M, N, l_score_id; double tol; { NRCTYPE score; DIRSET dirset; int iter; /* "SAVE" THE PASSED VARIABLES */ score_id = l_score_id; MM = M; NN = N; /* SETUP: INDEPENDENT OF ALPHA */ if (!vs_setup_alloc(MM, NN)) EXITBUG("1."); if (!vs_setup_ti(ti, NN)) EXITBUG("2."); if (!vs_setup_yi(yi, MM, NN)) EXITBUG("3."); if (!vs_setup_cov(covar, flag_iid)) EXITBUG("4."); if (!cv_setup_alloc()) EXITBUG("5."); if (!cv_setup_covar()) EXITBUG("6."); /* THE ROUTINES BELOW ALL DEPEND ON ALPHA */ dirset = SETUP_DIRSET(MM); powell(alpha, dirset, MM, (NRCTYPE) tol, &iter, &score, score_func); if ( !vs_setup_alpha(alpha) ) EXITBUG("6.2"); if ( !vs_smooth(ys, 1)) EXITBUG("6.5"); vs_cleanup(); cv_cleanup(); return (double) score; } /* * RETURN PARTICULAR SCORE */ static NRCTYPE score_func(alpha) ALPHA_SET alpha; { double ur, cv, gcv; double score; if ( !vs_setup_alpha(alpha) ) EXITBUG("8."); if ( !vs_smooth((VECT) 0, 0)) EXITBUG("9."); if ( !cross_val(&ur, &cv, &gcv, score_id) ) EXITBUG("10."); if (score_id == JOBUR) score = ur; else if (score_id == JOBCV) score = cv; else if (score_id == JOBGCV) score = gcv; else EXITBUG("Invalid score id."); message0("cv ", cv); PRINT_ALPHA(alpha, MM); return score; } End of smooth/best_score.c echo smooth/covar.c 1>&2 cat >smooth/covar.c <<'End of smooth/covar.c' /***************************** covar.c *********************************** * UTILITY SUBROUTINES FOR COVAR STRUCTURES * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" COVAR *alloc_covar(); void covar_free(); COVAR *inv_to_covar_set(); void copy_covar(), scale_covar(); extern void mat_to_sym_band(); extern int chol_sym_band(); extern void invert_sym_band(); extern double band_trace(); /* * ALLOCATE SPACE FOR A SET OF N COVARIANCE STRUCTURES * WITH COVARIANCE SIZE MxM */ COVAR *alloc_covar(M, N) int M, N; { COVAR *cov = (COVAR *) ALLOC(N, sizeof(COVAR)); if (!cov) return (COVAR *) 0; cov--; for (; N >= 1; N--) { if ( !(cov[N].cov = ALLOC_MATRIX(M, M)) ) return (COVAR *) 0; if ( !(cov[N].cov2 = ALLOC_MATRIX(M, M)) ) return (COVAR *) 0; if ( !(cov[N].inv = ALLOC_MATRIX(M, M)) ) return (COVAR *) 0; if ( !(cov[N].Linv = ALLOC_MATRIX(M, M)) ) return (COVAR *) 0; if ( !(cov[N].Dinv = ALLOC_VECTOR(M)) ) return (COVAR *) 0; } return cov; } void covar_free(cov, N) COVAR *cov; int N; { for (; N >= 1; N--) { FREE_MATRIX(cov[N].cov); FREE_MATRIX(cov[N].cov2); FREE_MATRIX(cov[N].inv); FREE_MATRIX(cov[N].Linv); FREE_VECTOR(cov[N].Dinv); } FREE(cov+1); } void copy_covar(src, dst, MM) COVAR *src, *dst; int MM; { COPY_MATRIX(src->cov, dst->cov, MM, MM); COPY_MATRIX(src->cov2, dst->cov2, MM, MM); COPY_MATRIX(src->inv, dst->inv, MM, MM); COPY_MATRIX(src->Linv, dst->Linv, MM, MM); COPY_VECTOR(src->Dinv, dst->Dinv, MM); } /* * MULTIPLY covariance by scale - correcting others */ void scale_covar(cov, scale, MM) COVAR *cov; double scale; int MM; { SCALE_MATRIX(cov->cov, scale, MM, MM); SCALE_MATRIX(cov->cov2, scale*scale, MM, MM); SCALE_MATRIX(cov->inv, 1. / scale, MM, MM); SCALE_VECTOR(cov->Dinv, scale, MM); } /****************************** inv_to_covar() *************************** * CONVERT INVERSE OF COVARAINCE MATRIX INTO COVAR STRUCTURE * INPUT: * inv: CMAT[1..MxM] * OUTPUT: * cov: * RETURNS: 1 if OK, 0 if singular *************************************************************************/ int inv_to_covar(cov, inv, M) COVAR *cov; CMAT inv; int M; { mat_to_sym_band(cov->inv, inv, M-1, M); if ( !chol_sym_band(cov->Linv, cov->Dinv, cov->inv, M-1, M) ) return 0; invert_sym_band(cov->cov, cov->Linv, cov->Dinv, M-1, M); cov->trace = band_trace(cov->cov, M-1, M); return 1; } /****************************** inv_to_covar_set() *********************** * CONVERT SET OF COVARIANCE MATRIX INVERSES INTO COVAR STRUCTURE * INPUT: * isigma[1..MxMx(Nc)] covariance inverses * flag_iid Nc = 1 or N ? * RETURN: * cov[1..(Nc)]: or 0 if problem *************************************************************************/ COVAR *inv_to_covar_set(isigma, M, N, flag_iid) CMAT isigma; int M, N, flag_iid; { int nn, Nc = flag_iid ? 1 : N; COVAR *covar; if ( !(covar = alloc_covar(M, Nc)) ) return 0; for (nn=Nc; nn >= 1; nn--) if (!inv_to_covar(&covar[nn], isigma+(nn-1)*M, M)) return 0; return covar; } End of smooth/covar.c echo smooth/cross_val.c 1>&2 cat >smooth/cross_val.c <<'End of smooth/cross_val.c' /*********************** cross_val.c ************************************* * SUBROUTINES FOR COMPUTING STATISTICS FOR * FIXED INTERVAL SMOOTHING OF VECTOR MEASUREMENTS WITH CUBIC SPLINES. * TO BE USED ONLY AFTER SMOOTHING WITH THE SUBROUTINES IN vs_smooth.c. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" int cv_setup_alloc(), cv_setup_covar(), cross_val(), trace_i_a(); void cv_cleanup(); extern double band_trace(); extern void invert_sym_band(), mult_full_band(); extern void select_nth_qbiqt(); extern int chol_sym_band(); extern double errn_term(), rss_term(), ur_term(), cv_term(), gcv_term(); /* VARIABLES DECLARED, ALLOCATED, AND PRECOMPUTED IN vs_smooth.c */ extern int NG, MG; /* 'G' SUFFIX FOR "GLOBAL" */ extern int NI, MI; /* 'I' SUFFIX FOR "INITIALIZED" */ extern int Flag_iid; extern BAND Qt, BL; extern VECT BD; extern OMATRIX sig_q_c; #define DELTA sig_q_c /* DELTA[1..N][1..M] = (yi-yhat)*/ /* POINTERS INITIALIZED IN vs_smooth.c */ extern COVAR *Gcovar; /* GLOBAL VARIABLES FOR THIS FILE */ static SYMBAND BI; /* BI(3M-1,M(N-2)) */ static SYMBAND Fn, FL; /* F, FL: SYMBAND(M-1,M) */ static VECT FD; /* FD[1..N] */ static double trace_pi = 0.; /* trace(Pi) */ /* * ALLOCATE MATRICES */ int cv_setup_alloc() { if ((BI = ALLOC_BAND(MI*(NI-2), 3*MI-1)) == 0) return 0; if ((Fn = ALLOC_MATRIX(MI, MI)) == 0) return 0; if ((FL = ALLOC_MATRIX(MI, MI)) == 0) return 0; if ((FD = ALLOC_VECTOR(MI)) == 0) return 0; return 1; } /* * PRECOMPUTE TERMS USED FOR CROSS VALIDATION AND OTHER SCORES * USES: * Gcovar */ int cv_setup_covar() { register int n; for (n = Flag_iid ? 1 : NG; n >= 1; n--) { COVAR *covar = &Gcovar[n]; covar->trace = band_trace(covar->cov, MG-1, MG); mult_full_band(covar->cov2, covar->cov, covar->cov, MG); trace_pi += covar->trace; } return 1; } /* * FREE ALL THE ALLOCATED MATRICES */ void cv_cleanup() { FREE_BAND(BI); FREE_MATRIX(Fn); FREE_MATRIX(FL); FREE_VECTOR(FD); } /* * COMPUTE CROSS VALIDATION SCORES AND OTHER STATISTICS * * CV = 1/N \sum_{n=1}^N | Sig_n^{-1/2} inv(I-A) (Y - Yhat) * UR = (1/N) [ | Y - Yhat |^2 - 2 trace(Sig (I-A)) + trace(Sig) ] * GCV = (1/N) RSS / [ (1/N) trace(I - A) ]^2 * * I - A = Pi F; F = (Q @ I) BI (Q' @ I); BI = inv(BL BD BL') * I - A_n = Sig_n F_n * e_n = Y_n - Yhat_n * CV_n = 1/N \sum b_n Sig_n^{-1/2} b_n, b_n = inv(Sig_n F_n) e_n * * INPUT: * job: which scores to compute * OUTPUT: * ur: (1) UNBIASED RISK SCORE * cv: (2) CROSS VALIDATION SCORE * gcv: (4) GENERALIZED CV SCORE * USES: * DELTA (sig_q_c) COMPUTED BY vect_smooth() * Gcovar error covariance structure(s) */ int cross_val(ur, cv, gcv, job) double *ur, *cv, *gcv; int job; { int nn; COVAR *covar = &Gcovar[1]; double rss = 0., errn = 0.; /* COMPUTE THE CENTRAL BANDS OF THE INVERSE MATRIX */ invert_sym_band(BI, BL, BD, 3*MG-1, MG*(NG-2)); *cv = *ur = *gcv = 0.; for (nn=1; nn <= NG; nn++) { VECT err = DELTA[nn]; select_nth_qbiqt(Fn, BI, Qt, MG, NG, nn); if ( !chol_sym_band(FL, FD, Fn, MG-1, MG) ) return 0; /* COMPUTE NORM OF RESIDUALS */ if (job & JOBUR) { errn += errn_term(err, MG); } /* COMPUTE RESIDUAL SUM OF SQUARES */ if (job & JOBGCV) { rss += rss_term(err, covar->inv, MG); } /* COMPUTE UNBIASED RISK ESTIMATOR */ if (job & JOBUR) { *ur += ur_term(Fn, covar->cov2, MG); } /* COMPUTE CV SCORE */ if (job & JOBCV) { *cv += cv_term(err, FL, FD, covar->inv, MG); } /* COMPUTE GENERALIZED CROSS VALIDATION SCORE */ if (job & JOBGCV) { *gcv += gcv_term(Fn, covar->cov, MG); } /* IF NOT IID NOISE, INCREMENT TO NEXT COVARIANCE */ if (!Flag_iid) covar++; } /* CORRECT SCORES FOR CONSTANT FACTORS (SAVES COMPUTATIONS) */ if (job & JOBUR) *ur = (errn - 2. * *ur + trace_pi) / NG; if (job & JOBCV) *cv /= NG; if (job & JOBGCV) *gcv = NG * rss / (*gcv * *gcv); return 1; } /* * COMPUTE TRACE (I - A) * INPUT: * alpha[1..M] smoothing parameters * OUTPUT: * trace: trace(I - A(alpha)) * USES: * Gcovar */ int trace_i_a(trace, alpha) double *trace; ALPHA_SET alpha; { int nn; COVAR *covar = &Gcovar[1]; if (!vs_setup_alpha(alpha)) return 0; if (!vs_transfer()) return 0; /* COMPUTE THE CENTRAL BANDS OF THE INVERSE MATRIX */ invert_sym_band(BI, BL, BD, 3*MG-1, MG*(NG-2)); *trace = 0.; for (nn=1; nn <= NG; nn++) { select_nth_qbiqt(Fn, BI, Qt, MG, NG, nn); if ( !chol_sym_band(FL, FD, Fn, MG-1, MG) ) return 0; *trace += gcv_term(Fn, covar->cov, MG); /* IF NOT IID NOISE, INCREMENT TO NEXT COVARIANCE */ if (!Flag_iid) covar++; } return 1; } End of smooth/cross_val.c echo smooth/cv_score.c 1>&2 cat >smooth/cv_score.c <<'End of smooth/cv_score.c' /***************************** cv_score.c ******************************** * SUBROUTINES FOR COMPUTING CROSS-VALIDATION LIKE SCORES. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" void select_nth_qbiqt(); double errn_term(), rss_term(), ur_term(), cv_term(), gcv_term(); extern void mult_full_vect(), solve_sym_band(); extern double trace_prod(); /***************************** select_nth_qbiqt() ************************ * COMPUTE THE n'TH CENTRAL M*M INFLUENCE MATRIX OF * Fnn = [(Q @ I_M) IB (Qt @ I_M)]nn * INPUT: * IB: SYMBAND(3M-1, M(N-2)) central inverse * Qt: see vs_init.c * n: WHICH ONE TO DO (1,..., N). * OUTPUT: * F: SYMBAND(M-1, M) FULL SYMMETRIC M*M MATRIX * FLOPS: M*(M+1)/2* 3 * (3*3) (1*1) for n=1, (2*2) for n=2 * all N calls: N(27/2M(M+1)) - 39/2M(M+1) *************************************************************************/ void select_nth_qbiqt(F, IB, Qt, M, N, n) SYMBAND F; SYMBAND IB, Qt; int n, M, N; { int i, j, m, k, n1, nN, i1, i2; double *sum; int J1 = 3*M-1 + 1; n1 = (n-2 <= 1) ? 1 : n-2; nN = (n <= N-2) ? n : N-2; for (i=1; i <= M; i++) /* LOWER TRIANGLE OF F */ { for (j=1; j <= i; j++) { sum = &F[0*M+i][j-i+M]; /* 0 instead of n-1 */ *sum = 0.; for (k=n1; k <= nN; k++) for (m=n1; m <= nN; m++) { i1 = (k-1)*M+i; i2 = (m-1)*M+j; *sum += Qt[m][3-(n-m)] * Qt[k][3-(n-k)] * ((i2 < i1) ? IB[i1][i2-i1+J1] : IB[i2][i1-i2+J1]); } } } } /*********************** errn_term() ************************************* * COMPUTE A SINGLE TERM IN THE SUM OF SQUARED ERROR * errn_n = | err_n |^2 * INPUT: * err[1..M] nth error term = yi^n - ys^n * RETURNS: * errn_n * FLOPS: *************************************************************************/ double errn_term(err, M) VECT err; int M; { return NORM_VECTOR(err, M); } /********************* rss_term() **************************************** * COMPUTE A SINGLE TERM IN THE RESIDUAL SUM OF SQUARES * rss_n = err_n' * inv(cov_n) * err_n * INPUT: * err[1..M] nth error term = yi^n - ys^n * inv: SYMBAND(M-1,M) inverse of covariance_n * RETURNS: * rss_n * FLOPS: *************************************************************************/ double rss_term(err, inv, M) VECT err; SYMBAND inv; int M; { register double rss; VECT werr; TMP_ALLOC_VECT(werr, M) mult_full_vect(werr, inv, err, M); rss = DOT_VECTOR(err, werr, M); TMP_FREE_VECT(werr); return rss; } /********************* ur_term() ***************************************** * COMPUTE A SINGLE TERM IN THE UNBAISED RISK SUMMATION. * UR = (1/N) [ || err ||^2 - 2 * trace(cov * (I - A)) + trace(cov) ] * ur_n = trace(cov_n * cov_n * Fn) * INPUT: * Fn: SYMBAND(M-1,M) nth qbiqt term * cov2: SYMBAND(M-1,M) covariance_n^2 * RETURNS: * ur_n: (part of) nth term of ur score * FLOPS: *************************************************************************/ double ur_term(Fn, cov2, M) SYMBAND Fn, cov2; int M; { return trace_prod(cov2, Fn, M); } /********************* cv_term() ***************************************** * COMPUTE A SINGLE TERM IN THE CROSS VALIDATION SUMMATION. * CV = (1/N) sum_1^N || sign^(-1/2) inv(I - An) err_n ||^2 * cv = werr' * inv(cov) * werr, Fn werr = inv(cov) * err * INPUT: * err[1..M] nth error term = yi^n - ys^n * FL: SYMBAND(M-1,M) * FD[1..M] Fn = FL' FD FL * inv: SYMBAND(M-1,M) inverse of covariance_n * RETURNS: * nth term of cv score * FLOPS: *************************************************************************/ double cv_term(err, FL, FD, inv, M) VECT err; SYMBAND FL, inv; VECT FD; int M; { double sum; VECT werr, invwerr; TMP_ALLOC_VECT(werr, M) TMP_ALLOC_VECT(invwerr, M) mult_full_vect(werr, inv, err, M); solve_sym_band((VECT) 0, FL, FD, werr, M-1, M, 1); mult_full_vect(invwerr, inv, werr, M); sum = DOT_VECTOR(invwerr, werr, M); TMP_FREE_VECT(werr); TMP_FREE_VECT(invwerr); return sum; } /********************* gcv_term() **************************************** * COMPUTE A SINGLE TERM IN THE GENERALIZED CROSS VALIDATION SUMMATION. * GCV = (1/N) RSS / [ (1/N) trace(I - A) ]^2 * den_gcv_n = trace(I - An) = trace(cov_n * Fn) * INPUT: * Fn: SYMBAND(M-1,M) * cov: SYMBAND(M-1,M) covariance_n * RETURNS: * den_gcv_n: nth term of denominator of gcv score * FLOPS: *************************************************************************/ double gcv_term(Fn, cov, M) SYMBAND Fn, cov; int M; { return trace_prod(cov, Fn, M); } End of smooth/cv_score.c echo smooth/example1.c 1>&2 cat >smooth/example1.c <<'End of smooth/example1.c' /*********************** example1.c ************************************** * EXAMPLE FRONT END ROUTINE FOR USING VECTOR-SPLINE SUBROUTINES TO * SMOOTH VECTOR MEASUREMENTS, COMPUTE CV SCORES, ETC. * THIS EXAMPLE DOES NOT REQUIRE THE numrecC LIBRARY. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" int vspline(), vect_smooth(); extern COVAR *inv_to_covar_set(); extern COVAR *alloc_covar(); extern void covar_free(); extern int vs_setup_alloc(), vs_setup_ti(), vs_setup_yi(); extern int vs_setup_alpha(), vs_setup_cov(); extern int cv_setup_alloc(), cv_setup_covar(); extern int vs_smooth(), cross_val(), penalty(); extern void vs_cleanup(), cv_cleanup(); extern double roughness(); /* * VECTOR SPLINE SMOOTH AND COMPUTE SCORES * INPUT: * alpha[1..M] smoothing parameters * ti[1..N] sample points * yi[1..NxM] measurements * isigma[1..MxMx(N)] inverse of measurement error covariances * flag_iid 0: one covariance, 1: N covariances * OUTPUT: * ys[1..NxM] smoothed estimates * RETURNS: 1 if successful, 0 if error */ int vspline(ys, alpha, ti, yi, isigma, flag_iid, MM, NN) CMAT ys; ALPHA_SET alpha; VECT ti; CMAT yi, isigma; int flag_iid, MM, NN; { COVAR *covar; /* convert the flattened covariance inverses to structure */ if ( !(covar = inv_to_covar_set(isigma, MM, NN, flag_iid)) ) return 0; if ( !vs_setup_alloc(MM, NN)) return 0; if ( !vs_setup_ti(ti, NN)) return 0; if ( !vs_setup_yi(yi, MM, NN)) return 0; if ( !vs_setup_cov(covar, flag_iid)) return 0; if ( !vs_setup_alpha(alpha)) return 0; if ( !vs_smooth(ys, 1) ) return 0; vs_cleanup(); covar_free(covar, flag_iid ? 1 : NN); return 1; } /* * VECTOR SPLINE SMOOTH AND COMPUTE SCORES * INPUT: * alpha[1..M] smoothing parameters * ti[1..N] sample points * yi[1..NxM] measurements * covar[1..?] measurement error covariances * flag_iid covariance varies? * OUTPUT: * ys[1..NxM] smoothed estimates * ur unbiased risk score * cv cross validation score * gcv generalized CV score * penin initial roughness penalty * penout final roughness penalty * rough final roughness penalty */ int vect_smooth(ys, ur, cv, gcv, penin, penout, rough, alpha, ti, yi, covar, flag_iid, MM, NN) CMAT ys; double *ur, *cv, *gcv, *penin, *penout, *rough; ALPHA_SET alpha; VECT ti; CMAT yi; COVAR *covar; int flag_iid, MM, NN; { /* SETUP: INDEPENDENT OF ALPHA */ if ( !vs_setup_alloc(MM, NN)) return 0; if ( !vs_setup_ti(ti, NN)) return 0; if ( !vs_setup_yi(yi, MM, NN)) return 0; if ( !vs_setup_cov(covar, flag_iid)) return 0; if ( !cv_setup_alloc()) return 0; if ( !cv_setup_covar()) return 0; /* THE ROUTINES BELOW ALL DEPEND ON ALPHA */ if ( !vs_setup_alpha(alpha)) return 0; if ( !vs_smooth(ys, 1)) return 0; if ( !cross_val(ur, cv, gcv, JOBUR | JOBCV | JOBGCV)) return 0; *rough = roughness(); if ( !penalty(penout, alpha, ys, MM, NN) ) return 0; if ( !penalty(penin, alpha, yi, MM, NN) ) return 0; vs_cleanup(); cv_cleanup(); return 1; } End of smooth/example1.c echo smooth/example2.c 1>&2 cat >smooth/example2.c <<'End of smooth/example2.c' /*********************** example2.c ************************************** * EXAMPLE FRONT END ROUTINE FOR USING VECTOR-SPLINE SUBROUTINES TO * SMOOTH VECTOR MEASUREMENTS, COMPUTE CV SCORES, ETC. * THIS EXAMPLE CHOOSES THE SMOOTHING PARAMETER AUTOMATICALLY, * AND RELIES ON THE numrecC LIBRARY FUNCTION powell() TO DO THE * MINIMIZATION. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" int bestcv(); extern double best_score(); extern COVAR *inv_to_covar_set(); extern void covar_free(); /* * USE CV SCORE TO DETERMINE SMOOTHING PARAMETER * alpha: initial value * all others as in vspline() above * OUTPUT: * alpha: optimal alpha * ys: best estimates * RETURNS: 1 if ok, 0 if not */ #define TOL 0.01 int bestcv(ys, alpha, ti, yi, isigma, flag_iid, MM, NN) CMAT ys; ALPHA_SET alpha; VECT ti; CMAT yi, isigma; int flag_iid, MM, NN; { COVAR *covar; if ( !(covar = inv_to_covar_set(isigma, MM, NN, flag_iid)) ) return 0; best_score(ys, alpha, ti, yi, covar, flag_iid, MM, NN, JOBCV, 0.01); covar_free(covar, flag_iid ? 1 : NN); return 1; } End of smooth/example2.c echo smooth/interp.c 1>&2 cat >smooth/interp.c <<'End of smooth/interp.c' /************************** interp.c ************************************* * CUBIC-SPLINE INTERPOLATION * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" int cubic_interp(); static int dumb_search(); static void simple_interp(); extern void spline_mat(); extern void form_qt_Y(); extern int chol_sym_band(); extern void solve_sym_band(); /* * INTERPOLATION * INPUT: * ti[1..N] knots * yi[1..N] function value at knots * xs[1..Nx] new sample points * OUTPUT: * ys[1..Nx] interpolated values */ int cubic_interp(ys, ti, yi, xs, NN, Nx) VECT ys, ti, yi, xs; int NN, Nx; { int ix; BAND Qt, R, RL; VECT RD, cc; if (NN == 2) { simple_interp(ys, ti, yi, xs, Nx); return 1; } TMP_ALLOC_BAND(Qt, NN-2, 2); TMP_ALLOC_BAND(R, NN-2, 1); TMP_ALLOC_BAND(RL, NN-2, 1); TMP_ALLOC_VECT(RD, NN-2); TMP_ALLOC_VECT(cc, NN-2); /* COMPUTE 'C' VECTOR */ spline_mat(R, Qt, ti, NN); form_qt_Y(cc, Qt, yi, 1, NN); if (!chol_sym_band(RL, RD, R, 1, NN-2)) return 0; solve_sym_band((VECT) 0, RL, RD, cc, 1, NN-2, 1); for (ix=1; ix <= Nx; ix++) { double xx = xs[ix]; /* FIX: WRONG END CONDITIONS */ if (xx >= ti[NN]) { double hn1 = ti[NN] - ti[NN-1]; double bn1= (yi[NN]-yi[NN-1])/hn1 + hn1/6.*cc[NN-2]; ys[ix] = yi[NN] + bn1 * (xx - ti[NN]); } else if (xx <= ti[1]) { double h1 = (ti[2]-ti[1]); double b1 = (yi[2]-yi[1])/h1 - h1/6.*cc[1]; ys[ix] = yi[1] + b1*(xx - ti[1]); } else { /* 1 <= it <= NN-1 */ int it = dumb_search(xx, ti); double dt = xx - ti[it]; double hi = ti[it+1] - ti[it]; double ai = yi[it]; double ci = (it > 1) ? cc[it-1] : 0; double ci1 = (it+1 < NN) ? cc[it] : 0; double di = (ci1 - ci) / hi; double bi = (yi[it+1]-ai)/hi - (ci1 + 2*ci) * hi / 6.; ys[ix] = ai + (bi + (ci/2. + di/6.*dt)*dt)*dt; } } TMP_FREE_BAND(Qt); TMP_FREE_BAND(R); TMP_FREE_BAND(RL); TMP_FREE_VECT(RD); TMP_FREE_VECT(cc); return 1; } /* * return i such that list[i] <= x < list[i+1] */ static int dumb_search(x, list) double x; VECT list; { int ii; for (ii=2; ;ii++) if (x < list[ii]) return ii-1; } /* * linear interpolation between (ti[1],yi[1]) and (ti[2],yi[2]). */ static void simple_interp(ys, ti, yi, xs, Nx) VECT ys, ti, yi, xs; int Nx; { int ix; double dx = ti[2] - ti[1]; double dy = yi[2] - yi[1]; for (ix=1; ix <= Nx; ix++) ys[ix] = yi[1] + (xs[ix] - ti[1])*dy/dx; } End of smooth/interp.c echo smooth/smooth.h 1>&2 cat >smooth/smooth.h <<'End of smooth/smooth.h' /***************************** smooth.h ********************************** * Declarations for vector-spline smoothing algorithm. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "defs.h" #include "matrix.h" /********************** BAND MATRIX FORMAT ******************************* * IF B IS A NxN SYMMETRIC BAND MATRIX WITH 2*J+1 NON-ZERO DIAGONALS * THEN I AM USING THE "TRANSPOSE" OF THE LINPACK BAND MATRIX FORMAT. * (IE, THE GOLUB - VAN LOAN FORMAT). * SUCH A MATRIX IS DENOTED BAND(J,N) IN COMMENTS. * * B(i,j) = B[i][j-i+J+1], i >= j * * B LOOKS LIKE: * | 0 0 b(1,1) | * | 0 b(2,1) b(2,2) | * | . : | * | . : | * | b(J+1,1) ... b(J+1,J+1) | * | : : | * | b(N,N-J) ... b(N,N) | * * SO THE COLUMNS OF B[][] ARE THE DIAGONALS OF B() * AND THE ROWS OF B[][] ARE THE (LOWER BAND) ROWS OF B(). *********************** BAND MATRIX FORMAT ******************************/ #define NRCTYPE float /* Num. Rec. in C subroutines */ #ifdef DoublePrec typedef double DATA; #define ALLOC_MATRIX ALLOC_MATRIX_D #define ALLOC_VECTOR ALLOC_VECTOR_D #define TMP_ALLOC_MATRIX TMP_ALLOC_MATRIX_D #define TMP_FREE_MATRIX TMP_FREE_MATRIX_D #define SCALE_VECTOR SCALE_VECTOR_D #define SCALE_MATRIX SCALE_MATRIX_D #define INC_SCALED_VECTOR INC_SCALED_VECTOR_D #define DOT_VECTOR DOT_VECTOR_D #define NORM_VECTOR NORM_VECTOR_D #define SUB_VECTOR SUB_VECTOR_D #define DEC_VECTOR DEC_VECTOR_D #define INC_VECTOR INC_VECTOR_D #define PRINT_VECTOR PRINT_VECTOR_D #define PRINT_MATRIX PRINT_MATRIX_D #define WRITE_VECTOR WRITE_VECTOR_D #define FPRINT_MATRIX FPRINT_MATRIX_D #define SCAN_VECTOR SCAN_VECTOR_D #define SCAN_MATRIX SCAN_MATRIX_D #define READ_VECTOR READ_VECTOR_D #define FSCAN_VECTOR FSCAN_VECTOR_D #else NOT IMPLEMENTED SINCE DOUBLE PRECISION SEEMS BEST! #endif #define ALLOC_BAND(n,j) ALLOC_MATRIX((n),(j+1)) #define FREE_BAND FREE_MATRIX #define COPY_BAND(src,dst,n,j) COPY_MATRIX(src,dst,n,(j+1)) #define ZERO_BAND(mat,n,j) ZERO_MATRIX(mat,(n),(j+1)) #define PRINT_BAND(m,n,j) PRINT_MATRIX(m,(n),(j+1)) #define TMP_ALLOC_VECT TMP_ALLOC_VECTOR_D #define TMP_FREE_VECT TMP_FREE_VECTOR_D #define TMP_ALLOC_BAND(m,n,j) TMP_ALLOC_MATRIX(m,n,(j+1)) #define TMP_FREE_BAND(m) TMP_FREE_MATRIX(m) #define JOBUR 1 /* CONTROL BITS FOR cross_val() */ #define JOBCV 2 #define JOBGCV 4 typedef DATA **OMATRIX; /* ORDINARY MATRIX [1..rows][1..cols] */ #define SYMBAND OMATRIX /* SYMMETRIC BAND MATRIX */ #define BAND SYMBAND /* BAND MATRIX */ typedef DATA *VECT; /* VECTOR[1..N] */ #define CMAT VECT /* MATRIX STORED COLUMN-WISE AS A VECT */ #define ALPHA NRCTYPE typedef ALPHA *ALPHA_SET; #define ALLOC_ALPHA ALLOC_VECTOR_F #define PRINT_ALPHA PRINT_VECTOR_F #define SCAN_ALPHA SCAN_VECTOR_F #define FSCAN_ALPHA FSCAN_VECTOR_F #define SETUP_DIRSET identity_f typedef NRCTYPE **DIRSET; /************************************************************************* * Structure for a MxM covariance matrix and its inverse. * (SYMMETRIC, POSITIVE DEFINITE) * inv(cov) = inv = Linv Dinv Linv' * For nonlinear problems, inv is typically the Hessian. * cov2 = cov * cov (used for unbiased risk score). *************************************************************************/ typedef struct { SYMBAND cov; /* SYMBAND(M-1,M) */ SYMBAND cov2; /* SYMBAND(M-1,M) */ SYMBAND inv; /* SYMBAND(M-1,M) */ SYMBAND Linv; /* SYMBAND(M-1,M) */ VECT Dinv; /* [1..M] */ DATA trace; /* trace(cov) */ } COVAR; End of smooth/smooth.h echo smooth/sym_band.c 1>&2 cat >smooth/sym_band.c <<'End of smooth/sym_band.c' /****************************** sym_band.c ******************************* * SUBROUTINES FOR SYMMETRIC BAND MATRIX MANIPULATIONS * SEE band.h FOR BAND MATRIX STORAGE FORMAT. * NOTE: BAND(N-1,N) MEANS THESE ROUTINES WORK FOR "FULL" MATRICES. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include "smooth.h" void kron_diag(); /* void kron_sym_band(); */ int chol_sym_band(); void solve_sym_band(); void invert_sym_band(); static int istiny(); static double ddot123(); static void dotdiff(); /**************************** kron_diag() ******************************* * TENSOR PRODUCT OF SYMMETRIC BAND MATRIX WITH DIAGONAL MATRIX * INPUT: * B SYMBAND(J,N) * D[1..M] DIAGONAL OF DIAGONAL MATRIX * OUTPUT: * K SYMBAND((J+1)M-1,NM) K = B @ Diag(D) * FLOPS: M[(N-J-1)(J+1) + (J+2)(J+1)/2] = M[N(J+1) - J(J+1)/2] * IF J=1, N<=N-2, THEN N(2M) - 5M ************************************************************************/ void kron_diag(K, B, J, N, D, M) SYMBAND K, B; VECT D; int J, N, M; { register int i, j, m, nk; nk = (J+1) * M; for (m=M; m >= 1; m--) /* UP THE DIAGONAL */ { for (i=J+1; i >= 1; i--) for (j=i; j >= 1; j--) K[(i-1)*M+m][(j-i)*M+nk] = B[i][j-i+J+1] * D[m]; for (i=J+2; i <= N; i++) for (j=i-J; j <= i; j++) K[(i-1)*M+m][(j-i)*M+nk] = B[i][j-i+J+1] * D[m]; } } #ifdef NOTUSED /**************************** kron_sym_band() *************************** * KRONECKER TENSOR PRODUCT OF SYMMETRIC BAND MATRIX WITH SYMMETRIC * FULL MATRIX. * INPUT: * B SYMBAND(J,N) * A[1..MM] M*M MATRIX STORED COLUMN-WISE. * OUTPUT: * K SYMBAND((J+1)M-1,NM) * FLOPS: MM[(N-J-1)J + J(J+1)/2] +NM((M-1)/2 = N[MMJ+M(M-1)/2]-MM/2(JJ+J) * OLD: M*M*[(J+1)*(J+2)/2 + (N-J)*(J+1)] * IF J=1, N<=N-2, THEN N(3/2MM+M/2)-4MM-M ************************************************************************/ void kron_sym_band(K, B, J, N, A, M) SYMBAND K, B; CMAT A; int J, N, M; { register int i, j, m, n, nk; nk = (J+1) * M; for (m=M; m >= 1; m--) /* DOWN THE MATRIX A */ { /* THE BELOW DIAGONAL BLOCKS */ for (n=M; n >= 1; n--) /* ACROSS THE MATRIX A */ { for (i=J+1; i >= 1; i--) for (j=i-1; j >= 1; j--) K[(i-1)*M+m][(j-i)*M+(n-m)+nk] = B[i][j-i+J+1] * A[(n-1)*M+m]; for (i=J+2; i <= N; i++) for (j=i-J; j < i; j++) K[(i-1)*M+m][(j-i)*M+(n-m)+nk] = B[i][j-i+J+1] * A[(n-1)*M+m]; } /* THE ON DIAGONAL BLOCKS */ for (n=m; n >= 1; n--) for (i=N; i >= 1; i--) K[(i-1)*M+m][ (n-m)+nk] = B[i][ J+1] * A[(n-1)*M+m]; } } #endif NOTUSED /****************************** chol_sym_band() ************************** * CHOLESKY DECOMPOSITION OF A NxN SYMMETRIC BAND MATRIX B (2*J+1) * J MUST BE AT LEAST 1 FOR THIS TO WORK * INPUT: * B: SYMBAND(J,N) -> L * diag(D) * trans(L) * OUTPUT: * L: BAND(J, N) UNIT LOWER TRIANGULAR * D[1..N] * RETURNS: * 1 IF SUCCESFUL, 0 IF THE MATRIX IS SINGULAR. * FLOPS: (N-1)(3/2JJ+5/2J) - [JJJ+JJ/2-3/2J] = N(3/2JJ+5/2J)-J(J+1)^2 * J==N-1 => N(N-1)(N-2)/2 * IF J=3M-1 AND N<=M(N-2), THEN N(27/2MMM -3/2MM -M)-54MMM+12MM+2M *************************************************************************/ int chol_sym_band(L, D, B, J, N) BAND L; VECT D; SYMBAND B; int J, N; { register int n, k, Jn, nJ; /* LOOP DOWN THROUGH ROWS OF BAND MATRIX STRUCTURE */ D[1] = B[1][J+1]; if (istiny(D[1])) return 0; if (N == 1) return 1; for (n=2; n <= J+1; n++) { Jn = J-n; L[n][Jn+2] = B[n][Jn+2] / D[1]; for (k=2; k < n; k++) L[n][k+Jn+1] = ddot123((double) B[n][k+Jn+1], &D[1], &L[n][Jn+2], &L[k][J-k+2], k-1) / D[ k ]; D[n] = ddot123((double) B[n][J+1], &D[1], &L[n][Jn+2], &L[n][Jn+2], n-1); if (istiny(D[n])) return 0; } for (; n <= N; n++) { nJ = n-J; L[n][1] = B[n][1] / D[ nJ ]; for (k=2; k <= J; k++) L[n][k] = ddot123((double) B[n][k], &D[nJ], &L[n][1], &L[nJ+k-1][J-k+2], k-1) / D[ nJ+k-1 ]; D[n] = ddot123((double) B[n][J+1], &D[nJ], &L[n][1], &L[n][1], J); if (istiny(D[n])) return 0; } return 1; } /* * SUBTRACT n POINT INNER PRODUCT OF p1, p2, p3 * FROM init * pi[0..n-1] */ static double ddot123(sum, p1, p2, p3, n) register double sum; register VECT p1, p2, p3; int n; { for (; n > 0; n--, p1++, p2++, p3++) sum -= *p1 * *p2 * *p3; return sum; } static int istiny(x) DATA x; { return (x == 0.); } /****************************** solve_sym_band() ************************* * Solve (L D L') x = z * SEE Golub and Van Loan Alg. 5.3-2,3 * INPUT: * L: BAND(J, N) * D[1..N] * z[1..N] * job: if 1, return in z, if 0 return in x * OUTPUT: * x[1..n] * FLOPS: N+2[J(J-1)+2J(N-J)] = N(1+4J) -2JJ -2J * J==N-1 => 2N^2 - N * IF J=3M-1, N<=M(N-2) THEN N(12MM-3M) -42MM +12M *************************************************************************/ void solve_sym_band(x, L, D, z, J, N, job) VECT x, D, z; BAND L; int J, N; { register int i, j; if (job) x = z; /* NOW x IS REALLY u (SOLVING Lu = z) */ else COPY_VECTOR(z, x, N); if (N == 1) { x[1] /= D[1]; return; } for (i=2; i <= J; i++) { /* x[i] = z[i]; */ dotdiff(&x[i], &L[i][J-i+2], &x[1], i-1); } for ( ; i <= N; i++) { /* x[i] = z[i]; */ dotdiff(&x[i], &L[i][1], &x[i-J], J); } /* NOW x IS REALLY v (SOLVING Dv = u) */ for (i=1; i <= N; i++) x[i] /= D[i]; /* NOW IS REALLY x ON LHS, v ON RHS (SOLVING L' x = v) */ /* x[N] = x[N]; */ for (i=N-1; i >= N-J; i--) for (j=i+1; j <= N; j++) x[i] -= L[j][i-j+J+1] * x[j]; for (; i >= 1; i--) for (j=i+1; j <= i+J; j++) x[i] -= L[j][i-j+J+1] * x[j]; } /* * SUBTRACT n POINT INNER PRODUCT OF p1, p2 FROM x * INPUT: * x * p1[0..n-1] * p2[0..n-1] * OUTPUT: * x x -= */ static void dotdiff(x, p1, p2, n) DATA *x, *p1, *p2; int n; { for (; n > 0; n--) *x -= *p1++ * *p2++; } /****************************** invert_sym_band() ************************ * FIND CENTRAL INVERSE OF A SYMMETRIC BAND MATRIX B = L * D * L' * See Hutchinson and deHoog, Numerische Mathematik, 1985. * IB = inv(B) * INPUT: * L: SYMBAND(J,N) * D[1..N] * OUTPUT: * IB: SYMBAND(J,N) central bands of inv(B) * FLOPS: 2JJN + ... * IF J=3M-1, N->M(N-2) THEN 18MMMN+... *************************************************************************/ void invert_sym_band(IB, L, D, J, N) SYMBAND IB, L; VECT D; int J, N; { register int n, l, k, nJ; IB[N][J+1] = 1. / D[N]; for (n=N-1 ; n >= 1; n--) { nJ = (J < N-n) ? J : N-n; for (l=1; l <= nJ; l++) { IB[n+l][-l+J+1] = 0.; /* WASTES A FLOP...*/ for (k=1; k <= l; k++) IB[n+l][-l+J+1] -= L[n+k][-k+J+1] * IB[n+l][k-l+J+1]; for ( ; k <= nJ; k++) IB[n+l][-l+J+1] -= L[n+k][-k+J+1] * IB[n+k][l-k+J+1]; } IB[n][J+1] = 1. / D[n]; for (l=1; l <= nJ; l++) IB[n][J+1] -= L[n+l][-l+J+1] * IB[n+l][-l+J+1]; } } End of smooth/sym_band.c echo smooth/vs_init.c 1>&2 cat >smooth/vs_init.c <<'End of smooth/vs_init.c' /****************************** vs_init.c ********************************* * SUBROUTINES FOR CUBIC SPLINE SMOOTHING ALGORITHM; * THESE ARE ALL INDEPENT OF THE SMOOTHING PARAMETERS. * SEE band.h FOR THE FORMAT OF THE BAND MATRICES. * SEE vs_smooth.c FOR DIMENSIONS OF * R, QtQ, QtSQ * * EXCEPTION: SINCE Qt IS UPPER BAND AND NOT SYMMETRIC, * THE FOLLOWING STORAGE FORMAT IS USED: * Qt(i,i+k) = Qt[i][3-k], for k=0,1,2 i = 1, ..., N-2 * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. **************************************************************************/ #include "smooth.h" void spline_mat(); void form_qt_Y(), form_qt_s_q(); static void form_qt_s_q_reg(), form_qt_s_q_iid(); /****************************** spline_mat() ***************************** * COMPUTE R, QtQ = Q'*Q, AND Qt = Q' GIVEN KNOTS ti * IF ti IS NULL, THEN THE SAMPLES ARE SPACED EQUALLY. * INPUT: * ti[1..N] sample points (increasing order) (or null) * N: # of points * OUTPUT: * R, QtQ, Qt * FLOPS: N-1 + 6(N-1) = 7N - 7 * PLUS (N-1) * 9 for QtQ - NOT NEEDED FOR NON-IID *************************************************************************/ void spline_mat(R, Qt, ti, N) SYMBAND R; BAND Qt; VECT ti; int N; { register int i, i1; VECT hi; TMP_ALLOC_VECT(hi, N-1) /* COMPUTE DIFFERENCES OF THE KNOTS */ if (ti) for (i=N-1; i >= 1; i--) { hi[i] = ti[i+1] - ti[i]; if (hi[i] <= 0.) ABORT("ti's not in increasing order"); } else for (i=N-1; i >= 1; i--) hi[i] = 1.; for (i=N-2, i1=N-1; i >= 1; i--, i1--) { R[i][1] = hi[i] / 6.; R[i][2] = (hi[i] + hi[i1]) / 3.; Qt[i][1] = 1. / hi[i1]; Qt[i][3] = 1. / hi[i]; Qt[i][2] = -Qt[i][3] - Qt[i][1]; } R[1][1] = 0.; /* CORRECTION */ #ifdef NOTUSED i = N-2; QtQ[i ][3] = Qt[i][3]*Qt[i ][3] + Qt[i][2]*Qt[i ][2]+ Qt[i][1]*Qt[i][1]; i1 = i--; QtQ[i1][2] = Qt[i][2]*Qt[i1][3] + Qt[i][1]*Qt[i1][2]; QtQ[i ][3] = Qt[i][3]*Qt[i ][3] + Qt[i][2]*Qt[i ][2]+ Qt[i][1]*Qt[i][1]; i2 = i1; i1 = i--; for ( ; i >= 1; i--, i1--, i2--) { QtQ[i2][1] = Qt[i][1]*Qt[i2][3]; QtQ[i1][2] = Qt[i][2]*Qt[i1][3] + Qt[i][1]*Qt[i1][2]; QtQ[i ][3] = Qt[i][3]*Qt[i ][3] + Qt[i][2]*Qt[i ][2]+ Qt[i][1]*Qt[i][1]; } #endif TMP_FREE_VECT(hi) } /****************************** form_qt_Y() ****************************** * COMPUTE (Q' @ I) Y * INPUT: * Qt: see above * Y[1..MN] measured data (Y1_1 .. YM_1, .., Y1_1 .. YM_N) * OUTPUT: * QtY[1..M(N-2)] (Qt @ I) * Y * FLOPS: M(5N-5) *************************************************************************/ void form_qt_Y(QtY, Qt, Y, M, N) VECT QtY; BAND Qt; VECT Y; int M, N; { register int m, i; for (m=M; m >= 1; m--) for (i=N-2; i >= 1; i--) QtY[(i-1)*M+m] = Y[(i-1)*M+m] * Qt[i][3] + Y[(i )*M+m] * Qt[i][2] + Y[(i+1)*M+m] * Qt[i][1]; } /****************************** form_qt_s_q() **************************** * COMPUTE (Qt @ I) Sigma (Q @ I) * INPUT: * Qt see above * covar[1..?] error covariance structure(s) * flag_iid: more than one covariance? * OUTPUT: * QtSQ *************************************************************************/ void form_qt_s_q(QtSQ, Qt, covar, M, N, flag_iid) SYMBAND QtSQ; BAND Qt; COVAR *covar; int M, N, flag_iid; { if (flag_iid) form_qt_s_q_iid(QtSQ, Qt, &covar[1], M, N); else form_qt_s_q_reg(QtSQ, Qt, covar, M, N); } /****************************** form_qt_s_q_reg() ************************ * COMPUTE QtSQ = (Q' @ I) diag(sigma) (Q @ I) * INPUT: * Qt: see above * covar[1..N] error covariance structures * OUTPUT: * QtSQ * FLOPS: (N-4)(MM+1)+(N-3)(3MM+2)+(N-2)(5MM+3) = N(9MM+6)-23MM-16 * OLD (N-2) M M 2 + (N-1) M M 5 + (N-1) M M 8 = 15MMN-17MM *************************************************************************/ static void form_qt_s_q_reg(QtSQ, Qt, cov, M, N) SYMBAND QtSQ; BAND Qt; COVAR *cov; int M, N; { register int m, k, n; for (n=N-2; n >= 3; n--) { register DATA t1 = Qt[n][3] * Qt[n-2][1]; for (m=M; m >= 1; m--) { for (k=M; k >= m; k--) QtSQ[(n-1)*M+m][k-m+M] = t1 * cov[n ].cov[k][m-k+M]; for (; k >= 1; k--) QtSQ[(n-1)*M+m][k-m+M] = t1 * cov[n ].cov[m][k-m+M]; } } for (n=N-2; n >= 2; n--) { register DATA t1 = Qt[n][3] * Qt[n-1][2]; register DATA t2 = Qt[n][2] * Qt[n-1][1]; for (m=M; m >= 1; m--) { for (k=M; k >= m; k--) QtSQ[(n-1)*M+m][k-m+2*M] = t1 * cov[n ].cov[k][m-k+M] + t2 * cov[n+1].cov[k][m-k+M]; for (; k >= 1; k--) QtSQ[(n-1)*M+m][k-m+2*M] = t1 * cov[n ].cov[m][k-m+M] + t2 * cov[n+1].cov[m][k-m+M]; } } for (n=N-2; n >= 1; n--) { register DATA t1 = Qt[n][3] * Qt[n][3]; register DATA t2 = Qt[n][2] * Qt[n][2]; register DATA t3 = Qt[n][1] * Qt[n][1]; for (m=M; m >= 1; m--) { for (k=m; k >= 1; k--) QtSQ[(n-1)*M+m][k-m+3*M] = t1 * cov[n ].cov[m][k-m+M] + t2 * cov[n+1].cov[m][k-m+M] + t3 * cov[n+2].cov[m][k-m+M]; } } } /****************************** form_qt_s_q_iid() ************************ * COMPUTE QtSQ = (Q' @ I) diag(sigma) (Q @ I) * INPUT: * Qt: see above * covar one covariance matrix structure * OUTPUT: * QtSQ * FLOPS: *************************************************************************/ static void form_qt_s_q_iid(QtSQ, Qt, covar, M, N) SYMBAND QtSQ; BAND Qt; COVAR *covar; int M, N; { register int m, k, n; SYMBAND cov = covar->cov; for (n=N-2; n >= 3; n--) { register DATA t1 = Qt[n][3] * Qt[n-2][1]; for (m=M; m >= 1; m--) { for (k=M; k >= m; k--) QtSQ[(n-1)*M+m][k-m+M] = t1 * cov[k][m-k+M]; for (; k >= 1; k--) QtSQ[(n-1)*M+m][k-m+M] = t1 * cov[m][k-m+M]; } } for (n=N-2; n >= 2; n--) { register DATA t1 = Qt[n][3] * Qt[n-1][2] + Qt[n][2] * Qt[n-1][1]; for (m=M; m >= 1; m--) { for (k=M; k >= m; k--) QtSQ[(n-1)*M+m][k-m+2*M] = t1 * cov[k][m-k+M]; for (; k >= 1; k--) QtSQ[(n-1)*M+m][k-m+2*M] = t1 * cov[m][k-m+M]; } } for (n=N-2; n >= 1; n--) { register DATA t1 = Qt[n][3] * Qt[n][3] + Qt[n][2] * Qt[n][2] + Qt[n][1] * Qt[n][1]; for (m=M; m >= 1; m--) { for (k=m; k >= 1; k--) QtSQ[(n-1)*M+m][k-m+3*M] = t1 * cov[m][k-m+M]; } } } End of smooth/vs_init.c echo smooth/vs_proc.c 1>&2 cat >smooth/vs_proc.c <<'End of smooth/vs_proc.c' /****************************** vs_proc.c ********************************* * SUBROUTINES FOR CUBIC SPLINE SMOOTHING ALGORITHM * THESE ALL DEPEND ON alpha - THE SMOOTHING PARAMETERS. * SEE vs_smooth.c and vs_init.c FOR THE DIMENSIONS OF Qt etc. * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. **************************************************************************/ #include "smooth.h" void form_band(), form_q_c(), form_sig_q_c(); extern void mult_full_vect(); /****************************** form_band() ****************************** * FORM THE BAND MATRIX USED IN THE REINSCH ALGORITHM * B = (R @ inv(alpha)) + (Qt @ I) Sigma (Q @ I) * INPUT: * QtSQ: SYMBAND(3*M-1, M*(N-2)) * R_A: SYMBAND(2*M-1, M*(N-2)) * OUTPUT: * sum: SYMBAND(3*M-1, M*(N-2)) R_A + QtQS * FLOPS: 2MM(N-2) *************************************************************************/ void form_band(sum, R_A, QtSQ, M, N) SYMBAND sum, R_A, QtSQ; int M, N; { register int n, m; for (m=3*M; m > M; m--) for (n = M*(N-2); n >= 1; n--) sum[n][m] = R_A[n][m-M] + QtSQ[n][m]; /* TRICKY: THIS ALSO FILLS IN DESIRABLE 0'S IN UPPER LEFT OF sum */ for (n = M*(N-2); n >= 1; n--) COPY_VECTOR(QtSQ[n], sum[n], M); } /****************************** form_q_c() ******************************* * COMPUTE (Q @ I) * c_hat * INPUT: * Qt: see above * c[1..M(N-2)] estimated coefficients * OUTPUT: * q_c[1..MN] (Q @ I) * c_hat * FLOPS: 6M(N-2) *************************************************************************/ void form_q_c(q_c, Qt, c, M, N) VECT q_c, c; BAND Qt; int M, N; { register int n; ZERO_VECTOR(q_c, M*N); INC_SCALED_VECTOR(q_c, (double) Qt[1][3], c, M); q_c += M; INC_SCALED_VECTOR(q_c, (double) Qt[1][2], c, M); INC_SCALED_VECTOR(q_c, (double) Qt[2][3], c + M, M); q_c += M; for (n=2; n <= N-3; n++, q_c += M) { INC_SCALED_VECTOR(q_c, (double) Qt[n-1][1], c + M*(n-2), M); INC_SCALED_VECTOR(q_c, (double) Qt[n ][2], c + M*(n-1), M); INC_SCALED_VECTOR(q_c, (double) Qt[n+1][3], c + M*n, M); } INC_SCALED_VECTOR(q_c, (double) Qt[N-3][1], c + M*(N-4), M); INC_SCALED_VECTOR(q_c, (double) Qt[N-2][2], c + M*(N-3), M); q_c += M; INC_SCALED_VECTOR(q_c, (double) Qt[N-2][1], c + M*(N-3), M); } /****************************** form_sig_q_c() *************************** * COMPUTE diag(sigma) * (Q' @ I) * c_hat * INPUT: * covar[1..?] error covariance structure(s) * q_c[1..MN] (Q' @ I) * c_hat * flag_iid is the noise iid? * OUTPUT: * sig_q_c[1..N][1..M] * FLOPS: NM(M+M-1) = N(2MM-M) *************************************************************************/ void form_sig_q_c(sig_q_c, covar, q_c, M, N, flag_iid) OMATRIX sig_q_c; CMAT q_c; COVAR *covar; int M, N; { register int n; for (n=1; n <= N; n++, q_c += M) { OMATRIX cov = covar[flag_iid ? 1 : n].cov; mult_full_vect(sig_q_c[n], cov, q_c, M); } } End of smooth/vs_proc.c echo smooth/vs_smooth.c 1>&2 cat >smooth/vs_smooth.c <<'End of smooth/vs_smooth.c' /************************* vs_smooth.c *********************************** * MAIN SUBROUTINES FOR * FIXED INTERVAL SMOOTHING OF VECTOR MEASUREMENTS WITH CUBIC SPLINES. * * REFERENCE: Non-parametric Fixed-Interval Smoothing with Vector Splines * Jeff Fessler, IEEE Trans-ASSP 1989 * PROBLEM: * minimize (Y - a)' inv(Sigma) (Y - a) * + \sum_{m=1}^M alpha_m c' R c * subject to (Q' @ I) a = (R @ I) c * (Sigma = blockdiag(sigma_1, ..., sigma_N)) * SOLUTION: * c_alpha = inv[(R@inv(alpha)) + (Q'@I) Sigma (Q@I)] (Q'@I) Y * ahat = Y - Sigma (Q@I) c_alpha * ("@" IS USED TO DENOTE MATRIX TENSOR (KRONECKER) PRODUCT.) * TO USE: * first call set_matrices() and set_covariances(), * then call vs_smooth() and compute_stats() * for each value of alpha. * THERE ARE ESSENTIALLY TWO VERSIONS CONTAINED WITHIN, DEPENDING * ON WHETHER THE NOISE IS IID OR NOT. IF SO, ONLY ONE COVARIANCE * MATRIX NEED BE PASSED, IE: * Flag_iid: * 1: noise is iid: Sigma[1..MxM] covar[1] * 0: noise varies: Sigma[1..MMxN] covar[1..N] * USES: exp() * Jeff Fessler * COPYRIGHT 1989, THE BOARD OF TRUSTEES OF THE * LELAND STANFORD JUNIOR UNIVERSITY. ALL RIGHTS RESERVED. *************************************************************************/ #include #include "smooth.h" int vs_setup_alloc(), vs_setup_ti(), vs_setup_alpha(); int vs_setup_yi(), vs_setup_cov(); int vs_smooth(), vs_transfer(); double roughness(), post_penalty(); int penalty(); void vs_cleanup(); extern void spline_mat(); extern void form_qt_Y(), form_qt_s_q(), kron_diag(), form_band(); extern int chol_sym_band(); extern void solve_sym_band(); extern void form_q_c(), form_sig_q_c(), mult_band_vect(); int NI=0, MI=0; /* 'I' SUFFIX: SIZES WHEN INITIALIZED */ int NG, MG; /* 'G' SUFFIX: MEANS "GLOBAL" */ int Flag_iid; int Flag_equispace; /* EQUALLY SPACED SAMPLE POINTS */ /* GLOBAL VARIABLES */ BAND Qt; /* SEE vs_init.c FOR FORMAT */ static SYMBAND R; /* R(1,N-2) */ static SYMBAND QtSQ; /* QtSQ(3M-1,M(N-2)) */ static VECT alpha_inv; /* [1..MG] diag(1/alpha) */ static SYMBAND R_A; /* R_A(2M-1,M(N-2)) */ BAND BL; /* BL(3M-1,M(N-2)) */ VECT BD; /* BD[1..M(N-2)] */ static VECT c_alpha; /* c_alpha[1..M(N-2)] */ OMATRIX sig_q_c; /* sig_q_c[1..N][1..M] */ static VECT QtY; /*, QtY[1..M(N-2)] */ /* GLOBAL PROINTERS */ CMAT Gyi; COVAR *Gcovar; /* * ALLOCATE MATRICES, USING UPPER BOUNDS FOR THE SIZES TO * BE USED IN SUBSEQUENT CALLS. * INPUT: * NI: # of sample points * MI: length of sample vectors * ALLOCATED: * R, Qt, QtSQ, QtY * alpha_inv, R_A, BL, BD, c_alpha, sig_q_c * RETURN: * 1 if OK, 0 if too big to initialize. */ int vs_setup_alloc(M, N) int M, N; { NI = N; MI = M; /* FIX: IF AN ALLOCATION IS UNSUCCESSFUL, FREE THE OTHERS */ if (!(Qt = ALLOC_BAND(NI-2, 2))) return 0; if (!(QtSQ = ALLOC_BAND(MI*(NI-2), 3*MI-1))) return 0; if (!(R = ALLOC_BAND(NI-2, 1))) return 0; if (!(R_A = ALLOC_BAND(MI*(NI-2), 2*MI-1))) return 0; if (!(BL = ALLOC_BAND(MI*(NI-2), 3*MI-1))) return 0; if (!(BD = ALLOC_VECTOR(MI*(NI-2)))) return 0; if (!(QtY = ALLOC_VECTOR(MI*(NI-2)))) return 0; if (!(c_alpha = ALLOC_VECTOR(MI*(NI-2)))) return 0; if (!(alpha_inv = ALLOC_VECTOR(MI))) return 0; if (!(sig_q_c = ALLOC_MATRIX(NI, MI))) return 0; return 1; } /* * SETUP MATRICES FOR SUBSEQUENT CALLS TO vs_smooth() * THIS NEED ONLY BE DONE ONCE FOR SUBSEQUENT ITERATION ON alpha * * INPUT: * ti[1..N] sample points (or null if equispaced) * NG: # of sample points * RETURN: * 1 if OK, 0 if error */ int vs_setup_ti(ti, N) VECT ti; int N; { NG = N; if (NG > NI) return 0; /* NO BIGGER THAN NI */ if (NG < 5) return 0; /* MUST HAVE AT LEAST 3 POINTS */ Flag_equispace = (ti == (VECT) 0); spline_mat(R, Qt, ti, NG); return 1; } /* * SETUP MATRICES FOR SUBSEQUENT CALLS TO vs_smooth() * THIS NEED ONLY BE DONE ONCE FOR SUBSEQUENT ITERATION ON alpha * * INPUT: * yi[1..MN] measurements * N: # of sample points * M: length of sample vectors * RETURN: * 1 if OK, 0 if error */ int vs_setup_yi(yi, M, N) VECT yi; int M, N; { MG = M; if (MG != MI) return 0; if (NG != N) return 0; Gyi = yi; form_qt_Y(QtY, Qt, yi, MG, NG); return 1; } /* * SETUP MATRICES FOR SUBSEQUENT CALLS TO vs_smooth() * THIS NEED ONLY BE DONE ONCE FOR SUBSEQUENT ITERATION ON alpha * * INPUT: * covar[1..?] noise covariances (see above) * flag_iid: see above * USES: * Qt * EFFECT: * QtSQ, Gcovar, Flag_iid */ int vs_setup_cov(covar, flag_iid) COVAR *covar; int flag_iid; { Flag_iid = flag_iid; Gcovar = covar; ZERO_BAND(QtSQ, 3*MI-1, MI*(NI-2)); form_qt_s_q(QtSQ, Qt, covar, MG, NG, Flag_iid); return 1; } /* * SETUP MATRICES FOR SUBSEQUENT CALLS TO vs_smooth() * THIS NEED ONLY BE DONE ONCE FOR SUBSEQUENT ITERATION ON alpha * * INPUT: * alpha[1..M] smoothing parameters */ int vs_setup_alpha(alpha) ALPHA_SET alpha; { int m; /* COMPUTE DIAGONAL OR MATRIX WITH 1/Alpha ON DIAGONAL */ for (m=MG; m >= 1; m--) alpha_inv[m] = (1./NG) * exp(-alpha[m]); ZERO_BAND(R_A, 2*MI-1, MI*(NI-2)); kron_diag(R_A, R, 1, NG-2, alpha_inv, MG); return 1; } /* * FREE ALL THE ALLOCATED MATRICES */ void vs_cleanup() { FREE_BAND(Qt); FREE_BAND(QtSQ); FREE_BAND(R); FREE_BAND(R_A); FREE_BAND(BL); FREE_VECTOR(BD); FREE_VECTOR(QtY); FREE_VECTOR(c_alpha); FREE_VECTOR(alpha_inv); FREE_MATRIX(sig_q_c); MI = NI = 0; } /* * ALGORITHM FOR SMOOTHING VECTOR MEASUREMENTS * job OPTIONS: * 0: ignore result, just compute sig_q_c (for cross_val()) * 1: result is yi * 2: result is (yi - yhat) on output * 4: result is yi on input -> yhat on output * INPUTS: * result[1..MxN] see above * job * OUTPUT: * result[1..MxN] see above * USES: * everything * EFFECT: * sig_q_c: set to (yi - yhat) */ int vs_smooth(result, job) CMAT result; int job; { VECT q_c; TMP_ALLOC_VECT(q_c,MG*NG) if (! vs_transfer() ) return 0; solve_sym_band(c_alpha, BL, BD, QtY, 3*MG-1, MG*(NG-2), 0); form_q_c(q_c, Qt, c_alpha, MG, NG); form_sig_q_c(sig_q_c, Gcovar, q_c, MG, NG, Flag_iid); if (job == 1) SUB_VECTOR(result, Gyi, sig_q_c[1], MG*NG); if (job == 2) COPY_VECTOR(sig_q_c[1], result, MG*NG); if (job == 4) DEC_VECTOR(result, sig_q_c[1], MG*NG); TMP_FREE_VECT(q_c); return 1; } /* * SETUP MATRIX FOR COMPUTING c_alpha * USES: * R_A, QtSQ * EFFECT: * BL, BD */ int vs_transfer() { SYMBAND sum; TMP_ALLOC_BAND(sum,MG*(NG-2),3*MG-1) form_band(sum, R_A, QtSQ, MG, NG); ZERO_BAND(BL, 3*MI-1, MI*(NI-2)); if (!chol_sym_band(BL, BD, sum, 3*MG-1, MG*(NG-2))) return 0; TMP_FREE_BAND(sum); return 1; } /* * COMPUTE ROUGHNESS PENALTY * vs_smooth() MUST BE CALLED FIRST. * USES: * c_alpha AND alpha_inv * RETURNS: * rough roughness term c' (R @ a) c */ double roughness() { int m, n; double rough = 0.; if (!Flag_equispace) return post_penalty(); for (m=MG; m >= 1; m--) { double ssq = 2. * c_alpha[(NG-3)*MG+m] * c_alpha[(NG-3)*MG+m]; for (n=NG-3; n >= 1; n--) { ssq += c_alpha[(n-1)*MG+m] * (2. * c_alpha[(n-1)*MG+m] + c_alpha[(n )*MG+m]); } rough += ssq * alpha_inv[m]; } return rough * 1./3.; } /* * COMPUTE ROUGHNESS PENALTY FOR SPLINE COEFFICIENTS * IMMEDIATELY AFTER SMOOTHING * RETURNS: \sum_m alpha_m c' R c = c_alpha' R_A c_alpha * USES: c_alpha, R_A, MG, NG */ double post_penalty() { double pen; VECT r_c; TMP_ALLOC_VECT(r_c,MG*(NG-2)) mult_band_vect(r_c, R_A, c_alpha, 2*MG-1, MG*(NG-2)); pen = DOT_VECTOR(r_c, c_alpha, MG*(NG-2)); TMP_FREE_VECT(r_c); return pen; } /* * COMPUTE ROUGHNESS PENALTY FOR SPLINE COEFFICIENTS * vs_setup_alloc() and vs_setup_ti() MUST BE CALLED FIRST. * INPUT: * alpha[1..M] smoothing parameters * yi[1..MxN] points to be interpolated * OUTPUT: * pen penalty term: c' (R @ a) c * = \sum_{m=1}^M alpha_m c_m' R c_m * RETURNS: 1 IF OK 0 if error * USES: QtY, c_alpha, alpha_inv, R_A, BL, BD, MG, NG */ int penalty(pen, alpha, yi, M, N) double *pen; ALPHA_SET alpha; VECT yi; int M, N; { VECT RAD, res; BAND RAL; if (!vs_setup_yi(yi, M, N)) return 0; if (!vs_setup_alpha(alpha)) return 0; TMP_ALLOC_VECT(RAD,MG*(NG-2)) TMP_ALLOC_VECT(res,MG*(NG-2)) TMP_ALLOC_BAND(RAL,MG*(NG-2),2*MG-1) if (!chol_sym_band(RAL, RAD, R_A, 2*MG-1, MG*(NG-2))) return 0; solve_sym_band(res, RAL, RAD, QtY, 2*MG-1, MG*(NG-2), 0); *pen = DOT_VECTOR(res, QtY, MG*(NG-2)); TMP_FREE_VECT(RAD); TMP_FREE_VECT(res); TMP_FREE_BAND(RAL); return 1; } End of smooth/vs_smooth.c echo Test/Makefile 1>&2 cat >Test/Makefile <<'End of Test/Makefile' # makefile for testing band matrix routines include ../List DEFS= -D$(PREC)Prec INCS= -I$(DEFDIR) -I$(MATRIX) -I$(SMOOTH) LIBS=-lm CFLAGS=-g $(FLOATFLAGS) $(DEFS) $(INCS) all: $(BIN)/test_ci $(BIN)/test_vs $(BIN)/test_ss include ../Vlist DOTH= $(SMOOTHH) data.h test_ci.o test_vs.o test_ss.o read_data.o: $(DOTH) TCC= \ test_ci.c \ $(SMOOTHC) $(IODEPC) $(MATRIXC) TCO= \ test_ci.o \ $(SMOOTHO) $(IODEPO) $(MATRIXO) $(BIN)/test_ci: $(TCO) $(CC) $(CFLAGS) -o $@ $(TCO) $(LIBS) ALLC= $(SMOOTHEXC) $(MATRIXC) $(IODEPC) $(MYNRC) #ALLO= $(SMOOTHEXO) $(MATRIXO) $(IODEPO) $(MYNRO) ALLO= $(DEFDIR)/libvspline.a $(NUMRECDIR)/libnumrecC.a TVC= test_vs.c read_data.c $(ALLC) TVO= test_vs.o read_data.o $(ALLO) $(BIN)/test_vs: $(TVO) $(CC) $(CFLAGS) -o $@ $(TVO) $(LIBS) TSC= test_ss.c read_data.c $(ALLC) TSO= test_ss.o read_data.o $(ALLO) $(BIN)/test_ss: $(TSO) $(CC) $(CFLAGS) -o $@ $(TSO) $(LIBS) LL: $(DOTH) $(TCC) $(TSC) $(TVC) rm -f LL; ( \ lint -u $(DEFS) $(INCS) $(TVC) $(LIBS); echo; \ lint -u $(DEFS) $(INCS) $(TCC) $(LIBS); echo; \ lint -u $(DEFS) $(INCS) $(TSC) $(LIBS); echo; \ ) > LL include /user/fessler/lib/Prmake End of Test/Makefile echo matrix/Makefile 1>&2 cat >matrix/Makefile <<'End of matrix/Makefile' # makefile for band matrices and linear non-parametric smoothing include ../List DEFS=-D$(PREC)Prec INCS=-I$(DEFDIR) LIBS=-lm CFLAGS=-O $(FLOATFLAGS) $(DEFS) $(INCS) LFLAGSREG=$(DEFS) $(INCS) #MEXDIR=/usr/local/lib/newmatlab/mex #LFLAGSMEX=$(DEFS) $(INCS) -DMATLAB -I$(MEXDIR) all: \ ofiles include List ofiles: \ $(IODEPO) $(IOMEXO) \ $(MATRIXO) LL: $(MATRIXH) $(MATRIXC) $(IODEPC) rm -f LL; ( \ lint -u $(LFLAGSREG) $(MATRIXC) $(LIBS); echo; \ ) > LL End of matrix/Makefile echo numrecC/Makefile 1>&2 cat >numrecC/Makefile <<'End of numrecC/Makefile' # double-precision versions of "Numerical Recipes in C" subroutines include ../List #DEFS= #INCS= LIBS=-lm #CFLAGS=-O $(FLOATFLAGS) $(DEFS) $(INCS) all: \ ofiles \ cfiles .SUFFIXES: .Z .Z.c: uncompress $@ include List cfiles: $(MYNRC) $(MYNRC): for i in $(CLIST) ; do uncompress < $(NRCDIR)/$$i.Z > $$i ; done ofiles: $(MYNRO) LL: $(MYNRH) $(MYNRC) rm -f LL; ( \ lint -v $(DEFS) $(INCS) $(MYNRC) $(LIBS); echo; \ ) > LL End of numrecC/Makefile echo smooth/Makefile 1>&2 cat >smooth/Makefile <<'End of smooth/Makefile' # makefile for band matrices and # for linear non-parametric vector-spline smoothing include ../List DEFS=-D$(PREC)Prec INCS=-I$(DEFDIR) -I$(MATRIX) LIBS=-lm CFLAGS=-O $(FLOATFLAGS) $(DEFS) $(INCS) all: \ ofiles include List ofiles: \ $(SMOOTHEXO) LL: $(SMOOTHH) $(SMOOTHEXC) rm -f LL; ( \ lint -u $(DEFS) $(INCS) $(SMOOTHEXC) $(LIBS); echo; \ ) > LL End of smooth/Makefile echo matrix/List 1>&2 cat >matrix/List <<'End of matrix/List' # THESE ARE NEEDED FOR MAKEFILES THAT USE THE SMOOTH LIBRARIES IODEPC= \ $(MATRIX)/debug.c \ $(MATRIX)/mexdep.c IODEPO= \ $(MATRIX)/debug.o \ $(MATRIX)/mexdep.o IOMEXO= \ $(MEXDEP)/debug.o \ $(MEXDEP)/mexdep.o MATRIXC= \ $(MATRIX)/alloc.c \ $(MATRIX)/ident.c \ $(MATRIX)/io_d.c \ $(MATRIX)/io_f.c \ $(MATRIX)/vector_f.c \ $(MATRIX)/vector_d.c # $(MATRIX)/sub_mat.c MATRIXO= \ $(MATRIX)/alloc.o \ $(MATRIX)/ident.o \ $(MATRIX)/io_d.o \ $(MATRIX)/io_f.o \ $(MATRIX)/vector_f.o \ $(MATRIX)/vector_d.o # $(MATRIX)/sub_mat.o MATRIXH= \ $(DEFDIR)/defs.h \ $(MATRIX)/matrix.h \ $(MATRIX)/mexdep.h $(IODEPO) $(IOMEXO) $(MATRIXO): $(MATRIXH) End of matrix/List echo numrecC/List 1>&2 cat >numrecC/List <<'End of numrecC/List' # THESE ARE NEEDED FOR MAKEFILES THAT USE THE numrecC LIBRARIES NRCDIR=/spin/usr/local/src/numrecC CLIST= \ brent.c \ f1dim.c \ linmin.c \ mnbrak.c \ powell.c \ nrutil.c # uncomment this for debugging only #MYNRC= \ # $(MYNR)/brent.c \ # $(MYNR)/f1dim.c \ # $(MYNR)/linmin.c \ # $(MYNR)/mnbrak.c \ # $(MYNR)/powell.c \ # $(MYNR)/nrutil.c MYNRO= \ $(NUMRECDIR)/libnumrecC.a End of numrecC/List echo smooth/List 1>&2 cat >smooth/List <<'End of smooth/List' # THESE ARE NEEDED FOR MAKEFILES THAT USE THE SMOOTH LIBRARIES BANDC= \ $(SMOOTH)/band_util.c \ $(SMOOTH)/covar.c \ $(SMOOTH)/sym_band.c BANDO= \ $(SMOOTH)/band_util.o \ $(SMOOTH)/covar.o \ $(SMOOTH)/sym_band.o SMOOTHC= \ $(BANDC) \ $(SMOOTH)/interp.c \ $(SMOOTH)/vs_init.c \ $(SMOOTH)/vs_smooth.c \ $(SMOOTH)/vs_proc.c SMOOTHO= \ $(BANDO) \ $(SMOOTH)/interp.o \ $(SMOOTH)/vs_init.o \ $(SMOOTH)/vs_smooth.o \ $(SMOOTH)/vs_proc.o SCOREC= \ $(SMOOTHC) \ $(SMOOTH)/cross_val.c \ $(SMOOTH)/cv_score.c SCOREO= \ $(SMOOTHO) \ $(SMOOTH)/cross_val.o \ $(SMOOTH)/cv_score.o BESTC= \ $(SCOREC) \ $(SMOOTH)/best_indiv.c \ $(SMOOTH)/best_mse.c \ $(SMOOTH)/best_score.c BESTO= \ $(SCOREO) \ $(SMOOTH)/best_indiv.o \ $(SMOOTH)/best_mse.o \ $(SMOOTH)/best_score.o SMOOTHEXC= \ $(BESTC) \ $(SMOOTH)/example1.c \ $(SMOOTH)/example2.c SMOOTHEXO= \ $(BESTO) \ $(SMOOTH)/example1.o \ $(SMOOTH)/example2.o SMOOTHH= \ $(DEFDIR)/defs.h \ $(MATRIX)/matrix.h \ $(SMOOTH)/smooth.h $(SMOOTHEXO): $(SMOOTHH) End of smooth/List echo Data/README 1>&2 cat >Data/README <<'End of Data/README' #This directory contains sample input data files and the resulting #output. Here is how they were generated: #test_vs 0 0 < data1 > data1_vs.out #test_vs 0 0 2 < data1 > data1_bs.out #test_vs 0 0 < data2 > data2_vs.out #test_vs 0 0 < data3i > data3i_vs.out #test_vs 0 1 < data3 > data3_vs.out #test_ss 2 < data4 > data4_ss.out #test_ci 7 241 < data5 > data5_ci.out End of Data/README echo Data/data1 1>&2 cat >Data/data1 <<'End of Data/data1' 7 2 -1 1 1 3 4 8 16 17 19 1 1 2 2 3 3 3 4 2 5 1 6 2 7 0.45454545454545 -0.18181818181818 -0.18181818181818 0.27272727272727 0.45454545454545 -0.18181818181818 -0.18181818181818 0.27272727272727 0.25000000000000 0 0 0.25000000000000 0.45454545454545 -0.18181818181818 -0.18181818181818 0.27272727272727 0.45454545454545 -0.18181818181818 -0.18181818181818 0.27272727272727 0.45454545454545 -0.18181818181818 -0.18181818181818 0.27272727272727 0.25000000000000 0 0 0.25000000000000 End of Data/data1 echo Data/data1_bs.out 1>&2 cat >Data/data1_bs.out <<'End of Data/data1_bs.out' -1 1 1 3 4 8 16 17 19 1 1 2 2 3 3 3 4 2 5 1 6 2 7 0 0.454545 -0.181818 0.272727 cv 0.415088 -1 1 cv 0.415088 -1 1 cv 0.389668 0 1 cv 0.431075 1.61803 1 cv 0.389668 0 1 cv 0.397547 0.618034 1 cv 0.393032 -0.381966 1 cv 0.389676 0.0132747 1 cv 0.389668 -0.00507164 1 cv 0.389668 -0.00584185 1 cv 0.389668 -0.00604069 1 cv 0.390144 -0.149631 1 cv 0.389737 -0.0608875 1 cv 0.389678 -0.0269903 1 cv 0.389669 -0.0140427 1 cv 0.389668 -0.00909722 1 cv 0.389668 -0.00720817 1 cv 0.389668 -0.00648665 1 cv 0.389668 -0.00623947 1 cv 0.389668 -0.00623947 1 cv 0.382795 -0.00623947 2 cv 0.376759 -0.00623947 3.61803 cv 0.375767 -0.00623947 4.36281 cv 0.37511 -0.00623947 5.56789 cv 0.375089 -0.00623947 5.64275 cv 0.375059 -0.00623947 5.76388 cv 0.374918 -0.00623947 6.66424 cv 0.374844 -0.00623947 8.12105 cv 0.374848 -0.00623947 7.96381 cv 0.374823 -0.00623947 10.4782 cv 0.374826 -0.00623947 9.6916 cv 0.374821 -0.00623947 14.2922 cv 0.374821 -0.00623947 12.5923 cv 0.374821 -0.00623947 20.4634 cv 0.374821 -0.00623947 17.4676 cv 0.374821 -0.00623947 30.4486 cv 0.374821 -0.00623947 20.4634 cv 0.374821 -0.00623947 24.2774 cv 0.374821 -0.00623947 26.6346 cv 0.374821 -0.00623947 28.0914 cv 0.374821 -0.00623947 28.9917 cv 0.374821 -0.00623947 29.5482 cv 0.374821 -0.00623947 29.8921 cv 0.374821 -0.00623947 30.1047 cv 0.374821 -0.00623947 30.236 cv 0.374821 -0.00623947 30.3172 cv 0.374821 -0.00623947 30.3674 cv 0.374821 -0.00623947 30.3984 cv 0.374821 -0.00623947 30.4175 cv 0.374821 -0.00623947 30.4294 cv 0.374821 -0.00623947 30.4367 cv 0.374821 -0.00623947 30.4426 cv 0.399042 0.987521 59.8852 cv 0.374821 -0.00623947 30.4426 cv 0.399292 0.993761 30.4426 cv 0.443869 -1.62427 30.4426 cv 0.374821 -0.00623947 30.4426 cv 0.382321 -0.624273 30.4426 cv 0.379714 0.375727 30.4426 cv 0.374664 -0.0720275 30.4426 cv 0.374659 -0.0903851 30.4426 cv 0.374658 -0.0872635 30.4426 cv 0.374658 -0.0872798 30.4426 cv 0.374658 -0.087296 30.4426 cv 0.374658 -0.0884759 30.4426 cv 0.374658 -0.0877467 30.4426 cv 0.374658 -0.0880252 30.4426 cv 0.374658 -0.0881974 30.4426 cv 0.374658 -0.0883038 30.4426 cv 0.374658 -0.0881316 30.4426 cv 0.374658 -0.088091 30.4426 cv 0.374658 -0.0880658 30.4426 cv 0.374658 -0.0880495 30.4426 cv 0.374658 -0.0880495 30.4426 cv 0.374658 -0.0880495 31.4426 cv 0.374658 -0.0880495 33.0606 cv 0.374658 -0.0880495 31.4426 cv 0.374658 -0.0880495 32.0606 cv 0.374658 -0.0880495 32.4426 cv 0.374658 -0.0880495 32.6787 cv 0.374658 -0.0880495 32.8246 cv 0.374658 -0.0880495 32.9147 cv 0.374658 -0.0880495 32.9705 cv 0.374658 -0.0880495 33.0049 cv 0.374658 -0.0880495 33.0262 cv 0.374658 -0.0880495 33.0393 cv 0.374658 -0.0880495 33.0475 cv 0.374658 -0.0880495 33.0525 cv 0.374658 -0.0880495 33.0556 cv 0.374658 -0.0880495 33.0575 cv 0.374658 -0.0880495 33.0587 cv 0.374658 -0.0880495 33.0594 cv 0.374658 -0.0880495 33.06 best score 0.374658 -0.0880495 33.06 1.5302 1.65711 2.12019 2.21907 2.3674 2.50006 2.81491 3.624 1.95067 5.87187 1.79186 6.15286 1.54557 6.71483 End of Data/data1_bs.out echo Data/data1_vs.out 1>&2 cat >Data/data1_vs.out <<'End of Data/data1_vs.out' -1 1 1 3 4 8 16 17 19 1 1 2 2 3 3 3 4 2 5 1 6 2 7 0 0.454545 -0.181818 0.272727 1.37196 1.53054 2.14842 2.21082 2.46427 2.54056 2.92836 3.72985 1.90505 5.81824 1.73788 6.11322 1.57102 6.72318 ur score -0.940108 cv score 0.415088 gcv score 0.0898478 roughness 0.136837 penalty yi (in) 33.6717 penalty ys (out) 0.136837 End of Data/data1_vs.out echo Data/data2 1>&2 cat >Data/data2 <<'End of Data/data2' 7 2 -1 1 1 3 4 8 16 17 19 10.1900 9.9168 19.6801 19.6389 30.2165 29.9494 28.9957 39.0335 30.4807 50.1602 19.5364 59.5550 9.9680 69.8931 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 100.0000 0 0 100.0000 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 100.0000 0 0 100.0000 End of Data/data2 echo Data/data2_vs.out 1>&2 cat >Data/data2_vs.out <<'End of Data/data2_vs.out' -1 1 1 3 4 8 16 17 19 10.19 9.9168 19.6801 19.6389 30.2165 29.9494 28.9957 39.0335 30.4807 50.1602 19.5364 59.555 9.968 69.8931 0 3 -1.5 2.25 12.2767 13.7389 24.8248 24.9624 30.1052 29.8181 31.2944 40.6251 29.5694 55.2449 22.8885 59.6666 9.95918 69.8455 ur score 17.0176 cv score 3636.94 gcv score 61.9313 roughness 206.21 penalty yi (in) 2954.03 penalty ys (out) 206.21 End of Data/data2_vs.out echo Data/data3 1>&2 cat >Data/data3 <<'End of Data/data3' 7 2 -1 1 10.1900 9.9168 19.6801 19.6389 30.2165 29.9494 28.9957 39.0335 30.4807 50.1602 19.5364 59.5550 9.9680 69.8931 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 100.0000 0 0 100.0000 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 100.0000 0 0 100.0000 End of Data/data3 echo Data/data3_vs.out 1>&2 cat >Data/data3_vs.out <<'End of Data/data3_vs.out' -1 1 10.19 9.9168 19.6801 19.6389 30.2165 29.9494 28.9957 39.0335 30.4807 50.1602 19.5364 59.555 9.968 69.8931 0 3 -1.5 2.25 12.5691 11.023 22.6026 20.4729 30.0222 30.0015 31.5261 39.7298 27.9329 49.6374 20.0683 59.7132 10.0219 69.8771 ur score 3.79109 cv score 911.019 gcv score 8.14922 roughness 250.887 penalty yi (in) 2159.87 penalty ys (out) 250.887 End of Data/data3_vs.out echo Data/data3i 1>&2 cat >Data/data3i <<'End of Data/data3i' 7 2 -1 1 1 2 3 4 5 6 7 10.1900 9.9168 19.6801 19.6389 30.2165 29.9494 28.9957 39.0335 30.4807 50.1602 19.5364 59.5550 9.9680 69.8931 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 100.0000 0 0 100.0000 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 3.0000 -1.5000 -1.5000 2.2500 100.0000 0 0 100.0000 End of Data/data3i echo Data/data3i_vs.out 1>&2 cat >Data/data3i_vs.out <<'End of Data/data3i_vs.out' -1 1 1 2 3 4 5 6 7 10.19 9.9168 19.6801 19.6389 30.2165 29.9494 28.9957 39.0335 30.4807 50.1602 19.5364 59.555 9.968 69.8931 0 3 -1.5 2.25 12.5691 11.023 22.6026 20.4729 30.0222 30.0015 31.5261 39.7298 27.9329 49.6374 20.0683 59.7132 10.0219 69.8771 ur score 3.79109 cv score 911.019 gcv score 8.14922 roughness 250.887 penalty yi (in) 2159.87 penalty ys (out) 250.887 End of Data/data3i_vs.out echo Data/data4 1>&2 cat >Data/data4 <<'End of Data/data4' 100 2 -14.14400000000000 -12.62400000000000 0.00872802734375 0.01776123046875 0.03652954101563 0.03778076171875 0.03804016113281 0.05033874511719 0.05046081542969 0.05799865722656 0.06211853027344 0.07255554199219 0.07278442382813 0.08236694335938 0.09878540039063 0.10624694824219 0.11671447753906 0.13525390625000 0.13577270507813 0.14022827148438 0.17289733886719 0.18214416503906 0.19433593750000 0.21131896972656 0.22599792480469 0.22840881347656 0.23481750488281 0.24028015136719 0.24079895019531 0.24317932128906 0.25314331054688 0.25602722167969 0.26106262207031 0.26126098632813 0.26403808593750 0.27488708496094 0.28308105468750 0.28486633300781 0.32711791992188 0.34045410156250 0.34812927246094 0.35504150390625 0.37854003906250 0.37930297851563 0.40206909179688 0.42366027832031 0.45243835449219 0.48316955566406 0.48445129394531 0.48989868164063 0.49433898925781 0.50692749023438 0.55101013183594 0.55981445312500 0.56024169921875 0.57824707031250 0.58763122558594 0.58871459960938 0.59646606445313 0.60411071777344 0.61349487304688 0.62188720703125 0.62496948242188 0.65248107910156 0.65376281738281 0.65525817871094 0.65875244140625 0.67666625976563 0.67686462402344 0.69070434570313 0.69076538085938 0.70008850097656 0.70338439941406 0.70542907714844 0.70745849609375 0.71948242187500 0.72564697265625 0.73609924316406 0.74693298339844 0.75987243652344 0.77017211914063 0.77409362792969 0.79241943359375 0.80749511718750 0.80964660644531 0.81486511230469 0.81591796875000 0.83602905273438 0.84744262695313 0.85531616210938 0.86418151855469 0.88067626953125 0.88789367675781 0.91096496582031 0.91664123535156 0.93675231933594 0.94479370117188 0.96264648437500 0.96868896484375 0.97003173828125 0.97642517089844 0.99330139160156 -1.12398497538699 -10.93178817794062 -3.17309365424752 -11.91723207564436 -4.27969385516744 -9.18600519088890 -5.10681925012676 -10.87013821377104 -2.05292317987267 -9.36133198594369 -6.31100010695271 -10.48891629832201 -8.98896849663828 -15.30358233616796 -8.23640607093421 -12.66339473143018 -5.67965665182554 -9.47015795730284 -7.22443958115417 -9.71984021134575 -7.95084001651953 -10.15225595454247 -8.99917958555789 -11.71945971075532 -8.31431689982978 -8.72415320136521 -8.76845119231010 -7.66350194719612 -9.13376938648027 -9.65290493278920 -6.26745661297852 -8.19053959358741 -8.94534543572409 -11.05767310381379 -11.88848559590654 -14.52345012895866 -7.18180278810573 -9.89277781490870 -6.72410924098239 -9.93054091245559 -5.36101552527228 -11.82760304437681 -4.54296780969114 -9.22863353341912 -1.83142749874624 -9.59526017505180 -3.68148699063362 -12.97103580634633 -1.56043751579175 -9.20820487995572 -0.74726910575254 -8.59612918852318 0.43579187877783 -8.72718092833845 -2.85793571731268 -11.65879549646217 -1.49253592736063 -12.43785371963662 2.91094827306969 -9.29753007303280 2.11659436320823 -8.13106997163606 0.36471226570417 -8.82210472257749 -0.24518282724765 -11.18500353006773 2.51772700505585 -10.28470111257078 3.95423775705867 -9.23131750685451 3.75726303965909 -8.57664926970630 7.82245608796869 -7.80190451752320 6.77795221925481 -8.38057519501939 6.97166427920412 -10.52922312090694 4.59273184493987 -12.08729399038268 7.06119135974869 -7.00462186564610 6.11879208809344 -7.50211461327799 3.94861482272061 -7.75520449817477 4.00540569985449 -5.35646651230471 2.67866417414087 -7.04138694487989 1.67002730869664 0.19817042256120 1.73414544348953 -1.30507292210797 0.17194446735798 -2.89237157977576 -0.44267173867388 -1.46614286363887 -2.27708671294323 -2.23615358864851 -4.96460170718095 1.98798911070409 -1.73359639073235 7.58565250861665 -5.17327378098329 3.92950780430648 -4.49051865894919 5.14418491970193 -1.87805069316382 7.68711640384167 -4.03140995488144 7.19279721847310 -2.49549389314198 9.26440488274794 -1.17195726725555 10.24894882825858 -4.55463111818626 7.59885874365692 -3.87089342733540 7.30866605324656 -4.09677360619903 10.43495045601807 -2.66827202225212 10.59632864254130 -3.95726114136361 6.54870806256987 -3.02159525246047 9.50258039220340 -1.58228166281845 10.57545071141379 -3.69342772319762 11.03657297754427 -2.13048660312614 11.13486924441208 -4.03099987094604 6.26836108589846 -0.42867403169730 11.79138545597352 -0.88867791027046 10.40096143363427 -2.89619641272182 6.83515227537886 -1.80410480914719 7.94162737372732 -3.01802917710480 10.14862808188024 2.93952429378248 14.89993516234490 -0.31692534824311 8.94967912549173 -1.49042932685630 7.98198800230913 0.08391972436066 9.80052380292071 0.14937734007512 9.02783111465752 1.87179134057138 11.85759376896628 -0.62431135534417 6.73695873726625 -0.70514841782816 8.82481668295866 2.50625632257919 12.28376603888035 0.37658379183970 8.74743489303171 3.39296270385137 11.24446985172180 -0.18326007699529 8.83463025359470 4.55458360024920 12.63549032058895 2.59134123432696 11.33391340551165 0.56363644487696 9.67596206940508 1.73211950273392 10.98938504718899 0.80282413258875 9.77336324549074 3.27689690999793 11.26607673763661 0.85076078049314 11.13773512597656 0.64804048943650 11.00631588502005 0.19018133990015 8.51645427916051 -0.64351814996506 9.58645257834105 -0.10426723894149 9.33816288748252 0.62740185183011 11.56160839947208 0.76121683766671 11.20751575708234 -1.36867229373008 8.94897208007103 0.31602611600448 11.01496086411661 1.23456790123457 -0.74074074074074 -0.74074074074074 0.69444444444444 End of Data/data4 echo Data/data4_ss.out 1>&2 cat >Data/data4_ss.out <<'End of Data/data4_ss.out' -8.48011 13.1262 96.072 98 End of Data/data4_ss.out echo Data/data5 1>&2 cat >Data/data5 <<'End of Data/data5' 1 3 4 8 16 17 19 1 2 3 3 2 1 2 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 2.7000 2.8000 2.9000 3.0000 3.1000 3.2000 3.3000 3.4000 3.5000 3.6000 3.7000 3.8000 3.9000 4.0000 4.1000 4.2000 4.3000 4.4000 4.5000 4.6000 4.7000 4.8000 4.9000 5.0000 5.1000 5.2000 5.3000 5.4000 5.5000 5.6000 5.7000 5.8000 5.9000 6.0000 6.1000 6.2000 6.3000 6.4000 6.5000 6.6000 6.7000 6.8000 6.9000 7.0000 7.1000 7.2000 7.3000 7.4000 7.5000 7.6000 7.7000 7.8000 7.9000 8.0000 8.1000 8.2000 8.3000 8.4000 8.5000 8.6000 8.7000 8.8000 8.9000 9.0000 9.1000 9.2000 9.3000 9.4000 9.5000 9.6000 9.7000 9.8000 9.9000 10.0000 10.1000 10.2000 10.3000 10.4000 10.5000 10.6000 10.7000 10.8000 10.9000 11.0000 11.1000 11.2000 11.3000 11.4000 11.5000 11.6000 11.7000 11.8000 11.9000 12.0000 12.1000 12.2000 12.3000 12.4000 12.5000 12.6000 12.7000 12.8000 12.9000 13.0000 13.1000 13.2000 13.3000 13.4000 13.5000 13.6000 13.7000 13.8000 13.9000 14.0000 14.1000 14.2000 14.3000 14.4000 14.5000 14.6000 14.7000 14.8000 14.9000 15.0000 15.1000 15.2000 15.3000 15.4000 15.5000 15.6000 15.7000 15.8000 15.9000 16.0000 16.1000 16.2000 16.3000 16.4000 16.5000 16.6000 16.7000 16.8000 16.9000 17.0000 17.1000 17.2000 17.3000 17.4000 17.5000 17.6000 17.7000 17.8000 17.9000 18.0000 18.1000 18.2000 18.3000 18.4000 18.5000 18.6000 18.7000 18.8000 18.9000 19.0000 19.1000 19.2000 19.3000 19.4000 19.5000 19.6000 19.7000 19.8000 19.9000 20.0000 20.1000 20.2000 20.3000 20.4000 20.5000 20.6000 20.7000 20.8000 20.9000 21.0000 21.1000 21.2000 21.3000 21.4000 21.5000 21.6000 21.7000 21.8000 21.9000 22.0000 22.1000 22.2000 22.3000 22.4000 22.5000 22.6000 22.7000 22.8000 22.9000 23.0000 23.1000 23.2000 23.3000 23.4000 23.5000 23.6000 23.7000 23.8000 23.9000 24.0000 End of Data/data5 echo Data/data5_ci.out 1>&2 cat >Data/data5_ci.out <<'End of Data/data5_ci.out' 0.709312 0.73838 0.767449 0.796518 0.825587 0.854656 0.883725 0.912793 0.941862 0.970931 1 1.02912 1.05856 1.08862 1.11962 1.15189 1.18572 1.22143 1.25934 1.29977 1.34302 1.38941 1.43925 1.49286 1.55055 1.61264 1.67944 1.75126 1.82842 1.91123 2 2.09477 2.19442 2.29756 2.4028 2.50873 2.61396 2.7171 2.81675 2.91152 3 3.08108 3.15475 3.22125 3.28085 3.33381 3.38037 3.4208 3.45536 3.4843 3.50788 3.52636 3.54 3.54905 3.55377 3.55441 3.55125 3.54452 3.53449 3.52143 3.50557 3.48719 3.46654 3.44387 3.41945 3.39353 3.36637 3.33823 3.30936 3.28002 3.25047 3.22097 3.19177 3.16313 3.13531 3.10857 3.08316 3.05934 3.03737 3.0175 3 2.98506 2.97264 2.96266 2.95501 2.9496 2.94634 2.94513 2.94588 2.94849 2.95287 2.95892 2.96655 2.97566 2.98616 2.99795 3.01094 3.02504 3.04015 3.05617 3.07301 3.09057 3.10877 3.1275 3.14668 3.16619 3.18596 3.20589 3.22588 3.24584 3.26566 3.28527 3.30456 3.32343 3.3418 3.35957 3.37664 3.39292 3.40832 3.42273 3.43607 3.44824 3.45914 3.46869 3.47678 3.48332 3.48822 3.49137 3.4927 3.4921 3.48947 3.48473 3.47777 3.46851 3.45684 3.44268 3.42592 3.40648 3.38426 3.35917 3.3311 3.29996 3.26567 3.22812 3.18723 3.14289 3.095 3.04349 2.98825 2.92918 2.86619 2.7992 2.72809 2.65278 2.57318 2.48918 2.4007 2.30763 2.20989 2.10738 2 1.88802 1.77316 1.6575 1.54312 1.43209 1.32649 1.2284 1.13991 1.06308 1 0.95227 0.919575 0.901124 0.896126 0.90379 0.923323 0.953935 0.994835 1.04523 1.10433 1.17134 1.24548 1.32595 1.41195 1.50271 1.59742 1.69529 1.79554 1.89738 2 2.10276 2.20551 2.30827 2.41102 2.51378 2.61654 2.71929 2.82205 2.9248 3.02756 3.13031 3.23307 3.33583 3.43858 3.54134 3.64409 3.74685 3.84961 3.95236 4.05512 4.15787 4.26063 4.36338 4.46614 4.5689 4.67165 4.77441 4.87716 4.97992 5.08268 5.18543 5.28819 5.39094 5.4937 5.59645 5.69921 5.80197 5.90472 6.00748 6.11023 6.21299 6.31575 6.4185 6.52126 6.62401 6.72677 6.82953 6.93228 7.03504 7.13779 End of Data/data5_ci.out echo numrecC/README 1>&2 cat >numrecC/README <<'End of numrecC/README' Sorry, we cannot redistribute the Numerical Recipes in C library. The only subroutines needed from it are: brent.c f1dim.c linmin.c mnbrak.c nrutil.c powell.c These subroutines simply do one- and multi-dimensional minimizations. See the Numerical Recipes in C book for the syntax. End of numrecC/README