tpclinopt.h File Reference

Header file for libtpclinopt. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "tpcextensions.h"

Go to the source code of this file.

Data Structures
struct	NNLSQDATA

Functions
int	fitLine (double *x, double *y, const int n, double *m, double *c)
int	fitLinePearson (double *x, double *y, const int n, double *m, double *msd, double *c, double *csd, double *d, double *dsd, double *r, double *ysd)
int	highestSlope (double *x, double *y, const int n, const int slope_n, double x_start, double *m, double *yi, double *xi, double *xh)
int	rootsCubic (const double a, const double b, const double c, double *x1, double *x2, double *x3)
int	rootsQuadratic (const double a, const double b, const double c, double *x1, double *x2)
int	qrLSQ (double **mat, double *rhs, double *sol, const unsigned int rows, const unsigned int cols, double *r2)
	QR least-squares solving routine.
int	qr (double **A, int m, int n, double *B, double *X, double *rnorm, double *tau, double *res, double **wws, double *ws)
int	qr_decomp (double **a, int M, int N, double *tau, double **cchain, double *chain)
int	qr_solve (double **QR, int M, int N, double *tau, double *b, double *x, double *residual, double *resNorm, double **cchain, double *chain)
int	qrWeight (int N, int M, double **A, double *b, double *weight, double *ws)
int	qrWeightRm (int N, int M, double **A, double *b, double *weight, double *ws)
int	qrSimpleLSQ (double **A, double *B, int M, int N, double *X, double *r2)
int	qrLH (const unsigned int m, const unsigned int n, double *a, double *b, double *x, double *r2)
	Solve over-determined least-squares problem A x ~ b using successive Householder rotations.
double	householder_transform (double *vector, int size)
int	householder_hm (double tau, double *vector, double **matrix, int rowNr, int columnNr)
int	householder_hv (double tau, int size, double *v, double *w)
double	householder_norm (double *v, int size)
int	nnls (double **a, int m, int n, double *b, double *x, double *rnorm, double *w, double *zz, int *index)
int	nnlsWght (int N, int M, double **A, double *b, double *weight)
int	nnlsWghtSquared (int N, int M, double **A, double *b, double *sweight)
int	nnlsq (NNLSQDATA *d, int verbose)
int	nnlsqWght (NNLSQDATA *d, double *weight)
int	nnlsqWghtSquared (NNLSQDATA *d, double *sweight)
void	nnlsqDataInit (NNLSQDATA *d)
void	nnlsqDataFree (NNLSQDATA *d)
int	nnlsqDataAllocate (NNLSQDATA *d, const int n, const int m)
int	bvls (int key, const int m, const int n, double *a, double *b, double *bl, double *bu, double *x, double *w, double *act, double *zz, int *istate, int *iter, int verbose)
	Bounded-value least-squares method to solve the linear problem A x ~ b , subject to limit1 <= x <= limit2.
int	llsqWght (int N, int M, double **A, double *a, double *b, double *weight)
int	llsqWghtSquared (int N, int M, double **A, double *a, double *b, double *weight)

Detailed Description

Header file for libtpclinopt.

Header file for template library libtpclinopt.

Author

Vesa Oikonen

Definition in file tpclinopt.h.

Function Documentation

bvls()

int bvls ( int key, const int m, const int n, double * a, double * b, double * bl, double * bu, double * x, double * w, double * act, double * zz, int * istate, int * iter, int verbose )

extern

Bounded-value least-squares method to solve the linear problem A x ~ b , subject to limit1 <= x <= limit2.

This routine is based on the text and Fortran code in C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, New Jersey, 1974, and Fortran codes by R.L. Parker and P.B. Stark, and by J. Burkardt.

See also

nnls, qrLH, qrLSQ, llsqWghtSquared, llsqWght

Returns

Returns 0 when successful, -1 if iteration count exceeded the limit, 1 in case of invalid data or settings, and 2 if problem cannot be solved.

Parameters

