#!/bin/sh
# This is a shell archive (produced by GNU sharutils 4.2).
# To extract the files from this archive, save it to some FILE, remove
# everything before the `!/bin/sh' line above, then type `sh FILE'.
#
# Made on 1999-08-02 17:27 WST by <handcock@jones>.
# Source directory was `/u/handcock/Projects/RDbook/WWW/library'.
#
# Existing files will *not* be overwritten unless `-c' is specified.
#
# This shar contains:
# length mode       name
# ------ ---------- ------------------------------------------
#    193 -rw-r-xr-x nlsd/Makefile
#   2121 -rw-r-xr-x nlsd/README
#    953 -rw-r-xr-x nlsd/.Data/.Help/rnlsds
#    953 -rw-r-xr-x nlsd/.Data/.Help/qnlsd
#    953 -rw-r-xr-x nlsd/.Data/.Help/pnlsd
#   1046 -rw-r-xr-x nlsd/.Data/.Help/nlsd.summary
#    783 -rw-r-xr-x nlsd/.Data/.Help/nlsd.plot
#   2877 -rw-r-xr-x nlsd/.Data/.Help/nlsd.fit
#    953 -rw-r-xr-x nlsd/.Data/.Help/dnlsd
#  84088 -rw-r-xr-x nlsd/ylsd.c
#   7130 -rw-r-xr-x nlsd/ylsd.S
#
save_IFS="${IFS}"
IFS="${IFS}:"
gettext_dir=FAILED
locale_dir=FAILED
first_param="$1"
for dir in $PATH
do
  if test "$gettext_dir" = FAILED && test -f $dir/gettext \
     && ($dir/gettext --version >/dev/null 2>&1)
  then
    set `$dir/gettext --version 2>&1`
    if test "$3" = GNU
    then
      gettext_dir=$dir
    fi
  fi
  if test "$locale_dir" = FAILED && test -f $dir/shar \
     && ($dir/shar --print-text-domain-dir >/dev/null 2>&1)
  then
    locale_dir=`$dir/shar --print-text-domain-dir`
  fi
done
IFS="$save_IFS"
if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED
then
  echo=echo
else
  TEXTDOMAINDIR=$locale_dir
  export TEXTDOMAINDIR
  TEXTDOMAIN=sharutils
  export TEXTDOMAIN
  echo="$gettext_dir/gettext -s"
fi
touch -am 1231235999 $$.touch >/dev/null 2>&1
if test ! -f 1231235999 && test -f $$.touch; then
  shar_touch=touch
else
  shar_touch=:
  echo
  $echo 'WARNING: not restoring timestamps.  Consider getting and'
  $echo "installing GNU \`touch', distributed in GNU File Utilities..."
  echo
fi
rm -f 1231235999 $$.touch
#
if mkdir _sh15024; then
  $echo 'x -' 'creating lock directory'
else
  $echo 'failed to create lock directory'
  exit 1
fi
# ============= nlsd/Makefile ==============
if test ! -d 'nlsd'; then
  $echo 'x -' 'creating directory' 'nlsd'
  mkdir 'nlsd'
fi
if test -f 'nlsd/Makefile' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/Makefile' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/Makefile' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/Makefile' &&
OBJECTS = ylsd.o
X
CFLAGS = -O
X
FFLAGS = -O
X
nlsd.o:	$(OBJECTS) 
X	ld -r -o nlsd.o $(OBJECTS) -lm 
X	Splus < ylsd.S
X
X
nlsd.so:	logsplineall.c mpack.f
X	Splus SHLIB -o nlsd.o logsplineall.c mpack.f
SHAR_EOF
  $shar_touch -am 1130072698 'nlsd/Makefile' &&
  chmod 0655 'nlsd/Makefile' ||
  $echo 'restore of' 'nlsd/Makefile' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/Makefile:' 'MD5 check failed'
3da7cb059eb1e729d5a27f2b06c43f97  nlsd/Makefile
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/Makefile'`"
    test 193 -eq "$shar_count" ||
    $echo 'nlsd/Makefile:' 'original size' '193,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/README ==============
