fitpack-curfit

subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp, wrk,lwrk,iwrk,ier)
c     iopt  :整数；不变
c           iopt=-1 通过给定的knots集和，计算加权最小二乘样条
c           iopt=0  光滑样条，初始化点集t(i)=xb, t(i+k+1)=xe, i=1,2,...k+1
c           iopt=1 样条曲线s(x)的knots个数和t(j),j=1,2,...,n的位置自动计算;
c                  通过最小化s(x)的k阶导数在结点t(j),j=k+2, k+3,…,n-k-1处的
c                  不连续跳变来实现s(x)的平滑性。平滑度的大小取决于条件
c                  f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, s给定的非负
c                  常数，称为平滑因子。拟合s(x)用b样条表示（b样条系数-
c                  c(j),j=1,2,...,n-k-1）可通过splev均值评估
c     m     :整数；不变
c           点个数；
c     x、y  :浮点数；不变
c           (x(i),y(i)), i=1,2,...,m; 点集；其中x(i)必须升序
c     w     :浮点数；不变
c           w(i),i=1,2,...,m; 正权重；
c     xb、xe:浮点数；不变
c           xb <= x <= xe点取值区间
c     k     :整数；不变；
c           1<=k<=5推荐；cubic splines (k=3)；k为奇数
c     s     :浮点；不变
c           正数；
c     nest  :整数；
c           nest >=2*k+2
c           nest=m/2一般足够了
c           如果s=0，则nest=m/2
c     n     :整数
c           iopt >=0;ier !=10，n将包含返回的样条近似的总节数。
c           iopt=1;不变。
c           iopt=-1，n的值必须在输入时指定。
c     t     :数组；
c           样条结点，即。内部结点的位置t(k+2),t(k+3)...,t(n-k-1);
c           以及附加结点的位置t(1)=t(2)=…=t(k+1)=xb
c           和t(n-k)=…=t(n)=xe需要b样条表示。
c           iopt=1，则t(1),c t(2)，…，t(n)不变。
c           iopt=-1，则值t(k+2)，…，t（n-k-1）必须在输入之前由用户提供。
c     c     :数组；
c           c(1),c(2),..,c(n-k-1)样条系数
c     fp    :浮点
c           样条近似的平方残差的加权和
c     wrk   :浮点数组
c           iopt=1，wrk(1)，…，wrk(n)不变。
c     lwrk  :整数;不变
c           lwrk必须数组wrk的维度。lwrk不能太小;
c     iwrk  :整数；
c           iopt=1，iwrk(1)，…，iwrk(n)不变。
c     ier   :整数；
c           ier=0       :正常；abs(fp-s)/s <= tol；tol=0.001;fp残差平和
c           ier=-1      :正常；fp=0插值样条
c           ier=-2      :正常；k次的加权最小二乘多项式;平滑因子s的上界f0
c           ier=1       :错误；存储空间超过可用的存储空间，nest太小。
c                       如果nest已经很大（比如nest b> m/2），它也可能表明s太小了。
c                       结点t(1),t(2)，…，t(n)（n=nest）加权最小二乘样条。参数fp
c                       给出相应的c平方残差（fp>s）的加权和。
c           ier=2       :错误;s太小。残差平方和不满足条件abs(fp-s)/s < tol
c           ier=3       :错误。fp=s的平滑样条;：s太小残差平方和不满足条件abs(fp-s)/s < tol。
c           ier=10      :错误。 -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
c                       xb<=x(1)<x(2)<...<x(m)<=xe, lwrk>=(k+1)*m+nest*(7+3*k)
c                       if iopt=-1: 2*k+2<=n<=min(nest,m+k+1) xb<t(k+2)<t(k+3)<...<t(n-k-1)<xe
c                       t(j) < xx(j) < t(j+k+1), j=1,2,...,n-k-1
c                       iopt>=0: s>=0 if s=0 : nest >= m+k+1
c
c  further comments:
c     通过参数s，用户可以控制近似拟合的紧密性和拟合的平滑性之间的权衡。
c     如果s太大，样条会太光滑，信号会丢失；如果s太小，样条将接收太多的噪声。
c     在极端情况下，如果s=0，程序将返回一个插值样条，如果s非常大，则k次的
c     加权最小二乘多项式。在这两个极端之间，适当选择的s将导致c在拟合的紧密
c     性和拟合的平滑性之间的良好折衷。为了确定的近似，对应于某一s是否满足
c     的要求，强烈建议用户以图形的方式检查是否符合。
c     s的推荐值取决于权重w(i)。如果将这些取为1/d(i), d(i)是y(i)的标准差估计，
c     则应该在范围(m-sqrt（2*m),m+c sqrt(2*m)）中找到一个良好的s值。
c     如果对y(i)c中的统计误差一无所知，则考虑到上面的注释，每个w(i)可以设置
c     为等于1和s，由试验和c误差确定。最好的方法是从一个非常大的s值开始（以确
c     定s的最小二乘多项式和相应的上限fp0），然后逐步减小s的值(例如在开始时减
c     少10倍，即s=fp0/10, fp0/100，…当逼近显示出更多的细节时，为了获得更接近
c     的拟合，为了节省寻找好的s值的时间，程序提供了不同的计算模式。在第一次调
c     用例程时，或者每当用户想用初始的节点集重新启动时，用户必须设置iopt=0。
c     如果iopt=1，程序将继续使用在例程最后调用时找到的结集。这将节省大量的计
c     算时间，如果对不同的s值重复调用曲线，返回的样条结的数量和它们的位置将
c     取决于s的值和数据底层函数形状的复杂性。但是，如果使用计算模式iopt=1，
c     则返回的节也可能取决于先前调用的s值（如果这些值较小）。因此，如果经过
c     不同s值和iopt=1的多次试验后，用户最终可以接受一个满意的拟合，那么他可
c     能值得再次调用curfit，使用s的选定值，但现在iopt=0.确实，curfit可能会
c     返回近似相同质量的拟合，但结更少，因此更好，如果数据缩减也是用户的一个重要目标。
c
c  other subroutines required:
c    fpback,fpbspl,fpchec,fpcurf,fpdisc,fpgivs,fpknot,fprati,fprota
c
c  references:
c   dierckx p. : an algorithm for smoothing, differentiation and integ-
c                ration of experimental data using spline functions,
c                j.comp.appl.maths 1 (1975) 165-184.
c   dierckx p. : a fast algorithm for smoothing data on a rectangular
c                grid while using spline functions, siam j.numer.anal.
c                19 (1982) 1286-1304.
c   dierckx p. : an improved algorithm for curve fitting with spline
c                functions, report tw54, dept. computer science,k.u.
c                leuven, 1981.
c   dierckx p. : curve and surface fitting with splines, monographs on
c                numerical analysis, oxford university press, 1993.
c
c  author:
c    p.dierckx
c    dept. computer science, k.u. leuven
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c  creation date : may 1979
c  latest update : march 1987
c
c  ..
c  ..scalar arguments..
      real*8 xb,xe,s,fp
      integer iopt,m,k,nest,n,lwrk,ier
c  ..array arguments..
      real*8 x(m),y(m),w(m),t(nest),c(nest),wrk(lwrk)
      integer iwrk(nest)
c  ..local scalars..
      real*8 tol
      integer i,ia,ib,ifp,ig,iq,iz,j,k1,k2,lwest,maxit,nmin
c  ..
c  we set up the parameters tol and maxit
      maxit = 20
      tol = 0.1d-02
c  before starting computations a data check is made. if the input data
c  are invalid, control is immediately repassed to the calling program.
      ier = 10
      if(k.le.0 .or. k.gt.5) go to 50
      k1 = k+1
      k2 = k1+1
      if(iopt.lt.(-1) .or. iopt.gt.1) go to 50
      nmin = 2*k1
      if(m.lt.k1 .or. nest.lt.nmin) go to 50
      lwest = m*k1+nest*(7+3*k)
      if(lwrk.lt.lwest) go to 50
      if(xb.gt.x(1) .or. xe.lt.x(m)) go to 50
      do 10 i=2,m
         if(x(i-1).gt.x(i)) go to 50
  10  continue
      if(iopt.ge.0) go to 30
      if(n.lt.nmin .or. n.gt.nest) go to 50
      j = n
      do 20 i=1,k1
         t(i) = xb
         t(j) = xe
         j = j-1
  20  continue
      call fpchec(x,m,t,n,k,ier)
      if (ier.eq.0) go to 40
      go to 50
  30  if(s.lt.0.) go to 50
      if(s.eq.0. .and. nest.lt.(m+k1)) go to 50
c we partition the working space and determine the spline approximation.
  40  ifp = 1
      iz = ifp+nest
      ia = iz+nest
      ib = ia+nest*k1
      ig = ib+nest*k2
      iq = ig+nest*k2
      call fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,
     * wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),iwrk,ier)
  50  return
      end

   void FitPack::curfit(int iopt, int m, const double* x, const double* y,
       const double* w, double xb, double xe, int k, double s,
       int nest, int* n, double* t, double* c, double* fp,
       double* wrk, int lwrk, int* iwrk, int* ier)
   {
       
        /*  C++ translation by Werner Stille. */
       int i, ia, ib, ifp, ig, iq, iz, j, k1, k2, lwest, maxit, nmin;
       double tol;
       double** a, ** b, ** g, ** q;
       /*  we set up the parameters tol and maxit */
       maxit = 20;
       tol = 0.001;
       /*  before starting computations a data check is made. if the input data
        are invalid, control is immediately repassed to the calling program. */
       *ier = 10;
       if (k <= 0 || k > 5)
           return;
       k1 = k + 1;
       k2 = k1 + 1;
       if (iopt < -1 || iopt > 1)
           return;
       nmin = k1 << 1;
       if (m < k1 || nest < nmin)
           return;
       lwest = m * k1 + nest * (k * 3 + 7);
       if (lwrk < lwest)
           return;
       if (xb > x[0] || xe < x[m - 1] || w[0] <= 0)
           return;
       for (i = 1; i < m; i++)
           if (x[i - 1] >= x[i] || w[i] <= 0)
               return;
       if (iopt >= 0) {
           if (s < 0)
               return;
           if (s == 0 && nest < m + k1)
               return;
           *ier = 0;
       }
       else {
           if (*n < nmin || *n > nest)
               return;
           j = *n;
           for (i = 0; i < k1; i++) {
               t[i] = xb;
               t[--j] = xe;
           }
           fpchec(x, m, t, *n, k, ier);
           if (*ier)
               return;
       }
       /* we partition the working space and determine the spline approximation. */
       ifp = 0;
       iz = ifp + nest;
       ia = iz + nest;
       ib = ia + nest * k1;
       ig = ib + nest * k2;
       iq = ig + nest * k2;
       a = new double* [nest];
       b = new double* [nest];
       g = new double* [nest];
       q = new double* [nest];
       for (i = 0; i < nest; i++) {
           a[i] = &wrk[ia + i * k1];
           b[i] = &wrk[ib + i * k2];
           g[i] = &wrk[ig + i * k2];
       }
       for (i = 0; i < m; i++)
           q[i] = &wrk[iq + i * k1];
       fpcurf(iopt, x, y, w, m, xb, xe, k, s, nest, tol, maxit, k1, k2, n, t, c,
           fp, &wrk[ifp], &wrk[iz], a, b, g, q, iwrk, ier);
       delete[] q;
       delete[] g;
       delete[] b;
       delete[] a;
       return;
   }