key	Enter 0 to solve the problem from scratch, or <>0 to initialize the routine using caller's guess about which parameters are active within their bounds, which are at their lower bounds, and which at their upper bounds, as supplied in the array istate ('warm start'). When key <> 0, the routine initially sets the active components to the averages of their upper and lower bounds.
m	Number of samples in matrix A and the length of vector b.
n	Number of parameters in matrix A and the length of vector x. The n must not be > m unless at least n-m variables are non-active (set to their bounds).
a	Pointer to matrix A; matrix must be given as an n*m array, containing n consecutive m-length vectors. Contents of A are modified in this routine.
b	Pointer to vector b of length m. Contents of b are modified in this routine.
bl	Array BL[0..n-1] of lower bounds for parameters.
bu	Array BU[0..n-1] of upper bounds for parameters.
x	Pointer to the result vector x of length n.
w	Pointer to an array of length n, which will be used as working memory, and the minimum 2-norm || a.x-b ||, (R^2), will be written in w[0].
act	Pointer to an array of length m*(n+2), or m*(m+2) if m<n, which will be used as working memory.
zz	Pointer to an array of length m, to be used as working memory.
istate	Pointer to an integer array of length n+1. If parameter key <>0, then the last position istate[n] must contain the total number of components at their bounds (nbound, the bound variables'). The absolute values of the first nbound entries of istate[] are the indices of these bound' components of x[]. The sign of istate[0.. nbound-1] entries indicates whether x[|istate[i]|] is at its upper (positive) or lower (negative) bound. Entries istate[nbound..n-1] contain the indices of the active components.
iter	Number of performed iterations is returned in this variable. Maximum nr of iterations is set with this variable, or set it to 0 to use the default maximum, 3*n.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 32 of file bvls.c.

  {
  if(verbose>0) {printf("bvls(%d, %d, %d, ...)\n", key, m, n); fflush(stdout);}
  /* Check the input */
  if(a==NULL || b==NULL || bl==NULL || bu==NULL || x==NULL || w==NULL ||
     act==NULL || zz==NULL || istate==NULL || iter==NULL || n<1 || m<1) 
  {
    if(verbose>0) fprintf(stderr, "Error: invalid input to BVLS.\n");
    return(1);
  }
 
  int maxIter=*iter; if(maxIter<3) maxIter=3*n;
  *iter=0; 
  const double eps=1.0E-13; // stopping rule
 
  /* Step 1. Initialize everything -- active and bound sets, initial values, etc. */
  if(verbose>1) {printf("step 1\n"); fflush(stdout);}
 
  /* Set mm to the smaller of matrix dimensions n and m. */
  int mm; if(m<n) mm=m; else mm=n;
  int aindx=0; // one-based index of X[] that determined the alpha; zero means not set.
  int aindxsign=0; // +1 if zz[aindx]> bu[aindx], -1 if zz[aindx]<bl[aindx]
  /* istateFromStep5 is the one-based index of the parameter that most wants to be active;
     zero value means its is currently not set. */
  int istateFromStep5 = 0;
 
  /* Check the consistency of given bounds bl[] and bu[]. */
  {
    double maxrange=0.0;
    for(int ni=0; ni<n; ni++) {
      double d=bu[ni]-bl[ni];
      if(verbose>3) printf("  bounds[%d]: %g %g\n", 1+ni, bl[ni], bu[ni]);
      if(d<0.0) {
        if(verbose>0) fprintf(stderr, "Error: inconsistent bounds in BVLS.\n"); 
        return(1);
      }
      maxrange=fmax(maxrange, d);
    }
    if(verbose>2) printf("  maxrange := %g\n", maxrange);
    if(maxrange<1.0E-10) {
      if(verbose>0) fprintf(stderr, "Error: no free variables in BVLS.\n");
      return(1);
    }
  }
 
  /* In a fresh initialization (key=0), bind all variables at their lower bounds. 
     If key<>0, use the supplied istate[] array to initialize the variables. */
  int nbound, nact; // number of bound and active parameters, respectively.
  if(key==0) {
    nbound=n;
    /* Write the indices; negative sign indicates lower limit */
    /* Note that these indices must start from 1, because 0 can not have sign */
    for(int ni=0; ni<nbound; ni++) istate[ni]=-(1+ni);
  } else {
    nbound=istate[n];
  }
  nact=n-nbound;
  if(nact>mm) {
    if(verbose>0) fprintf(stderr, "Error: too many active variables in BVLS starting solution.\n");
    return(2);
  }
  for(int ni=0; ni<nbound; ni++) {
    int i=abs(istate[ni])-1;
    if(istate[ni]<0) x[i]=bl[i]; else x[i]=bu[i];
  }
 
  /* In a warm start (key<>0, and nbound<n) initialize the active variables to the mean of
     their bounds. This is needed in case the initial QR results in active variables 
     out-of-bounds and Steps 8-11 get executed the first time through. */
  for(int ni=nbound; ni<n; ni++) {
    int i=abs(istate[ni])-1; // indices of active variables should not have signs, but let's be sure
    x[i]=0.5*(bl[i]+bu[i]);
  }
 
  /* Compute bnorm, the norm of the data vector b, for reference. */
  double bnorm=0.0;
  for(int mi=0; mi<m; mi++) bnorm+=b[mi]*b[mi];
  bnorm=sqrt(bnorm); if(verbose>2) printf("  initial_bnorm := %g\n", bnorm);
 
 
  /*
   *  Main loop
   */
  int skipStep2=0;
  double obj=0.0;
  int iact=0; // The component x[iact] is the one that most wants to become active
 
  for(*iter=1; *iter<=maxIter; ++(*iter)) {
 
    if(verbose>1) {printf("iteration %d\n", *iter); fflush(stdout);}
 
    if(!skipStep2) {
      /* Step 2. */ if(verbose>1) {printf("  step 2\n"); fflush(stdout);}
 
      /*  Initialize the negative gradient vector w(*). */
      for(int ni=0; ni<n; ni++) w[ni]=0.0;
 
      /* Compute the residual vector b-a.x , the negative gradient vector w[*], and 
         the current objective value obj = || a.x - b ||.
         The residual vector is stored in the mm+1'st column of act[*,*]. */
      obj=0.0;
      for(int mi=0; mi<m; mi++) {
        double ri=b[mi];
        for(int ni=0; ni<n; ni++) ri-=a[mi+ni*m]*x[ni];
        obj+=ri*ri;
        for(int ni=0; ni<n; ni++) w[ni]+=a[mi+ni*m]*ri;
        act[mi+mm*m]=ri;
      }
      if(verbose>3) {printf("    obj := %g\n", obj); fflush(stdout);}
 
      /* Converged?  Stop if the misfit << || b ||, or if all components are active 
        (unless this is the first iteration from a 'warm start'). */
      if((sqrt(obj)<=bnorm*eps) || ((*iter)>1 && nbound==0)) {
        if(verbose>1) {printf("bvls converged.\n"); fflush(stdout);}
        istate[n]=nbound;
        w[0]=sqrt(obj);
        return(0);
      }
 
      /* Add the contribution of the active components back into the residual. */
      for(int ni=nbound; ni<n; ni++) {
        int i=abs(istate[ni])-1;
        for(int mi=0; mi<m; mi++) act[mi+mm*m] += a[mi+i*m]*x[i];
      }
      if(verbose>9) {
        printf("Residual vector:\n");
        for(int mi=0; mi<m; mi++) printf("\t%g", act[mi+mm*m]);
        printf("\n");
      }
 
    }
 
 
 
    /* The first iteration in a 'warm start' requires immediate QR in Step 6
       but mostly we want go through Steps 3-5 first . */
    if(key!=0 && (*iter)==1) {
      if(verbose>1) printf("  'warm start' requires immediate QR in Step 6\n");
    } else {
 
      int it; // variable indicating the element in istate that most wants to be active
 
      do {
        /* Steps 3, 4. */ if(verbose>1) {printf("  steps 3 and 4\n"); fflush(stdout);}
        /* Find the bound element that most wants to be active. */
        double worst=0.0;
        it=1;
        for(int ni=0; ni<nbound; ni++) {
          int i=abs(istate[ni])-1;
          double bad; if(istate[ni] < 0) bad=-w[i]; else bad=+w[i];
          if(bad < worst) { it=ni+1; worst=bad; iact=i; }
        }
 
        /* Test whether the Kuhn-Tucker condition is met. */
        if(worst>=0.0) {
          if(verbose>1) {printf("Kuhn-Tucker condition is met.\n"); fflush(stdout);}
          istate[n]=nbound;
          w[0]=sqrt(obj);
          return(0);
        }
 
        /* The component x[iact] is the one that most wants to become active.
           If the last successful change in the active set was to move x[iact] to a bound, 
           don't let x[iact] in now: set the derivative of the misfit with respect to x[iact] 
           to zero and return to the Kuhn-Tucker test. */
        if(iact==(aindx-1)) w[aindx-1]=0.0;
      } while(iact==(aindx-1)); // Step 3 again
 
      /* Step 5. */ if(verbose>1) {printf("  step 5\n"); fflush(stdout);}
 
      /* Undo the effect of the new (potentially) active variable on the residual vector. */
      if(istate[it-1]==0) { // remove this if never happening
        if(verbose>0) fprintf(stderr, "Error: BVLS istate is zero!\n"); 
        return(1);
      }
      {
        double bnd;
        if(istate[it-1]>0) bnd=bu[iact]; else bnd=bl[iact];
        for(int mi=0; mi<m; mi++) act[mi+mm*m]+=bnd*a[mi+iact*m];
      }
 
      /* Set flag istateFromStep5, indicating that Step 6 was entered from Step 5.
         This forms the basis of a test for instability: the gradient calculation shows that x[iact] 
         wants to join the active set; if QR puts x[iact] beyond the bound from which it came, 
         the gradient calculation was in error and the variable should not have been introduced. */
      istateFromStep5=istate[it-1]; 
 
      /* Swap the indices (in istate) of the new active variable and the rightmost bound variable; 
         `unbind' that location by decrementing nbound. */
      istate[it-1]=istate[nbound-1];
      nbound--;  nact++;
      istate[nbound]=1+iact;
      if(mm<nact) {
        if(verbose>0) fprintf(stderr, "Error: too many free variables in BVLS.\n");
        return(2);
      }
 
    } // finalized steps 3-5
 
    do {
 
      skipStep2=0;
 
      /* Step 6. */ if(verbose>1) {printf("  step 6\n"); fflush(stdout);}
 
      /* Load array act with the appropriate columns of A for QR. For added stability, reverse 
         the column ordering so that the most recent addition to the active set is in the last 
         column. Also copy the residual vector from act[., mm] into act[., mm+1]. */
      for(int mi=0; mi<m; mi++) {
        act[mi+(mm+1)*m]=act[mi+mm*m]; // vector b for QR
        for(int ni=nbound; ni<n; ni++) {
          int i=abs(istate[ni])-1;
          act[mi+(nact+nbound-ni-1)*m]=a[mi+i*m];
        }
      }
      if(verbose>9) {
        printf("Matrix A for QR:\n");
        for(int ni=0; ni<nact; ni++) {
          for(int mi=0; mi<m; mi++) printf("\t%g", act[mi+ni*nact]);
          printf("\n");
        }
        printf("Vector B for QR:\n");
        for(int mi=0; mi<m; mi++) printf("\t%g", act[(mm+1)*m + mi]);
        printf("\n");
      }
 
      /* Test for linear dependence in QR, and for an instability that moves the variable 
         just introduced away from the feasible region (rather than into the region or 
         all the way through it). In either case, remove the latest vector introduced from 
         the active set and adjust the residual vector accordingly.
         Set the gradient component (w[iact]) to zero and return to the Kuhn-Tucker test. */
      double r2;
      if(qrLH(m, nact, act, &act[(mm+1)*m], zz, &r2) !=0 || 
         (istateFromStep5>0 && zz[nact-1]>bu[iact]) || 
         (istateFromStep5<0 && zz[nact-1]<bl[iact]) )
      {
        nbound++;
        if(bu[iact]>x[iact]) istate[nbound-1]=-istate[nbound-1];
        nact--;
        for(int mi=0; mi<m; mi++) act[mi+mm*m]-=x[iact]*a[mi+iact*m];
        istateFromStep5 = 0; // not from step 5
        w[iact]=0.0;
        skipStep2=1; // we want to skip Step 2 and go directly to Step 3
        if(verbose>3) {printf("    going from step 6 to step 3\n"); fflush(stdout);}
        break; // go to step 3
      }
 
      /* If Step 6 was entered from Step 5 and we are here, a new variable has been successfully 
         introduced into the active set; the last variable that was fixed at a bound is again 
         permitted to become active. */
      if(istateFromStep5!=0) aindx=0;
      istateFromStep5=0;
 
      /* Step 7. */ if(verbose>1) {printf("  step 7\n"); fflush(stdout);}
      /* Check for strict feasibility of the new QR solution. */
      int qr_solution_feasible=1;
      int indexHolder=0;
      if(verbose>8) printf("    nact=%d  nbound=%d\n", nact, nbound);
      for(int ni=0; ni<nact; ni++) {
        indexHolder=ni; // Loop in step 8 will start from this
        int i=abs(istate[ni+nbound])-1;
        if(verbose>8) {
          printf("      istate[%d]=%d\n", ni+nbound, 1+i);
          printf("      zz[%d]=%g  bl[%d]=%g  bu[%d]=%g\n", nact-ni-1, zz[nact-ni-1], i, bl[i], i, bu[i]);
        }
        if(zz[nact-ni-1]<bl[i] || zz[nact-ni-1]>bu[i]) {
          if(verbose>3) {printf("    new iterate is not feasible\n"); fflush(stdout);}
          qr_solution_feasible=0; break; // go to Step 8
        }
      }
      if(verbose>8) printf("    indexHolder=%d\n", indexHolder);
      if(qr_solution_feasible) {
        if(verbose>3) {printf("    new iterate is feasible\n"); fflush(stdout);}
        for(int ni=0; ni<nact; ni++) {
          int i=abs(istate[ni+nbound])-1;
          x[i]=zz[nact-ni-1];
        }
        /* New iterate is feasible; back to the top. */
        break; // Back to the start of the main loop
      }
 
      { // keep local variables alpha and alf local
        double alpha=2.0;
        /* Steps 8 and 9 */ if(verbose>1) {printf("  steps 8 and 9\n"); fflush(stdout);}
        double alf=alpha;
        for(int ni=indexHolder; ni<nact; ni++) {
          int i=abs(istate[ni+nbound])-1;
          if(zz[nact-ni-1] > bu[i]) alf=(bu[i]-x[i])/(zz[nact-ni-1]-x[i]);
          if(zz[nact-ni-1] < bl[i]) alf=(bl[i]-x[i])/(zz[nact-ni-1]-x[i]);
          if(alf<alpha) {
            alpha=alf;
            aindx=1+i;
            if((zz[nact-ni-1]-bl[i])<0.0) aindxsign=-1; else aindxsign=+1;
          }
        }
        /* Step 10 */ if(verbose>1) {printf("  step 10\n"); fflush(stdout);}
        for(int ni=0; ni<nact; ni++) {
          int i=abs(istate[ni+nbound])-1;
          x[i]+=alpha*(zz[nact-ni-1]-x[i]);
        }
      }
 
      /* Step 11 */ if(verbose>1) {printf("  step 11\n"); fflush(stdout);}
      /* Move the variable that determined alpha to the appropriate bound
         (aindx is its index; sj is + if zz[aindx]> bu[aindx], - if zz[aindx]<bl[aindx] ).
         If any other component of  x  is infeasible at this stage, it must be due to round-off.
         Bind every infeasible component and every component at a bound to the appropriate bound.
         Correct the residual vector for any variables moved to bounds. Since at least one 
         variable is removed from the active set in this step, Loop B
         (Steps 6-11) terminates after at most nact steps. */
      {
        int noldb=nbound;
        for(int ni=0; ni<nact; ni++) {
          int i=abs(istate[ni+noldb])-1;
          if((bu[i]-x[i]<=0.0) || (i==(aindx-1) && aindxsign>0)) {
            /* Move x[i] to its upper bound. */
            x[i]=bu[i];
            istate[ni+noldb]=istate[nbound]; istate[nbound]=+(1+i); nbound++;
            for(int mi=0; mi<m; mi++) act[mi+mm*m]-=bu[i]*a[mi+i*m];
          } else if( ((x[i]-bl[i])<=0.0) || (i==(aindx-1) && aindxsign<0)) {
            /* Move x(j) to its lower bound. */
            x[i]=bl[i];
            istate[ni+noldb]=istate[nbound]; istate[nbound]=-(1+i); nbound++;
            for(int mi=0; mi<m; mi++) act[mi+mm*m]-=bl[i]*a[mi+i*m];
          }
        }
        nact=n-nbound;
      }
      /* If there are still active variables left, repeat the QR; if not, go back to step 6. */
    } while(nact>0);
 
  } // main loop
 
  /* iterMax reached */
  if(verbose>0) fprintf(stderr, "Error: BVLS fails to converge.\n");
  return(-1);
}

fitLine()

int fitLine ( double * x, double * y, const int n, double * m, double * c )

extern

Calculate regression line slope (m) and y axis intercept (c).

See also

fitLinePearson, statMeanSD, doubleMean, statSortDouble, highestSlope, linRegr3p

Returns

Returns the number of samples used in the regression, 0 in case of an error.

Parameters

x	Pointer to an array of x axis values. NaN's and infinite values are ignored.
y	Pointer to an array of y axis values. NaN's and infinite values are ignored.
n	The number of samples (length of x[] and y[]).
m	Pointer where calculated slope is written; enter NULL if not needed.
c	Pointer where calculated y axis intercept is written; enter NULL if not needed.

Definition at line 23 of file regression.c.

  {
  /* Check the data */
  if(x==NULL || y==NULL) return(0);
  if(m!=NULL) *m=nan("");
  if(c!=NULL) *c=nan("");
  if(n<1) return(0);
 
  double xsum=0.0, ysum=0.0, x2sum=0.0, xysum=0.0;
  int nn=0;
  for(int i=0; i<n; i++) if(isfinite(x[i]) && isfinite(y[i])) {
    xsum+=x[i]; ysum+=y[i]; x2sum+=x[i]*x[i]; xysum+=x[i]*y[i]; nn++;
  } //printf("nn=%d\n", nn);
  // printf("xsum=%g ysum=%g x2sum=%g xysum=%g\n", xsum, ysum, x2sum, xysum);
  if(nn<1) return(0);
  if(nn==1) {
    double mm=ysum/xsum;
    if(isfinite(mm)) {
      if(m!=NULL) *m=mm;
      if(c!=NULL) *c=0.0;
      return(nn);
    } else
      return(0);
  }
  double delta=(double)nn*x2sum - xsum*xsum; //printf("delta=%g\n", delta);
  if(delta==0.0) return(0);
  if(m!=NULL) *m=((double)nn*xysum - xsum*ysum)/delta;
  if(c!=NULL) *c=(x2sum*ysum - xsum*xysum)/delta;
  return(nn);
}

Referenced by highestSlope().

fitLinePearson()

int fitLinePearson ( double * x, double * y, const int n, double * m, double * msd, double * c, double * csd, double * d, double * dsd, double * r, double * ysd )

extern

Calculate regression line slope (m), y axis intercept (c), their SDs, and Pearson's correlation coefficient (r).

See also

fitLine, statMeanSD, doubleMean, statSortDouble

Returns

Returns the number of samples used in the regression, 0 in case of an error.

Todo

Check SD of x axis intercept against some other program.

Parameters

x	Pointer to an array of x axis values. NaN's and infinite values are ignored.
y	Pointer to an array of y axis values. NaN's and infinite values are ignored.
n	The number of samples (length of x[] and y[]).
m	Pointer where calculated slope is written; enter NULL if not needed.
msd	Pointer where SD of slope is written; enter NULL if not needed.
c	Pointer where calculated y axis intercept is written; enter NULL if not needed.
csd	Pointer where SD of y axis intercept is written; enter NULL if not needed.
d	Pointer where calculated x axis intercept is written; enter NULL if not needed.
dsd	Pointer where SD of x axis intercept is written; enter NULL if not needed.
r	Pointer where Pearson's correlation coefficient is written; enter NULL if not needed.
ysd	Pointer where residual variance of y values is written; enter NULL if not needed.

Definition at line 72 of file regression.c.

  {
  /* Check the data */
  if(x==NULL || y==NULL) return(0);
  if(m!=NULL) *m=nan("");
  if(c!=NULL) *c=nan("");
  if(d!=NULL) *d=nan("");
  if(msd!=NULL) *msd=nan("");
  if(csd!=NULL) *csd=nan("");
  if(dsd!=NULL) *dsd=nan("");
  if(r!=NULL) *r=nan("");
  if(ysd!=NULL) *ysd=nan("");
  if(n<1) return(0);
 
  double xsum=0.0, ysum=0.0, x2sum=0.0, y2sum=0.0, xysum=0.0;
  int nn=0;
  for(int i=0; i<n; i++) if(isfinite(x[i]) && isfinite(y[i])) {
    xsum+=x[i]; ysum+=y[i]; x2sum+=x[i]*x[i]; y2sum+=y[i]*y[i]; xysum+=x[i]*y[i]; nn++;
  }
  if(nn<1) return(0);
  if(nn==1) {
    double mm=ysum/xsum;
    if(isfinite(mm)) {
      if(m!=NULL) *m=mm;
      if(c!=NULL) *c=0.0;
      if(d!=NULL) {*d=0.0; if(fabs(mm)<1.0E-10) *d=nan("");}
      if(msd!=NULL) *msd=0.0;
      if(csd!=NULL) *csd=0.0;
      if(dsd!=NULL) {*dsd=0.0; if(fabs(mm)<1.0E-10) *dsd=nan("");}
      if(r!=NULL) *r=1.0;
      if(ysd!=NULL) *ysd=0.0;
      return(nn);
    } else
      return(0);
  }
  double xmean=xsum/(double)nn;
  double ymean=xsum/(double)nn;
  double dx2sum=0.0, /*dy2sum=0.0,*/ dxdysum=0.0;
  for(int i=0; i<n; i++) if(isfinite(x[i]) && isfinite(y[i])) {
    double f=x[i]-xmean; dx2sum+=f*f;
    double g=y[i]-ymean; /*dy2sum+=g*g;*/
    dxdysum+=f*g;
  }
  // slope and intercept
  double slope=dxdysum/dx2sum;
  double ic=(dx2sum*ysum - xsum*dxdysum)/((double)nn*dx2sum);
  if(!isfinite(slope) || !isfinite(ic)) return(0);
  if(m!=NULL) *m=slope;
  if(c!=NULL) *c=ic;
  double xic=-ic/slope; if(!isfinite(xic)) xic=nan("");
  if(d!=NULL) *d=xic;
  // Errors
  double ry2sum=0.0;
  for(int i=0; i<n; i++) if(isfinite(x[i]) && isfinite(y[i])) {
    double f=slope*x[i]+ic-y[i];
    ry2sum+=f*f;
  }
  double ye=sqrt(ry2sum/(double)(nn-2)); if(!isfinite(ye)) ye=0.0;
  if(ysd!=NULL) *ysd=ye;
  double slopee=ye/sqrt(dx2sum); if(!isfinite(slopee)) slopee=0.0;
  double ice=ye/sqrt((double)nn-xsum*xsum/x2sum); if(!isfinite(ice)) ice=0.0;
  if(msd!=NULL) *msd=slopee;
  if(csd!=NULL) *csd=ice;
  if(dsd!=NULL) { // SD of x axis intercept
    double xice=fabs(ye/slope)*sqrt((1.0/(double)nn) + ymean*ymean/(slope*slope*dx2sum));
    if(!isfinite(xice)) xice=nan("");
    *dsd=xice;
  }
  // Pearson's correlation coefficient
  double rr=(xysum-((xsum*ysum)/(double)nn)) /
            sqrt((x2sum-xsum*xsum/(double)nn)*(y2sum-ysum*ysum/(double)nn));
  // and correct for small sample size
  if(nn>4) rr*=1.0+(1.0-rr*rr)/(double)(2*(nn-4));
  if(!isfinite(rr)) rr=0.0; else rr=fabs(rr);
  if(r!=NULL) *r=rr;
 
  return(nn);
}

highestSlope()

int highestSlope ( double * x, double * y, const int n, const int slope_n, double x_start, double * m, double * yi, double * xi, double * xh )

extern

Find the regression line with the highest slope for x,y data.

See also

fitLine

Returns

Return 0 if no errors were found.

Parameters

x	Pointer to an array of x axis values. NaN's and infinite values are ignored. Data must be sorted by increasing x, and overlapping x values may cause problem.
y	Pointer to an array of y axis values. NaN's and infinite values are ignored. Data is not modified.
n	The number of samples (length of x[] and y[]).
slope_n	The number of samples used to fit the line; must be at least 2.
x_start	Estimation start x value, samples with smaller x are ignored; can usually be set to zero.
m	Pointer where calculated max slope is written; NULL if not needed.
yi	Pointer where calculated y axis intercept is written; NULL if not needed.
xi	Pointer where calculated x axis intercept is written; NULL if not needed. This could be used as an estimate of radioactivity appearance time in TAC data, but you must then check that the max slope m is positive.
xh	Pointer where the place (x) of the highest slope is written; NULL if not needed.

Definition at line 179 of file regression.c.

  {
  /* Check the data */
  if(x==NULL || y==NULL || n<2 || slope_n<2) return(1);
  if(m!=NULL) *m=nan("");
  if(yi!=NULL) *yi=nan("");
  if(xi!=NULL) *xi=nan("");
  if(xh!=NULL) *xh=nan("");
 
  /* Make copy of original data */
  double xx[n], yy[n]; int nn=0;
  for(int i=0; i<n; i++)
    if(isfinite(x[i]) && isfinite(y[i])) {
      if(isfinite(x_start) && x[i]<x_start) continue;
      xx[nn]=x[i]; yy[nn]=y[i]; nn++;
    }
  if(nn<slope_n) return(2);
 
  /* Compute regression lines */
  double max_slope=nan(""), ic_at_max=nan(""), slope, ic;
  int i_at_max=0;
  for(int i=0; i<=nn-slope_n; i++) {
    if(fitLine(xx+i, yy+i, slope_n, &slope, &ic)<slope_n) continue;
    if(!isfinite(slope) || !isfinite(ic)) continue;
    if(isnan(max_slope) || slope>max_slope) {
      max_slope=slope; ic_at_max=ic; i_at_max=i;
    }
  }
  if(!isfinite(max_slope)) return(3);
  if(m!=NULL) *m=max_slope;
  if(yi!=NULL) *yi=ic_at_max;
  if(xi!=NULL) {if(max_slope!=0.0) *xi=-ic_at_max/max_slope;}
  if(xh!=NULL) {
    *xh=0.0; for(int i=i_at_max; i<i_at_max+slope_n; i++) *xh+=xx[i];
    *xh/=(double)slope_n;
  }
 
  return(0);
}

householder_hm()

int householder_hm ( double tau, double * vector, double ** matrix, int rowNr, int columnNr )

extern

Applies a householder transformation defined by vector "vector" and scalar tau to the left-hand side of the matrix. (I - tau vector vector^T)*matrix The result of the transform is stored in matrix.

Returns

Returns 0 if ok.

Parameters

tau	Coefficient defining householder transform.
vector	Vector defining householder transform (of size rowNr).
matrix	the matrix that is to be transformed.
rowNr	Nr of rows in matrix.
columnNr	Nr of columns in matrix.

Definition at line 87 of file hholder.c.

  {
  int i, j;
  double wj;
 
  if(tau==0.0) return(0); // success
  if(rowNr<1 || columnNr<1) return(1);
  for(j=0; j<columnNr; j++) {
    /* Compute wj = vk Akj */
    wj=matrix[0][j];
    for(i=1; i<rowNr; i++)  /* note, computed for v(0) = 1 above */
      wj += vector[i]*matrix[i][j];
    /* Aij = Aij - tau vi wj */
    /* i = 0 */
    matrix[0][j]-=tau*wj;
    /* i = 1 .. M-1 */
    for(i=1; i<rowNr; i++) matrix[i][j]-=tau*vector[i]*wj;
  }
  return(0);
}

Referenced by qr_decomp().

householder_hv()

int householder_hv ( double tau, int size, double * v, double * w )

extern

Applies a householder transformation defined by vector v and coefficient tau to vector w w = (I - tau v v^T) w.

Returns

Returns 0 if ok.

Parameters

tau	Coefficient defining householder transform.
size	Size of vectors v and w.
v	Vector v.
w	Vector w.

Definition at line 126 of file hholder.c.

  {
  int i;
  double d;
 
  if(tau==0) return(0); // success
  if(size<1) return(1);
  /* d = v'w */
  d=w[0]; for(i=1; i<size; i++) d+=v[i]*w[i];
  /* w = w - tau (v) (v'w) */
  w[0]-=tau*d;
  for(i=1; i<size; i++) w[i]-=tau*v[i]*d;
  return(0);
}

householder_norm()

double householder_norm ( double * v, int size )

extern

Calculates the euclidean norm of vector v[].

Returns

Returns the euclidean norm of a vector.

Parameters

v	Vector v.
size	Size of vector v[].

Definition at line 155 of file hholder.c.

  {
  double ss=0.0;
  int i;
 
  for(i=0; i<size; i++) ss+=v[i]*v[i];
  return sqrt(ss);
}

householder_transform()

double householder_transform ( double * v, int N )

extern

This function prepares a Householder transformation P = I - tau h h^T which can be used to zero all the elements of the input vector except the first one that will get value beta. On output the elements 1 - size-1 of the vector h are stored in locations vector[1] - vector[size-1] of the input vector and value of beta is stored in location vector[0].

Returns

The scalar tau is returned.

Parameters

v	The N-vector to be transformed.
N	size of the vector.

Definition at line 33 of file hholder.c.

  {
  double vnorm, alpha, beta, tau;
  int n;
 
  if(N<1) return 0.0; // tau = 0
  /* Euclidean norm of the vector starting from the second value */
  vnorm=0.0; for(n=1; n<N; n++) vnorm+=v[n]*v[n];
  vnorm=sqrt(vnorm); if(isnan(vnorm) || vnorm==0.0) return 0.0; // tau = 0
 
  /* Computing the coefficient tau */
  alpha=v[0];
  beta= - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, vnorm);
  tau=(beta-alpha)/beta ;
 
  /* Scale the Householder vector so that the first element will be 1.
   * (Scaling is also affecting the coefficient tau).
   * Without scaling, the first element would have value (alpha - beta). */
#if(0) // Kaisa's version
  {
    double os=1.0/(alpha-beta);
    for(n=1; n<N; n++) v[n]*=os;
    v[0]=beta;
  }
#else // Vesa's version
  {
    double s=alpha-beta;
    if(fabs(s)>DBL_MIN) {
      v[0]=beta;
      for(n=1; n<N; n++) v[n]*=(1.0/s);
    } else {
      v[0]=beta;
      for(n=1; n<N; n++) v[n]*=(doubleMachEps()/s);
      for(n=1; n<N; n++) v[n]*=(1.0/doubleMachEps());
    }
  }
#endif
 
  return tau;
}

