fitpack-parcur

      subroutine parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,
     * nc,c,fp,wrk,lwrk,iwrk,ier)
c     在idim维空间中给定m个有序点集x(i),并且给定一个对应的严格递增的u(i)和
c     正数w(i)的集合，i=1,2，…m，子程序parcur确定光滑近似样条曲线s(u)，即。
c      x1 = s1(u)
c      x2 = s2(u)       ub <= u <= ue
c      .........
c      xidim = sidim(u)
c     sj(u),j=1,2,...,idim是knots t(j),j=1,2,...,n的k阶普通样条函数
c     如果ipar=1，则值ub，ue和u(i),i=1,2，…，m必须由用户设置。
c     如果ipar=0，这些值将由parcur自动选择
c      v(1) = 0
c      v(i) = v(i-1) + dist(x(i),x(i-1)) ,i=2,3,...,m
c      u(i) = v(i)/v(m) ,i=1,2,...,m
c      ub = u(1) = 0, ue = u(m) = 1.
c     如果iopt=-1，parcur根据一组给定的knots。计算加权最小二乘样条曲线
c     如果iopt>=0，则样条的结点数sj(u)和位置t(j),j=1,2,...,n由例程自动
c     选择。通过最小化s(u)的k阶导数在结点t(j),j=k+2,k+3,...,n-k-1处最小
c     化不连续来实现光滑。平滑度的大小取决于条件
c     f(p)=sum((w(i)*dist(x(i),s(u(i))))**2) be <= s,s给定的非负常数，
c     称为平滑因子。
c     拟合s(u)用b样条表示，可以是用子程序曲线求值。
c
c  calling sequence:
c     call parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,
c    * fp,wrk,lwrk,iwrk,ier)
c
c  parameters:
c   iopt  : 整数。iopt必须指定是否加权最小二乘样条曲线（iopt=-1）或平滑样条
c           曲线（iopt=0或1）必须确定。如果iopt=0，则例程从一组初始结点开始
c           ，t(i)=ub,t(i+k+1)=ue,i=1,2,...,k+1.如果iopt=1，例程将继续执
c           行在最后一次表演时发现结。
c           注意：iopt=1的调用必须是立即的,在iopt=1或iopt=0的另一个调用之前。
c           退出时不变。
c   ipar  : 整数。在条目ipar必须指定是否（ipar=1）用户将提供参数值u(i)，ub和
c           ue或者（ipar=0）这些值是否由parcur。退出时不变。
c   idim  : 整数。在条目中必须指定curve的维度0 < idim < 11。退出时不变。
c   m     : 整数。在条目m上必须指定数据点的个数。m > k，退出时不变。
c   u     : 数组。在ipar=1的情况下，尺寸至少为(m)的实数组。表项u(i)必须设置为
c           该参数的第i个值变量u =1,2，…，m。这些值必须是严格按升序提供，在退
c           出不会改变。如果ipar=0，在退出时，数组u将包含值u(i)由parcur决定。
c   mx    : 整数。在条目mx上必须指定的在调用（子）程序中声明的数组x的实际尺寸不
c           要太小（见x）。退出时不变。
c   x     : 数组。尺寸至少为idim*m的实数组。在输入前，x(idim*(i-1)+j)必须包含
c           第i个数据点第j个坐标i=1,2,...,m and j=1,2,...,idim退出时不变。
c   w     : 数组。维数至少为(m)的实数组。在进入前，w(i)必须设置为权重集合中的第
c           i个值。w(i)必须是严格正的。退出时不变。
c   ub,ue : 浮点数。在输入时（在ipar=1的情况下），ub、ue必须包含参数u的上界和下
c           界。 ub <=u(1), ue>= u(m)。如果ipar = 0，这些值将自动设置为0和1。
c   k     : 整数。必须指定样条的阶k。建议使用三次样条k=3（1<=k<=5）；强烈建议用
c           户不要同时选择k偶数,同时s值很小。退出时不变。
c   s     : 浮点。在输入上（在iopt>=0的情况下）必须指定平滑因子s >=0。退出不变。
c   nest  : 整数。在输入时必须包含一个高估计返回的样条结点总数nest，表示例程可用
c           的存储空间nest >=2*k+2, 在大多数实际情况下，nest=m/2就足够了。足够
c           大的nest=m+k+1，结点的数量需要插值（s=0）。退出时不变。
c   n     : 整数。除非ier = 10（在iopt >=0的情况下），否则n将包含返回的光滑样条
c           曲线的结点总数。如果使用计算模式iopt=1，则此值为n在随后的调用之间应该
c           保持不变。如果iopt=-1，则必须在条目时指定n的值。
c   t     : 数组。维度至少nest。成功退出后，该数组将包含样条曲线。内部结点的位置
c           t(k+2),t(k+3),..,t(n-k-1)以及附加的位置t(1)=t(2)=...=t(k+1)=ub 
c           and t(n-k)=...=t(n)=ue b样条表示。
c           若采用计算模式iopt=1，则t(1),t(2),...,t(n)的值应保持不变.如果计算模
c           式为iopt=-1，则t(k+2),...,t(n-k-1)必须由用户提供。
c   nc    : 整数。在输入必须指定的，在调用（子）程序中声明的数组c的实际尺寸nc。nc
c           不能太小（见c）。退出时不变。
c   c     : 数组。实际数组的维度至少（nest*idim）。成功退出后，该数组将包含系数在
c           样条曲线s(u)的b样条表示中，即将给出样条sj(u)的b样条系数在 
c           c(n*(j-1)+i),i=1,2,...,n-k-1 对于j=1,2,...,idim.
c   fp    : 浮点。除非ier = 10，否则fp包含返回样条曲线的残差的平方和。
c   wrk   : 维数至少为m*(k+1)+nest*(6+idim+3*k)的实数组。用作工作空间。如果计算
c           模式iopt=1时，应保持值wrk(1),...,wrk(n)不变
c   lwrk  : 整数。在输入时，lwrk必须指定的在调用（子）程序中声明的数组wrk的实际尺
c           寸，lwrk不能太小（见wrk）。退出时不变。
c   iwrk  : 正数数组，维数至少为（nest）的整数数组。用作工作空间。如果计算模式
c           iopt=1为时，值iwrk(1),...,iwrk(n) 应保持不变。
c   ier   : 整数。除非例程检测到错误，否则它包含退出时的非正值，即
c     ier=0  : 正常返回。返回曲线的残差和为对fp，使abs(fp-s)/s <= tol，程
c           序将tol的相对公差设置为0.001。
c     ier=-1 : 正常返回。返回的曲线是插值样条曲线（fp=0）。
c     ier=-2 : 正常返回。返回的曲线k次的加权最小平方多项式曲线,是极端情况下fp
c              给出平滑因子s的上界fp0。
c     ier=1  : 错误。存储空间需求超过可用存储空间容量，由参数nest指定。可能原因：
c              nest太小。如果nest已经足够大（例如，nest > m/2），它也可能表示s是
c              太小了。返回近似最小二乘样条knots t(1),t(2),...,t(n). (n=nest)
c              参数fp给出了相应的加权残差的平方和（fp>s）
c     ier=2  : 错误。寻找光滑样条曲线的迭代过程中,发现了一个理论上不可能的结果fp=s，
c              可能原因：s太小。近似返回，但是对应的残差的加权平方和不满足条件
c              abs(fp-s)/s < tol。
c     ier=3  : 错误。最大迭代次数maxit（设置为20）允许找到平滑曲线fp=s；可能原因：
c             s太小近似返回，但是对应的残差的加权平方和不满足条件abs(fp-s)/s<tol。
c     ier=10 : 错误。对输入数据的有效性进行控制必须满足以下限制条件。
c              -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
c              0<=ipar<=1, 0<idim<=10, lwrk>=(k+1)*m+nest*(6+idim+3*k),
c              nc>=nest*idim
c              ipar=0: sum j=1,idim (x(idim*i+j)-x(idim*(i-1)+j))**2>0
c                      i=1,2,...,m-1.
c              ipar=1: ub<=u(1)<u(2)<...<u(m)<=ue
c              iopt=-1: 2*k+2<=n<=min(nest,m+k+1)  
c                       ub<t(k+2)<t(k+3)<...<t(n-k-1)<ue
c                       (ub=0 and ue=1 in case ipar=0)
c                       舍恩伯格-惠特尼条件必须是数据点uu(j)的子集
c                       t(j) < uu(j) < t(j+k+1), j=1,2,...,n-k-1
c              iopt>=0: s>=0 
c                           s=0 : nest >= m+k+1
c              发现违反下列条件之一的，予以控制立即被重新传递给调用程序。没有返回近似.
c  further comments:
c   通过参数s，用户可以权衡在拟合的接近度和拟合的平滑度。如果s太大，曲线就会太平滑，
c   信号就会丢失;如果s太小，曲线会吸收太多的噪声。在极端情况下s=0，程序将返回插值曲
c   线，s非常大, k阶的最小二乘多项式曲线。在这两个极端之间，会有一个正确选择的结果,
c   拟合的紧密度和拟合的平滑度之间很好的折衷。强烈建议配合图形化检查.
c   s的推荐值取决于权重w(i)。如果这些是取1/d(i), d(i)为标准差的估计值x(i)，一个好
c   的s值应该在(m-sqrt(2*m),m+sqrt(2*m))。如果对x(i)中的统计误差一无所知,每个w(i)可以
c   设为1，s由实验和误差确定，考虑到上面的注释。最好的办法就是从一个非常大的s（来确定最小
c   二乘,对于s多项式曲线和上界fp0）开始，然后逐步减小s的值(比如在开始，即s=fp0/10, 
c   fp0/100，…而且更小心近似曲线显示更多细节)以获得更接近的拟合。为了节省寻找一个好的s值
c   的程序提供不同的计算模式。在例程的第一次调用时，或者当他想用初始的knots重新开始时，用
c   户必须设置iopt=0。
c   如果iopt=1，程序将继续使用在例程的最后一次调用处找到的knots集。这将节省大量对于不同的
c   s值重复调用parcur时的计算时间。返回数据背后的曲线的样条结点的数目和它们的位置将取决于s
c   的值和形状的复杂程度。但是，如果计算模式iopt=1是时，返回的结点也可能取决于先前的s值（如
c   果这些较小）。因此，如果用不同的s值和iopt=1计算模式进行试验，用户最终可以接受一个的结果
c   ，也许值得调用parcur再次对s的选定值进行并行处理，但现在使用iopt=0。实际上，parcur可以
c   返回近似相同质量的结果但knots更少，因此更好的数据减少也是用户的一个重要目标。
c   近似曲线的形式会受到以下因素的强烈影响参数值u(i)的选择。如果没有物理选择某一特定参数u的
c   原因，将通过选择parpar（在ipar=0的情况下）往往获得效果良好，即
c        v(1)=0, v(i)=v(i-1)+q(i), i=2,...,m, u(i)=v(i)/v(m), i=1,..,m
c   where
c        q(i)= sqrt(sum j=1,idim (xj(i)-xj(i-1))**2 )
c   other possibilities for q(i) are
c        q(i)= sum j=1,idim (xj(i)-xj(i-1))**2
c        q(i)= sum j=1,idim abs(xj(i)-xj(i-1))
c        q(i)= max j=1,idim abs(xj(i)-xj(i-1))
c        q(i)= 1
c
c  other subroutines required:
c    fpback,fpbspl,fpchec,fppara,fpdisc,fpgivs,fpknot,fprati,fprota
c
c  references:
c   dierckx p. : algorithms for smoothing data with periodic and
c                parametric splines, computer graphics and image
c                processing 20 (1982) 171-184.
c   dierckx p. : algorithms for smoothing data with periodic and param-
c                etric splines, report tw55, dept. computer science,
c                k.u.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 ub,ue,s,fp
      integer iopt,ipar,idim,m,mx,k,nest,n,nc,lwrk,ier