if test -f 'nlsd/README' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/README' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/README' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/README' &&
This section contains programs for implementing
logspline density estimation with case weights.
The functions were written by Charles Kooperberg,
and are documented more fully in "Logspline density
estimation for censored data" by Charles Kooperberg
and Charles J. Stone (Journal of Computational and
Graphical Statistics, 1992, 1, 301-328), and "A study
of logspline density estimation" by the same authors
(Computational Statistics and Data Analysis, 1991,
12, 327-348.
X
They are based on ideas in Section 3 and 4 of "The
use of polynomial splines and their tensor products
in extended linear modeling," by Charles J. Stone,
Mark Hansen, Charles Kooperberg, and Young K. Truong,
Annals of Statistics, 25 (1997), 1371-1470.  The code
was adapted from "Logspline density estimation under
censoring and truncation," by Ja-Yong Koo, Charles
Kooperberg and Jinho Park, Scandinavian Journal of
Statistics, 26 (1999), no. 1, 87--105.  For additional
information see http://lynx.fhcrc.org/~clk.
X
This software was created by Charles Kooperberg,
who hereby gives permission to use the software for
noncommercial purposes and to freely distribute it
as long as this README file is intact.
X
Contains:
X
ylsd.c   ylsd.S   7 helpfiles
X
To use load the S-file, compile the c-file and
dyn.load() or dyn.load.shared().
X
This README written by  Mark S. Handcock to be used
as part of the software to the implement the methods
described in the book
"Relative Distribution Methods in the Social Sciences" 
by Mark S. Handcock and Martina Morris, 
Springer-Verlag, 1999, ISBN 0387987789 
X
Other software is available from  the Relative Distribution Website:
X
http://www.stat.psu.edu/RelDist
X
A link to the website is maintained by the publisher at:
X
http://www.springer-ny.com/stats
X
under the heading "Author/Editor Home Pages."
The website contains descriptions of the variables
and data file formats.
X
The authors can be reached via
X
Mark S. Handcock
Department of Statistics
The Pennsylvania State University
311 Joab L. Thomas Building
University Park, PA 16802-2111
ph:  814-865-3993
fax: 814-863-7114
internet: handcock@stat.psu.edu
SHAR_EOF
  $shar_touch -am 0802164999 'nlsd/README' &&
  chmod 0655 'nlsd/README' ||
  $echo 'restore of' 'nlsd/README' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/README:' 'MD5 check failed'
ef250f7ee3d43721322d1ac5371695f9  nlsd/README
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/README'`"
    test 2121 -eq "$shar_count" ||
    $echo 'nlsd/README:' 'original size' '2121,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/rnlsds ==============
if test ! -d 'nlsd/.Data'; then
  $echo 'x -' 'creating directory' 'nlsd/.Data'
  mkdir 'nlsd/.Data'
fi
if test ! -d 'nlsd/.Data/.Help'; then
  $echo 'x -' 'creating directory' 'nlsd/.Data/.Help'
  mkdir 'nlsd/.Data/.Help'
fi
if test -f 'nlsd/.Data/.Help/rnlsds' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/rnlsds' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/rnlsds' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/rnlsds' &&
X.BG
X.SH NAME
Logspline Density Estimation
X.SH DESCRIPTION
Density, cumulative probability, quantiles and  random samples from
a nlsd (new logspline) density.
X.SH USAGE
dnlsd(q, fit)
X.br
pnlsd(q, fit)
X.br
qnlsd(p, fit)
X.br
rnlsd(n, fit)
X.PP
X.AG q
vector of quantiles. Missing values (NAs) are allowed.
X.AG p
vector of probabilities. Missing values (NAs) are allowed.
X.AG n
sample size. If length(n) is larger than 1, then
length(n) random values are returned.
X.AG fit
a list like the output from nlsd.fit.
X.RT
densities (dnlsd), probabilities (pnlsd), quantiles (qnlsd),
or a random sample (rnlsd) from a nlsd density.
X.SH SIDE EFFECTS
The  function  rnlsd causes creation of the  
dataset .Random.seed if it does not already exist, otherwise its
value is updated.
X.SH DETAILS:
Elements of q or p that are missing will cause the
corresponding elements of the result to be missing.
X.SH SEE ALSO
nlsd.fit, nlsd.plot, nlsd.summary,
pnlsd, qnlsd, rnlsd.
X.WR
SHAR_EOF
  $shar_touch -am 1023032397 'nlsd/.Data/.Help/rnlsds' &&
  chmod 0655 'nlsd/.Data/.Help/rnlsds' ||
  $echo 'restore of' 'nlsd/.Data/.Help/rnlsds' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/rnlsds:' 'MD5 check failed'
ec7442d7833b973f6aa7f7726768b537  nlsd/.Data/.Help/rnlsds
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/rnlsds'`"
    test 953 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/rnlsds:' 'original size' '953,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/qnlsd ==============
if test -f 'nlsd/.Data/.Help/qnlsd' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/qnlsd' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/qnlsd' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/qnlsd' &&
X.BG
X.SH NAME
Logspline Density Estimation
X.SH DESCRIPTION
Density, cumulative probability, quantiles and  random samples from
a nlsd (new logspline) density.
X.SH USAGE
dnlsd(q, fit)
X.br
pnlsd(q, fit)
X.br
qnlsd(p, fit)
X.br
rnlsd(n, fit)
X.PP
X.AG q
vector of quantiles. Missing values (NAs) are allowed.
X.AG p
vector of probabilities. Missing values (NAs) are allowed.
X.AG n
sample size. If length(n) is larger than 1, then
length(n) random values are returned.
X.AG fit
a list like the output from nlsd.fit.
X.RT
densities (dnlsd), probabilities (pnlsd), quantiles (qnlsd),
or a random sample (rnlsd) from a nlsd density.
X.SH SIDE EFFECTS
The  function  rnlsd causes creation of the  
dataset .Random.seed if it does not already exist, otherwise its
value is updated.
X.SH DETAILS:
Elements of q or p that are missing will cause the
corresponding elements of the result to be missing.
X.SH SEE ALSO
nlsd.fit, nlsd.plot, nlsd.summary,
pnlsd, qnlsd, rnlsd.
X.WR
SHAR_EOF
  $shar_touch -am 1023032397 'nlsd/.Data/.Help/qnlsd' &&
  chmod 0655 'nlsd/.Data/.Help/qnlsd' ||
  $echo 'restore of' 'nlsd/.Data/.Help/qnlsd' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/qnlsd:' 'MD5 check failed'
ec7442d7833b973f6aa7f7726768b537  nlsd/.Data/.Help/qnlsd
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/qnlsd'`"
    test 953 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/qnlsd:' 'original size' '953,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/pnlsd ==============
if test -f 'nlsd/.Data/.Help/pnlsd' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/pnlsd' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/pnlsd' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/pnlsd' &&
X.BG
X.SH NAME
Logspline Density Estimation
X.SH DESCRIPTION
Density, cumulative probability, quantiles and  random samples from
a nlsd (new logspline) density.
X.SH USAGE
dnlsd(q, fit)
X.br
pnlsd(q, fit)
X.br
qnlsd(p, fit)
X.br
rnlsd(n, fit)
X.PP
X.AG q
vector of quantiles. Missing values (NAs) are allowed.
X.AG p
vector of probabilities. Missing values (NAs) are allowed.
X.AG n
sample size. If length(n) is larger than 1, then
length(n) random values are returned.
X.AG fit
a list like the output from nlsd.fit.
X.RT
densities (dnlsd), probabilities (pnlsd), quantiles (qnlsd),
or a random sample (rnlsd) from a nlsd density.
X.SH SIDE EFFECTS
The  function  rnlsd causes creation of the  
dataset .Random.seed if it does not already exist, otherwise its
value is updated.
X.SH DETAILS:
Elements of q or p that are missing will cause the
corresponding elements of the result to be missing.
X.SH SEE ALSO
nlsd.fit, nlsd.plot, nlsd.summary,
pnlsd, qnlsd, rnlsd.
X.WR
SHAR_EOF
  $shar_touch -am 1023032397 'nlsd/.Data/.Help/pnlsd' &&
  chmod 0655 'nlsd/.Data/.Help/pnlsd' ||
  $echo 'restore of' 'nlsd/.Data/.Help/pnlsd' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/pnlsd:' 'MD5 check failed'
ec7442d7833b973f6aa7f7726768b537  nlsd/.Data/.Help/pnlsd
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/pnlsd'`"
    test 953 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/pnlsd:' 'original size' '953,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/nlsd.summary ==============
if test -f 'nlsd/.Data/.Help/nlsd.summary' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/nlsd.summary' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/nlsd.summary' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/nlsd.summary' &&
X.BG
X.SH NAME
nlsd.summary
X.SH DESCRIPTION
nlsd.summary: summarizes various nlsd (new logspline) density estimates
X.SH USAGE
nlsd.summary(fit)
X.PP
X.AG fit
a list like the output from nlsd.fit
X.RT
This function produces only printed output. The main body
is a table with five columns: the first column is a possible number
of knots for the fitted model;
X.br
the second column is the log-likelihood for the fit;
X.br
the third column is -2*loglikelihood + penalty*(number of knots-1),
which is the AIC criterion - nlsd.fit selected the model with
the smallest value of AIC;
X.br
the fourth and fifth columns give the
endpoints of the interval of values of penalty that would yield the
model with the indicated number of knots. (NAs imply that the model is
not optimal for any choice of penalty.) At the bottom of the table the
number of knots corresponding to the selected model is reported, as is
the value of penalty that was used.
X.EX
fit <- nlsd.fit(y)      
X.br
nlsd.summary(fit)
X.SH SEE ALSO
nlsd.fit, nlsd.plot, dnlsd, pnlsd, qnlsd, rnlsd.
X.WR
SHAR_EOF
  $shar_touch -am 1023032297 'nlsd/.Data/.Help/nlsd.summary' &&
  chmod 0655 'nlsd/.Data/.Help/nlsd.summary' ||
  $echo 'restore of' 'nlsd/.Data/.Help/nlsd.summary' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/nlsd.summary:' 'MD5 check failed'
651326fc07bf6aee5bcca21db0daac0e  nlsd/.Data/.Help/nlsd.summary
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/nlsd.summary'`"
    test 1046 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/nlsd.summary:' 'original size' '1046,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/nlsd.plot ==============
if test -f 'nlsd/.Data/.Help/nlsd.plot' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/nlsd.plot' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/nlsd.plot' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/nlsd.plot' &&
X.BG
X.SH NAME
nlsd.plot
X.SH DESCRIPTION
nlsd.plot: plots a nlsd (new logspline) density, distribution function, hazard
function or survival function.
X.SH USAGE
nlsd.plot(fit, n=100, what="d", ...)
X.PP
X.AG fit
a list like the output from nlsd.fit.
X.AG n
the number of equally spaced points at which to plot the density.
X.AG what
what should be plotted: d (density), p (distribution function), s (survival
function) or h (hazard function).
X.AG ...
all regular plotting options as desired.
X.RT
This function produces a plot of a nlsd fit at n equally
spaced points roughly covering the support of the density. (Use
xlim=c(from,to) to change the range of these points.)
X.EX
fit <- nlsd.fit(y)      
X.br
nlsd.plot(fit)
X.SH SEE ALSO
nlsd.fit, nlsd.summary, dnlsd, pnlsd, qnlsd, 
rnlsd.
X.WR
SHAR_EOF
  $shar_touch -am 1023032097 'nlsd/.Data/.Help/nlsd.plot' &&
  chmod 0655 'nlsd/.Data/.Help/nlsd.plot' ||
  $echo 'restore of' 'nlsd/.Data/.Help/nlsd.plot' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/nlsd.plot:' 'MD5 check failed'
20339e9c5c6dc1c4626bb2064e0db0a2  nlsd/.Data/.Help/nlsd.plot
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/nlsd.plot'`"
    test 783 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/nlsd.plot:' 'original size' '783,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/nlsd.fit ==============
if test -f 'nlsd/.Data/.Help/nlsd.fit' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/nlsd.fit' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/nlsd.fit' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/nlsd.fit' &&
X.BG
X.SH NAME
nlsd.fit - nlsd density estimation
X.SH DESCRIPTION
nlsd.fit: density estimation using splines to
approximate the log-density.
X.SH USAGE
nlsd.fit(data, wgt, lbound, ubound, maxknots = 0, knots, nknots = 0,
penalty = -1, silent = T, mind = -1, wttype = 1)
X.EX
The simplest example is:
X.br
fit <- nlsd.fit(y)
X.br
nlsd.plot(fit)
X.PP
X.AG data
data vector
X.AG wgt
vector with case weights. The program does not check whether these
numbers are positive. Warning: if some weights are <=0 the behavior
of the program is unpredictable, in particular it will likely crash if
used in conjunction with wttype = 1.
X.AG lbound,ubound
lower/upper bound for the support of the density. For example, if there
is a priori knowledge that the density equals zero to the left of 0,
and has a discontinuity at 0,
the user could specify lbound=0. However, if the density is 
essentially zero near 0, one does not need to specify lbound.
X.AG maxknots
the maximum nuber of knots. The routine stops adding knots
when this number of knots is reached.
The method has an automatic rule
for selecting maxknots if this parameter is not specified.
X.AG knots 
ordered vector of values (that should cover the complete range of the
observations), which forces the method to start with these knots.
Overrules knots.	
If knots is not specified, a default knot-placement rule is employed.
X.AG nknots
forces the method to start with nknots knots.
The method has an automatic rule
for selecting nknots if this parameter is not specified.
X.AG penalty
the parameter to be used in the AIC criterion. The method chooses
the number of knots that minimizes -2*loglikelihood+penalty*(number of knots-1).
The default (-1)
is to use penalty=log(samplesize) as in BIC. The effect of
this parameter is summarized in nlsd.summary().
X.AG silent
should diagnostic output be printed?
X.AG mind
minimum distance, in order statistics, between knots.
X.AG wttype
if wttype=0 Rao and Wald statistics are computed as usual, and model selection
is carried out using BIC. Note that when weights can be rescaled, BIC
is not uniquely defined, and another model would be selected if the weights
were rescaled. If wttype=1 (the default) the Rao and Wald statistic are
divided by the caseweight of the observation at which the knot is located.
The model selection now always stops during knot deletion, at the moment
that the Wald statistic is larger than penalty (the default of penalty
is still log n). Note that if there are multiple observations at the
same location, with different weights, pretty much one of the weights
is taken at random.
X.RT
The output is organized to serve as input for nlsd.plot,
nlsd.summary, dnlsd, pnlsd, qnlsd and rnlsd.
X.SH EXAMPLES
estimate the density of a positive random variable,
and generate 50 random numbers from
the estimated density:
X.br
fit <- nlsd.fit(y,  lbound=0)
X.br
rnlsd(50, fit)
X.WR
SHAR_EOF
  $shar_touch -am 1023032097 'nlsd/.Data/.Help/nlsd.fit' &&
  chmod 0655 'nlsd/.Data/.Help/nlsd.fit' ||
  $echo 'restore of' 'nlsd/.Data/.Help/nlsd.fit' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/nlsd.fit:' 'MD5 check failed'
812149952f4b3f3d1766c41b516d827f  nlsd/.Data/.Help/nlsd.fit
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/nlsd.fit'`"
    test 2877 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/nlsd.fit:' 'original size' '2877,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/.Data/.Help/dnlsd ==============
if test -f 'nlsd/.Data/.Help/dnlsd' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/.Data/.Help/dnlsd' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/.Data/.Help/dnlsd' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/.Data/.Help/dnlsd' &&
X.BG
X.SH NAME
Logspline Density Estimation
X.SH DESCRIPTION
Density, cumulative probability, quantiles and  random samples from
a nlsd (new logspline) density.
X.SH USAGE
dnlsd(q, fit)
X.br
pnlsd(q, fit)
X.br
qnlsd(p, fit)
X.br
rnlsd(n, fit)
X.PP
X.AG q
vector of quantiles. Missing values (NAs) are allowed.
X.AG p
vector of probabilities. Missing values (NAs) are allowed.
X.AG n
sample size. If length(n) is larger than 1, then
length(n) random values are returned.
X.AG fit
a list like the output from nlsd.fit.
X.RT
densities (dnlsd), probabilities (pnlsd), quantiles (qnlsd),
or a random sample (rnlsd) from a nlsd density.
X.SH SIDE EFFECTS
The  function  rnlsd causes creation of the  
dataset .Random.seed if it does not already exist, otherwise its
value is updated.
X.SH DETAILS:
Elements of q or p that are missing will cause the
corresponding elements of the result to be missing.
X.SH SEE ALSO
nlsd.fit, nlsd.plot, nlsd.summary,
pnlsd, qnlsd, rnlsd.
X.WR
SHAR_EOF
  $shar_touch -am 1023031997 'nlsd/.Data/.Help/dnlsd' &&
  chmod 0655 'nlsd/.Data/.Help/dnlsd' ||
  $echo 'restore of' 'nlsd/.Data/.Help/dnlsd' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/.Data/.Help/dnlsd:' 'MD5 check failed'
ec7442d7833b973f6aa7f7726768b537  nlsd/.Data/.Help/dnlsd
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/.Data/.Help/dnlsd'`"
    test 953 -eq "$shar_count" ||
    $echo 'nlsd/.Data/.Help/dnlsd:' 'original size' '953,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/ylsd.c ==============
if test -f 'nlsd/ylsd.c' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/ylsd.c' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/ylsd.c' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/ylsd.c' &&
/* #define S_alloc calloc  */
#include <math.h>
#include <malloc.h>
#include <stdio.h>
#define MAXKNOTS 60
char *S_alloc(); 
X
struct datastruct { 
int ndata; 
double *data,*weights,dsw,wttype;
X   int *idata;
X   short *same;
};
X
/* ndata;
X   ndat - # number of cases
X   dat  - data
X   idat - the ips are the integration points: idat indicates
X          what the integration point immediately to the left of a datapoint is
X   same - is the observation the same as the previous in the same category? */
X
struct space {
X   int ndim,nk,nip,*iknots,ilow,iupp;
X   double *knots,aic,**info,*score,*ips,low,upp,cth;
X   struct basisfunct *basis;
};
X
/* ndim    - dimension
X   nk      - number of knots (=ndim+1)
X   nip     - number of integration points
X   iknots  - datapoint at or just left of knot
X   ilow    - is the lower bound -infinity? (1=yes)
X   iupp    - is the upper bound +infinity? (1=yes)
X   knots   - the knots
X   aic     - present value of aic
X   info    - the hessian
X   score   - score function
X   ips     - integration points
X   low     - lower integration boundary
X   upp     - upper integration boundary 
X   cth     - ctheta */
X
struct basisfunct {
X   double beta,*c1,**c2,sumunc;
X   int c3[2],iks[5];
};
X
/* beta   - coefficient
X   c1     - to translate the basis function in the truncated power basis
X   c2     - to translate the basisfunction at an integration point in a
X            polynomial
X   c3     - first and last integration point for which this function is nonzero
X   iks    - which knots are involved with this basisfunction - integrationpt 
X   sumunc - sum_i B(x_i) over the uncensored data */
X
X
int *isvector();
short *issvector();
double *dsvector(),**dsmatrix();
/* allocation */
X
int lusolve2();
void luinverse();
/* matrix inversion, solve a system */
X
double ctheta,*betaaddsorted;
double **kints,*cuu;
/* see piter - partial integrals and so on, which we want to keep */
struct basisfunct *bbx;
/* storage */
double ww6[7],yy6[7],ww7[33],yy7[33],*pompalcy,**pompalcyy;
int *rearix,*getiips,*luwi;
double *fiveee,*fiveh1,*fiverr,*betaaddv1,*betaremr1,*raoss,*luw,*luw2,**luww;
double *itertmp1,*itertmp2,*rearsorted,**solc1,**solc2,**solc3;
double **itertmp3,**pompcoef,**betaaddt1,**raoii,**raoc2,**betaremm2;
/******************************************************************************/
void nlogcensor(intpars,data0,wgt0,dpars,logs,ad,kts)
int *intpars,*ad;
double *data0,*wgt0,*logs,*kts,*dpars;
X
/* data0  - uncensored data; coefs on exit
X   intpars- integer parameters
X   dpars  - double parameters
X   ad     - is a model added during addition (1), deletion (2) or not at all
X   kts    - knots */
{
X   struct datastruct *data;
X   struct space *spc,*definespace();
X   void getsame(),quadalloc(),five();
X   int i,j,nlsd(),strt,mind,ndmax,silent;
X   double x,y,alpha,mylog();
X
/* data  - datastructure for all the data
X   spc   - datastructure for a model
X   definespace() - allocation for a model 
X   i,j   - counters
X   k     - one line for kdata
X   nlsd()- does the work
X   strt  - where starting knots provided
X   mind  - minimum distance between knots
X   ndmax - maximum dimension, the sign indicates whether it should be attained
X   silent- print diagnostic output? (1=yes)
X   alpha - penalty parameter
X   x,y   - utility */
X
/* we only want parameters and leave... */
X   if(intpars[0]<-10){
X      intpars[0]=MAXKNOTS+5;
X      return;
X   }
X
/* define the data */
X   data=(struct datastruct *)S_alloc((unsigned)1,sizeof(struct datastruct));
X   (*data).ndata=intpars[0];
X
X   (*data).data=data0;
X   (*data).weights=wgt0;
X   (*data).dsw=0.;
X   (*data).wttype=intpars[7];
X   for(i=0;i<(*data).ndata;i++)((*data).dsw)+=(*data).weights[i];
X   (*data).same=issvector(intpars[0]+1);
/* get the "same" vectors */
X   getsame(data0,intpars[0],(*data).same);
X   (*data).idata=isvector(intpars[0]);
/* allocate the space */
X   spc=definespace((*data).ndata);
X
X   getiips=isvector((*spc).nip+10);
X   luwi=isvector(2*MAXKNOTS+20);
X   rearix=isvector((*data).ndata);
X   fiverr=dsvector((*data).ndata+2*MAXKNOTS);
X   fiveee=dsvector((*data).ndata+MAXKNOTS+5);
X   fiveh1=dsvector((*data).ndata+MAXKNOTS+5);
X   betaaddv1=dsvector((*data).ndata+MAXKNOTS+5);
X   betaremr1=dsvector((*data).ndata+MAXKNOTS+5);
X   raoss=dsvector((*data).ndata+MAXKNOTS+5);
X   itertmp1=dsvector((*data).ndata+MAXKNOTS+5);
X   itertmp2=dsvector((*data).ndata+MAXKNOTS+5);
X   rearsorted=dsvector((*data).ndata+MAXKNOTS+5);
X   luw=dsvector(2*MAXKNOTS+20);
X   luw2=dsvector(2*MAXKNOTS+20);
X   itertmp3=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   solc1=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   solc2=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   solc3=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   luww=dsmatrix(2*MAXKNOTS+20,2*MAXKNOTS+20);
X   pompcoef=dsmatrix((*spc).nip+2,4);
X   betaaddt1=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   betaremm2=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   raoii=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   raoc2=dsmatrix((*spc).nip+10,(*spc).nip+10);
X
X   pompalcy=dsvector(2*MAXKNOTS+10);
X   betaaddsorted=dsvector((*data).ndata);
X   pompalcyy=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X
X   bbx=
X   (struct basisfunct *)S_alloc((unsigned)(MAXKNOTS),sizeof(struct basisfunct));
/* get the integer and double parameters */
X   ndmax=intpars[1];
X   mind=intpars[6];
X   if(mind<1){
X   mind=2.5*pow((double)(*data).ndata,(double)0.2)+0.5;
X   if((*data).ndata/mind<10)mind=(*data).ndata/10;
X   if(mind<3)mind=3;
X   }
X   intpars[6]=mind;
X   strt=intpars[2];
X   silent=intpars[3];
X   alpha=dpars[0];
X   if(strt==547) {
X       strt= floor(2.5*pow((double)intpars[0],(double)0.2));
X       if(strt>25)strt=25;
X       if(strt>intpars[0]/4)strt=intpars[0]/4;
X   }
X   if(alpha<0) alpha= mylog((double)intpars[0]);
X   (*spc).ilow=intpars[4];
X   (*spc).iupp=intpars[5];
X   (*spc).low=dpars[1];
X   (*spc).upp=dpars[2];
X
X   i=0;
/* starting knots */
X   if(ndmax==0)
X      ndmax = - floor(4.*pow((double)intpars[0],(double)0.2)+1);
X   if(ndmax>MAXKNOTS)ndmax=MAXKNOTS;
X   if(strt<=0){
X      if(intpars[2]<0)intpars[2]= -intpars[2];
X      else intpars[2]=floor(2.5*pow((double)intpars[0],(double)0.2)+1);
X      if(intpars[2]<0)intpars[2]= -intpars[2];
X      if(intpars[2]<3)intpars[2]=3;
X      five(data0,kts,intpars,(*data).same);
X      strt= intpars[2];
X   }
/* they were user provided */
X   if(strt>0){
X      (*spc).nk=strt;
X      (*spc).ndim=strt-1;
X      for(i=0;i<strt;i++)(*spc).knots[i]=kts[i];
X      strt=1;
/* find the iknots */
X      j=0;
X      y= -pow((double)10.,(double)100.);
X      for(i=0;i<(*data).ndata;i++){
X         x=(*data).data[i];
X         if(y<=(*spc).knots[j]&&x>(*spc).knots[j]){
X            (*spc).iknots[j]=i;
X            j++;
X            i--;
X            if(j==(*spc).nk)i=(*data).ndata+10;
X         }
X         else y=x;
X      }
X      if(j==((*spc).nk)-1)
X         (*spc).iknots[(*spc).nk-1]=(*data).ndata-1;
/* two knots outside the range of the data is not allowed */
X      if(j<((*spc).nk)-1){
X         intpars[0]=17;
X         return;
X      }
X      if((*spc).iknots[1]==0){
X         intpars[0]=18;
X         return;
X      }
X   }
/* allocations */
X   cuu = dsvector(MAXKNOTS+5);
X   kints = dsmatrix((*spc).nip+10,7);
X   quadalloc();
X
/* do the work */
X   intpars[0]=nlsd(spc,data,alpha,ndmax,mind,strt,silent,logs,ad);
X   dpars[0]= alpha;
/* output */
X   if(intpars[0]>0 && intpars[0]<100)return;
X   intpars[1]=(*spc).nk;
X   intpars[2]=(*spc).ndim;
X   for(i=0;i<((*spc).nk)+2;i++){
X      data0[i] = 0.;
X      for(j=0;j<(*spc).nk-1;j++)
X         data0[i]+=(*spc).basis[j].beta*(*spc).basis[j].c1[i];
X   }
X   data0[0]+=mylog((*spc).cth);
X   for(i=0;i<((*spc).nk);i++)kts[i]=(*spc).knots[i];
X   return;
}
/******************************************************************************/
/* the work */
int nlsd(best,data,alpha,ndmax,mind,strt,silent,logs,ad)
struct space *best;
struct datastruct *data;
double alpha,*logs;
int ndmax,mind,strt,silent,*ad;
X
/* best  - best space up to now
X   data  - the data
X   alpha - penalty parameter
X   ndmax - maximum dimension size: negative: does not have to be attained
X   mind  - minimum distance between knots
X   strt  - were starting knots provided (1=yes)
X   silent- should diagnostic info be printed? (1=yes)
X   logs  - log-likelihood of models 
X   ad    - fit during addition (1), deletion (2), not at all (0) */
{
X   struct space *current,*trynew,*definespace();
X   int add=1,adddim(),i,oops=0,ndm2,iter(),oops2=0,oops3=0,startspace(),j,icon;
X   int wttype=(*data).wttype; double wxx=0.,lold=0.,lold2=0.;
X   double xxa=0;
X   void swapspace(); double remdim(),wwd;
X
/* current - present space
X   trynew  - needed during addition and deletion
X   definespace()- allocates a space
X   add     - adding=1, deleting=something else
X   i       - counter
X   oops    - error status
X   ndm2    - sign of ndmax
X   iter()  - fits a model
X   swapspace() - copies a model
X   adddim()- adds a dimension
X   remdim()- removes a dimension
X   startspace()- the starting model */
X
/* allocates storage for spaces */
X   trynew=definespace((*data).ndata);
X   current=definespace((*data).ndata);
X
/* starting */
X   swapspace(current,best);
X   i=startspace(current,data,strt,silent);
X   if(i==0)return 39;
X
/* initialization */
X   ndm2=ndmax;
X   if(ndmax<0)ndmax= -ndmax;
X   (*best).aic=pow((double)10.,(double)150.);
X   for(i=0;i<MAXKNOTS;i++)logs[i]= -pow((double)10.,(double)150.);
X
/* we start in adding mode */
X   do{
X
/* fits the model */
X      if(oops3==0)oops=iter(current,data,silent,&xxa);
X      if((oops>0 || oops3>0)&& ndmax > (*current).ndim ){
/* problems. Jonge vriend, verzin een list! */
X         do{
X            for(i= -1;i> -4; i--){
/* begin opnieuw */
X               xxa=0.;
X               j=startspace(current,data,i,silent);
X               if(j==0)return 39;
X               oops=iter(current,data,silent,&xxa);
X               oops2++;
X               if(oops==0)i= -10;
X            }
X         }while(oops!=0 && ndmax > (*current).ndim);
X      }
X      if(oops2>2)oops2--;
X      if(oops!=0){ 
X         if((*best).aic< -1.0e149)return 40;
X         else swapspace(current,best);
X         add=0;
X      }
X      
X      if(oops==0){
/* compute aic */
X         if(wttype!=2) logs[(*current).ndim-1]=(*current).aic;
X         else logs[(*current).ndim-1]=lold2+wxx*((*current).aic-lold);
X         lold=(*current).aic;
X         lold2=logs[(*current).ndim-1];
X         ad[(*current).ndim-1]=1;
X         (*current).aic=(*current).ndim*alpha-2*logs[(*current).ndim-1];
X         if((*current).ndim==ndmax)add=0;
X         /* printf("%f %f\n",wxx,(*current).aic); */
/* did we improve */
X         if((*current).aic<(*best).aic) swapspace(best,current);
X      }
X
/* continue */
/*    if(add==1 && ndm2<0){   */
/* was there any recent improvement? */
/*       for(i=2;i<(*current).ndim-2;i++){
X            if(logs[(*current).ndim-1]-logs[i-1]<((*current).ndim-i)/2.-0.5){
X               add=0;
X               ndmax=(*current).ndim;
X            }
X         }
X      } */
/* adds dimensions, computes new starting values */
X      if(add==1){
X         add=adddim(current,trynew,data,mind,silent,&wxx);
X         if(add!=1 && oops2<2) ndmax=(*current).ndim;
X         if(add!=1 && oops2>=2){
X            oops3=1;
X            add=1;
X         }
X      }
X
/* keep on adding? */
X   }while(add==1);
X
/* start deleting */
X   if((*current).ndim>2)do{
X
/* removes dimensions, computes new starting values */
X      if(ndmax>2)wwd=remdim(current,data,trynew,silent,&wxx);
X
/* fits the model */
X      oops=iter(current,data,silent,&xxa);
X      if(oops!=0){
X         oops=oops+100;
X         (*best).ndim=ndmax-1;
X         return oops;
X      }
X      lold2+=wxx*((*current).aic-lold);
X      lold=(*current).aic;
X      if(wttype==2)(*current).aic=lold2;
X      
/* compute aic */
X      if((*current).aic>logs[(*current).ndim-1]){
X         logs[(*current).ndim-1]=(*current).aic;
X         ad[(*current).ndim-1]=2;
X      }
X      (*current).aic=(*current).ndim*alpha-2*(*current).aic;
X         /* printf("%f %f\n",wxx,(*current).aic); */
X
/* did we improve */
X      if(wttype!=1){if((*current).aic<(*best).aic) swapspace(best,current); }
X      else{if(wwd<alpha) swapspace(best,current); }
/* does further deleting make sense */
X      icon=1;
X      if(wttype==1 && wwd>=alpha)icon=0.;
X      if((*current).ndim==2)icon=0.;
X      if(wttype!=1&&
X         (*current).aic-(*best).aic>=alpha*((*current).ndim-2))icon=0;
X   }while(icon==1);
X
X   (*best).ndim=ndmax-1;
X   if(oops2>0)return 100;
X   return 0;
}
/******************************************************************************/
/* allocates storage for a space, and initializes elements */
struct space *definespace(nd)
int nd;
{
X   int i,j,k;
X   struct space *spc;
X   spc=(struct space *)S_alloc((unsigned)1,sizeof(struct space));
X   (*spc).aic=pow(10.,100.);
X   (*spc).ndim=0;
X   (*spc).nk=0;
X   (*spc).nip=0;
X   (*spc).iknots=isvector(MAXKNOTS+5);
X   (*spc).knots=dsvector(MAXKNOTS+5);
X   (*spc).score=dsvector(MAXKNOTS+5);
X   (*spc).info=dsmatrix(MAXKNOTS+5,MAXKNOTS+5);
X   k=MAXKNOTS+10+nd/100+300;
X   (*spc).ips=dsvector(k);
X   (*spc).basis=
X   (struct basisfunct *)S_alloc((unsigned)(MAXKNOTS),sizeof(struct basisfunct));
X   for(i=0;i<MAXKNOTS;i++){
X      (*spc).basis[i].beta=0.;
X      (*spc).basis[i].sumunc=0.;
X      for(j=0;j<2;j++)(*spc).basis[i].c3[j]=0;
X      (*spc).basis[i].c1=dsvector(MAXKNOTS+5);
X      (*spc).basis[i].c2=dsmatrix(k,4);
X   }
X   for(j=0;j<5;j++)(*spc).basis[i].iks[j]=0;
X   (*spc).nip=k;
X   return spc;
}
/******************************************************************************/
/* figure out which datapoints are the same.... */
void getsame(x,n,s)
double *x;
int n;
short *s;
{
X   int i;
X   double r;
X   s[0]=0;
X   for(i=1;i<n;i++){
/* exactly the same */
X      if(x[i]==x[i-1])s[i]=1;
X      else {
X         if(x[i]!=0){
/* almost exactly the same */
X            r=fabs(x[i-1]/x[i]-1.);
X            if(r<0.0000000000001)s[i]=1;
X            else s[i]=0;
X         }
X         else s[i]=0;
X      }
X   }
}
/******************************************************************************/
main()
{
X   int *intpars,*ad,i,maxn;
X   double *data,*logs,*kts,*dpars,rr,*wgt;
X   FILE *f1;
X   void nlogcensor(),exit();
X   maxn=500;
X   data=dsvector(8000); logs=dsvector(maxn); wgt=dsvector(8000);
X   kts=dsvector(maxn); intpars=isvector(maxn); ad=isvector(maxn);
X   dpars=dsvector(maxn);
X   for(i=0;i<maxn;i++){
X      data[i]=0.; logs[i]=0.; ad[i]=0.; kts[i]=0.; intpars[i]=0.; dpars[i]=0.;
X   }
X   intpars[0]=100; /* ndata */
X   intpars[1]=0;   /* maxdim - becomes maxmcmc */
X   intpars[2]=0;   /* start */
X   intpars[3]=1;   /* silent */
X   intpars[4]=1;   /* ilow */
X   intpars[5]=1;   /* iupp */
X   maxn=intpars[0]; 
X   dpars[0]= -1.; /* alpha */
X   dpars[1]=0.;   /* lbound */
X   dpars[2]=1.;   /* ubound */
X   f1=fopen("uu","r");
X   for(i=0;i<intpars[0];i++){ (void)fscanf(f1,"%lf",&rr); data[i]=rr; }
X   for(i=0;i<intpars[0];i++) wgt[i]=1.;
X   nlogcensor(intpars,data,wgt,dpars,logs,ad,kts);
X   exit(1);
}
void five(data,kts,intpars,same)
double *data,*kts;
int *intpars;
short *same;
{
X   int n1,k,l,i,j,g1,g2,md;
X   double *rr,*ee,*h1,h2,h3;
X   void five00(),five01();
X   n1 = intpars[0];
X   k=intpars[2];
X   md=intpars[6];
X   rr=fiverr;
X   if(intpars[4]+intpars[5]>=1)five00(rr,k,n1,md);
X   if(intpars[4]+intpars[5]==5)five01(rr,k,n1,intpars[5],md);
X   if(intpars[4]+intpars[5]==0)for(i=0;i<k;i++)rr[i]=(double)i/(double)(k-1);
X   ee=fiveee;
X   h1=fiveh1;
X   ee[0]=1.;
X   h1[0]=0.;
X   h1[1]=1.;
X   l=0;
X   for(i=1;i<n1;i++){
X      if(same[i]==1){
X         h1[l+1]++;
X      }
X      else{
X         l++;
X         ee[l]=data[i];
X         h1[l+1]=h1[l]+1;
X     }
X   }
X   for(i=0;i<l;i++) h1[i]=0.5*(h1[i]+h1[i+1]);   
X   h2=h1[0];
X   h3=0.;
X   for(i=0;i<l;i++) {
X      h1[i]-=h2;
X      if(h1[i]>h3)h3=h1[i];
X   }
X   for(i=0;i<l;i++) h1[i]/=h3;
X   kts[0]=data[0];
X   kts[k-1]=data[n1-1];
X   for(i=1;i<k-1;i++){
X      g1=0;
X      g2=l;
X      for(j=0;j<l;j++){
X         if(h1[j]<=rr[i] && j>g1)g1=j;
X         if(h1[j]>rr[i] && j<g2)g2=j;
X      }
X      h2=h1[g1];
X      h3=h1[g2];
X      h3=(rr[i]-h2)/(h3-h2);
/*    kts[i]=h3*ee[g2] + (1 - h3) * ee[g1]; */
X      if(h3>0.5)kts[i]=ee[g2];
X      else kts[i]=ee[g1];
X   }
}     
void five01(rr,k,n,il,md)
int k,n,il,md;
double *rr;
{
X   void five00();
X   int i;
X   five00(rr,2*k-1,2*n,md);
X   for(i=0;i<k;i++)rr[i]=2*rr[i];
X   if(il==0) for(i=0;i<k;i++)rr[i]=1-rr[2*k-2-i];
} 
void five00(rr,k,n,md)
int k,n,md;
double *rr;
{
X   double fi=4.,eps1,eps2,eps,s,w,v;
X   int i,i1,j,j2;
X   if((double)md-1.>fi)fi=md-1.;
X   j=floor((double)((k-1.)/2.+0.99));
X   j2=floor((double)((k-1.)/2.));
X   eps1=fi - pow((double)((n-1.)/fi),(double)(1./(j-1.)));
X   if(eps1>0) eps1 = 0;
X   eps2=fi-1;
X   for(i=0;i<k;i++)rr[i]=0;
X   rr[0]=1.;
X   rr[k]=n;
X   for(i1=0;i1<100;i1++)if(i1==0||eps2-eps1>0.0001){
X      eps = (eps1+eps2)/2.;
X      s=1;
X      w=fi;
X      for(i=1;i<=j2;i++){
X         v=i;
X         s+=w;
X         rr[i]=s;
X         rr[k-i-1]=n+1.-s;
X         v=fi-v*eps;
X         if(v<1)v=1;
X         w*=v;
X      }
X      if(2*j==k)s+=w/2.;
X      else rr[j]=(n+1)/2.;
X      if(2.*s>=n+1)eps1=eps;
X      else eps2=eps;
X   }
X   else i1=100;
X   for(i=0;i<k;i++)rr[i]=(rr[i]-1)/(n-1.);
} 
/******************************************************************************/
void luinverse(a,n)
double **a;
int n;
{
X   int k,j,*i,ludcmp();
X   void lubksb();
X   double *d,**c,r;
X   c=luww;
X   d=luw2;
X   i=luwi;
X   for(j=1;j<=n;j++) for(k=1;k<=n;k++) c[j][k]=a[j-1][k-1];
X   if(ludcmp(c,n,i,&r)==0)(void)printf("singular matrix in routine LUDCMP\n");
X   for(j=1;j<=n;j++){
X      for(k=1;k<=n;k++) d[k]=0.;
X      d[j]=1.;
X      lubksb(c,n,i,d);
X      for(k=1;k<=n;k++) a[k-1][j-1]=d[k];
X   }
}
/******************************************************************************/int lusolve2(a,n,b)
double **a,*b;
int n;
{
X   int *i,j,k,ludcmp();
X   void lubksb();
X   double **c,r;
X   c=luww;
X   i=luwi;
X   for(j=0;j<(n+1);j++)i[j]=0;
X   for(j=0;j<n;j++) for(k=0;k<n;k++) c[j+1][k+1]=a[j][k];
X   if(ludcmp(c,n,i,&r)==0)return 0;
X   lubksb(c,n,i,b-1);
X   return 1;
}
/******************************************************************************/
void lubksb(a,n,indx,b)
double **a,b[];
int n,*indx;
{
X   int i,ii=0,ip,j;
X   double sum;
X   for (i=1;i<=n;i++) {
X      ip=indx[i];
X      sum=b[ip];
X      b[ip]=b[i];
X      if (ii) for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
X      else if (sum) ii=i;
X      b[i]=sum;
X   }
X   for (i=n;i>=1;i--) {
X      sum=b[i];
X      for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
X      b[i]=sum/a[i][i];
X   }
}
/******************************************************************************/
#define TINY 1.0e-20;
int ludcmp(a,n,indx,d)
int n,*indx;
double **a,*d;
{
X   int i,imax,j,k;
X   double big,dum,sum,temp,*vv;
X   vv=luw;
X   for(i=0;i<=n+1;i++)vv[i]=0.;
X   *d=1.0;
X   for (i=1;i<=n;i++) {
X      big=0.0;
X      for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp;
X      if (big == 0.0) return 0;
X      vv[i]=1.0/big;
X   }
X   for (j=1;j<=n;j++) {
X      for (i=1;i<j;i++) {
X         sum=a[i][j];
X         for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
X         a[i][j]=sum;
X      }
X      big=0.0;
X      for (i=j;i<=n;i++) {
X         sum=a[i][j];
X         for (k=1;k<j;k++) sum -= a[i][k]*a[k][j];
X         a[i][j]=sum;
X         if ( (dum=vv[i]*fabs(sum)) >= big) {
X            big=dum;
X            imax=i;
X         }
X      }
X      if (j != imax) {
X         for (k=1;k<=n;k++) {
X            dum=a[imax][k];
X            a[imax][k]=a[j][k];
X            a[j][k]=dum;
X         }
X         *d = -(*d);
X         vv[imax]=vv[j];
X      }
X      indx[j]=imax;
X      if (a[j][j] == 0.0) a[j][j]=TINY;
X      if (j != n) {
X         dum=1.0/(a[j][j]);
X         for (i=j+1;i<=n;i++) a[i][j] *= dum;
X      }
X   }
X   return 1;
}
#undef TINY
/******************************************************************************/
int adddim(spc,spc2,data,mind,silent,wxx)
struct space *spc,*spc2;
struct datastruct *data;
int mind,silent;
double *wxx;
{
X   int i,nx,uu,ll,nowloc1,findr(),findl(),findm(),loloc,uploc,bestloc= -1;
X   int besti= -1,findyl(),findyr(),nowloc2;
X   double *sorted,nowrao1,rao(),bestrao= -1.,nowrao2;
X   void betaadd(),setupspace(),swapspace(),save3coden();
X   sorted=betaaddsorted;
X   swapspace(spc2,spc);
X   for(i=0;i<(*data).ndata;i++)
X      sorted[i]=(*data).data[i];
X   nx=(*data).ndata;
/* find the interval */
X   for(i=0;i<=(*spc).nk;i++){
/* before first knot */
X      if(i==0) nowloc1=findl(&ll,&uu,mind,sorted,nx,(*spc).knots[0]);
/* after last knot */
X      if(i==(*spc).nk) 
X         nowloc1=findr(&ll,&uu,mind,sorted,nx,(*spc).knots[(*spc).nk-1]);
/* in between knots */
X      if(i>0 && i<(*spc).nk)nowloc1=
X         findm(&ll,&uu,mind,sorted,nx,(*spc).knots[i-1],(*spc).knots[i]);
/* possible location */
X      if(nowloc1>=0){
X         nowrao1=rao(spc,data,sorted[nowloc1],nowloc1);
X         if(nowrao1>bestrao){
X            loloc=ll;
X            uploc=uu;
X            bestloc=nowloc1;
X            bestrao=nowrao1;
X            besti=i;
X         }
X      }
X   }
X   if(bestloc<0)return -1;
/* as long as the locations are different, do interval halving */
X   do{
X      if(sorted[uploc]>sorted[loloc]){
X         nowloc2=findyr(uploc,bestloc,sorted);
/* two search points, the upper one */
X         if(nowloc2>=0) nowrao2=rao(spc,data,sorted[nowloc2],nowloc2);
X         else nowrao2=bestrao;
/* two search points, the lower one */
X         nowloc1=findyl(bestloc,loloc,sorted);
X         if(nowloc1>=0) nowrao1=rao(spc,data,sorted[nowloc1],nowloc1);
X         else nowrao1=bestrao;
/* the middle one is the best, we call it quits */
X         if(bestrao>=nowrao2 && bestrao>=nowrao1) loloc=uploc;
X         else{
/* the lower search point is the best */
X            if(nowrao1>bestrao){
X               uploc=bestloc;
X               bestloc=nowloc1;
X               bestrao=nowrao1;
X            }
/* the upper search point is the best */
X            else{
X               loloc=bestloc;
X               bestloc=nowloc2;
X               bestrao=nowrao2;
X            }
X         }
X      }
X   }while(sorted[uploc]>sorted[loloc]);
/* failure */
X   if(bestloc<0)return bestloc;
X   (*wxx)=(*data).weights[bestloc];
/* done record the new knot in its correct position */
X   if(besti==(*spc).nk){
X       (*spc).knots[(*spc).nk]=sorted[bestloc];
X       (*spc).iknots[(*spc).nk]=bestloc;
X   }
X   else{
X      for(i=(*spc).nk;i>besti;i=i-1){
X         (*spc).knots[i]=(*spc).knots[i-1];
X         (*spc).iknots[i]=(*spc).iknots[i-1];
X      }
X      (*spc).knots[besti]=sorted[bestloc];
X      (*spc).iknots[besti]=bestloc;
X   }
X   ((*spc).nk)++;
X   ((*spc).ndim)++;
X   if(silent==1) (void)printf("add(%.2f), rao=%.2f ",sorted[bestloc],bestrao);
/* get (*spc).ips (*spc).nip (*data).idatx */
/* get (*spc).basis.c1 (*spc).basis.c2 (*spc).basis.c3 (*spc).basis.sumunc */
X   setupspace(spc,data);
/* get (*spc).basis.beta */
X   betaadd(spc,spc2,besti);
X   return 1;
}
/******************************************************************************/
/* finds location in an interval (l,b) - l might not have been tested yet */
int findyl(u,l,x)
int l,u;
double *x;
{
X   int i;
X   if(x[l]==x[u])return -1;
X   i=(u+l-1)/2;
X   if(x[i]!=x[u])return i;
X   i=(i+l)/2;
X   if(x[i]!=x[u])return i;
X   return l;
}
/******************************************************************************/
/* finds location in an interval (b,u) - u might not have been tested yet */
int findyr(u,l,x)
int l,u;
double *x;
{
X   int i;
X   if(x[l]==x[u])return -1;
X   i=(u+l+1)/2;
X   if(x[i]!=x[l])return i;
X   i=(i+u)/2;
X   if(x[i]!=x[l])return i;
X   return u;
}
/******************************************************************************/
/* Finds a possible location for a knot on the interval (0,knot1)
X   ll - lowest number we can search on in the future
X   uu - highest number we can search on in the future
X   mind minimum distance between knots
X   x  - data
X   nx - length of data
X   knt- knot */
int findl(ll,uu,mind,x,nx,knt)
double *x,knt;
int nx,*ll,*uu,mind;
{
/* dlocation - finds uu */
X   int i,dlocation();
X   (*uu)=dlocation(0,x,nx,knt);
X   if((*uu)<mind)return -1;
X   i=((*uu)-1)/2;
X   if((*uu)-i<mind+1)i=(*uu)-mind-1;
X   *ll=0;
X   *uu=(*uu)-mind-1;
X   return i;
}
/******************************************************************************/
/* Finds a possible location for a knot on the interval (knot-last,nx-1)
X   ll - lowest number we can search on in the future
X   uu - highest number we can search on in the future
X   mind minimum distance between knots
X   x  - data
X   nx - length of data
X   knt- knot */
int findr(ll,uu,mind,x,nx,knt)
double *x,knt;
int nx,*ll,*uu,mind;
{
/* dlocation - finds ll */
X   int i,dlocation();
X   (*ll)=dlocation(1,x,nx,knt);
X   if(nx-1-(*ll)<mind)return -1;
X   i=(nx+(*ll))/2;
X   if(i-(*ll)<mind+1)i=(*ll)+mind+1;
X   *uu=nx-1;
X   *ll=(*ll)+mind+1;
X   return i;
}
/******************************************************************************/
/* Finds a possible location for a knot on the interval (k0,k1)
X   ll - lowest number we can search on in the future
X   uu - highest number we can search on in the future
X   mind minimum distance between knots
X   x  - data
X   nx - length of data
X   k0 - knot
X   k1 - knot */
int findm(ll,uu,mind,x,nx,k0,k1)
double *x,k0,k1;
int nx,*ll,*uu,mind;
{
/* dlocation - finds ll */
X   int dlocation();
X   (*ll)=dlocation(1,x,nx,k0);
X   (*uu)=dlocation(0,x,nx,k1);
X   if((*uu)-(*ll)<2*mind+1)return -1;
X   *uu=(*uu)-mind-1;
X   *ll=(*ll)+mind+1;
X   return ((*uu)+(*ll))/2;
}
/******************************************************************************/
/* finds the lowest (if what = 0) or the highest (if what = 1) index of x for
X   which x==k
X   what - see above
X   x    - data
X   nx   - length data
X   k    - see above */
int dlocation(what,x,nx,k)
int nx,what;
double k,*x;
{
X   int i;
X   if(what==1){
X      if(x[0]>k)return 0;
X      if(x[nx-1]<=k)return nx-1;
X      for(i=0;i<nx-1;i++) if(x[i+1]>k && x[i]<=k) return i;
X   }
X   if(x[nx-1]<k)return nx-1;
X   if(x[0]>=k)return 0;
X   for(i=1;i<nx;i++) if(x[i]>=k && x[i-1]<k)return i;
X   return nx;
}
/******************************************************************************/
/* gets the new starting values, solve one coef matrix w.r.t. the other */
void betaadd(spc,spc2,besti)
struct space *spc,*spc2;
int besti;
{
X   int i,j,k;
X   double **t1,*v1;
X   k=(*spc).ndim+3;
X   t1=betaaddt1;
X   v1=betaaddv1;
X   for(i=0;i<((*spc2).nk)+2;i++){
X      v1[i] = 0.;
X      for(j=0;j<(*spc2).ndim;j++)
X         v1[i]+=(*spc2).basis[j].beta*(*spc2).basis[j].c1[i];
X   }
X   for(i=(*spc2).nk;i>besti;i=i-1) v1[i+2]=v1[i+1];
X   v1[besti+2]=0.;
X   for(i=0;i<k;i++)for(j=0;j<k;j++)t1[i][j]=0.;
X   for(j=0;j<k;j++)for(i=0;i<k-3;i++)t1[i][j]=(*spc).basis[i].c1[j];
X   for(i=k-2;i<k;i++)t1[i][i]=1.;
X   t1[k-3][0]=1.;
X   luinverse(t1,k);
X   for(i=0;i<k-3;i++){
X      (*spc).basis[i].beta=0.;
X      for(j=0;j<k;j++)(*spc).basis[i].beta+=t1[j][i]*v1[j];
X   }
X   return;
}
/******************************************************************************/
/* top iteration - governs bounds */
int iter(spc,data,silent,xxa)
int silent;
double *xxa;
struct space *spc;
struct datastruct *data;
{
X   double xpp=(*spc).upp,lpp=(*spc).upp,lll=(*spc).low,xll=(*spc).low;
X   int n=0,oops,iterx(),mll=(*spc).ilow,muu=(*spc).iupp;
/* bounded intervals */
X   if((*spc).ilow==0 && (*spc).iupp==0) return iterx(spc,data,silent,xxa);
X   do{
/* try unbounded intervals */
X      n++;
X      if(((*spc).basis[1].beta<0||muu==0)&&((*spc).basis[0].beta<0||mll==0)){
X         (*spc).low=xll;
X         (*spc).ilow=mll;
X         (*spc).upp=xpp;
X         (*spc).iupp=muu;
X         oops=iterx(spc,data,silent,xxa);
X         if(oops==0 || n==6)return oops;
X      }
X      (*spc).iupp=0;
X      (*spc).ilow=0;
/* widen bounded intervals */
X      if(muu==1)(*spc).upp=4*lpp-3*(*spc).low;
X      if(mll==1)(*spc).low=4*lll-3*(*spc).upp;
X      lpp=(*spc).upp;
X      lll=(*spc).low;
X      oops=iterx(spc,data,silent,xxa);
X      (*spc).iupp=muu;
X      (*spc).upp=xpp;
X      (*spc).ilow=mll;
X      (*spc).low=xll;
X      if(oops!=0)return oops;
X   }while(n<6);
X   return 9999;
}
/******************************************************************************/
/* the works */
int iterx(spc,data,silent,xxa)
int silent;
double *xxa;
struct space *spc;
struct datastruct *data;
X
/* spc   - model
X   data  - data
X   silent- print info? (1=yes) */
X
{
X   double error=0.000001,factor,logl,lnew,pompall(),mylog(),myexp();
X   double *tmp1,*tmp2,**tmp3;
X   int iter,maxiter=100,j,i,i1,i2,kk;
X   void startnow();
X
/* error - convergence criterion
X   j     - counter
X   logl  - log-likelihood
X   lnew  - new loglikelihood
X   iter  - iteration counter
X   factor- stepsize for the NR algorithm
X   evmin/evmax - minimum and maximum eigenvalue
X   pompall()   - compute logl/score/hessian
X   maxiter     - maximum number of iterations */
X
X   for(iter=0;iter<maxiter;iter++){
X      tmp1=itertmp1;
X      tmp2=itertmp2;
X      tmp3=itertmp3;
/* compute logl/score/hessian */
X      logl=pompall(spc,data,2,&i1);
X      if(iter==0 && fabs((*xxa))>0.01 && (*xxa)-logl > 100){
/* try alternate starting values */
X         lnew=logl;
X         for(j=0;j<(*spc).ndim;j++){
X            tmp1[j]=(*spc).score[j];
X            tmp2[j]=(*spc).basis[j].beta;
X            (*spc).basis[j].beta=0.;;
X            for(i=0;i<(*spc).ndim;i++)tmp3[i][j]=(*spc).info[i][j];
X         }
X         startnow(spc,data);
X         logl=pompall(spc,data,2,&i2);
X         if(lnew>logl ){
X            logl=lnew;
X            for(j=0;j<(*spc).ndim;j++){
X               (*spc).score[j]=tmp1[j];
X               (*spc).basis[j].beta=tmp2[j];
X               for(i=0;i<(*spc).ndim;i++)(*spc).info[i][j]=tmp3[i][j];
X            }
X         }
X         else i1=i2;
X      }
/* serious ctheta problems */
X      if(i1==1)return 7;
/* solve the system */
X      j=lusolve2((*spc).info,(*spc).ndim,(*spc).score);
/*    return 2 - something wrong with system */
X      if(j==0) return 2;
/* adjust the tail shifts */
X      if((*spc).ilow==1){
X         (*spc).score[0]= -(*spc).score[0]/(*spc).basis[0].beta;
X         if((*spc).score[0]<-100)(*spc).score[0]=-100;
X      }
X      if((*spc).iupp==1){
X         (*spc).score[1]= -(*spc).score[1]/(*spc).basis[1].beta;
X         if((*spc).score[1]<-100)(*spc).score[1]=-100;
X      }
/* find the right step size */
X      factor= -1.;
/* tail check */
X      if((*spc).ilow==1 && (*spc).iupp==1 &&
X         (*spc).basis[0].beta==0 && (*spc).basis[1].beta==0)return 6;
X      if((*spc).ilow==1 && (*spc).basis[0].beta>=0)return 4;
X      if((*spc).iupp==1 && (*spc).basis[1].beta>=0)return 3;
/* adjust beta */
X      if((*spc).ilow==0)(*spc).basis[0].beta-=factor*(*spc).score[0];
X      else (*spc).basis[0].beta= 
X                 -myexp(factor*(*spc).score[0]+mylog(-(*spc).basis[0].beta));
X      if((*spc).iupp==0)(*spc).basis[1].beta-=factor*(*spc).score[1];
X      else (*spc).basis[1].beta= 
X                 -myexp(factor*(*spc).score[1]+mylog(-(*spc).basis[1].beta));
X      for(j=2;j<(*spc).ndim;j++)(*spc).basis[j].beta-=factor*(*spc).score[j];
X      do{
/* new logl */
X         if((*spc).ilow==1 && (*spc).iupp==1 &&
X            (*spc).basis[0].beta==0 && (*spc).basis[1].beta==0)return 6;
X         if((*spc).ilow==1 && (*spc).basis[0].beta>=0)return 4;
X         if((*spc).iupp==1 && (*spc).basis[1].beta>=0)return 3;
X         lnew=pompall(spc,data,0,&i);
/* did we win? */
X         kk=0;
X         if((lnew-logl)< -error)kk=1;
X         if((lnew-logl)< -error * 100 && (*spc).ilow==1 &&
X               (*spc).basis[0].beta> -1.e8 )kk=1;
X         if((lnew-logl)< -error * 100 && (*spc).iupp==1 &&
X            (*spc).basis[1].beta> -1.e8 )kk=1;
X         if(kk==1 || (i==1 && fabs(factor)>0.1)){
/* adjust the stepsize */
X            i=0;
X            factor=factor/2.;
X            if((*spc).ilow==0)(*spc).basis[0].beta+=factor*(*spc).score[0];
X            else (*spc).basis[0].beta= 
X                   -myexp(-factor*(*spc).score[0]+mylog(-(*spc).basis[0].beta));
X            if((*spc).iupp==0)(*spc).basis[1].beta+=factor*(*spc).score[1];
X            else (*spc).basis[1].beta= 
X                   -myexp(-factor*(*spc).score[1]+mylog(-(*spc).basis[1].beta));
X            for(j=2;j<(*spc).ndim;j++)
X               (*spc).basis[j].beta+=factor*(*spc).score[j];
X            if(fabs(factor)< 0.00001 && 
X              (((*spc).iupp==1 && (*spc).basis[1].beta> -1.e8) ||
X               ((*spc).ilow==1 && (*spc).basis[0].beta> -1.e8))) return 5;
X            if(fabs(factor)< 0.00001) return 8;
/*             return 5/8 - too much step-halving */
X         }
X         if(i==1)return 7;
X      } while(kk==1);
/* convergence */
X      if(fabs(lnew-logl)<error && fabs(factor) > 0.96 )iter=maxiter+1000;
X      if(fabs(lnew-logl)<error * 100 && (*spc).ilow==1 &&
X            (*spc).basis[0].beta> -1.e8 )iter=maxiter+1000;
X      if(fabs(lnew-logl)<error * 100 && (*spc).iupp==1 &&
X            (*spc).basis[1].beta> -1.e8 )iter=maxiter+1000;
X   }
X   if(iter<maxiter+500) return 1;
/* return 1 - no convergence */
X   logl=pompall(spc,data,1,&i);
X   (*xxa)=logl;
X   (*spc).aic=logl;
X   if(silent==1){
X      (void)printf("|| logl= %.2f (nd=%d)\n",logl,(*spc).ndim);
X      (void)fflush(stdout);
X   }
X   (*spc).cth=ctheta;
X   return 0;
}
/******************************************************************************/
double pompall(spc,data,what,xp)
int what,*xp;
struct space *spc;
struct datastruct *data;
{
X   double *ips,f,mylog(),myexp(),logl,pol3(),inp3(),mat3();
X   double *cy,**cyy,**coef;
X   void l1int(),l2int(),m1int(),savecode1(),initk();
X   int i,j,k,nip=(*spc).nip,ndim=(*spc).ndim,savecoden();
X   cy=pompalcy;
X   cyy=pompalcyy;
X   coef=pompcoef;
X   ips=(*spc).ips;
/* numerical integration */
X   for(i=0;i<nip-1;i++) for(j=0;j<4;j++){
X      coef[i][j]=0.;
X      for(k=0;k<ndim;k++)
X         coef[i][j]+=(*spc).basis[k].beta*(*spc).basis[k].c2[i][j];
X   }
X   if((*spc).ilow==1) l1int(kints[0],ips[1],coef[0],1,what);
X   else l2int(kints[0],(*spc).low,ips[1],coef[0],what);
X   for(i=1;i<nip-2;i++)
X      m1int(kints[i],ips[i],ips[i+1],what,coef[i],0);
X   if((*spc).iupp==1) l1int(kints[nip-2],ips[nip-2],coef[nip-2],-1,what);
X   else l2int(kints[nip-2],ips[nip-2],(*spc).upp,coef[nip-2],what);
/* ctheta */
X   ctheta=0.;
X   for(i=0;i<nip-1;i++) ctheta+=kints[i][0];
X   if(ctheta>0)(*xp)=0;
X   else{
X      (*xp)=1.;
X      return 0.;
X   }
X   ctheta=mylog(ctheta);
/* logl - uncensored */
X   logl=0.;
X   for(i=0;i<(*data).ndata;i++){
X      if((*data).same[i]==0) 
X         f=pol3(coef[(*data).idata[i]],(*data).data[i])-ctheta;
X      logl+=(*data).weights[i]*f;
X   }
X   ctheta=myexp(-ctheta);
X   if(what==0){
X      return logl;
X   }
/* get ctheta-j and ctheta-jk */
X   initk(0,ndim,(*spc).score,(*spc).info,cy,cyy);
X   (void)savecoden(spc,0,nip-1,(*spc).score,(*spc).info);
X
/* score and hessian - basic */
X   for(i=0;i<ndim;i++){
X      for(j=0;j<ndim;j++){
X         (*spc).info[i][j]= (*data).dsw*(*spc).info[i][j]*ctheta;
X      }
X      (*spc).score[i]= -(*data).dsw*(*spc).score[i]*ctheta;
X   }
X   for(i=0;i<ndim;i++) for(j=0;j<ndim;j++)
X      (*spc).info[i][j]-= (*spc).score[i]*(*spc).score[j]/(*data).dsw;
X   if(what==1)for(i=0;i<ndim;i++)cuu[i]=(*spc).score[i];
/* uncensored - score */
X   for(i=0;i<ndim;i++)(*spc).score[i]+=(*spc).basis[i].sumunc;
/* symmatrize */
X   for(i=0;i<ndim;i++) for(j=i+1;j<ndim;j++)(*spc).info[i][j]=(*spc).info[j][i];
X   return logl;
}
/******************************************************************************/
void savecode1(spc,j,cz,czz,what)
int j;
struct space *spc;
double *cz,**czz,*what;
{
X   int k,j2;
X   double inp3(),mat3();
X   for(k=0;k<(*spc).ndim;k++){
X      if(j>=(*spc).basis[k].c3[0]&&j<=(*spc).basis[k].c3[1]){
X         cz[k]+=inp3(what,(*spc).basis[k].c2[j]);
X         for(j2=0;j2<=k;j2++){
X            if(j>=(*spc).basis[j2].c3[0]&&j<=(*spc).basis[j2].c3[1])
X               czz[k][j2]+=
X                  mat3(what,(*spc).basis[k].c2[j],(*spc).basis[j2].c2[j]);
X         }
X      }
X   }
X   return;
}
/******************************************************************************/
int savecoden(spc,i0,i1,cz,czz)
int i0,i1;
struct space *spc;
double *cz,**czz;
{
X   int j;
X   for(j=i0;j<i1;j++) savecode1(spc,j,cz,czz,kints[j]);
X   return i1;
}
/******************************************************************************/
void initk(i,ndim,v,mm,v2,mm2)
int ndim,i;
double *v,*v2,**mm,**mm2;
{
X   int j,k;
X   if(i==0){
X      for(j=0;j<ndim;j++){
X         for(k=0;k<ndim;k++)mm[j][k]=0.;
X         v[j]=0.;
X      }
X   }
X   else{
X      for(j=0;j<ndim;j++){
X         for(k=0;k<ndim;k++)mm[j][k]=mm2[j][k];
X         v[j]=v2[j];
X      }
X   }
X   return;
}
/******************************************************************************/
/* gets the rao statistic */
double rao(spc,data,loc,iloc)
struct space *spc;
struct datastruct *data;
double loc;
int iloc;
{
X   double **ii,*ss,r,praox(),iext[7],c2ext[4];
X   int i,j,getnewc2(),j0ext,ndim=(*spc).ndim;
X   ii=raoii;
X   ss=raoss;
X   for(i=0;i<ndim;i++){
X      for(j=0;j<ndim;j++)ii[i][j]=(*spc).info[i][j];
X      ss[i]=0.;
X   }
X   (*bbx).c2=raoc2;
X   j0ext=getnewc2(spc,data,loc,bbx,iext,c2ext);
X   ss[ndim]=praox(spc,data,bbx,ii[ndim],iext,c2ext,j0ext);
X   for(i=0;i<ndim;i++)ii[i][ndim]=ii[ndim][i];
X   r=ss[ndim];
X   i=lusolve2(ii,ndim+1,ss);
X   if(i<0)r=0.;
X   if((*data).wttype>0) return ss[ndim]*r/(*data).weights[iloc];
X   return ss[ndim]*r;
}
/******************************************************************************/
/* computes the new parts of score and hessian */
double praox(spc,data,bb,iext,intext,c2ext,j0ext)
struct space *spc;
struct datastruct *data;
struct basisfunct *bb;
double *iext,intext[7],c2ext[4];
{
X   double pol3(),inp3(),mat3(),sext;
X   int i,j,ndim=(*spc).ndim;
X   double int2ext[7],save22coden();
X   if(j0ext>=0)for(i=0;i<7;i++)int2ext[i]=kints[j0ext][i]-intext[i];
/* get ctheta-j and ctheta-jk */
X   for(j=0;j<=ndim;j++) iext[j]=0.;
X   sext=save22coden(spc,iext,bb,int2ext,j0ext,c2ext);
X   for(j=0;j<=ndim;j++) iext[j]= (*data).dsw*iext[j]*ctheta;
X   sext= -(*data).dsw*sext*ctheta;
X   for(j=0;j<ndim;j++) iext[j]-= sext*cuu[j]/(*data).dsw;
X   iext[ndim]-= sext*sext/(*data).dsw;
X   sext+=(*bb).sumunc;
X   return sext;
}
/******************************************************************************/
/* coefficients for a test-basis function */
int getnewc2(spc,data,loc,bb,intext,c2ext)
struct space *spc;
struct datastruct *data;
struct basisfunct *bb;
double loc,intext[7],c2ext[4];
{
X   int i,j,k,j0,j1,ii[3],j0ext;
X   double pol3(),coef[10],rrr[10],t[3],cc[3];
X   void l1int(),l2int(),m1int();
/* get (*spc).basis.c3 */
X   t[1]=(*spc).knots[(*spc).nk-2];
X   t[2]=(*spc).knots[(*spc).nk-1];
X   t[0]=loc;
X   for(j=0;j<3;j++){
X      ii[j]=(*spc).nip-1;
X      for(i=1;i<(*spc).nip;i++) if((*spc).ips[i]>=t[j]){
X         ii[j]=i;
X         i=(*spc).nip;
X      }
X   }
X   (*bb).c3[0]=ii[0]-1;
X   if(ii[1]<ii[0]-1)(*bb).c3[0]=ii[1];
X   (*bb).c3[1]=(*spc).nip+1;
X   if((*bb).c3[0]<0)(*bb).c3[0]=0;
X   cc[0]=1.;
X   cc[1]=(t[0]-t[2])/(t[2]-t[1]);
X   cc[2]=(t[1]-t[0])/(t[2]-t[1]);
/* get (*spc).basis.c2 */
X   for(j=0;j<(*spc).nip;j++) for(j0=0;j0<4;j0++)(*bb).c2[j][j0]=0.;
X   for(j=(*bb).c3[0];j<=(*bb).c3[1];j++) for(j1=0;j1<3;j1++) if(j>=ii[j1]){
X      (*bb).c2[j][3]+=cc[j1];
X      (*bb).c2[j][2]-=3.*cc[j1]*t[j1];
X      (*bb).c2[j][1]+=3.*cc[j1]*t[j1]*t[j1];
X      (*bb).c2[j][0]-=cc[j1]*t[j1]*t[j1]*t[j1];
X   }
/* get j0ext */
X   j0ext=(*spc).nip+100;
X   if(t[0]<(*spc).ips[1])j0ext=0;
X   else for(i=1;i<(*spc).nip-2;i++){
X      if(t[0]==(*spc).ips[i])j0ext= -1;
X      else if(t[0]<(*spc).ips[i+1])j0ext=i;
X      if(j0ext<(*spc).nip+50)i=(*spc).nip;
X   }
X   if(j0ext>(*spc).nip+50)j0ext=(*spc).nip-2;
/* get c2ext */
X   for(i=0;i<4;i++)c2ext[i]=0.;
X   for(j1=0;j1<3;j1++) if(t[j1]<=t[0]){
X      c2ext[3]+=cc[j1];
X      c2ext[2]-=3.*cc[j1]*t[j1];
X      c2ext[1]+=3.*cc[j1]*t[j1]*t[j1];
X      c2ext[0]-=cc[j1]*t[j1]*t[j1]*t[j1];
X   }
/* get intext */
X   if(j0ext>=0){
X      for(j=0;j<4;j++){
X         coef[j]=0.;
X         for(k=0;k<(*spc).ndim;k++)
X            coef[j]+=(*spc).basis[k].beta*(*spc).basis[k].c2[j0ext][j];
X      }
X      if(j0ext==0){
X         if((*spc).ilow==1) l1int(rrr,t[0],coef,1,1);
X         else l2int(rrr,(*spc).low,t[0],coef,1);
X      }
X      if(j0ext==((*spc).nip)-2) l2int(rrr,(*spc).ips[j0ext],t[0],coef,1);
X      if(j0ext>0&&j0ext<((*spc).nip)-2)
X         m1int(rrr,(*spc).ips[j0ext],t[0],1,coef,0);
X     for(i=0;i<7;i++)intext[i]=rrr[i];
X   }
/* get (*spc).basis.sumunrc */
X   (*bb).sumunc=0.;
X   for(j1=0;j1<(*data).ndata;j1++){
X      j0=(*data).idata[j1];
X      if(j0>=(*bb).c3[0]&&j0<=(*bb).c3[1]){
X        if(j0!=j0ext || t[0]>(*data).data[j1])
X           (*bb).sumunc+=(*data).weights[j1]*pol3((*bb).c2[j0],(*data).data[j1]);
X        else (*bb).sumunc+=(*data).weights[j1]*pol3(c2ext,(*data).data[j1]);
X      }
X   }
X   return j0ext;
}
/******************************************************************************/
/* integrates all steps for score and hessian */
double save22coden(spc,czz,bb,int2ext,j0ext,c2ext)
int j0ext;
struct basisfunct *bb;
struct space *spc;
double *czz,int2ext[7],c2ext[4];
{
X   int j,k,i1=((*spc).nip)-1;
X   double inp3(),mat3(),cz=0;
/* correct the new one */
X   if(j0ext>=0 && j0ext<i1)for(k=0;k<7;k++)kints[j0ext][k]-=int2ext[k];
/* regular stuff */
X   for(j=0;j<i1;j++){
X      if(j>=(*bb).c3[0]&&j<=(*bb).c3[1]){
X         for(k=0;k<(*spc).ndim;k++)
X            if(j>=(*spc).basis[k].c3[0]&&j<=(*spc).basis[k].c3[1])
X               czz[k]+=mat3(kints[j],(*spc).basis[k].c2[j],(*bb).c2[j]);
X         cz+=inp3(kints[j],(*bb).c2[j]);
X         czz[(*spc).ndim]+=mat3(kints[j],(*bb).c2[j],(*bb).c2[j]);
X      }
X   }
/* correct the new one  part II */
X   if(j0ext>=0 && j0ext<i1){
X      for(k=0;k<(*spc).ndim;k++)
X         if(j0ext>=(*spc).basis[k].c3[0]&&j0ext<=(*spc).basis[k].c3[1])
X            if(j0ext>=(*bb).c3[0]&&j0ext<=(*bb).c3[1])
X            czz[k]+=mat3(int2ext,(*spc).basis[k].c2[j0ext],c2ext);
X      cz+=inp3(int2ext,c2ext);
X      czz[(*spc).ndim]+=mat3(int2ext,c2ext,c2ext);
X      for(k=0;k<7;k++)kints[j0ext][k]+=int2ext[k];
X   }
X   return cz;
}
/******************************************************************************/
double remdim(spc,data,spc2,silent,wxx)
struct space *spc,*spc2;
struct datastruct *data;
int silent;
double *wxx;
X
/* spc   - model to be worked on
X   spc2  - temporary copy of the space
X   data  - data
X   silent- should info be printed? (1=yes) */
X
{
X   double ratmax=0.,se,phi;
X   int i,j,k,irmax=1,ndim=(*spc).ndim;
X   double wbest=0,lbest=0;int l;
X   void betarem(),setupspace(),swapspace();
X
/* ratmax - largest phi/se ratio
X   phi    - coefficient in power basis
X   se     - standard errors
X   i,j,k  - counters
X   irmax  - for which coefficient is ratmax attained
X   getip  - gets the
X   setupspace() - sets up a new space
X   swapspace()  - copies a space
X   betarem()    - new starting values */
X
/* invert the Hessian */
X   luinverse((*spc).info,ndim);
/* copy for later use */
X   swapspace(spc2,spc);
X
X   for(i=0;i<(*spc).nk;i++){
/* compute the coefficient */
X      phi = 0.;
X      for(j=0;j<ndim;j++) phi+=(*spc).basis[j].beta*(*spc).basis[j].c1[i+2];
X      phi=fabs(phi);
/* the standard error */
X      se = 0.;
X      for(j=0;j<ndim;j++)for(k=0;k<ndim;k++)
X         se-=(*spc).basis[j].c1[i+2]*(*spc).basis[k].c1[i+2]*(*spc).info[j][k];
X      l=(*spc).iknots[i];
X      se = sqrt(fabs(se));
X      if((*data).wttype>0) se=se*sqrt((*data).weights[l]);
/* Select for which knot se/phi takes it maximal value */
X      if(se > phi * ratmax){
X         ratmax = se / phi;
X         irmax = i;
X         wbest=1./(ratmax*ratmax);
X         lbest=wbest*(*data).weights[l]/2;
X         (*wxx)=(*data).weights[l];
X      }
X   }
X   if(silent==1) (void)printf("rem(%.2f), wald=%.2f pll=%.2f ",
X                       (*spc).knots[irmax],wbest,lbest);
X
/* get (*spc).nk and (spc).ndim */
X   (*spc).nk -= 1;
X   (*spc).ndim -= 1;
/* remove the knot */
X   for(i=irmax;i<(*spc).nk;i++){
X       (*spc).iknots[i]=(*spc).iknots[i+1];
X       (*spc).knots[i]=(*spc).knots[i+1];
X   }
/* get (*spc).ips (*spc).nip (*data).idatx and (*spc).basis.iks */
/* get (*spc).basis.c1  (*spc).basis.c2 (*spc).basis.c3 (*spc).basis.sumunc */
X   setupspace(spc,data);
/* get (*spc).basis.beta */
X   betarem(spc2,spc,irmax);
X   return 1./(ratmax*ratmax);
}
/******************************************************************************/
void betarem(spc2,spc,irmax)
struct space *spc,*spc2;
int irmax;
{
X   int i,j,k;
X   double **mm2,*r1,x,y;
X   void solver(),redo1(),redo2();
X   k=(*spc2).ndim;
X   mm2=betaremm2;
X   r1=betaremr1;
X
/* find A, the restriction */
X   for(i=0;i<k;i++)mm2[0][i]=(*spc2).basis[i].c1[irmax+2];
/* solve the quadratic problem */
X   solver(mm2,k,1,r1,spc2);
/* problems, a coefficient needing to be <0 is >=0 */
X   if(((*spc).ilow==1 && r1[0]>=0) ||( (*spc).iupp==1 && r1[1]>=0)){ 
/* only restrictions on the lower tail */
X      if((*spc).ilow==1 && (*spc).iupp==0){ 
X         if(irmax<=2){ 
X            for(i=0;i<((*spc2).nk)+2;i++){
X               r1[i] = 0.;
X               for(j=0;j<(*spc2).ndim;j++)
X                  r1[i]+=(*spc2).basis[j].beta*(*spc2).basis[j].c1[i];
X            }
X            redo1(spc2,irmax,k);
X            for(i=0;i<k+3;i++)for(j=0;j<k+3;j++)mm2[i][j]=0.;
X            for(j=0;j<k+3;j++)for(i=0;i<k;i++)mm2[i][j]=(*spc2).basis[i].c1[j];
X            for(i=k+1;i<k+3;i++)mm2[i][i]=1.;
X            mm2[k][0]=1.;
X            luinverse(mm2,k+3);
X            for(i=0;i<k;i++){
X               (*spc2).basis[i].beta=0.;
X               for(j=0;j<k+3;j++)(*spc2).basis[i].beta+=mm2[j][i]*r1[j];
X            }
X         }
X         for(i=0;i<k;i++)mm2[0][i]=(*spc2).basis[i].c1[irmax+2];
X         for(i=0;i<k;i++)mm2[1][i]=0.;
X         mm2[1][0]=1.;
X         x=(*spc2).basis[0].beta;
X         (*spc2).basis[0].beta=0.;
X         solver(mm2,k,2,r1,spc2);
X         r1[0]+=x;
X      } 
/* only restrictions on the upper tail */
X      if((*spc).ilow==0 && (*spc).iupp==1){  
X         if(irmax>=k-2){
X            for(i=0;i<((*spc2).nk)+2;i++){
X               r1[i] = 0.;
X               for(j=0;j<(*spc2).ndim;j++)
X                  r1[i]+=(*spc2).basis[j].beta*(*spc2).basis[j].c1[i];
X            }
X            redo2(spc2,irmax,(*spc2).nk); 
X            for(i=0;i<k+3;i++)for(j=0;j<k+3;j++)mm2[i][j]=0.;
X            for(j=0;j<k+3;j++)for(i=0;i<k;i++)mm2[i][j]=(*spc2).basis[i].c1[j];
X            for(i=k+1;i<k+3;i++)mm2[i][i]=1.;
X            mm2[k][0]=1.;
X            luinverse(mm2,k+3);
X            for(i=0;i<k;i++){
X               (*spc2).basis[i].beta=0.;
X               for(j=0;j<k+3;j++)(*spc2).basis[i].beta+=mm2[j][i]*r1[j];
X            }
X            for(i=0;i<((*spc2).nk)+2;i++){
X               r1[i] = 0.;
X               for(j=0;j<(*spc2).ndim;j++)
X                  r1[i]+=(*spc2).basis[j].beta*(*spc2).basis[j].c1[i];
X            }
X         }
X         for(i=0;i<k;i++)mm2[0][i]=(*spc2).basis[i].c1[irmax+2];
X         for(i=0;i<k;i++)mm2[1][i]=0.;
X         mm2[1][1]=1.;
X         x=(*spc2).basis[1].beta;
X         (*spc2).basis[1].beta=0.;
X         solver(mm2,k,2,r1,spc2);
X         r1[1]+=x;
X      } 
/* restrictions on both tails */
X      if((*spc).ilow==1 && (*spc).iupp==1){ 
X         if(irmax>k-3 || irmax<=2){
X            for(i=0;i<((*spc2).nk)+2;i++){
X               r1[i] = 0.;
X               for(j=0;j<(*spc2).ndim;j++)
X                  r1[i]+=(*spc2).basis[j].beta*(*spc2).basis[j].c1[i];
X            }
X            if(irmax<=2)redo1(spc2,irmax,k);
X            if(irmax>k-3) redo2(spc2,irmax,(*spc2).nk);
X            for(i=0;i<k+3;i++)for(j=0;j<k+3;j++)mm2[i][j]=0.;
X            for(j=0;j<k+3;j++)for(i=0;i<k;i++)mm2[i][j]=(*spc2).basis[i].c1[j];
X            for(i=k+1;i<k+3;i++)mm2[i][i]=1.;
X            mm2[k][0]=1.;
X            luinverse(mm2,k+3);
X            for(i=0;i<k;i++){
X               (*spc2).basis[i].beta=0.;
X               for(j=0;j<k+3;j++)(*spc2).basis[i].beta+=mm2[j][i]*r1[j];
X            }
X         }
X         for(i=0;i<k;i++)mm2[0][i]=(*spc2).basis[i].c1[irmax+2];
X         for(i=0;i<k;i++)mm2[1][i]=0.;
X         mm2[1][0]=1.;
X         for(i=0;i<k;i++)mm2[2][i]=0.;
X         mm2[2][1]=1.;
X         x=(*spc2).basis[0].beta;
X         (*spc2).basis[0].beta=0.;
X         y=(*spc2).basis[1].beta;
X         (*spc2).basis[1].beta=0.;
X         solver(mm2,k,3,r1,spc2);
X         r1[0]+=x;
X         r1[1]+=y;
X      }  
X   } 
/* record beta */
X   for(i=0;i<(*spc2).ndim;i++)(*spc2).basis[i].beta=r1[i];
/* get it to be a power basis */
X   for(i=0;i<((*spc2).nk)+2;i++){
X      r1[i] = 0.;
X      for(j=0;j<(*spc2).ndim;j++)
X         r1[i]+=(*spc2).basis[j].beta*(*spc2).basis[j].c1[i];
X   } 
X   for(i=irmax;i<(*spc2).nk;i++) r1[i+2]=r1[i+3];
X   k=k-1;
X   for(i=0;i<k+3;i++)for(j=0;j<k+3;j++)mm2[i][j]=0.;
X   for(j=0;j<k+3;j++)for(i=0;i<k;i++)mm2[i][j]=(*spc).basis[i].c1[j];
X   for(i=k+1;i<k+3;i++)mm2[i][i]=1.;
X   mm2[k][0]=1.;
X   luinverse(mm2,k+3);
X   for(i=0;i<k;i++){
X      (*spc).basis[i].beta=0.;
X      for(j=0;j<k+3;j++)(*spc).basis[i].beta+=mm2[j][i]*r1[j];
X   }
X   return;
}
/******************************************************************************/
void redo1(spc,irmax,k)
int k,irmax;
struct space *spc;
{
X   int i,i0=0,i1=2,i2=3;
X   double a,b,*t,*c;
X   c=(*spc).basis[0].c1;
X   t=(*spc).knots;
X   if(irmax==2)i1=1;
X   if(irmax==0)i0=1;
X   for(i=0;i<=k+1;i++)c[i]=0.;
X   a=t[i2]-t[i0];
X   b=t[i2]-t[i1];
X   c[i0+2]=1.;
X   c[i1+2]= -a/b;
X   c[i2+2]= -1.-c[i1+2];
X   c[1]= -3.*(t[i0]*t[i0]+c[i1+2]*t[i1]*t[i1]+c[i2+2]*t[i2]*t[i2]);
X   c[0]= -t[i2]*c[1]-c[i0+2]*a*a*a-c[i1+2]*b*b*b;
}
/******************************************************************************/
void redo2(spc,irmax,k)
int k,irmax;
struct space *spc;
{
X   int i,i0=k-4,i1=k-3,i2=k-1;
X   double *t,*c;
X   c=(*spc).basis[1].c1;
X   t=(*spc).knots;
X   if(irmax==k-3)i1=k-2;
X   if(irmax==k-1)i2=k-2;
X   for(i=0;i<=k+1;i++)c[i]=0.;
X   c[i0+2]=1.;
X   c[i1+2]=(t[i0]-t[i2])/(t[i2]-t[i1]);
X   c[i2+2]= -1.-c[i1+2];
}
/******************************************************************************/
/* This routine computes  Inv(H)%*%t(A)%*%Inv(A%*%Inv(H)%*%t(A))%*%A          
X   mm2 = A jxi, mm1 = Inv (H) ixi */
void solver (mm2,i,j,r1,spc)
double **mm2,*r1;
int i,j;
struct space *spc;
/* r1   - new beta
X   i    - (*spc).ndim
X   j    - number of restrictions */
{
X   int k0,k1,k2;
X   double **c1,**c2,**c3,**mm1;
X   c1=solc1;
X   c2=solc2;
X   c3=solc3;
/* k0,k1,k2 - counters
X   c1,c2,c3 - half-products, as indicated below
X   mm1      - hessian */
X   if(i==j)for(k0=0;k0<i;k0++)r1[k0]=0.;
X   mm1=(*spc).info;
/* c1=Inv(H)%*%t(A) ixi * ixj = ixj */
X   for(k0=0;k0<i;k0++)for(k1=0;k1<j;k1++){
X      c1[k0][k1]=0.;
X      for(k2=0;k2<i;k2++)c1[k0][k1]+=mm1[k0][k2]*mm2[k1][k2];
X   }
/* c2=A%*%c1 jxi * ixj = jxj */
X   for(k0=0;k0<j;k0++)for(k1=0;k1<j;k1++){
X      c2[k0][k1]=0.;
X      for(k2=0;k2<i;k2++)c2[k0][k1]+=mm2[k0][k2]*c1[k2][k1];
X   }
/* c2=Inv(c2)=Inv(A%*%Inv(H)%*%t(A)) */
X   luinverse(c2,j);
/* c3=c1%*%c2=Inv(H)%*%t(A)%*%Inv(A%*%Inv(H)%*%t(A)) ixj * jxj = ixj */
X   for(k0=0;k0<i;k0++)for(k1=0;k1<j;k1++){
X      c3[k0][k1]=0.;
X      for(k2=0;k2<j;k2++)c3[k0][k1]+=c1[k0][k2]*c2[k2][k1];
X   }
/* c1=c3%*%A=Inv(H)%*%t(A)%*%Inv(A%*%Inv(H)%*%t(A))%*%A ixj * jxi = ixi */
X   for(k0=0;k0<i;k0++)for(k1=0;k1<i;k1++){
X      c1[k0][k1]=0.;
X      for(k2=0;k2<j;k2++)c1[k0][k1]+=c3[k0][k2]*mm2[k2][k1];
X   }
/* shift beta */
X   for(k0=0;k0<i;k0++){
X      r1[k0]=(*spc).basis[k0].beta;
X      for(k1=0;k1<i;k1++)r1[k0]-=c1[k0][k1]*(*spc).basis[k1].beta;
X   }
X   return;
}
/******************************************************************************/
/* get the c2, c3 and sumunc elements of a bsis function */
void getc2(spc,data)
struct space *spc;
struct datastruct *data;
X
/* spc  - model
X   data - data */
{
X   int i,j,k,l,m,n;
X   double pol3(),a,b;
X   struct basisfunct *bn;
/* no deep meanings */
/* get (*spc).basis.c3 */
X   for(i=0;i<(*spc).ndim;i++){
X      bn=&((*spc).basis[i]);
X      (*bn).c3[0]=(*bn).iks[0]-1;
X      (*bn).c3[1]=(*bn).iks[2]+1;
X      if(i==0) (*bn).c3[0]=0;
X      if(i==1 || i==2) (*bn).c3[1]=(*spc).nip;
X      if(i>2) (*bn).c3[1]=(*bn).iks[4]+1;
/* get (*spc).basis.c2 */
X      for(j=0;j<(*spc).nip;j++) for(k=0;k<4;k++)(*bn).c2[j][k]=0.;
X      for(j=(*bn).c3[0];j<=(*bn).c3[1];j++){
X         l=5;
X         if(i==0||i==1)l=3;
X         if(i==2)l=4;
X         if(i==0){
X            (*bn).c2[j][0]+=(*bn).c1[0];
X            (*bn).c2[j][1]+=(*bn).c1[1];
X         }
X         for(n=0;n<l;n++){
X            if(i==0)m=n;
X            if(i==1)m=(*spc).nk-3+n;
X            if(i==2)m=(*spc).nk-4+n;
X            if(i>2)m=i+n-3;
X            a=(*spc).knots[m];
X            b=(*bn).c1[2+m];
X            if(j>=(*bn).iks[n]){
X               (*bn).c2[j][3]+=b;
X               b=a*b;
X               (*bn).c2[j][2]-=3.*b;
X               b=a*b;
X               (*bn).c2[j][1]+=3.*b;
X               b=a*b;
X               (*bn).c2[j][0]-=b;
X            }
X         }
X      }
/* get (*spc).basis.sumunc */
X      (*bn).sumunc=0.;
X      for(m=0;m<(*data).ndata;m++){
X         l=(*data).idata[m];
X         if(l>=(*bn).c3[0]&&l<=(*bn).c3[1])
X           (*bn).sumunc+=(*data).weights[m]*pol3((*bn).c2[l],(*data).data[m]);
X      } 
X   }
X   return;
}
/******************************************************************************/
/* get c1 - the power basis representation - for a basisfunction */
void getc1(t,c,i,k)
double *c,*t;
int i,k;
{
/* get (*spc).basis.c1 */
X   double a,b,d[10],r;
X   void getonec1();
X   int j;
X   for(j=0;j<=k+1;j++)c[j]=0.;
X   if(i==0){
X      a=t[2]-t[0];
X      b=t[2]-t[1];
X      c[2]= 1.;
X      c[3]= -a/b;
X      c[4]= -c[2]-c[3];
X      c[1]= -3.*(t[0]*t[0]+c[3]*t[1]*t[1]+c[4]*t[2]*t[2]);
X      c[0]= -t[2]*c[1]-c[2]*a*a*a-c[3]*b*b*b;
X   }
X   if(i==1){
X      c[k-1]=1.;
X      c[k]=(t[k-3]-t[k-1])/(t[k-1]-t[k-2]);
X      c[k+1]= -c[k]-c[k-1];
X   }
X   if(i==2) getonec1(c,k-2,t,k-4);
X   if(i>2){
X      getonec1(c,i-1,t,i-3);
X      getonec1(d,0,t,i-2);
X      a=0.;
X      b=0.;
X      for(j=0;j<4;j++){
X         r=(t[k-1]-t[i+j-2]);
X         a+=d[j]*r*r*r;
X         r=(t[k-1]-t[i+j-3]);
X         b+=c[i+j-1]*r*r*r;
X      }
X      for(j=0;j<4;j++)c[i+j]-=d[j]*b/a;
X   } 
}
/******************************************************************************/
void getonec1(c,i,t,j)
double *c,*t;
int i,j;
{
X   c[i]=1.;
X   c[i+3]=(t[j+2]-t[j])*(t[j]-t[j+1])/((t[j+2]-t[j+3])*(t[j+1]-t[j+3]));
X   c[i+2]=(c[i+3]*(t[j+1]-t[j+3])+t[j+1]-t[j])/(t[j+2]-t[j+1]);
X   c[i+1]=-1.-c[i+3]-c[i+2];
X   return;
}
X   
/******************************************************************************/
void setupspace(spc,data)
struct space *spc;
struct datastruct *data;
{
X   int i;
X   void getip();
/* get (*spc).ips (*spc).nip (*data).idatx and (*spc).basis.iks */
X   getip(spc,data);
/* get (*spc).basis.c1 */
X   for(i=0;i<(*spc).ndim;i++)getc1((*spc).knots,(*spc).basis[i].c1,i,(*spc).nk);
/* get (*spc).basis.c2 (*spc).basis.c3 and (*spc).basis.sumunc */
X   getc2(spc,data);
X   return;
}
/******************************************************************************/
int startspace(spc,data,strt,silent)
int silent;
struct space *spc;
struct datastruct *data;
{
X   int i,k,l,ok,rearrange();
X   double r,s;
X   void setupspace(),startnow();
/* place the knots */
X   k=(*data).ndata;
X   ok=1;
X   if(strt==0){
X      (*spc).iknots[0]=0;
X      (*spc).iknots[1]=(int)(k/2);
X      (*spc).iknots[2]=k-1;
X      for(i=0;i<3;i++)
X         (*spc).knots[i]=(*data).data[(*spc).iknots[i]];
X      (*spc).nk=3;
X      if(silent==1)(void)printf("Starting knots at %.2f, %.2f and %.2f ",
X               (*spc).knots[0],(*spc).knots[1],(*spc).knots[2]);
X      (*spc).ndim=2;
X      if(silent==1)(void)fflush(stdout);
X   }
X   if(strt<0){
X      if(strt== -1){
X         l=((*spc).nk);
X         r=l+2.;
X         if(l>3){
X         s=(double)l/((double)l-3.);
X         for(i=1;i<l-1;i++)  (*spc).iknots[i]=(((i-1)*s+1.)/r)*k;
X         }
X         else (*spc).iknots[1]=(k-1)/2.;
X         (*spc).iknots[0]=0;
X         (*spc).iknots[l-1]=k-1;
X      }
X      if(strt== -2){
X         l=((*spc).nk);
X         r=l+4.;
X         if(l>3){
X         s=(double)(l+2)/((double)l-3.); 
X         for(i=1;i<l-1;i++)  (*spc).iknots[i]=(((i-1)*s+1.)/r)*k;
X         }
X         else (*spc).iknots[1]=(k-1)/2.;
X         (*spc).iknots[0]=1;
X         (*spc).iknots[l-1]=k-2;
X      }
X      if(strt == -3){
X         l=(*spc).nk-2;
X         r=l+2.;
X         s=0;
X         if(l>3) s=(double)(l)/((double)l-3.);
X         for(i=1;i<l-1;i++)  (*spc).iknots[i+1]=(((i-1)*s+1.)/r)*(k-8.)+4;
X         (*spc).iknots[1]=4;
X         (*spc).iknots[l]=k-5;
X         (*spc).iknots[0]=0;
X         (*spc).iknots[l+1]=k-1;
X         l=l+2;
X      }
X      for(i=0;i<l;i++){
X         (*spc).knots[i]=(*data).data[(*spc).iknots[i]];
X         if(i>0)if((*spc).knots[i]<=(*spc).knots[i-1])ok=0;
X      }
X      (*spc).nk=l;
X      if(ok==0)ok=rearrange(spc,data);
X      if(ok==0)return ok;
X      (*spc).ndim=l-1;
X      if(silent==1){
X         (void)printf("\nRestart: knots at ");
X         for(i=0;i<l;i++)(void)printf("%.2f ",(*spc).knots[i]);
X         (void)fflush(stdout);
X      }
X   }
X   if((*spc).ilow==1) (*spc).low=(*data).data[0];
X   if((*spc).iupp==1) (*spc).upp=(*data).data[k-1];
/* compute c1, c2, c3 and so on */
X   setupspace(spc,data);
/* starting values */
X   startnow(spc,data);
X   return ok;
}
/******************************************************************************/
void startnow(spc,data)
struct space *spc;
struct datastruct *data;
{
X   int i,j0,j1;
X   double r0,r1,s0,s1;
X   for(i=0;i<(*spc).ndim;i++)(*spc).basis[i].beta=0.;
X   r0=((*spc).knots[0]+(*spc).knots[1])/2.;
X   r1=((*spc).knots[(*spc).nk-2]+(*spc).knots[(*spc).nk-1])/2.;
X   j0=0;
X   j1=0;
X   s0=0.;
X   s1=0.;
X   for(i=0;i<(*data).ndata;i++){
X      if((*data).data[i]<r0){
X         s0+= r0-(*data).data[i];
X         j0+=2;
X      }
X      if((*data).data[i]>r1){
X         s1+= (*data).data[i]-r1;
X         j1+=2;
X      }
X   }
X   s0=2.*s0/(double)j0;
X   s1=2.*s1/(double)j1;
X   if((*spc).ilow==1) (*spc).basis[0].beta= -1./fabs(s0*(*spc).basis[0].c1[1]);
X   if((*spc).iupp==1)(*spc).basis[1].beta=
X      -1./fabs(s1*(*spc).basis[1].c2[(*spc).nip][1]);
X   return;
}
/******************************************************************************/
int rearrange(spc,data)
struct space *spc;
struct datastruct *data;
{
X   int i,k,*ix,jx[500],nk=(*spc).nk,is,j,l;
X   double *sorted;
X   sorted=rearsorted;
X   ix=rearix;
X   for(i=0;i<(*data).ndata;i++){
X      sorted[i]=(*data).data[i];
X      ix[i]=i;
X   }
X   k=1;
X   for(i=1;i<(*data).ndata;i++){
X      if(sorted[i]>sorted[k-1]){
X         sorted[k]=sorted[i];
X         ix[k]=ix[i];
X         k++;
X      }
X   }
X   is=0;
X   for(i=0;i<nk;i++){
X      for(j=is;j<k;j++){
X         if((*spc).knots[i]<=sorted[j]){
X            jx[i]=j;
X            is=j;
X            j=k;
X         }
X      }
X   }
X   for(is=0;is<10;is++){
X      for(j=1;j<nk-1;j++) if(jx[j]==jx[j-1]) if(jx[j+1]>jx[j])jx[j]++;
X      for(j=nk-2;j>0;j--) if(jx[j]==jx[j+1]) if(jx[j-1]<jx[j])jx[j]--;
X   }
X   l=1;
X   for(j=1;j<nk;j++)if(jx[j]==jx[j-1])l=0;
X   if(l==0)return 0;
X   for(i=0;i<nk;i++){
X      (*spc).iknots[i]=ix[jx[i]];
X      (*spc).knots[i]=sorted[jx[i]];
X   }
X   return 1;
}
/******************************************************************************/
/* selects integration points */
void getip(spc,data)
struct space *spc;
struct datastruct *data;
X
/* spc - model
X   data- data */
{
X   int i,nip=1,*iips,j,k,*kips;
X   double *ips;
X
/* i,j,k - counters
X   kips  - integration points for knots
X   iips  - index of integration points
X   ips   - integration points
X   nip   - number of integration points */
X
/* allocation */
X   iips=getiips;
X   ips=(*spc).ips;
X
/* the first two points */
X   iips[0]= -4;
X   iips[1]=(*spc).iknots[0];
X   ips[0]=(*spc).low;
X   ips[1]=(*spc).knots[0];
/* get (*spc).ips - the actual integration points */
X   for(i=1;i<(*spc).nk;i++){
X      j=(*spc).iknots[i]-(*spc).iknots[i-1];
X      j=floor((double)(j/100.))+1;
/* integration points in between knots */
X      for(k=1;k<j;k++){
X         nip++;
X         iips[nip]=(*spc).iknots[i-1]+k*((*spc).iknots[i]-(*spc).iknots[i-1])/j;
X         ips[nip]=(*data).data[iips[nip]];
X      }
/* every knot is an integration point */
X      nip++;
X      ips[nip]=(*spc).knots[i];
X      iips[nip]=(*spc).iknots[i];
X   }
/* the last integration point */
X   nip++;
X   ips[nip]=(*spc).upp;
X   iips[nip]=(*data).ndata+4;
/* get (*spc).nip */
X   nip++;
X   (*spc).nip=nip;
/* get (*data).idatx */
X   for(i=0;i<nip-1;i++) for(j=iips[i];j<=iips[i+1]-1;j++){
X      if(j>=0 && j<(*data).ndata) (*data).idata[j]=i;
X   }
/* get (*spc).basis.iks */
X   kips=iips;
X   for(i=0;i<(*spc).nk;i++){
X      j=(*spc).iknots[i];
X      kips[i]=(*data).idata[j];
X   }
X   for(i=0;i<3;i++)(*spc).basis[0].iks[i]=kips[i];
X   for(i=0;i<3;i++)(*spc).basis[1].iks[i]=kips[i-3+(*spc).nk];
X   if((*spc).ndim>2)for(i=0;i<4;i++)(*spc).basis[2].iks[i]=kips[i-4+(*spc).nk];
X   for(j=3;j<(*spc).ndim;j++)for(i=0;i<5;i++)(*spc).basis[j].iks[i]=kips[j+i-3];
X   return;
}
void rpqlsd(coef,knots,bnd,ipq,pq,lk,lp)
double *coef,*knots,*pq,*bnd;
int *ipq,*lk,*lp;
{
X   double *kpl,**cpl,*ppl,*dsvector(),**dsmatrix(),*pqx,r;
X   double z1int(),z2int(),z3int(),*zz,cor;
X   int i,j,nk,fst,lst;
X   void getp0(),getp1(),getp2(),getq0(),getq1(),getq2();
X   /* Gaussian quadrature coefficients */
X   ww6[1 ]= 0.467913934572691; yy6[1 ]= 0.238619186083197;
X   ww6[2 ]= 0.360761573048139; yy6[2 ]= 0.661209386466265;
X   ww6[3 ]= 0.171324429379170; yy6[3 ]= 0.932469514203152;
X   ww7[1 ]=  0.00178328072169643; yy7[1 ]=  0.99930504173577217;
X   ww7[2 ]=  0.00414703326056247; yy7[2 ]=  0.99634011677195533;
X   ww7[3 ]=  0.00650445796897836; yy7[3 ]=  0.99101337147674429;
X   ww7[4 ]=  0.00884675982636395; yy7[4 ]=  0.98333625388462598;
X   ww7[5 ]=  0.01116813946013113; yy7[5 ]=  0.97332682778991098;
X   ww7[6 ]=  0.01346304789671864; yy7[6 ]=  0.96100879965205377;
X   ww7[7 ]=  0.01572603047602472; yy7[7 ]=  0.94641137485840277;
X   ww7[8 ]=  0.01795171577569734; yy7[8 ]=  0.92956917213193957;
X   ww7[9 ]=  0.02013482315353021; yy7[9 ]=  0.91052213707850282;
X   ww7[10]=  0.02227017380838325; yy7[10]=  0.88931544599511414;
X   ww7[11]=  0.02435270256871087; yy7[11]=  0.86599939815409277;
X   ww7[12]=  0.02637746971505466; yy7[12]=  0.84062929625258032;
X   ww7[13]=  0.02833967261425948; yy7[13]=  0.81326531512279754;
X   ww7[14]=  0.03023465707240248; yy7[14]=  0.78397235894334139;
X   ww7[15]=  0.03205792835485155; yy7[15]=  0.75281990726053194;
X   ww7[16]=  0.03380516183714161; yy7[16]=  0.71988185017161088;
X   ww7[17]=  0.03547221325688239; yy7[17]=  0.68523631305423327;
X   ww7[18]=  0.03705512854024005; yy7[18]=  0.64896547125465731;
X   ww7[19]=  0.03855015317861563; yy7[19]=  0.61115535517239328;
X   ww7[20]=  0.03995374113272034; yy7[20]=  0.57189564620263400;
X   ww7[21]=  0.04126256324262353; yy7[21]=  0.53127946401989457;
X   ww7[22]=  0.04247351512365359; yy7[22]=  0.48940314570705296;
X   ww7[23]=  0.04358372452932345; yy7[23]=  0.44636601725346409;
X   ww7[24]=  0.04459055816375657; yy7[24]=  0.40227015796399163;
X   ww7[25]=  0.04549162792741814; yy7[25]=  0.35722015833766813;
X   ww7[26]=  0.04628479658131442; yy7[26]=  0.31132287199021097;
X   ww7[27]=  0.04696818281621002; yy7[27]=  0.26468716220876742;
X   ww7[28]=  0.04754016571483031; yy7[28]=  0.21742364374000708;
X   ww7[29]=  0.04799938859645831; yy7[29]=  0.16964442042399283;
X   ww7[30]=  0.04834476223480295; yy7[30]=  0.12146281929612056;
X   ww7[31]=  0.04857546744150343; yy7[31]=  0.07299312178779904;
X   ww7[32]=  0.04869095700913972; yy7[32]=  0.02435029266342443;
X
/* allocation */
X   kpl=dsvector((*lk)*4);
X   ppl=dsvector((*lk)*4);
X   cpl=dsmatrix((*lk)*4,4);
X   pqx=dsvector((*lp));
/* get the integration points: the knots */
X   nk=(*lk)+1;
X   for(i=0;i<=nk;i++){
X      if(i<nk && i>0)kpl[i]=knots[i-1];
X      if(i==0){
X         kpl[0]=knots[0]-1;
X         if(bnd[0]>0.5)kpl[0]=bnd[1];
X         else if(pq[0]<kpl[0])kpl[0]=pq[0]-1.;
X      }
X      if(i==nk){
X         kpl[nk]=knots[nk-2]+1;
X         if(bnd[2]>0.5)kpl[nk]=bnd[3];
X         else if(pq[(*lp)-1]+1.>kpl[nk])kpl[nk]=pq[(*lp)-1]+1.;
X      }
/* get the coeffiecients */
X      cpl[i][0]=coef[0];
X      cpl[i][1]=coef[1];
X      cpl[i][2]=0.;
X      cpl[i][3]=0.;
X      for(j=0;i>j && j<(*lk);j++){
X         cpl[i][3]+=coef[j+2];
X         cpl[i][2]-=3.*coef[j+2]*knots[j];
X         cpl[i][1]+=3.*coef[j+2]*knots[j]*knots[j];
X         cpl[i][0]-=coef[j+2]*knots[j]*knots[j]*knots[j];
X      }
X      if(i>=nk-1){
X         cpl[i][3]=0.;
X         cpl[i][2]=0.;
X      }
X   }
/* compute the density */
X   ppl[0]=0.;
X   if(bnd[0]>0.5)ppl[1]=z2int(kpl[0],kpl[1],cpl[0]);
X   else ppl[1]=z1int(kpl[1],cpl[0],1);
X   for(i=1;i<nk-1;i++)ppl[i+1]=z3int(kpl[i],kpl[i+1],cpl[i],0)+ppl[i];
X   if(bnd[2]>0.5) ppl[nk]=z2int(kpl[nk-1],kpl[nk],cpl[nk-1])+ppl[nk-1];
X   else ppl[nk]=z1int(kpl[nk-1],cpl[nk-1],-1)+ppl[nk-1];
/* higher precision needed */
X   if(ppl[nk]<0.99999 || ppl[nk]>1.00001){
/* integration points: knots times four */
X      nk=4*(*lk)-2;
X      for(i=0;i<=nk;i++){
X         if(i<nk && i>0){
X            j=floor((double)(i-1)/4.);
X            r=(double)i/4.-j-0.25;
X            kpl[i]=(1.-r)*knots[j]+r*knots[j+1];
X         }
X         if(i==0){
X            kpl[0]=knots[0]-1;
X            if(bnd[0]>0.5)kpl[0]=bnd[1];
X            else if(pq[0]<kpl[0])kpl[0]=pq[0]-1.;
X         }
X         if(i==nk){
X            kpl[nk]=knots[(*lk)-1]+1;
X            if(bnd[2]>0.5)kpl[nk]=bnd[3];
X            else if(pq[(*lp)-1]+1>kpl[nk])kpl[nk]=pq[(*lp)-1]+1.;
X         }
/* get the coeffiecients */
X         cpl[i][0]=coef[0];
X         cpl[i][1]=coef[1];
X         cpl[i][2]=0.;
X         cpl[i][3]=0.;
X         for(j=0;i>j*4 && j<(*lk);j++){
X            cpl[i][3]+=coef[j+2];
X            cpl[i][2]-=3.*coef[j+2]*knots[j];
X            cpl[i][1]+=3.*coef[j+2]*knots[j]*knots[j];
X            cpl[i][0]-=coef[j+2]*knots[j]*knots[j]*knots[j];
X         }
X         if(i>=nk-1){
X            cpl[i][3]=0.;
X            cpl[i][2]=0.;
X         }
X      }
/* compute the density */
X      if(bnd[0]>0.5)ppl[1]=z2int(kpl[0],kpl[1],cpl[0]);
X      else ppl[1]=z1int(kpl[1],cpl[0],1);
X      for(i=1;i<nk-1;i++) ppl[i+1]=z3int(kpl[i],kpl[i+1],cpl[i],0)+ppl[i];
X      if(bnd[2]>0.5) ppl[nk]=z2int(kpl[nk-1],kpl[nk],cpl[nk-1])+ppl[nk-1];
X      else ppl[nk]=z1int(kpl[nk-1],cpl[nk-1],-1)+ppl[nk-1];
X   }
/* correction factor */
X   cor=ppl[nk];
/* correct the density */
X   for(i=0;i<=nk;i++)ppl[i]=ppl[i]/ppl[nk];
X   j=0;
/* initialize */
X   if((*ipq)==0)zz=ppl;  
X   else zz=kpl;
/* before the first point */
X   for(j=0;j<(*lp) && pq[j]<=zz[0];j++){
X      if((*ipq)==0){
X         if(bnd[0]>0.5)pqx[j]=kpl[0];
X         else pqx[j]= -1.0e100;
X      }
X      else pqx[j]=0.;
X   }
/* before the first knot */
X   fst=j;
X   lst=j-1;
X   for(j=j;j<(*lp) && pq[j]<=zz[1];j++) lst=j;
X   if(lst>=fst){
X      if((*ipq)==0) getq0(pq,pqx,fst,lst,cpl[0],kpl[0],bnd[0],cor);
X      else getp0(pq,pqx,fst,lst,cpl[0],kpl[0],bnd[0],cor);
X   }
/* per interval between integration points */
X   for(i=1;i<nk-1;i++){
X      fst=j;
X      lst=j-1;
X      for(j=j;j<(*lp) && pq[j]<=zz[i+1];j++) lst=j;
X      if(lst>=fst){
X         if((*ipq)==0)
X            getq1(pq,pqx,fst,lst,cpl[i],kpl[i],kpl[i+1],ppl[i],ppl[i+1]);
X         else getp1(pq,pqx,fst,lst,cpl[i],kpl[i],kpl[i+1],ppl[i],ppl[i+1]);
X      }
X   }
/* beyond the larst knot */
X   fst=j;
X   lst=j-1;
X   for(j=j;j<(*lp) && pq[j]<zz[nk];j++) lst=j;
X   if(lst>=fst){
X      if((*ipq)==0) getq2(pq,pqx,fst,lst,cpl[nk-1],kpl[nk],bnd[2],cor);
X      else getp2(pq,pqx,fst,lst,cpl[nk-1],kpl[nk],bnd[2],cor);
X   }
/* outside the range */
X   for(j=j;j<(*lp);j++){
X      if((*ipq)==0){
X         if(bnd[2]>0.5)pqx[j]=kpl[nk];
X         else pqx[j]= 1.0e100;
X      }
X      else pqx[j]=1.;
X   }
X  for(j=0;j<(*lp);j++)pq[j]=pqx[j];
}
/******************************************************************************/
void getp0(q,p,f,l,cf,k,b,cr)
double *q,*p,*cf,k,b,cr;
int f,l;
{
X   int i; 
X   double z2int(),z1int();
X   if(b>0.5) for(i=f;i<=l;i++) p[i]=z2int(k,q[i],cf)/cr;
X   else for(i=f;i<=l;i++) p[i]=z1int(q[i],cf,1)/cr;
}
/******************************************************************************/
void getq0(p,q,f,l,cf,k,b,cr)
double *q,*p,*cf,k,b,cr;
int f,l;
{
X   int i; 
X   double pqexpi();
X   if(b>0.5)for(i=f;i<=l;i++) q[i]=pqexpi(2,k,p[i]/cr,cf);
X   else for(i=f;i<=l;i++) q[i]=pqexpi(1,k,p[i]/cr,cf);
}
/******************************************************************************/
void getp2(q,p,f,l,cf,k,b,cr)
double *q,*p,*cf,k,b,cr;
int f,l;
{
X   int i; 
X   double z2int(),z1int();
X   if(b>0.5) for(i=f;i<=l;i++) p[i]=1.-z2int(q[i],k,cf)/cr;
X   else for(i=f;i<=l;i++) p[i]=1.-z1int(q[i],cf,-1)/cr;
}
/******************************************************************************/
void getq2(p,q,f,l,cf,k,b,cr)
double *q,*p,*cf,k,b,cr;
int f,l;
{
X   int i; 
X   double pqexpi();
X   if(b>0.5)for(i=f;i<=l;i++) q[i]=pqexpi(4,k,1.-p[i]/cr,cf);
X   else for(i=f;i<=l;i++) q[i]=pqexpi(3,k,1.-p[i]/cr,cf);
}
/******************************************************************************/
void getp1(q,p,f,l,cf,k0,k1,p0,p1)
double *p,*q,*cf,k0,k1,p0,p1;
int f,l;
{
X   int i,j=0;
X   double z3int(),r;
X   if(l-f>5)j=1;
X   p[f]=z3int(k0,q[f],cf,j);
X   r=p[f]+z3int(q[l],k1,cf,j);
X   for(i=f+1;i<=l;i++) p[i]=z3int(q[i-1],q[i],cf,j)+p[i-1];
X   r=p[l]+z3int(q[l],k1,cf,j);
X   r=(p1-p0)/r;
X   for(i=f;i<=l;i++)p[i]=p0+p[i]*r;
}
/******************************************************************************/
void getq1(p,q,f,l,cf,k0,k1,p0,p1)
double *p,*q,*cf,k0,k1,p0,p1;
int f,l;
{
X   int i,j;
X   double y[51],f1[101],r,s,getf();
X   r=(k1-k0)/100.;
X   for(i=0;i<101;i++) f1[i]=getf(cf,(double)(k0+r*i));
X   y[0]=0.;
X   for(i=1;i<=50;i++)y[i]=y[i-1]+r*(f1[2*(i-1)]+4*f1[2*i-1]+f1[2*i])/3.;
X   s=(p1-p0)/y[50];
X   for(i=0;i<=50;i++)y[i]=p0+y[i]*s;
X   i=0;
X   s=2.*r;
X   for(j=f;j<=l;j++){
X      q[j]=k0-1.;
X      do{
X         if(p[j]>=y[i] && p[j]<=y[i+1]) q[j]=k0+s*i+s*(p[j]-y[i])/(y[i+1]-y[i]);
X         else i++;
X      }
X      while (q[j]<k0);
X   }
}
/******************************************************************************/
/* computes integrals from -inf to t1 (j=1) or from t1 to inf (j1= -1)
X   of exp(polynomial(c0)) */
X
double z1int(t1,c0,j)
double t1,*c0;
int j;
{
X   double f1,mylog(),myexp();
X   if(c0[1]<0) j = -j;
X   f1 = mylog(fabs(1./c0[1])) + c0[1]*t1+c0[0];
X   if(f1>600.) f1=600.;
X   return (double)(j*myexp(f1));
}
/******************************************************************************/
/* computes integrals from t1 to t2 of exp(polynomial(c0)) */
X
double z2int(t1,t2,c0)
double t1,t2,*c0;
{
X   int i1=1;
X   double f1,f2,mylog(),myexp();
X   if(t2==t1)return 0.;
X   if(c0[1]!=0){
X      if(c0[1]<0) i1 = -1;
X      f1 = mylog(fabs(1./c0[1])) + c0[1]*t1+c0[0];
X      f2 = f1 + c0[1]*(t2-t1);
X      if(f1>600.) f1=600.;
X      if(f2>600.) f2=600.;
X      return (double)(i1*myexp(f2)-i1*myexp(f1));
X   }
X   else return (t2-t1)*myexp(c0[0]);
}
/******************************************************************************/
/* computes integrals from t1 to t2 of exp(polynomial(coef)) */
X
double z3int(k1,k2,coef,accuracy)
double k1,k2,*coef;
int accuracy;
{
X   double r1,r2,x,y,v,vv=0.,myexp();
X   int i1;
X   if(k2==k1)return 0.;
X   r1 = ((k2 - k1) / 2);
X   r2 = ((k2 + k1) / 2);
X   if(accuracy==1){
X      for(i1=1;i1<4;i1++){
X         y=yy6[i1]*r1;
X         v=r1*ww6[i1];
X         x=r2-y;
X         vv+=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X         x=r2+y;
X         vv+=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X      }
X   }
X   else{
X      for(i1=1;i1<33;i1++){
X         y=yy7[i1]*r1;
X         v=r1*ww7[i1];
X         x=r2-y;
X         vv+=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X         x=r2+y;
X         vv+=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X      }
X   }
X   return vv;
}
/******************************************************************************/
/* 1: -inf -> x / 2: t -> x / 3: x -> inf / 4: x -> t  */
double pqexpi(version,t,p,cf)
int version;
double t,p,*cf;
{
X   double mylog(),myexp();
X   if(cf[1]!=0. || version == 1 || version == 3){
X      p=p*cf[1];
X      if(version == 1 && p < 0)return myexp((double)600.);
X      if(version == 3 && p > 0)return -myexp((double)600.);
X      if(version==2 || version ==4)t=myexp(t*cf[1]+cf[0]);
X      if(version == 2 && t+p < 0)return myexp((double)600.);
X      if(version == 4 && t-p < 0)return -myexp((double)600.);
X      if(version==1)return (mylog(p)-cf[0])/cf[1];
X      if(version==2)return (mylog(t+p)-cf[0])/cf[1];
X      if(version==3)return (mylog(-p)-cf[0])/cf[1];
X      return (mylog(t-p)-cf[0])/cf[1];
X   }
X   if(version==2)return t+p/myexp(cf[0]);
X   return t-p/myexp(cf[0]);
}
/******************************************************************************/
double *dsvector(l)
int l;
/* allocate a double vector with subscript range v[0...l] */
{
X   double *v;
X   int i;
X   char *S_alloc();
X   v=(double *)S_alloc((unsigned)(l+1),sizeof(double));
X   for(i=0;i<=l;i++)v[i]=0.;
X   return v;
}
/******************************************************************************/
double getf(c,x)
double *c,x;
{
X   return exp(c[0]+x*(c[1]+x*(c[2]+x*c[3])));
}
/******************************************************************************/
double mylog(x)
double x;
{
if(x < 10.e-250)return (double)(-575.64627);
else return log(x);
}
/******************************************************************************/
double myexp(x)
double x;
{
if(x > 576.)return exp((double)576.);
else return exp(x);
}
/******************************************************************************/
/* allocate an short vector with subscript range v[0...l] */
short *issvector(l)
int l;
{
X   int i;
X   short *v;
X   v=(short *)S_alloc((unsigned)(l+1),sizeof(short));
X   for(i=0;i<=l;i++)v[i]=0;
X   return v;
}
/******************************************************************************/
/* allocate an int vector with subscript range v[0...l] */
int *isvector(l)
int l;
{
X   int *v,i;
X   v=(int *)S_alloc((unsigned)(l+1),sizeof(int));
X   for(i=0;i<=l;i++)v[i]=0;
X   return v;
}
/******************************************************************************/
/* computes integrals from t1 to t2  (numerically)
X   of x^i (i<1, what=0; i<7 o.w.) times exp(polynomial(coef)) */
void m1int(vv,k1,k2,what,coef,accuracy)
X
/* accuracy  - accuracy
X   r1 and r2 - from (k1,k2) to (-1,1)         */
X
double k1,k2,*vv,*coef;
int accuracy,what;
{
X   double r1,r2,x,y,z,v,myexp();
X   int i1,i2,j;
X   r1 = ((k2 - k1) / 2);
X   r2 = ((k2 + k1) / 2);
X   for(i1=0;i1<7;i1++)vv[i1]=0.;
X   if(k2==k1)return;
X   j=7;
X   if(what==0)j=1;
X   if(accuracy==1){
X      for(i1=1;i1<4;i1++){
X         y=yy6[i1]*r1;
X         v=r1*ww6[i1];
X         x=r2-y;
X         z=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X         vv[0]+=z;
X         for(i2=1;i2<j;i2++){
X            z=z*x;
X            vv[i2]+=z;
X         }
X         x=r2+y;
X         z=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X         vv[0]+=z;
X         for(i2=1;i2<j;i2++){
X            z=z*x;
X            vv[i2]+=z;
X         }
X      }
X   }
X   else{
X      for(i1=1;i1<33;i1++){
X         y=yy7[i1]*r1;
X         v=r1*ww7[i1];
X         x=r2-y;
X         z=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X         vv[0]+=z;
X         for(i2=1;i2<j;i2++){
X            z=z*x;
X            vv[i2]+=z;
X         }
X         x=r2+y;
X         z=v*myexp(coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3])));
X         vv[0]+=z;
X         for(i2=1;i2<j;i2++){
X            z=z*x;
X            vv[i2]+=z;
X         }
X      }
X   }
}
/******************************************************************************/
/* computes integrals from -inf to t1 (j=1) or from t1 to inf (j1= -1)
X   of x^i (i<1, what=0; i<7 o.w.) times exp(polynomial(coef)) */
X
void l1int(results,t1,coef,j,what)
int j,what;
double t1,*coef,*results;
X
{
X   double u6,u5,u4,u3,u2,u1,u0,f1,b1=coef[1],b0=coef[0],b3[7],fctf1();
X   b3[0]=(double)1./b1;
X
X   f1 = b3[0];
X   results[0]=fctf1(b0,b1,t1,f1,j);
X   if(what==0)return;
X   b3[1]=b3[0]*b3[0];
X   b3[2]=b3[1]*b3[0];
X   b3[3]=b3[2]*b3[0];
X   b3[4]=b3[3]*b3[0];
X   b3[5]=b3[4]*b3[0];
X   b3[6]=b3[5]*b3[0];
X
X   u1 = b3[0];
X   u0 = -b3[1];
X   f1 = u1*t1+u0;
X   results[1]=fctf1(b0,b1,t1,f1,j);
X
X   u2 = b3[0];
X   u1 = -2*b3[1];
X   u0 = 2*b3[2];
X   f1 = (u2*t1+u1)*t1+u0;
X   results[2]=fctf1(b0,b1,t1,f1,j);
X
X   u3 = b3[0];
X   u2 = -3*b3[1];
X   u1 = 6*b3[2];
X   u0 = -6*b3[3];
X   f1 = ((u3*t1+u2)*t1+u1)*t1+u0;
X   results[3]=fctf1(b0,b1,t1,f1,j);
X
X   u4 = b3[0];
X   u3 = -4*b3[1];
X   u2 = 12*b3[2];
X   u1 = -24*b3[3];
X   u0 = 24*b3[4];
X   f1 = (((u4*t1+u3)*t1+u2)*t1+u1)*t1+u0;
X   results[4]=fctf1(b0,b1,t1,f1,j);
X
X   u5 = b3[0];
X   u4 = -5*b3[1];
X   u3 = 20*b3[2];
X   u2 = -60*b3[3];
X   u1 = 120*b3[4];
X   u0 = -120*b3[5];
X   f1 = ((((u5*t1+u4)*t1+u3)*t1+u2)*t1+u1)*t1+u0;
X   results[5]=fctf1(b0,b1,t1,f1,j);
X
X   u6 = b3[0];
X   u5 = -6*b3[1];
X   u4 = 30*b3[2];
X   u3 = -120*b3[3];
X   u2 = 360*b3[4];
X   u1 = -720*b3[5];
X   u0 = 720*b3[6];
X   f1 = (((((u6*t1+u5)*t1+u4)*t1+u3)*t1+u2)*t1+u1)*t1+u0;
X   results[6]=fctf1(b0,b1,t1,f1,j);
X
X   return;
}
/******************************************************************************/
/* computes integrals from t1 to t2  (exactly)
X   of x^i (i<1, what=0; i<7 o.w.) times exp(polynomial(coef)) */
void l2int(results,t1,t2,coef,what)
int what;
double t1,t2,*coef,*results;
X
{
X   double u6,u5,u4,u3,u2,u1,u0,f1,b1=coef[1],b0=coef[0],b3[7],f2,fctf2();
X   double myexp();
X   int i;
X   if(b1!=0.){
X      b3[0]=(double)1./b1;
X 
X      f1 = b3[0];
X      f2 = b3[0];
X      results[0]=fctf2(b0,b1,t1,t2,f1,f2);
X      if(what==0)return;
X      b3[1]=b3[0]*b3[0];
X      b3[2]=b3[1]*b3[0];
X      b3[3]=b3[2]*b3[0];
X      b3[4]=b3[3]*b3[0];
X      b3[5]=b3[4]*b3[0];
X      b3[6]=b3[5]*b3[0];
X
X      u1 = b3[0];
X      u0 = -b3[1];
X      f1 = u1*t1+u0;
X      f2 = u1*t2+u0;
X      results[1]=fctf2(b0,b1,t1,t2,f1,f2);
X
X      u2 = b3[0];
X      u1 = -2*b3[1];
X      u0 = 2*b3[2];
X      f1 = (u2*t1+u1)*t1+u0;
X      f2 = (u2*t2+u1)*t2+u0;
X      results[2]=fctf2(b0,b1,t1,t2,f1,f2);
X
X      u3 = b3[0];
X      u2 = -3*b3[1];
X      u1 = 6*b3[2];
X      u0 = -6*b3[3];
X      f1 = ((u3*t1+u2)*t1+u1)*t1+u0;
X      f2 = ((u3*t2+u2)*t2+u1)*t2+u0;
X      results[3]=fctf2(b0,b1,t1,t2,f1,f2);
X
X      u4 = b3[0];
X      u3 = -4*b3[1];
X      u2 = 12*b3[2];
X      u1 = -24*b3[3];
X      u0 = 24*b3[4];
X      f1 = (((u4*t1+u3)*t1+u2)*t1+u1)*t1+u0;
X      f2 = (((u4*t2+u3)*t2+u2)*t2+u1)*t2+u0;
X      results[4]=fctf2(b0,b1,t1,t2,f1,f2);
X
X      u5 = b3[0];
X      u4 = -5*b3[1];
X      u3 = 20*b3[2];
X      u2 = -60*b3[3];
X      u1 = 120*b3[4];
X      u0 = -120*b3[5];
X      f1 = ((((u5*t1+u4)*t1+u3)*t1+u2)*t1+u1)*t1+u0;
X      f2 = ((((u5*t2+u4)*t2+u3)*t2+u2)*t2+u1)*t2+u0;
X      results[5]=fctf2(b0,b1,t1,t2,f1,f2);
X
X      u6 = b3[0];
X      u5 = -6*b3[1];
X      u4 = 30*b3[2];
X      u3 = -120*b3[3];
X      u2 = 360*b3[4];
X      u1 = -720*b3[5];
X      u0 = 720*b3[6];
X      f1 = (((((u6*t1+u5)*t1+u4)*t1+u3)*t1+u2)*t1+u1)*t1+u0;
X      f2 = (((((u6*t2+u5)*t2+u4)*t2+u3)*t2+u2)*t2+u1)*t2+u0;
X      results[6]=fctf2(b0,b1,t1,t2,f1,f2);
X      return;
X   }
X   b0=myexp(b0);
X   results[0]=(t2-t1)*b0;
X   if(what==0)return;
X   f2=t2;
X   f1=t1;
X   for(i=1;i<7;i++){
X     f2=f2*t2;
X     f1=f1*t1;
X     results[i]=b0*(f2-f1)/(double)(i+1);
X   }
X   return;
}
/******************************************************************************/double fctf1(b0,b1,t1,f1,j)
double b0,b1,t1,f1;
int j;
{
X   double mylog(),myexp();
X   if(f1<0) j = -j;
X   f1 = mylog(fabs(f1)) + b1*t1+b0;
X   if(f1>600.) f1=600.;
X   return (double)(j*myexp(f1));
}
/******************************************************************************/
double fctf2(b0,b1,t1,t2,f1,f2)
double b0,b1,t1,t2,f1,f2;
{
X   int i1=1,i2=1;
X   double mylog(),myexp();
X   if(f1<0) i1 = -1;
X   f1 = mylog(fabs(f1)) + b1*t1+b0;
X   if(f1>600.) f1=600.;
X   if(f2<0) i2 = -1;
X   f2 = mylog(fabs(f2)) + b1*t2+b0;
X   if(f2>600.) f2=600.;
X   return (double)(i2*myexp(f2)-i1*myexp(f1));
}
/******************************************************************************/
double pol3(coef,x)
double *coef,x;
{
X   return coef[0]+x*(coef[1]+x*(coef[2]+x*coef[3]));
}
/******************************************************************************/
double inp3(c1,c2)
double *c1,*c2;
{
X   return c1[0]*c2[0]+c1[1]*c2[1]+c1[2]*c2[2]+c1[3]*c2[3];
}
/******************************************************************************/
double mat3(c1,c2,c3)
double *c1,*c2,*c3;
{
X   double x=0.;
X   int i,j;
X   for(i=0;i<4;i++)for(j=0;j<4;j++)x+=c1[i+j]*c2[i]*c3[j];
X   return x;
}
/******************************************************************************/
/* copies one space into another space */
void swapspace(s1,s2)
struct space *s1,*s2;
{
X   int i,j,k;
X   (*s1).ndim=(*s2).ndim;
X   (*s1).nk=(*s2).nk;
X   (*s1).cth=(*s2).cth;
X   (*s1).nip=(*s2).nip;
X   (*s1).aic=(*s2).aic;
X   (*s1).low=(*s2).low;
X   (*s1).upp=(*s2).upp;
X   (*s1).ilow=(*s2).ilow;
X   (*s1).iupp=(*s2).iupp;
X   for(i=0;i<(*s1).nip;i++) (*s1).ips[i]=(*s2).ips[i];
X   for(i=0;i<(*s1).nk;i++){
X      (*s1).knots[i]=(*s2).knots[i];
X      (*s1).iknots[i]=(*s2).iknots[i];
X   }
X   for(i=0;i<(*s1).ndim;i++){
X      for(j=0;j<5;j++)(*s1).basis[i].iks[j]=(*s2).basis[i].iks[j];
X      (*s1).score[i]=(*s2).score[i];
X      for(j=0;j<(*s1).ndim;j++) (*s1).info[i][j]=(*s2).info[i][j];
X      (*s1).basis[i].beta=(*s2).basis[i].beta;
X      for(j=0;j<2;j++)(*s1).basis[i].c3[j]=(*s2).basis[i].c3[j];
X      (*s1).basis[i].sumunc=(*s2).basis[i].sumunc;
X      for(j=0;j<(*s1).nk+2;j++)(*s1).basis[i].c1[j]=(*s2).basis[i].c1[j];
X      for(j=0;j<4;j++)for(k=0;k<(*s1).nip;k++)
X         (*s1).basis[i].c2[k][j]=(*s2).basis[i].c2[k][j];
X   }
X   return;
}
/******************************************************************************/
void quadalloc()
{
/* Gaussian quadrature coefficients */
X   ww6[1 ]= 0.467913934572691; yy6[1 ]= 0.238619186083197;
X   ww6[2 ]= 0.360761573048139; yy6[2 ]= 0.661209386466265;
X   ww6[3 ]= 0.171324429379170; yy6[3 ]= 0.932469514203152;
X   ww7[1 ]=  0.00178328072169643; yy7[1 ]=  0.99930504173577217;
X   ww7[2 ]=  0.00414703326056247; yy7[2 ]=  0.99634011677195533;
X   ww7[3 ]=  0.00650445796897836; yy7[3 ]=  0.99101337147674429;
X   ww7[4 ]=  0.00884675982636395; yy7[4 ]=  0.98333625388462598;
X   ww7[5 ]=  0.01116813946013113; yy7[5 ]=  0.97332682778991098;
X   ww7[6 ]=  0.01346304789671864; yy7[6 ]=  0.96100879965205377;
X   ww7[7 ]=  0.01572603047602472; yy7[7 ]=  0.94641137485840277;
X   ww7[8 ]=  0.01795171577569734; yy7[8 ]=  0.92956917213193957;
X   ww7[9 ]=  0.02013482315353021; yy7[9 ]=  0.91052213707850282;
X   ww7[10]=  0.02227017380838325; yy7[10]=  0.88931544599511414;
X   ww7[11]=  0.02435270256871087; yy7[11]=  0.86599939815409277;
X   ww7[12]=  0.02637746971505466; yy7[12]=  0.84062929625258032;
X   ww7[13]=  0.02833967261425948; yy7[13]=  0.81326531512279754;
X   ww7[14]=  0.03023465707240248; yy7[14]=  0.78397235894334139;
X   ww7[15]=  0.03205792835485155; yy7[15]=  0.75281990726053194;
X   ww7[16]=  0.03380516183714161; yy7[16]=  0.71988185017161088;
X   ww7[17]=  0.03547221325688239; yy7[17]=  0.68523631305423327;
X   ww7[18]=  0.03705512854024005; yy7[18]=  0.64896547125465731;
X   ww7[19]=  0.03855015317861563; yy7[19]=  0.61115535517239328;
X   ww7[20]=  0.03995374113272034; yy7[20]=  0.57189564620263400;
X   ww7[21]=  0.04126256324262353; yy7[21]=  0.53127946401989457;
X   ww7[22]=  0.04247351512365359; yy7[22]=  0.48940314570705296;
X   ww7[23]=  0.04358372452932345; yy7[23]=  0.44636601725346409;
X   ww7[24]=  0.04459055816375657; yy7[24]=  0.40227015796399163;
X   ww7[25]=  0.04549162792741814; yy7[25]=  0.35722015833766813;
X   ww7[26]=  0.04628479658131442; yy7[26]=  0.31132287199021097;
X   ww7[27]=  0.04696818281621002; yy7[27]=  0.26468716220876742;
X   ww7[28]=  0.04754016571483031; yy7[28]=  0.21742364374000708;
X   ww7[29]=  0.04799938859645831; yy7[29]=  0.16964442042399283;
X   ww7[30]=  0.04834476223480295; yy7[30]=  0.12146281929612056;
X   ww7[31]=  0.04857546744150343; yy7[31]=  0.07299312178779904;
X   ww7[32]=  0.04869095700913972; yy7[32]=  0.02435029266342443;
}
/******************************************************************************/
double **dsmatrix(r,c)
int r,c;
/* allocate a double matrix with subscript range m[0..r][0..c] */
{
X   int i;
X   double **m,*dsvector();
X   m=(double **) S_alloc((unsigned)(r+1),sizeof(double*));
X   for(i=0;i<=r;i++) m[i]=dsvector(c);
X   return m;
}
SHAR_EOF
  $shar_touch -am 1023032197 'nlsd/ylsd.c' &&
  chmod 0655 'nlsd/ylsd.c' ||
  $echo 'restore of' 'nlsd/ylsd.c' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/ylsd.c:' 'MD5 check failed'