Referenced by qr_decomp().

llsqWght()

int llsqWght ( int N, int M, double ** A, double * a, double * b, double * weight )

extern

Algorithm for weighting the problem that is given to a LLSQ algorithm.

Square roots of weights are used because in LLSQ algorithms the difference w*A-w*b is squared.

See also

bvls, llsqWghtSquared, nnlsWghtSquared

Returns

Algorithm returns zero if successful, 1 if arguments are inappropriate.

Parameters

N	Matrix A dimension N (nr of parameters).
M	Matrix A dimension M (nr of samples).
A	Pointer to matrix A[N][M]; enter NULL to use following matrix format instead.
a	Pointer to matrix A[N*M]; enter NULL to use previous matrix format instead.
b	Vector B of length M.
weight	Weights for each sample (array of length M).

Definition at line 425 of file bvls.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || (A==NULL && a==NULL) || b==NULL || weight==NULL) return(1);
 
  /* Allocate memory */
  w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
 
  /* Check that weights are not zero and get the square roots of them to w[] */
  for(m=0; m<M; m++) {
    if(weight[m]<=1.0e-20) w[m]=0.0;
    else w[m]=sqrt(weight[m]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) {
      if(A!=NULL) A[n][m]*=w[m];
      if(a!=NULL) a[m+n*M]*=w[m];
    }
    b[m]*=w[m];
  }
 
  free(w);
  return(0);
}