c  ..array arguments..
      real*8 u(m),x(mx),w(m),t(nest),c(nc),wrk(lwrk)
      integer iwrk(nest)
c  ..local scalars..
      real*8 tol,dist
      integer i,ia,ib,ifp,ig,iq,iz,i1,i2,j,k1,k2,lwest,maxit,nmin,ncc
c ..function references
      real*8 sqrt
c  ..
c  we set up the parameters tol and maxit
      maxit = 20
      tol = 0.1e-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(iopt.lt.(-1) .or. iopt.gt.1) go to 90
      if(ipar.lt.0 .or. ipar.gt.1) go to 90
      if(idim.le.0 .or. idim.gt.10) go to 90
      if(k.le.0 .or. k.gt.5) go to 90
      k1 = k+1
      k2 = k1+1
      nmin = 2*k1
      if(m.lt.k1 .or. nest.lt.nmin) go to 90
      ncc = nest*idim
      if(mx.lt.m*idim .or. nc.lt.ncc) go to 90
      lwest = m*k1+nest*(6+idim+3*k)
      if(lwrk.lt.lwest) go to 90
      if(ipar.ne.0 .or. iopt.gt.0) go to 40
      i1 = 0
      i2 = idim
      u(1) = 0.
      do 20 i=2,m
         dist = 0.
         do 10 j=1,idim
            i1 = i1+1
            i2 = i2+1
            dist = dist+(x(i2)-x(i1))**2
  10     continue
         u(i) = u(i-1)+sqrt(dist)
  20  continue
      if(u(m).le.0.) go to 90
      do 30 i=2,m
         u(i) = u(i)/u(m)
  30  continue
      ub = 0.
      ue = 1.
      u(m) = ue
  40  if(ub.gt.u(1) .or. ue.lt.u(m) .or. w(1).le.0.) go to 90
      do 50 i=2,m
         if(u(i-1).ge.u(i) .or. w(i).le.0.) go to 90
  50  continue
      if(iopt.ge.0) go to 70
      if(n.lt.nmin .or. n.gt.nest) go to 90
      j = n
      do 60 i=1,k1
         t(i) = ub
         t(j) = ue
         j = j-1
  60  continue
      call fpchec(u,m,t,n,k,ier)
      if (ier.eq.0) go to 80
      go to 90
  70  if(s.lt.0.) go to 90
      if(s.eq.0. .and. nest.lt.(m+k1)) go to 90
      ier = 0