50881089ada478f62d7911c5e2624be7  nlsd/ylsd.c
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/ylsd.c'`"
    test 84088 -eq "$shar_count" ||
    $echo 'nlsd/ylsd.c:' 'original size' '84088,' 'current size' "$shar_count!"
  fi
fi
# ============= nlsd/ylsd.S ==============
if test -f 'nlsd/ylsd.S' && test "$first_param" != -c; then
  $echo 'x -' SKIPPING 'nlsd/ylsd.S' '(file already exists)'
else
  $echo 'x -' extracting 'nlsd/ylsd.S' '(text)'
  sed 's/^X//' << 'SHAR_EOF' > 'nlsd/ylsd.S' &&
dnlsd <- function(x, fit)
{
X	y <- fit$coef.pol[1] + x * fit$coef.pol[2]
X	for(i in 1:length(fit$knots))
X		y <- y + fit$coef.kts[i] * ((abs(x - fit$knots[i]) + x - fit$
X			knots[i])/2)^3
X	y <- exp(y)
X	if(fit$bound[1] > 0)
X		y[x < fit$bound[2]] <- 0
X	if(fit$bound[3] > 0)
X		y[x > fit$bound[4]] <- 0
X	y
}
nlsd.plot <- function(fit, n = 100, what = "d", add = F, xlim, xlab, ylab, type, ...)
{
X	if(add) {
X		plim <- (par()$usr)[1:2]
X		u4 <- plim[1]
X		u3 <- plim[2]
X		if(!missing(xlim)) {
X			u4 <- max(xlim[1], plim[1])
X			u3 <- min(xlim[2], plim[2])
X		}
X	}
X	else {
X		if(missing(xlim)) {
X			u1 <- qnlsd(0.01, fit)
X			u2 <- qnlsd(0.99, fit)
X			u3 <- 1.1 * u1 - 0.1 * u2
X			u4 <- 1.1 * u2 - 0.1 * u1
X		}
X		else {
X			u3 <- xlim[1]
X			u4 <- xlim[2]
X		}
X	}
X	xx <- (0:(n - 1))/(n - 1) * (u4 - u3) + u3
X	if(what == "d" || what == "D")
X		yy <- dnlsd(xx, fit)
X	if(what == "f" || what == "F" || what == "p" || what == "P")
X		yy <- pnlsd(xx, fit)
X	if(what == "s" || what == "S")
X		yy <- 1 - pnlsd(xx, fit)
X	if(what == "h" || what == "H")
X		yy <- dnlsd(xx, fit)/(1 - pnlsd(xx, fit))
X	if(missing(xlab))
X		xlab <- ""
X	if(missing(ylab))
X		ylab <- ""
X	if(missing(type))
X		type <- "l"
X	if(add)
X		lines(xx, yy, type=type,  ...)
X	else plot(xx, yy, xlab = xlab, ylab = ylab, type = type, ...)
X	invisible()
}
nlsd.summary <- function(fit)
{
X	ul <- fit$penalty
X	um <- fit$samples[1]
X	if(length(fit$samples) > 1)
X		um <- fit$samples[1] + fit$samples[4]
X	else um <- fit$samples
X	kk <- fit$logl[fit$logl[, 2] != 0, 1]
X	ad <- fit$logl[fit$logl[, 2] != 0, 2]
X	ll <- fit$logl[fit$logl[, 2] != 0, 3]
X	bb <- -2 * ll + ul * (kk - 1)
X	cc1 <- bb
X	cc2 <- bb
X	cc2[1] <- Inf
X	cc1[length(bb)] <- 0
X	if(length(bb) > 1) {
X		for(i in 1:(length(bb) - 1)) {
X			cc1[i] <- max((ll[(i + 1):(length(bb))] - ll[i])/(kk[(i +
X				1):(length(bb))] - kk[i]))
X			cc2[i + 1] <- min((ll[1:i] - ll[i + 1])/(kk[1:i] - kk[i +
X				1]))
X		}
X	}
X	c3 <- cc2 - cc1
X	cc1[c3 < 0] <- NA
X	cc2[c3 < 0] <- NA
X	uu <- cbind(kk, ad, ll, bb, 2 * cc1, 2 * cc2)
X	ww <- rep("", length(bb))
X	dimnames(uu) <- list(ww, c("knots", "A(1)/D(2)", "loglik", "AIC", 
X		"minimum penalty", "maximum penalty"))
X	print(round(uu, 2))
X	cat(paste("the present optimal number of knots is ", kk[bb == min(bb)], 
X		"\n"))
X	if(ul == log(um))
X		cat(paste("penalty(AIC) was the default: BIC=log(samplesize): log(",
X			um, ")=", round(ul, 2), "\n"))
X	else cat(paste("penalty(AIC) was ", round(ul, 2), 
X			", the default (BIC) ", "would have been", round(log(um
X			), 2), "\n"))
X	invisible()
}
pnlsd <- function(q, fit)
{
X	sq <- rank(q)
X	q <- sort(q)
X	z <- .C("rpqlsd",
X		as.double(c(fit$coef.pol, fit$coef.kts)),
X		as.double(fit$knots),
X		as.double(fit$bound),
X		as.integer(1),
X		pp = as.double(q),
X		as.integer(length(fit$knots)),
X		as.integer(length(q)))
X	zz <- z$pp[sq]
X	if(fit$bound[1] > 0)
X		zz[q < fit$bound[2]] <- 0
X	if(fit$bound[3] > 0)
X		zz[q > fit$bound[4]] <- 1
X	zz
}
qnlsd <- function(p, fit)
{
X	sp <- rank(p)
X	p <- sort(p)
X	z <- .C("rpqlsd",
X		as.double(c(fit$coef.pol, fit$coef.kts)),
X		as.double(fit$knots),
X		as.double(fit$bound),
X		as.integer(0),
X		qq = as.double(p),
X		as.integer(length(fit$knots)),
X		as.integer(length(p)))
X	zz <- z$qq[sp]
X	zz[p < 0] <- NA
X	zz[p > 1] <- NA
X	zz
}
rnlsd <- function(n, fit)
{
X	pp <- runif(n)
X	qnlsd(pp, fit)
}
nlsd.fit <- function(data, wgt, lbound, ubound, maxknots = 0, knots, nknots = 0,
X	penalty = -1, silent = T, mind = -1, wttype = 1)
{
X	call <- match.call()
X	u2 <- length(data)
X	data <- data[!is.na(data)]
X	nsample <- length(data)
X	if(nsample < 10)
X		stop("not enough data")
X	if(u2 != nsample)
X		print(paste("***", u2 - nsample, " NAs ignored in data"))
X	data <- sort(data)
X	if(missing(wgt)) wgt <- rep(1, nsample)	
X	# data can not be beyond the boundaries of the density
X	if(!missing(lbound))
X		if(data[1] < lbound)
X			stop("data below lbound")
X	if(!missing(ubound))
X		if(data[nsample] > ubound)
X			stop("data above ubound")
X	mm <- range(data)
X	if(!missing(lbound))
X		mm <- range(c(mm, lbound))
X	if(!missing(ubound))
X		mm <- range(c(mm, ubound))	# boundaries
X	ilow <- (!missing(lbound)) * 1
X	iupp <- (!missing(ubound)) * 1
X	low <- 0
X	upp <- 0
X	if(ilow == 1)
X		low <- lbound
X	if(iupp == 1) upp <- ubound	# get the maximal dimension
X	intpars <- c(-100, rep(0, 9))
#
X	if(!is.loaded(symbol.C("nlogcensor")))
X		dyn.load("nlsd.o")
#
X	z <- .C("nlogcensor",
X		z = as.integer(intpars))
X	maxp <- z$z[1]	# organize knots
X	kts <- vector(mode = "double", length = max(maxp))
X	if(maxknots > maxp - 5)
X		warning(paste("maxknots reduced to", maxp))
X	nknots <-  - nknots
X	if(!missing(knots)) {
X		nknots <- length(knots)
X		knots <- sort(knots)
X		if(!missing(lbound))
X			if(min(knots) < lbound)
X				stop("data (knots) below lbound")
X		if(!missing(ubound))
X			if(max(knots) > ubound)
X				stop("data (knots) above ubound")
X		if(nknots < 3)
X			stop("need at least three starting knots")
X		if(nknots > maxp - 5)
X			stop(paste("at most", maxp - 5, "knots possible"))
X		kts[1:nknots] <- knots
X	}
X	silent <- silent == F	# group parameters
X	intpars <- c(nsample, maxknots, nknots, silent, 1 - ilow, 1 - iupp, 
X		mind, wttype)
X	dpars <- c(penalty, low, upp)
X	data <- c(data, rep(0, maxp))	# do it
X	z <- .C("nlogcensor",
X		ip = as.integer(intpars),
X		coef = as.double(data),
X		as.double(wgt),
X		dp = as.double(dpars),
X		logl = as.double(rep(0, maxp)),
X		ad = as.integer(rep(0, maxp)),
X		kts = as.double(kts))	# error messages
X	if(z$ip[1] != 0 && z$ip[1] < 100) {
X		if(z$ip[1] == 17)
X			stop("too many knots beyond data")
X		if(z$ip[1] == 18)
X			stop("too many knots before data")
X		if(z$ip[1] == 39)
X			stop("too much data close together")
X		if(z$ip[1] == 40)
X			stop("no model could be fitted")
X		if(z$ip[1] == 2)
X			stop("error while solving system")
X		if(z$ip[1] == 8)
X			stop("too much step-halving")
X		if(z$ip[1] == 5)
X			stop("too much step-halving")
X		if(z$ip[1] == 7)
X			stop("numerical problems, likely tail related. Try lbound/ubound"
X				)
X		if(z$ip[1] == 1)
X			stop("no convergence")
X		i <- 0
X		if(missing(knots))
X			i <- 1
X		if(z$ip[1] == 3 && i == 1)
X			stop("right tail extremely heavy, try running with ubound"
X				)
X		if(z$ip[1] == 4 && i == 1)
X			stop("left tail extremely heavy, try running with lbound"
X				)
X		if(z$ip[1] == 6 && i == 1)
X			stop("both tails extremely heavy, try running with lbound and ubound"
X				)
X		if(z$ip[1] == 3 && i == 0)
X			stop("right tail too heavy or not enough knots in right tail"
X				)
X		if(z$ip[1] == 4 && i == 0)
X			stop("left tail too heavy or not enough knots in left tail"
X				)
X		if(z$ip[1] == 6 && i == 0)
X			stop("both tails too heavy or not enough knots in both tail"
X				)
X	}
X	if(z$ip[1] > 100) {
X		warning(" Not all models could be fitted")
X	}
# organize logl
X	logl <- cbind(z$ad, z$logl)
X	logl <- cbind(2 + (1:z$ip[3]), logl[1 + (1:z$ip[3]),  ])
X	kk <- (1:length(logl[, 1]))
X	kk <- kk[logl[, 2] == 0]
X	if(length(kk) > 0) logl <- logl[ - kk,  ]	# bye bye
X	list(call = call, nknots = z$ip[2], coef.pol = z$coef[1:2], coef.kts = 
X		z$coef[2 + (1:z$ip[2])], knots = z$kts[1:z$ip[2]], maxknots = z$
X		ip[3] + 2, penalty = z$dp[1], bound = c(ilow, low, iupp, upp), 
X		samples = nsample, logl = logl, range = mm, mind = z$ip[7], 
X		wttype = wttype)
}
SHAR_EOF
  $shar_touch -am 1222065098 'nlsd/ylsd.S' &&
  chmod 0655 'nlsd/ylsd.S' ||
  $echo 'restore of' 'nlsd/ylsd.S' 'failed'
  if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' ) >/dev/null 2>&1 \
  && ( md5sum --version 2>&1 | grep -v 'textutils 1.12' ) >/dev/null; then
    md5sum -c << SHAR_EOF >/dev/null 2>&1 \
    || $echo 'nlsd/ylsd.S:' 'MD5 check failed'
c23f1a0f354e3656dd08305c74e8ccd6  nlsd/ylsd.S
SHAR_EOF
  else
    shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c < 'nlsd/ylsd.S'`"
    test 7130 -eq "$shar_count" ||
    $echo 'nlsd/ylsd.S:' 'original size' '7130,' 'current size' "$shar_count!"
  fi
fi
rm -fr _sh15024
exit 0