llsqWghtSquared()

int llsqWghtSquared ( int N, int M, double ** A, double * a, double * b, double * sweight )

extern

Algorithm for weighting the problem that is given to a LLSQ algorithm.

Square roots of weights are used because in LLSQ algorithms the difference w*A-w*b is squared.

Here user must give squared weights; this makes calculation faster, when this function needs to be called many times.

See also

llsqWght, bvls

Returns

Algorithm returns zero if successful, 1 if arguments are inappropriate.

Parameters

N	Matrix A dimension N (nr of parameters).
M	Matrix A dimension M (nr of samples).
A	Pointer to matrix A[N][M]; enter NULL to use following matrix format instead.
a	Pointer to matrix A[N*M]; enter NULL to use previous matrix format instead.
b	Vector B of length M.
sweight	Squared weights for each sample (array of length M).

Definition at line 480 of file bvls.c.

  {
  int n, m;
 
  /* Check the arguments */
  if(N<1 || M<1 || (A==NULL && a==NULL) || b==NULL || sweight==NULL) return(1);
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) {
      if(A!=NULL) A[n][m]*=sweight[m];
      if(a!=NULL) a[m+n*M]*=sweight[m];
    }
    b[m]*=sweight[m];
  }
 
  return(0);
}

nnls()

int nnls ( double ** a, int m, int n, double * b, double * x, double * rnorm, double * wp, double * zzp, int * indexp )

extern

Algorithm NNLS (Non-negative least-squares).

Given an m by n matrix A, and an m-vector B, computes an n-vector X, that solves the least squares problem A * X = B , subject to X>=0

Instead of pointers for working space, NULL can be given to let this function to allocate and free the required memory.

See also

qrLSQ, nnlsWghtSquared

Returns

Function returns 0 if successful, 1, if iteration count exceeded 3*N, or 2 in case of invalid problem dimensions or memory allocation error.

Parameters

a	On entry, a[ 0... N ][ 0 ... M ] contains the M by N matrix A. On exit, a[][] contains the product matrix Q*A, where Q is an m by n orthogonal matrix generated implicitly by this function.
m	Matrix dimension m.
n	Matrix dimension n.
b	On entry, b[] must contain the m-vector B. On exit, b[] contains Q*B.
x	On exit, x[] will contain the solution vector.
rnorm	On exit, rnorm contains the squared Euclidean norm of the residual vector, R^2 (Sum-of-Squares). If NULL is given, no rnorm is calculated.
wp	An n-array of working space, wp[]. On exit, wp[] will contain the dual solution vector. wp[i]=0.0 for all i in set p and wp[i]<=0.0 for all i in set z. Can be NULL, which causes this algorithm to allocate memory for it.
zzp	An m-array of working space, zz[]. Can be NULL, which causes this algorithm to allocate memory for it.
indexp	An n-array of working space, index[]. Can be NULL, which causes this algorithm to allocate memory for it.