c we partition the working space and determine the spline curve.
  80  ifp = 1
      iz = ifp+nest
      ia = iz+ncc
      ib = ia+nest*k1
      ig = ib+nest*k2
      iq = ig+nest*k2
      call fppara(iopt,idim,m,u,mx,x,w,ub,ue,k,s,nest,tol,maxit,k1,k2,
     * n,t,ncc,c,fp,wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),
     * iwrk,ier)
  90  return
      end

 void FitPack::parcur(int iopt, int ipar, int idim, int m, double* u, int mx, double* x, double* w, double* ub, double* ue, int k, double* s, int nest, int n,
     int nc, double* c, double* fp, double* wrk, int lwrk, int* iwrk, int* ier) {

    double tol, dist;
    int i, ia, ib, ifp, ig, iq, iz, i1, i2, j, k1, k2, lwest, maxit, nmin, ncc;
    double** a, ** b, ** g, ** q;
    double* t, * fpint, * z;
    //we set up the parameters tol and maxit
    maxit = 20;
    tol = 0.1e-02;
    /*
    * 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 (iopt < -1 || iopt > 1) return;
    if (ipar < 0 || ipar > 1) return;
    if (idim <= 0 || idim > 10) return;
    if (k <= 0 || k > 5) return;
    k1 = k + 1;
    k2 = k1 + 1;
    nmin = 2 * k1;
    if (m < k1 || nest < nmin) return;
    ncc = nest * idim;
    if (mx < m * idim || nc < ncc) return;
    lwest = m * k1 + nest * (6 + idim + 3 * k);
    if (lwrk < lwest) return;
    if (ipar != 0 || iopt > 0)
    {
        if (*ub > u[0] || *ue < u[m - 1] || w[0] <= 0.0) return;
        for (i = 1; i < m; i++) {
            if (u[i - 1] >= u[i] || w[i] <= 0.0) return;
        }
    }

    i1 = 0;
    i2 = idim;
    u[0] = 0.0;
    for (i = 1; i < m; i++) {
        dist = 0.0;
        for (j = 0; j < idim; j++) {
            i1++;
            i2++;
            dist += (x[i2] - x[i1]) * (x[i2] - x[i1]);
        }
        u[i] = u[i - 1] + sqrt(dist);
    }
    if (u[m - 1] <= 0.0) return;
    for (i = 1; i < m; i++) {
        u[i] = u[i] / u[m - 1];
    }
    *ub = 0.0;
    *ue = 1.0;
    u[m - 1] = *ue;
    if (iopt >= 0)
    {
        if (*s < 0) return;
        if (*s == 0.0 && nest < (m + k1)) return;
        *ier = 0;
    }
    if (n < nmin || n > nest) return;
    t = new double[nest];
    j = n;
    for (i = 0; i < k1; i++) {
        t[i] = *ub;
        t[j] = *ue;
        j--;
    }
    fpchec(u, m, t, n, k, ier);
    if (*ier == 0)
    {
        /*we partition the working space and determine the spline curve.*/
        ifp = 0;
        iz = ifp + nest;
        ia = iz + ncc;
        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];
        fpint = new double[nest];
        z = 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];
            q[i] = &wrk[iq + i * k2];
            fpint[i] = 0.0;
            z[i] = 0.0;
        }
        fppara(iopt, idim, m, u, mx, x, w, *ub, *ue, k, *s, nest, tol, maxit,
            k1, k2, n, t, ncc, c, fp, fpint, z, a, b, g, q, iwrk, ier);
        delete[] z;
        delete[] fpint;
        delete[] q;
        delete[] g;
        delete[] b;
        delete[] a;
        delete[] t;
    }
    return;
 }