Definition at line 43 of file nnls.c.

  {
  /* Check the parameters and data */
  if(m<=0 || n<=0 || a==NULL || b==NULL || x==NULL) return(2);
  if(rnorm!=NULL) *rnorm=nan("");
  /* Allocate memory for working space, if required */
  int *index=NULL;
  double *w=NULL, *zz=NULL;
  if(wp!=NULL) w=wp; else w=(double*)calloc(n, sizeof(double));
  if(zzp!=NULL) zz=zzp; else zz=(double*)calloc(m, sizeof(double));
  if(indexp!=NULL) index=indexp; else index=(int*)calloc(n, sizeof(int));
  if(w==NULL || zz==NULL || index==NULL) return(2);
 
  /* Initialize the arrays INDEX[] and X[] */
  for(int ni=0; ni<n; ni++) {x[ni]=0.; index[ni]=ni;}
  int iz1=0;
  int iz2=n-1;
  int nsetp=0;
  int npp1=0;
 
  /* Main loop; quit if all coefficients are already in the solution or
     if M cols of A have been triangulated */
  double up=0.0;
  int itmax; if(n<3) itmax=n*3; else itmax=n*n;
  int iter=0; 
  int k, j=0, jj=0;
  while(iz1<=iz2 && nsetp<m) {
    /* Compute components of the dual (negative gradient) vector W[] */
    for(int iz=iz1; iz<=iz2; iz++) {
      int ni=index[iz];
      double sm=0.;
      for(int mi=npp1; mi<m; mi++) sm+=a[ni][mi]*b[mi];
      w[ni]=sm;
    }
 
    double wmax;
    int izmax=0;
    while(1) {
 
      /* Find largest positive W[j] */
      wmax=0.0;
      for(int iz=iz1; iz<=iz2; iz++) {
        int i=index[iz];
        if(w[i]>wmax) {wmax=w[i]; izmax=iz;}
      }
 
      /* Terminate if wmax<=0.; it indicates satisfaction of the Kuhn-Tucker conditions */
      if(wmax<=0.0) break;
      j=index[izmax];
 
      /* The sign of W[j] is ok for j to be moved to set P.
         Begin the transformation and check new diagonal element to avoid
         near linear dependence. */
      double asave=a[j][npp1];
      up=0.0;
      nnls_lss_h12(1, npp1, npp1+1, m, a[j], &up, NULL);
      double unorm=0.0;
      if(nsetp!=0) for(int mi=0; mi<nsetp; mi++) unorm+=a[j][mi]*a[j][mi];
      unorm=sqrt(unorm);
      double d=unorm+fabs(a[j][npp1])*0.01;
      if((d-unorm)>0.0) {
        /* Col j is sufficiently independent. Copy B into ZZ, update ZZ
           and solve for ztest ( = proposed new value for X[j] ) */
        for(int mi=0; mi<m; mi++) zz[mi]=b[mi];
        nnls_lss_h12(2, npp1, npp1+1, m, a[j], &up, zz);
        double ztest=zz[npp1]/a[j][npp1];
        /* See if ztest is positive */
        if(ztest>0.) break;
      }
 
      /* Reject j as a candidate to be moved from set Z to set P. Restore
         A[npp1,j], set W[j]=0., and loop back to test dual coefficients again */
      a[j][npp1]=asave; w[j]=0.;
    } /* while(1) */
    if(wmax<=0.0) break;
 
    /* Index j=INDEX[izmax] has been selected to be moved from set Z to set P.
       Update B and indices, apply householder transformations to cols in
       new set Z, zero sub-diagonal elements in col j, set W[j]=0. */
    for(int mi=0; mi<m; mi++) b[mi]=zz[mi];
    index[izmax]=index[iz1]; index[iz1]=j; iz1++; nsetp=npp1+1; npp1++;
    if(iz1<=iz2)
      for(int jz=iz1; jz<=iz2; jz++) {
        jj=index[jz];
        nnls_lss_h12(2, nsetp-1, npp1, m, &a[j][0], &up, a[jj]);
      }
    if(nsetp!=m) for(int mi=npp1; mi<m; mi++) a[j][mi]=0.;
    w[j]=0.;
    
    /* Solve the triangular system; store the solution temporarily in Z[] */
    for(int mi=0; mi<nsetp; mi++) {
      int ip=nsetp-(mi+1);
      if(mi!=0) for(int ii=0; ii<=ip; ii++) zz[ii]-=a[jj][ii]*zz[ip+1];
      jj=index[ip]; zz[ip]/=a[jj][ip];
    }
 
    /* Secondary loop begins here */
    while(++iter<itmax) {
      /* See if all new constrained coefficients are feasible; if not, compute alpha */
      double alpha=2.0;
      for(int ip=0; ip<nsetp; ip++) {
        int ni=index[ip];
        if(zz[ip]<=0.) {
          double t=-x[ni]/(zz[ip]-x[ni]);
          if(alpha>t) {alpha=t; jj=ip-1;}
        }
      }
 
      /* If all new constrained coefficients are feasible then still alpha==2.
         If so, then exit from the secondary loop to main loop */
      if(alpha==2.0) break;
 
      /* Use alpha (0.<alpha<1.) to interpolate between old X and new ZZ */
      for(int ip=0; ip<nsetp; ip++) {
        int ni=index[ip]; x[ni]+=alpha*(zz[ip]-x[ni]);
      }
 
      /* Modify A and B and the INDEX arrays to move coefficient i from set P to set Z. */
      int pfeas=1;
      k=index[jj+1];
      do {
        x[k]=0.;
        if(jj!=(nsetp-1)) {
          jj++;
          for(int ni=jj+1; ni<nsetp; ni++) {
            int ii=index[ni]; index[ni-1]=ii;
            double ss, cc;
            nnls_lss_g1(a[ii][ni-1], a[ii][ni], &cc, &ss, &a[ii][ni-1]);
            a[ii][ni]=0.0;
            for(int nj=0; nj<n; nj++) if(nj!=ii) {
              /* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
              double temp=a[nj][ni-1];
              a[nj][ni-1]=cc*temp+ss*a[nj][ni];
              a[nj][ni]=-ss*temp+cc*a[nj][ni];
            }
            /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
            double temp=b[ni-1]; b[ni-1]=cc*temp+ss*b[ni]; b[ni]=-ss*temp+cc*b[ni];
          }
        }
        npp1=nsetp-1; nsetp--; iz1--; index[iz1]=k;
 
        /* See if the remaining coefficients in set P are feasible; they should be because of 
           the way alpha was determined. If any are infeasible it is due to round-off error. 
           Any that are non-positive will be set to zero and moved from set P to set Z. */
        for(jj=0, pfeas=1; jj<nsetp; jj++) {
          k=index[jj]; if(x[k]<=0.) {pfeas=0; break;}
        }
      } while(pfeas==0);
 
      /* Copy B[] into zz[], then solve again and loop back */
      for(int mi=0; mi<m; mi++) zz[mi]=b[mi];
      for(int mi=0; mi<nsetp; mi++) {
        int ip=nsetp-(mi+1);
        if(mi!=0) for(int ii=0; ii<=ip; ii++) zz[ii]-=a[jj][ii]*zz[ip+1];
        jj=index[ip]; zz[ip]/=a[jj][ip];
      }
    } /* end of secondary loop */
 
    if(iter>=itmax) break;
    for(int ip=0; ip<nsetp; ip++) {k=index[ip]; x[k]=zz[ip];}
  } /* end of main loop */
 
  if(wp != NULL) {
    if(npp1>=m) for(int ni=0; ni<n; ni++) w[ni]=0.;
  }
 
  /* Compute the norm of the final residual vector (sum-of-squares) */
  if(rnorm != NULL) {
    double ss=0.0;
    if(npp1<m) for(int mi=npp1; mi<m; mi++) ss+=(b[mi]*b[mi]);
    *rnorm=ss;
  }
 
  /* Free working space, if it was allocated here */
  if(wp==NULL) free(w); 
  if(zzp==NULL) free(zz); 
  if(indexp==NULL) free(index);
  if(iter>=itmax) return(1);
  return(0);
} /* nnls */

Referenced by spectralDExp(), spectralDMSurge(), and tacDelayCMFit().

nnlsq()

int nnlsq ( NNLSQDATA * d, int verbose )

extern

Algorithm NNLS (Non-negative least-squares).

The same algorithm as nnls(), but this one gets its data in a struct, and max iteration number and dependence factor (previously fixed to 0.01) can be set as function parameters in the struct.

Given an m by n matrix A, and an m-vector B, computes an n-vector X, that solves the least squares problem A * X = B , subject to X>=0

Instead of pointers for working space, NULL can be given to let this function to allocate and free the required memory.

See also

nnlsqDataInit, nnlsqDataAllocate, nnlsqDataFree, nnlsqWghtSquared, nnls, qrLSQ

Returns

Function returns 0 if successful, 1, if iteration count exceeded specified limit, 2 in case of invalid problem dimensions, and >2 in case of other errors.

Parameters

d	Pointer to structure containing the data and place for the results.

Precondition

Initiate the structure, allocate memory for the data, and possibly weight the data.

Postcondition

Free the allocated memory.

See also

nnlsqDataInit, nnlsqDataAllocate, nnlsqDataFree

Parameters

verbose

Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 55 of file nnlsq.c.

  {
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
 
  /* Check the data */
  if(d==NULL || d->m<1 || d->n<1 || d->a==NULL || d->b==NULL || 
     d->x==NULL || d->w==NULL || d->zz==NULL || d->index==NULL)
       return(2);
  if(d->depf<1.0E-08 || d->depf>0.5) return(3);
 
  /* Initialize */
  for(int ni=0; ni<d->n; ni++) d->x[ni]=0.0;
  for(int ni=0; ni<d->n; ni++) d->index[ni]=ni;
  d->rnorm=nan("");
  int iz1=0;
  int iz2=d->n-1;
  int nsetp=0;
  int npp1=0;
 
 
  /* Main loop; quit if all coefficients are already in the solution or
     if M cols of A have been triangulated */
  double up=0.0;
  int itermax, iter=0; 
  if(d->iternr>=3) itermax=d->iternr; else itermax=3*d->n;
  int j=0, jj=0;
  while(iz1<=iz2 && nsetp<d->m) {
    /* Compute components of the dual (negative gradient) vector W[] */
    for(int iz=iz1; iz<=iz2; iz++) {
      int ni=d->index[iz];
      double sm=0.;
      for(int mi=npp1; mi<d->m; mi++) sm+=d->a[ni][mi]*d->b[mi];
      d->w[ni]=sm;
    }
 
    double wmax;
    int izmax=0;
    while(1) {
 
      /* Find largest positive W */
      wmax=0.0;
      for(int iz=iz1; iz<=iz2; iz++) {
        int i=d->index[iz];
        if(d->w[i]>wmax) {wmax=d->w[i]; izmax=iz;}
      }
 
      /* Terminate if wmax<=0.; it indicates satisfaction of the Kuhn-Tucker conditions */
      if(wmax<=0.0) break;
 
      /* The sign of W[j] is ok for j to be moved to set P.
         Begin the transformation and check new diagonal element to avoid near linear dependence. */
      j=d->index[izmax];
      double asave=d->a[j][npp1];
//      up=0.0;
      if(nnlsq_lss_h12(1, npp1, npp1+1, d->m, d->a[j], &up, NULL)) return(2);
      double unorm=0.0;
      if(nsetp!=0) for(int mi=0; mi<nsetp; mi++) unorm+=d->a[j][mi]*d->a[j][mi];
      unorm=sqrt(unorm);
      double e=unorm+fabs(d->a[j][npp1])*d->depf;
      if((e-unorm)>0.0) {
        /* Col j is sufficiently independent. Copy B into ZZ, update ZZ
           and solve for ztest ( = proposed new value for X[j] ) */
        for(int mi=0; mi<d->m; mi++) d->zz[mi]=d->b[mi];
        nnlsq_lss_h12(2, npp1, npp1+1, d->m, d->a[j], &up, d->zz);
        double ztest=d->zz[npp1]/d->a[j][npp1];
        /* See if ztest is positive */
        if(ztest>0.) break;
      }
 
      /* Reject j as a candidate to be moved from set Z to set P. Restore
         A[npp1,j], set W[j]=0., and loop back to test dual coefficients again */
      d->a[j][npp1]=asave; d->w[j]=0.;
    } /* while(1) */
    if(wmax<=0.0) break;
 
    /* Index j=INDEX[izmax] has been selected to be moved from set Z to set P.
       Update B and indices, apply householder transformations to cols in
       new set Z, zero sub-diagonal elements in col j, set W[j]=0. */
    for(int mi=0; mi<d->m; mi++) d->b[mi]=d->zz[mi];
    d->index[izmax]=d->index[iz1]; d->index[iz1]=j; iz1++; nsetp=npp1+1; npp1++;
    if(iz1<=iz2)
      for(int jz=iz1; jz<=iz2; jz++) {
        jj=d->index[jz];
        nnlsq_lss_h12(2, nsetp-1, npp1, d->m, d->a[j], &up, d->a[jj]);
      }
    if(nsetp!=d->m) for(int mi=npp1; mi<d->m; mi++) d->a[j][mi]=0.;
    d->w[j]=0.;
    
    /* Solve the triangular system; store the solution temporarily in Z[] */
    for(int mi=0; mi<nsetp; mi++) {
      int ip=nsetp-(mi+1);
      if(mi!=0) for(int ii=0; ii<=ip; ii++) d->zz[ii]-=d->a[jj][ii]*d->zz[ip+1];
      jj=d->index[ip]; d->zz[ip]/=d->a[jj][ip];
    }
 
    /* Secondary loop begins here */
    while(++iter<itermax) {
      /* See if all new constrained coefficients are feasible; if not, compute alpha */
      double alpha=2.0;
      for(int ip=0; ip<nsetp; ip++) {
        int ni=d->index[ip];
        if(d->zz[ip]<=0.) {
          double t=-d->x[ni]/(d->zz[ip]-d->x[ni]);
          if(alpha>t) {alpha=t; jj=ip-1;}
        }
      }
 
      /* If all new constrained coefficients are feasible then still alpha==2.
         If so, then exit from the secondary loop to main loop */
      if(alpha==2.0) break;
 
      /* Use alpha (0.<alpha<1.) to interpolate between old X and new ZZ */
      for(int ip=0; ip<nsetp; ip++) {
        int ni=d->index[ip]; d->x[ni]+=alpha*(d->zz[ip]-d->x[ni]);
      }
 
      /* Modify A and B and the INDEX arrays to move coefficient i from set P to set Z. */
      int pfeas=1;
      int k=d->index[jj+1];
      do {
        d->x[k]=0.;
        if(jj!=(nsetp-1)) {
          jj++;
          for(int ni=jj+1; ni<nsetp; ni++) {
            int ii=d->index[ni]; d->index[ni-1]=ii;
            double ss, cc;
            nnlsq_lss_g1(d->a[ii][ni-1], d->a[ii][ni], &cc, &ss, &d->a[ii][ni-1]);
            d->a[ii][ni]=0.0;
            for(int nj=0; nj<d->n; nj++) if(nj!=ii) {
              /* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
              double temp=d->a[nj][ni-1];
              d->a[nj][ni-1]=cc*temp+ss*d->a[nj][ni];
              d->a[nj][ni]=-ss*temp+cc*d->a[nj][ni];
            }
            /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
            double temp=d->b[ni-1];
            d->b[ni-1]=cc*temp+ss*d->b[ni];
            d->b[ni]=-ss*temp+cc*d->b[ni];
          }
        }
        npp1=nsetp-1; nsetp--; iz1--; d->index[iz1]=k;
 
        /* See if the remaining coefficients in set P are feasible; they should be because of 
           the way alpha was determined. If any are infeasible it is due to round-off error. 
           Any that are non-positive will be set to zero and moved from set P to set Z. */
        pfeas=1;
        for(jj=0; jj<nsetp; jj++) {
          k=d->index[jj]; if(d->x[k]<=0.) {pfeas=0; break;}
        }
      } while(pfeas==0);
 
      /* Copy B[] into zz[], then solve again and loop back */
      for(int mi=0; mi<d->m; mi++) d->zz[mi]=d->b[mi];
      for(int mi=0; mi<nsetp; mi++) {
        int ip=nsetp-(mi+1);
        if(mi!=0) for(int ii=0; ii<=ip; ii++) d->zz[ii]-=d->a[jj][ii]*d->zz[ip+1];
        jj=d->index[ip]; d->zz[ip]/=d->a[jj][ip];
      }
    } /* end of secondary loop */
 
    if(iter>=itermax) break;
    for(int ip=0; ip<nsetp; ip++) {int k=d->index[ip]; d->x[k]=d->zz[ip];}
  } /* end of main loop */
 
  if(npp1>=d->m) for(int ni=0; ni<d->n; ni++) d->w[ni]=0.;
 
  /* Compute the norm of the final residual vector (sum-of-squares) */
  d->rnorm=0.0;
  for(int mi=npp1; mi<d->m; mi++) d->rnorm+=(d->b[mi]*d->b[mi]);
 
  d->iternr=iter;
  if(verbose>2) printf("  %d iterations.\n", iter);
  if(iter>=itermax) {
    if(verbose>1) printf("  max iterations reached.\n");
    return(1);
  }
  return(0);
} /* nnlsq */

nnlsqDataAllocate()

int nnlsqDataAllocate ( NNLSQDATA * d, const int n, const int m )

extern

Allocate memory for NNLSQDATA (and set data pointers inside the structure). Any previous contents are deleted. Contents are set to NaN or zero.

See also

nnlsDataInit, nnlsqDataFree, nnlsq, nnlsqWght

Returns

Returns TPCERROR status (0 when successful).

Parameters

d	Pointer to initiated structure; any old contents are deleted.
n	Number of parameters (columns of matrix A).
m	Number of samples (rows of matrix A).

Definition at line 411 of file nnlsq.c.

  {
  if(d==NULL) return(TPCERROR_FAIL);
  /* Delete any previous contents */
  nnlsqDataFree(d);
  /* If no memory is requested, then return fail */
  if(n<1 || m<1) return(TPCERROR_FAIL);
 
  /* Allocate memory for all double arrays and matrix */
  int s = n*m + m + n + n + m;
#ifdef TEST_BUFOVERFLOW
  s+=5;
#endif
  d->_data=(double*)malloc(sizeof(double)*s);
  if(d->_data==NULL) return(TPCERROR_OUT_OF_MEMORY);
  for(int i=0; i<s; i++) d->_data[i]=nan("");
  /* Allocate memory for matrix pointers */
  d->a=(double**)malloc(sizeof(double*)*n);
  if(d->a==NULL) {free(d->_data); return(TPCERROR_OUT_OF_MEMORY);}
  /* Set up matrix a and double vectors */
  for(int i=0; i<n; i++) d->a[i] = d->_data + i*m;
  d->b =  d->_data + n*m;
  d->x =  d->_data + n*m + m;
  d->w =  d->_data + n*m + m + n;
  d->zz = d->_data + n*m + m + n + n;
#ifdef TEST_BUFOVERFLOW
  d->b++; d->x+=2; d->w+=3; d->zz+=4;
#endif
 
  /* Allocate memory for integer array */
#ifdef TEST_BUFOVERFLOW
  d->index=(int*)malloc(sizeof(int)*(1+n));
  d->index[n] = 999;
#else
  d->index=(int*)malloc(sizeof(int)*n);
#endif
  if(d->index==NULL) {free(d->_data); free(d->a); return(TPCERROR_OUT_OF_MEMORY);}
 
  /* Set data sizes */
  d->n=n;
  d->m=m;
 
  return(TPCERROR_OK);
}

nnlsqDataFree()

void nnlsqDataFree ( NNLSQDATA * d )

extern

Free memory allocated for NNLSQDATA data. All contents are destroyed.

Precondition

Before first use initialize the structure with nnlsqDataInit().

See also

nnlsqDataInit, nnlsqDataAllocate

Author

Vesa Oikonen

Parameters

d	Pointer to structure.

Definition at line 378 of file nnlsq.c.

  {
  if(d==NULL) return;
  if(d->n<1) return;
  if(d->m<1) return;
 
#ifdef TEST_BUFOVERFLOW
  fprintf(stderr, "testing if buffer overflow happened\n"); fflush(stderr);
  if(d->index[d->n]!=999) fprintf(stderr, "BUFFER OVERFLOW 1\n"); 
  if(!isnan(d->_data[d->m*d->n])) fprintf(stderr, "BUFFER OVERFLOW 2\n"); 
  if(!isnan(d->b[d->m])) fprintf(stderr, "BUFFER OVERFLOW 3\n"); 
  if(!isnan(d->x[d->n])) fprintf(stderr, "BUFFER OVERFLOW 4\n"); 
  if(!isnan(d->w[d->n])) fprintf(stderr, "BUFFER OVERFLOW 5\n"); 
  if(!isnan(d->zz[d->m])) fprintf(stderr, "BUFFER OVERFLOW 6\n"); 
  fprintf(stderr, "tested buffer overflow\n"); fflush(stderr);
#endif
 
  free(d->a);
  free(d->index);
  free(d->_data);
  // then set everything to zero or NULL again
  nnlsqDataInit(d);
}

Referenced by nnlsqDataAllocate().

nnlsqDataInit()

void nnlsqDataInit ( NNLSQDATA * d )

extern

Initiate the NNLSQDATA structure before any use.

See also

nnlsqDataFree, nnlsqDataAllocate, nnlsq, nnlsqWght

Author

Vesa Oikonen

Parameters

d	Pointer to structure to initiate before use.

Definition at line 357 of file nnlsq.c.

  {
  if(d==NULL) return;
  d->n = d->m = 0;
  d->a = NULL;
  d->b = d->x = d->w = d->zz = d->_data = NULL;
  d->index = NULL;
  d->rnorm = 0.0;
  d->iternr = 0;
  d->depf = 0.01; // default
}

Referenced by nnlsqDataFree().

nnlsqWght()

int nnlsqWght ( NNLSQDATA * d, double * weight )

extern

Algorithm for weighting the problem that is given to NNLS-algorithm.

Square roots of weights are used because in NNLS the difference w*A-w*b is squared.

Returns

Algorithm returns zero if successful, 1 if arguments are inappropriate.

See also

nnlsqWghtSquared, nnlsq, llsqWght

Parameters

d	Pointer to filled data structure.
weight	Weights for each sample (array of length M).

Definition at line 470 of file nnlsq.c.

  {
  if(d==NULL || d->n<1 || d->m<1 || d->a==NULL || d->b==NULL || weight==NULL) return(1);
 
  /* Check that weights are not zero and get the square roots of them to w[] */
  double w[d->m];
  for(int mi=0; mi<d->m; mi++) {
    if(weight[mi]<=1.0e-20) w[mi]=0.0;
    else w[mi]=sqrt(weight[mi]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(int mi=0; mi<d->m; mi++) {
    for(int ni=0; ni<d->n; ni++) {
      d->a[ni][mi]*=w[mi];
    }
    d->b[mi]*=w[mi];
  }
 
  return(0);
}

nnlsqWghtSquared()

int nnlsqWghtSquared ( NNLSQDATA * d, double * sweight )

extern

Algorithm for weighting the problem that is given to NNLS-algorithm.

Square roots of weights are used because in NNLS the difference w*A-w*b is squared. Here user must give squared weights; this makes calculation faster, when this function needs to be called many times.

Returns

Algorithm returns zero if successful, 1 if arguments are inappropriate.

See also

nnlsqWght, nnlsq

Parameters

d	Pointer to filled data structure.
sweight	Squared weights for each sample (array of length d->m).

Definition at line 506 of file nnlsq.c.

  {
  if(d==NULL || d->n<1 || d->m<1 || d->a==NULL || d->b==NULL || sweight==NULL) return(1);
 
  /* Multiply rows of matrix A and elements of vector b with weights */
  for(int mi=0; mi<d->m; mi++) {
    for(int ni=0; ni<d->n; ni++) {
      d->a[ni][mi]*=sweight[mi];
    }
    d->b[mi]*=sweight[mi];
  }
 
  return(0);
}

nnlsWght()

int nnlsWght ( int N, int M, double ** A, double * b, double * weight )

extern

Algorithm for weighting the problem that is given to nnls-algorithm.

Square roots of weights are used because in nnls the difference w*A-w*b is squared.

Returns

Algorithm returns zero if successful, 1 if arguments are inappropriate.

See also

nnlsWghtSquared, nnls, qrLSQ, llsqWght

Parameters

N	NNLS dimension N (nr of parameters).
M	NNLS dimension M (nr of samples).
A	NNLS matrix A.
b	NNLS vector B.
weight	Weights for each sample (array of length M).

Definition at line 259 of file nnls.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || A==NULL || b==NULL || weight==NULL) return(1);
 
  /* Allocate memory */
  w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
 
  /* Check that weights are not zero and get the square roots of them to w[] */
  for(m=0; m<M; m++) {
    if(weight[m]<=1.0e-20) w[m]=0.0;
    else w[m]=sqrt(weight[m]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) {
      A[n][m]*=w[m];
    }
    b[m]*=w[m];
  }
 
  free(w);
  return(0);
}

Referenced by spectralDExp(), and spectralDMSurge().

nnlsWghtSquared()

int nnlsWghtSquared ( int N, int M, double ** A, double * b, double * sweight )

extern

Algorithm for weighting the problem that is given to nnls-algorithm.

Square roots of weights are used because in nnls the difference w*A-w*b is squared. Here user must give squared weights; this makes calculation faster, when this function needs to be called many times.

Returns

Algorithm returns zero if successful, 1 if arguments are inappropriate.

See also

nnlsWghtSquared, qrLSQ, nnls

Parameters

N	NNLS dimension N (nr of parameters).
M	NNLS dimension M (nr of samples).
A	NNLS matrix A.
b	NNLS vector B.
sweight	Squared weights for each sample (array of length M).

Definition at line 308 of file nnls.c.

  {
  int n, m;
 
  /* Check the arguments */
  if(N<1 || M<1 || A==NULL || b==NULL || sweight==NULL) return(1);
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) {
      A[n][m]*=sweight[m];
    }
    b[m]*=sweight[m];
  }
 
  return(0);
}

qr()

int qr ( double ** A, int m, int n, double * B, double * X, double * rnorm, double * tau, double * res, double ** wws, double * ws )

extern

Algorithm QR.

Solves a matrix form least square problem min||A x - b|| => A x =~ b (A is m*n matrix, m>=n) using the QR decomposition for overdetermined systems. Based on GNU Scientific Library, edited by Kaisa Liukko.

Instead of pointers for working space, NULL can be given to let this function to allocate and free the required memory.

See also

qrLSQ, nnls, qr_decomp

Returns

Function returns 0 if successful and 1 in case of invalid problem dimensions or memory allocation error.

Parameters

A	On entry, a[m][n] contains the m by n matrix A. On exit, a[][] contains the QR factorization.
m	Dimensions of matrix A are a[m][n].
n	Dimensions of matrix A are a[m][n].
B	B[] is an m-size vector containing the right-hand side vector b.
X	On exit, x[] will contain the solution vector x (size of n).
rnorm	On exit, rnorm (pointer to double) contains the squared Euclidean norm of the residual vector (R^2); enter NULL if not needed.
tau	On exit, tau[] will contain the householder coefficients (size of n); enter NULL, if not needed.
res	An m-size array of working space, res[]. On output contains residual b - Ax. Enter NULL to let qr() to handle it.
wws	m*n array of working space. Enter NULL to let qr() to handle it.
ws	2m-array of working space. Enter NULL to let qr() to handle it.

Definition at line 411 of file qrlsq.c.

  {
  int i;
  double *qrRes, *qrTau, **qrWws, *qrWs, *chain;
 
  /* Check the parameters and data */
  if(m<=0 || n<=0 || A==NULL || B==NULL || X==NULL) return(1);
  if(m<n) return(1);
 
  /* Allocate memory for working space, if required */
  if(tau!=NULL) qrTau=tau; else qrTau=(double*)calloc(n, sizeof(double));
  if(res!=NULL) qrRes=res; else qrRes=(double*)calloc(m, sizeof(double));
  if(wws!=NULL) {
    qrWws=wws; chain=(double*)NULL;
  } else {
    qrWws=(double**)malloc(m * sizeof(double*));
    chain=(double*)malloc(m*n * sizeof(double));
    for(i=0; i<m; i++) qrWws[i]=chain + i*n;
  }
  if(ws!=NULL) qrWs=ws; else qrWs=(double*)calloc(2*m, sizeof(double));
  if(qrTau==NULL || qrRes==NULL || qrWws==NULL || qrWs==NULL) return(1);
 
  /* Form the householder decomposition and solve the least square problem */
  if(qr_decomp(A, m, n, qrTau, qrWws, qrWs)) return(2);
  if(qr_solve(A, m, n, qrTau, B, X, qrRes, rnorm, qrWws, qrWs)) return(3);
 
  /* Free working space, if it was allocated here */
  //if(tau==NULL || res==NULL || wws==NULL || ws==NULL) printf("free in qr()\n");
  if(tau==NULL) free(qrTau);
  if(res==NULL) free(qrRes);
  if(wws==NULL) {free(qrWws); free(chain);}
  if(ws==NULL) free(qrWs);
  for(i=0; i<n; i++) if(isnan(X[i])) return(4);
  return(0);
} /* qr */

qr_decomp()

int qr_decomp ( double ** a, int M, int N, double * tau, double ** cchain, double * chain )

extern

Factorise a general M x N matrix A into A = Q R , where Q is orthogonal (M x M) and R is upper triangular (M x N).

Q is stored as a packed set of Householder vectors in the strict lower triangular part of the input matrix A and a set of coefficients in vector tau.

R is stored in the diagonal and upper triangle of the input matrix.

The full matrix for Q can be obtained as the product Q = Q_1 Q_2 .. Q_k and it's transform as the product Q^T = Q_k .. Q_2 Q_1 , where k = min(M,N) and Q_i = (I - tau_i * h_i * h_i^T) and where h_i is a Householder vector h_i = [1, A(i+1,i), A(i+2,i), ... , A(M,i)]. This storage scheme is the same as in LAPACK.

NOTICE! The calling program must take care that pointer tau is of size N or M, whichever is smaller.

See also

qr

Returns

Function returns 0 if ok.

Parameters

a	contains coefficient matrix A (m*n) as input and factorisation QR as output.
M	nr of rows in matrix A.
N	nr of columns in matrix A.
tau	Vector for householder coefficients, of length N or M, whichever is smaller.
cchain	m*n array of working space.
chain	m size array of working space.

Definition at line 493 of file qrlsq.c.

  {
  //printf("qr_decomp()\n");
  int i, m, n, MNmin;
  double *subvector, **submatrix;
 
  /* Local variables */
  if(M<N) MNmin=M; else MNmin=N;
  if(MNmin<1 || a==NULL || tau==NULL || cchain==NULL || chain==NULL) return(1);
  subvector=chain;
  submatrix=cchain;
 
  for(i=0; i<MNmin; i++) {
    //printf("i=%d (MNmin=%d)\n", i, MNmin);
    /* Compute the Householder transformation to reduce the j-th column of the matrix A to 
       a multiple of the j-th unit vector. Householder vector h_i is saved in the lower triangular 
       part of the column and Householder coefficient tau_i in the vector tau. */
    for(m=i; m<M; m++) {
      //printf("subvector[%d]=a[%d][%d]\n", m-1, m, i);
      subvector[m-i]=a[m][i];
    }
    tau[i] = householder_transform(subvector, M-i);
    for(m=i; m<M; m++) a[m][i]=subvector[m-i];
 
    /* Apply the transformation to the remaining columns to get upper triangular part of matrix R  */
    if(i+1 < N) {
      printf("     apply transformation\n");
      for(m=i; m<M; m++) for(n=i+1; n<N; n++) {
        if((m-i)<0 || (m-i)>=M || (n-i-1)<0 || (n-i-1)>=N) printf("OVERFLOW!\n"); 
        submatrix[m-i][n-i-1]=a[m][n];
      }
      if(householder_hm(tau[i], subvector, submatrix, M-i, N-i)) return(2);
      for(m=i; m<M; m++) for(n=i+1; n<N; n++)
        a[m][n]=submatrix[m-i][n-i-1];
    }
  }
 
  return(0);
}

Referenced by qr().

qr_solve()

int qr_solve ( double ** QR, int M, int N, double * tau, double * b, double * x, double * residual, double * resNorm, double ** cchain, double * chain )

extern

Find the least squares solution to the overdetermined system

A x = b

for m >= n using the QR factorisation A = Q R. qr_decomp() must be used prior to this function in order to form the QR factorisation of A. Solution is formed in the following order: QR x = b => R x = Q^T b => x = R^-1 (Q^T b)

NOTICE! The calling program must take care that pointers b, x and residual are of the right size.

See also

qr, qr_decomp

Returns

Function returns 0 if ok.

Parameters

QR	m*n matrix containing householder vectors of A.
M	nr of rows in matrix A.
N	nr of columns in matrix A.
tau	vector containing householder coefficients tau; length is M or N, whichever is smaller.
b	Contains m-size vector b of A x = b.
x	solution vector x of length n.
residual	residual vector of length m.
resNorm	norm^2 of the residual vector; enter NULL if not needed.
cchain	m*n array of the working space.
chain	2m length array of the working space.

Definition at line 563 of file qrlsq.c.

  {
  //printf("qr_solve()\n");
  if(QR==NULL || tau==NULL || b==NULL || x==NULL || residual==NULL) return(1);
  if(cchain==NULL || chain==NULL) return(1);
  if(M<1 || N<1) return(2);
 
  int m, n;
  double **Rmatrix;
 
  /* Local variable */
  Rmatrix=cchain;
 
  /* Get matrix R from the upper triangular part of QR
     First the rows N - M-1 are eliminated from R matrix*/
  for(m=0; m<N; m++) for(n=0; n<N; n++) Rmatrix[m][n]=QR[m][n];
  for(m=0; m<M; m++) residual[m]=b[m];
 
  /* Compute b = Q^T b */
  if(qr_QTvec(M, N, QR, tau, residual, chain)) return(3);
 
  /* Solve R x = b by computing x = R^-1 b */
  for(n=0; n<N; n++) x[n]=residual[n];
  if(qr_invRvec(N, Rmatrix, x)) return(4);
 
  /* Compute residual = b - A x = Q (Q^T b - R x) */
  for(n=0; n<N; n++) residual[n]=0.0;
  /* Compute residual= Q*residual */
  if(qr_Qvec(M, N, QR, tau, residual, chain)) return(4);
 
  /* Compute norm^2 of the residual vector, if needed */
  if(resNorm!=NULL)
    for(m=0, *resNorm=0.0; m<M; m++) *resNorm +=residual[m]*residual[m];
 
  return(0);
}

Referenced by qr().

qrLH()

int qrLH ( const unsigned int m, const unsigned int n, double * a, double * b, double * x, double * r2 )

extern

Solve over-determined least-squares problem A x ~ b using successive Householder rotations.

This routine is based on the text and Fortran code in C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, New Jersey, 1974, and Fortran code by R.L. Parker and P.B. Stark.

Returns

Returns 0 when successful, 1 if system is singular, and 2 in case of other errors, including that system is under-determined.

Parameters

m	Number of samples in matrix A and the length of vector b.
n	Number of parameters in matrix A and the length of vector x. The n must be smaller or equal to m.
a	Pointer to matrix A; matrix must be given as an n*m array, containing n consecutive m-length vectors. Contents of A are modified in this routine.
b	Pointer to vector b of length n. Contents of b are modified in this routine.
x	Pointer to the result vector x of length n.
r2	Pointer to a double value, in where the sum of squared residuals is written.

Definition at line 946 of file qrlsq.c.

  {
  /* Check the input */
  if(a==NULL || b==NULL || x==NULL || r2==NULL) return(2);
  if(n<1 || m<n) {*r2=nan(""); return(2);}
 
  /* Initiate output to zeroes, in case of exit because of singularity */
  for(unsigned int ni=0; ni<n; ni++) x[ni]=0.0;
  *r2=0.0;
 
  /* Rotates matrix A into upper triangular form */
  for(unsigned int ni=0; ni<n; ni++) {
    /* Find constants for rotation and diagonal entry */
    double sq=0.0;
    for(unsigned int mi=ni; mi<m; mi++) sq+=a[mi + ni*m]*a[mi + ni*m];
    if(sq==0.0) return(1);
    double qv1=-copysign(sqrt(sq), a[ni + ni*m]);
    double u1=a[ni + ni*m] - qv1;
    a[ni + ni*m]=qv1;
    unsigned int ni1=ni+1;
    /*  Rotate the remaining columns of sub-matrix. */
    for(unsigned int nj=ni1; nj<n; nj++) {
      double dot=u1*a[ni + nj*m];
      for(unsigned int mi=ni1; mi<m; mi++)
        dot+=a[mi + nj*m] * a[mi + ni*m];
      double c=dot/fabs(qv1*u1);
      for(unsigned int mi=ni1; mi<m; mi++)
        a[mi + nj*m]-=c*a[mi + ni*m];
      a[ni + nj*m]-=c*u1;
    }
    /* Rotate vector B */
    double dot=u1*b[ni];
    for(unsigned int mi=ni1; mi<m; mi++)
      dot+=b[mi]*a[mi + ni*m];
    double c=dot/fabs(qv1*u1);
    b[ni]-=c*u1;
    for(unsigned int mi=ni1; mi<m; mi++)
      b[mi]-=c*a[mi + ni*m];
  } // end of rotation loop
 
  /* Solve triangular system by back-substitution. */
  for(unsigned int ni=0; ni<n; ni++) {
    int k=n-ni-1;
    double s=b[k];
    for(unsigned int nj=k+1; nj<n; nj++) 
      s-=a[k + nj*m] * x[nj];
    if(a[k + k*m]==0.0) return(1);
    x[k]=s/a[k + k*m];
  }
 
  /* Calculate the sum of squared residuals. */
  *r2=0.0;
  for(unsigned int mi=n; mi<m; mi++)
    *r2 += b[mi]*b[mi];
 
  return(0);
}

Referenced by bvls().

qrLSQ()

int qrLSQ ( double ** mat, double * rhs, double * sol, const unsigned int rows, const unsigned int cols, double * r2 )

extern

QR least-squares solving routine.

Solves a matrix form least square problem min||A x - b|| => A x =~ b (A is m*n matrix, m>=n) using the QR decomposition for overdetermined systems

Author

: Michael Mazack License: Public Domain. Redistribution and modification without restriction is granted. If you find this code helpful, please let the author know (https://mazack.org). Small editions by Vesa Oikonen when added to tpcclib.

See also

nnls, qr

Todo

Check that modification in one subroutine is ok, and add possibility to provide working space for the main and sub-functions.

Returns

Returns 0 if ok.

Parameters

mat	Pointer to the row-major matrix (2D array), A[rows][cols]; not modified.
rhs	Vector b[rows]; modified to contain the computed (fitted) b[], that is, matrix-vector product A*x.
sol	Solution vector x[cols]; solution x corresponds to the solution of both the modified and original systems A and B.
rows	Number of rows (samples) in matrix A and length of vector B.
cols	Number of cols (parameters) in matrix A and length of vector X.
r2	Pointer to value where R^2, the difference between original and computed right hand side (b); enter NULL if not needed.

Definition at line 72 of file qrlsq.c.

  {
  if(rows<1 || cols<1) return(1);
  if(mat==NULL || rhs==NULL || sol==NULL) return(2);
 
  unsigned int i, j, max_loc;
  unsigned int *p=NULL;
  double *buf=NULL, *ptr, *v=NULL, *orig_b=NULL, **A=NULL;
 
  /* Allocate memory for index vector and Householder transform vector */
  /* and original data matrix and vector. */
  p=malloc(sizeof(unsigned int)*cols);
  A=malloc(sizeof(double*)*rows);
  buf=malloc(sizeof(double)*(rows*cols+rows+rows));
  if(p==NULL || A==NULL || buf==NULL) {
    free(p); free(A); free(buf); 
    return(3);
  }
  ptr=buf; 
  for(i=0; i<rows; i++) {A[i]=ptr; ptr+=cols;}
  v=ptr; ptr+=rows; orig_b=ptr;
  /* copy the original data */
  for(i=0; i<rows; i++) {
    orig_b[i]=rhs[i];
    for(j=0; j<cols; j++) A[i][j]=mat[i][j];
  }
 
  /* Initial permutation vector. */
  for(i=0; i<cols; i++) p[i]=i;
  
  /* Apply rotators to make R and Q'*b */
  for(i=0; i<cols; i++) {
    if(qr_get_next_col(A, rows, cols, i, p, &max_loc)==0)
      qr_swap_cols(p, i, max_loc);
    qr_householder(A, rows, cols, i, p[i], v);
    qr_apply_householder(A, rhs, rows, cols, v, i, p);
  }
 
  /* Back solve Rx = Q'*b */
  qr_back_solve(A, rhs, rows, cols, sol, p);
 
  /* Compute fitted b[] (matrix-vector product A*x) using non-modified A */
  for(i=0; i<rows; i++) {
    rhs[i]=0.0;
    for(j=0; j<cols; j++) rhs[i]+=mat[i][j]*sol[j];
  }
 
  /* Compute R^2, if requested */
  if(r2!=NULL) {
    double ss=0, d;
    for(i=0; i<rows; i++) {
      d=orig_b[i]-rhs[i];
      ss+=d*d;
    }
    *r2=ss;
  }
 
  /* Collect garbage. */
  free(p); free(A); free(buf);
  return(0);
}

qrSimpleLSQ()

int qrSimpleLSQ ( double ** A, double * B, int M, int N, double * X, double * r2 )

extern

Simplistic QR least-squares solving routine.

Reference: Máté A. Introduction to Numerical Analysis with C programs. Chapter 38. Overdetermined systems of linear equations. Brooklyn College of the City University of New York, 2014.

See also

qrLSQ

Returns

Function returns 0 when successful.

Parameters

A	Matrix A[M][N].
B	Vector B[M].
M	nr of rows in matrix A, must be >=N.
N	nr of columns in matrix A, must be <=M.
X	Vector X of length N for the results.
r2	Pointer to value where R^2, the difference between original and computed right hand side (b); enter NULL if not needed.

Definition at line 854 of file qrlsq.c.

  {
  //printf("qrSimpleLSQ()\n");
  if(A==NULL || B==NULL || X==NULL || N<1 || M<N) return(1);
 
  const double closeZero=1.0E-20;
 
  double cc[M], orig_B[M];
  double p[M][N], c[M][N];
 
  /* Householder transformation of matrix A into matrix P */
  for(int n=0; n<N; n++) {
    for(int nn=0; nn<n; nn++) p[nn][n]=0.0;
    p[n][n]=0.0;
    for(int m=n; m<M; m++) p[n][n]+=A[m][n]*A[m][n];
    p[n][n]=copysign(sqrt(p[n][n]), A[n][n]);
    p[n][n]+=A[n][n];
    if(fabs(p[n][n])<closeZero) return(2);
    for(int m=n+1; m<M; m++) p[m][n]=A[m][n];
    double norm=0.0;
    for(int m=n; m<M; m++) norm+=p[m][n]*p[m][n];
    norm=sqrt(norm);
    for(int m=n; m<M; m++) p[m][n]/=norm;
    for(int m=n; m<M; m++) {
      for(int nn=n; nn<N; nn++) {
        c[m][nn]=A[m][nn];
        for(int mm=n; mm<M; mm++) c[m][nn]-=2.0*p[m][n]*p[mm][n]*A[mm][nn];
      }
    }
    for(int m=n; m<M; m++)
      for(int nn=n; nn<N; nn++) A[m][nn]=c[m][nn];
  }
 
  /* Keep the original B[] for R^2 calculation */
  for(int m=0; m<M; m++) orig_B[m]=B[m];
 
  /* Compute P'b */
  for(int n=0; n<N; n++) {
    for(int m=n; m<M; m++) {
      cc[m]=B[m];
      for(int mm=n; mm<M; mm++) cc[m]-=2.0*p[m][n]*p[mm][n]*B[mm];
    }
    for(int m=n; m<M; m++) B[m]=cc[m];
  }
 
  /* Solve the linear system with backward substitution */
  X[N-1]=B[N-1]/A[N-1][N-1];
  for(int n=N-2; n>=0; n--) {
    X[n]=B[n];
    for(int nn=n+1; nn<N; nn++) X[n]-=A[n][nn]*X[nn];
    X[n]/=A[n][n];
  }
 
  /* Compute R^2, if requested */
  if(r2!=NULL) {
    double ss=0, d;
    for(int m=0; m<M; m++) {
      d=orig_B[m]-B[m];
      ss+=d*d;
    }
    *r2=ss;
  }
 
  return(0);
}

qrWeight()

int qrWeight ( int N, int M, double ** A, double * b, double * weight, double * ws )

extern

Algorithm for weighting the problem that is given to QR algorithm. Square roots of weights are used because in QR the difference w*A-w*b is squared.

See also

qrWeightRm, qrLSQ, qr

Todo

Add tests.

Returns

Algorithm returns zero if successful, otherwise <>0.

Parameters

N	Dimensions of matrix A are a[m][n].
M	Dimensions of matrix A are a[m][n]; size of vector B is m.
A	Matrix a[m][n] for QR, contents will be weighted here.
b	B[] is an m-size vector for QR, contents will be weighted here.
weight	Pointer to array of size m, which contains sample weights, used here to weight matrix A and vector b.
ws	m-sized vector for working space; enter NULL to allocate locally.

Definition at line 745 of file qrlsq.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || A==NULL || b==NULL || weight==NULL) return(1);
 
  /* Allocate memory, if necessary */
  if(ws==NULL) {
    w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
  } else {
    w=ws;
  }
 
  /* Check that weights are not zero and get the square roots of them to w[]. */
  for(m=0; m<M; m++) {
    if(weight[m]<=1.0e-100) w[m]=1.0e-50;
    else w[m]=sqrt(weight[m]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights. */
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) A[m][n]*=w[m];
    b[m]*=w[m];
  }
 
  if(ws==NULL) free(w);
  return(0);
}

qrWeightRm()

int qrWeightRm ( int N, int M, double ** A, double * b, double * weight, double * ws )

extern

Remove the weights from data provided to QR LSQ algorithms.

See also

qrWeight, qrLSQ, qr

Todo

Add tests.

Returns

Algorithm returns zero if successful, otherwise <>0.

Parameters

N	Dimensions of matrix A are a[m][n].
M	Dimensions of matrix A are a[m][n]; size of vector B is m.
A	Matrix a[m][n] for QR, weights will be removed here; enter NULL if not needed.
b	B[] is an m-size vector for QR, weights will be removed here; enter NULL if not needed here.
weight	Pointer to array of size m, which contains sample weights, used here to weight matrix A and vector b.
ws	m-sized vector for working space; enter NULL to allocate locally

Definition at line 796 of file qrlsq.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || weight==NULL) return(1);
  if(A==NULL && b==NULL) return(0); // nothing to do
 
  /* Allocate memory, if necessary */
  if(ws==NULL) {
    w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
  } else {
    w=ws;
  }
 
  /* Check that weights are not zero and get the square roots of them to w[] */
  for(m=0; m<M; m++) {
    if(weight[m]<=1.0e-100) w[m]=1.0e-50;
    else w[m]=sqrt(weight[m]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights */
  for(m=0; m<M; m++) {
    if(A!=NULL) for(n=0; n<N; n++) A[m][n]/=w[m];
    if(b!=NULL) b[m]/=w[m];
  }
 
  if(ws==NULL) free(w);
  return(0);
}

rootsCubic()

int rootsCubic ( const double a, const double b, const double c, double * x1, double * x2, double * x3 )

extern

Find the real roots of cubic equation x^3 + a x^2 + b x + c = 0.

The number of real roots can be 1 or 3, but coincident roots are possible. Root values, excluding duplicates and triplicates, will be placed in pointers x1, x2, and x3 in ascending order, or set to NaN.

Based on GNU Scientific Library function gsl_poly_solve_cubic.

See also

rootsQuadratic, fitLine

Returns

Returns the number of real individual roots, and zero in case of an error.

Parameters

a	Input: Parameter a of the cubic equation.
b	Input: Parameter b of the cubic equation.
c	Input: Parameter c of the cubic equation.
x1	Output: Pointer for root 1 value. Will be set to NaN if there are no real roots.
x2	Output: Pointer for root 2 value. Will be set to NaN if there is only one individual real root.
x3	Output: Pointer for root 3 value. Will be set to NaN if there are only two individual real roots.

Definition at line 30 of file roots.c.

  {
  if(x1==NULL || x2==NULL || x3==NULL) return(0);
  *x1=*x2=*x3=nan("");
 
  if(a==0 && b==0 && c==0) return(0);
  if(!isfinite(a) || !isfinite(b) || !isfinite(c)) return(0);
 
  long double q=(a*a - 3*b);
  long double r=(2*a*a*a - 9*a*b + 27*c);
 
  long double Q=q/9;
  long double R=r/54;
 
  long double Q3=Q*Q*Q;
  long double R2=R*R;
 
  long double CR2=729*r*r;
  long double CQ3=2916*q*q*q;
 
  if(R==0 && Q==0) {
    *x1=-a/3; // *x2=-a/3; *x3=-a/3;
    return(1);
  } else if(CR2==CQ3) {
    long double sqrtQ=sqrt(Q);
    if(R>0) {
      *x1=-2*sqrtQ - a/3;
      *x2=sqrtQ - a/3; // *x3=sqrtQ - a/3;
    } else {
      *x1=-sqrtQ - a/3; // *x2=-sqrtQ - a/3; *x3=2*sqrtQ - a/3;
      *x2=2*sqrtQ - a/3;
    }
    return(2);
  } else if(R2<Q3) {
    double ratio=copysign(sqrt(R2/Q3), R);
    //if(!isfinite(ratio)) return(0);
    double theta=acos(ratio);
    double norm=-2*sqrt(Q);
    //if(!isfinite(norm)) return(0);
    *x1=norm*cos(theta/3) - a/3;
    *x2=norm*cos((theta + 2.0*M_PI)/3) - a/3;
    *x3=norm*cos((theta - 2.0*M_PI)/3) - a/3;
    if(*x1>*x2) {double s=*x1; *x1=*x2; *x2=s;}
    if(*x2>*x3) {
      double s=*x2; *x2=*x3; *x3=s;
      if(*x1>*x2) {double s=*x1; *x1=*x2; *x2=s;}
    }
    return(3);
  } else {
    long double A = -copysign(pow(fabsl(R) + sqrt(R2-Q3), 1.0/3.0), R);
    long double B=Q/A;
    *x1=A + B - a/3;
    return(1);
  }
}

rootsQuadratic()

int rootsQuadratic ( const double a, const double b, const double c, double * x1, double * x2 )

extern

Find the real roots of quadratic equation a x^2 + b x + c = 0.

The number of real roots can be 0, 1, or 2, but coincident roots are possible. Root values, excluding duplicate value, will be placed in pointers x1 and x2 in ascending order, or set to NaN.

Based on GNU Scientific Library function gsl_poly_solve_quadratic.

See also

rootsCubic, fitLine

Returns

Returns the number of real roots, excluding duplicate solutions.

Parameters

a	Input: Parameter a of the quadratic equation.
b	Input: Parameter b of the quadratic equation.
c	Input: Parameter c of the quadratic equation.
x1	Output: Pointer for root 1 value. Will be set to NaN if there are no real roots.
x2	Output: Pointer for root 2 value. Will be set to NaN if there is only one individual real root.

Definition at line 111 of file roots.c.

  {
  if(x1==NULL || x2==NULL) return(0);
  *x1=*x2=nan("");
 
  if(!isfinite(a) || !isfinite(b) || !isfinite(c)) return(0);
 
  if(a==0) {
    if(b==0) return(0);
    else {*x1=-c/b; /* *x2=-c/b;*/ return(1);}
  }
  long double d=b*b - 4*a*c;
  if(d>0) {
    if(b==0) {
      long double r=fabs(sqrt(-c/a));
      *x1=-r; *x2=r;
    } else {
      long double t=-0.5*(b + copysign(sqrt(d), b));
      long double r1=t/a; 
      long double r2=c/t;
      if(r1<r2) {*x1=r1; *x2=r2;} else {*x1=r2; *x2=r1;}
    }
    return(2);
  } else if(d==0) {
    *x1=-0.5*b/a; // *x2=-0.5*b/a;
    return(1);
  } else {
    return(0);
  }
}