関連情報
ホーム > 製品 > NAG数値計算ライブラリ > サンプルソースコード集 > 最小二乗双3次スプライン曲面フィット (C言語/C++)

最小二乗双3次スプライン曲面フィット

C言語によるサンプルソースコード
使用関数名:nag_2d_spline_fit_panel (e02dac)

Keyword: 最小二乗, 双3次, スプライン, フィット

概要

本サンプルは最小二乗双3次スプライン曲面フィットを行うC言語によるサンプルプログラムです。 本サンプルは以下に示されるデータについて最小二乗双3次スプライン曲面フィットを行い、スプラインの値を求めて出力します。

双3次スプラインのデータ 

※本サンプルはNAG Cライブラリに含まれる関数 nag_2d_spline_fit_panel() のExampleコードです。本サンプル及び関数の詳細情報は nag_2d_spline_fit_panel のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

(本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照)

このデータをダウンロード
nag_2d_spline_fit_panel (e02dac) Example Program Data
0.000001
  30
   8
  10
   -0.52    0.60    0.93   10.
   -0.61   -0.95   -1.79   10.
    0.93    0.87    0.36   10.
    0.09    0.84    0.52   10.
    0.88    0.17    0.49   10.
   -0.70   -0.87   -1.76   10.
    1.00    1.00    0.33    1.
    1.00    0.10    0.48    1.
    0.30    0.24    0.65    1.
   -0.77   -0.77   -1.82    1.
   -0.23    0.32    0.92    1.
   -1.00    1.00    1.00    1.
   -0.26   -0.63    8.88    1.
   -0.83   -0.66   -2.01    1.
    0.22    0.93    0.47    1.
    0.89    0.15    0.49    1.
   -0.80    0.99    0.84    1.
   -0.88   -0.54   -2.42    1.
    0.68    0.44    0.47    1.
   -0.14   -0.72    7.15    1.
    0.67    0.63    0.44    1.
   -0.90   -0.40   -3.34    1.
   -0.84    0.20    2.78    1.
    0.84    0.43    0.44    1.
    0.15    0.28    0.70    1.
   -0.91   -0.24   -6.52    1.
   -0.35    0.86    0.66    1.
   -0.16   -0.41    2.32    1.
   -0.35   -0.05    1.66    1.
   -1.00   -1.00   -1.00    1.
  -0.5
   0.0

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に一次方程式の有効ランクを決定するための閾値(eps)を指定しています。
  • 3行目にデータ点の数(m)を指定しています。
  • 4行目に変数xのノット数の合計(px)を指定しています。
  • 5行目に変数yのノット数の合計(py)を指定しています。
  • 6〜35行目にデータ点yの座標、xの座標、fの座標、重み(w)を指定しています。
  • 36〜37行目に変数yに関する内部ノット(mu)を指定しています。

出力結果

(本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照)

この出力例をダウンロード
nag_2d_spline_fit_panel (e02dac) Example Program Results

Interior y -knots
    -0.5000
     0.0000

Interior x -knots
None

Sum of squares of residual RHS       1.467e+01

Rank   22

x and y have been interchanged

      X          Y          Data       Fit      Residual
    -0.9500    -0.6100    -1.7900    -1.7931  -3.126e-03
    -0.8700    -0.7000    -1.7600    -1.7521   7.893e-03
    -0.7700    -0.7700    -1.8200    -2.4301  -6.101e-01
    -0.6300    -0.2600     8.8800     7.6346  -1.245e+00
    -0.6600    -0.8300    -2.0100    -1.5815   4.285e-01
    -0.5400    -0.8800    -2.4200    -2.6795  -2.595e-01
    -0.7200    -0.1400     7.1500     7.5708   4.208e-01
    -1.0000    -1.0000    -1.0000    -1.0228  -2.277e-02
    -0.4000    -0.9000    -3.3400    -4.6955  -1.356e+00
    -0.2400    -0.9100    -6.5200    -4.7072   1.813e+00
    -0.4100    -0.1600     2.3200     2.7039   3.839e-01
    -0.0500    -0.3500     1.6600     2.2865   6.265e-01
     0.6000    -0.5200     0.9300     0.9441   1.407e-02
     0.8700     0.9300     0.3600     0.3529  -7.097e-03
     0.8400     0.0900     0.5200     0.5024  -1.761e-02
     0.1700     0.8800     0.4900     0.4705  -1.951e-02
     1.0000     1.0000     0.3300     0.6315   3.015e-01
     0.1000     1.0000     0.4800     1.4910   1.011e+00
     0.2400     0.3000     0.6500     0.9241   2.741e-01
     0.3200    -0.2300     0.9200    -0.3692  -1.289e+00
     1.0000    -1.0000     1.0000     1.0835   8.352e-02
     0.9300     0.2200     0.4700     1.4912   1.021e+00
     0.1500     0.8900     0.4900     0.4414  -4.861e-02
     0.9900    -0.8000     0.8400     0.5495  -2.905e-01
     0.4400     0.6800     0.4700     1.5862   1.116e+00
     0.6300     0.6700     0.4400     0.6288   1.888e-01
     0.2000    -0.8400     2.7800     1.7123  -1.068e+00
     0.4300     0.8400     0.4400     0.6888   2.488e-01
     0.2800     0.1500     0.7000     0.7713   7.134e-02
     0.8600    -0.3500     0.6600     0.9347   2.747e-01

Sum of squared residuals       1.467e+01

Spline coefficients
    -1.0228   115.4668  -433.5558   -68.1973
    24.8426  -140.1485   258.5042    15.6756
   -29.4878   132.2933  -173.5103    20.0983
     9.9575   -51.6200    67.6666    -5.8765
    10.0577     4.7543   -15.3533    -0.3260
     1.0835    -2.7932     7.7708     0.6315

  • 3〜5行目に変数yの内部ノットが出力されています。
  • 7〜8行目に変数xの内部ノットがないことが示されています。
  • 10行目に残差平方和が出力されています。
  • 12行目にランクが出力されています。
  • 14行目にxとyが入れ替えられたことが示されています。
  • 16〜46行目にx、y、f、フィットされた値、残差が出力されています。
  • 48行目に残差平方和が出力されています。
  • 50〜56行目にスプライン係数が出力されています。

ソースコード

(本関数の詳細はnag_2d_spline_fit_panel のマニュアルページを参照)

※本サンプルソースコードはNAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法


このソースコードをダウンロード
/* nag_2d_spline_fit_panel  (e02dac) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage02.h>

int main(void)
{
  /* Initialized data */
  char label[] = "xy";

  /* Scalars */
  double d, eps, sigma, sum, temp;
  Integer exit_status = 0, i, iadres, itemp, j, m, nc, np, npoint, px, py,
         rank;

  /* Arrays */
  double *dl = 0, *f = 0, *ff = 0, *lamda = 0, *mu = 0, *w = 0, *x = 0;
  double *y = 0;
  Integer *point = 0;

  /* Nag Types */
  Nag_2dSpline spline;
  NagError fail;

  exit_status = 0;
  INIT_FAIL(fail);

  /* Initialize spline */
  spline.lamda = 0;
  spline.mu = 0;
  spline.c = 0;

  printf("nag_2d_spline_fit_panel (e02dac) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  while (scanf("%lf", &eps) != EOF && exit_status == 0)
  {
    /* Read data, interchanging X and Y axes if PX.LT.PY */
    scanf("%ld%*[^\n] ", &m);
    if (m > 1) {
      /* Allocate memory */
      if (!(f = NAG_ALLOC(m, double)) ||
          !(ff = NAG_ALLOC(m, double)) ||
          !(w = NAG_ALLOC(m, double)) ||
          !(x = NAG_ALLOC(m, double)) || !(y = NAG_ALLOC(m, double))
             )
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
    }
    else {
      printf("Invalid m.\n");
      exit_status = 1;
      goto END;
    }
    scanf("%ld%ld%*[^\n] ", &px, &py);
    if (px < 8 && py < 8) {
      printf("px or py is too small.\n");
      exit_status = 1;
      goto END;
    }
    nc = (px - 4) * (py - 4);
    np = (px - 7) * (py - 7);
    npoint = m + (px - 7) * (py - 7);

    /* Allocate memory */
    if (!(dl = NAG_ALLOC(nc, double)) ||
        !(point = NAG_ALLOC(npoint, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    if (px < py) {
      itemp = px;
      px = py;
      py = itemp;
      itemp = 1;
      /* Allocate memory */
      if (!(lamda = NAG_ALLOC(px, double)) || !(mu = NAG_ALLOC(py, double))
             )
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }

      for (i = 0; i < m; ++i)
        scanf("%lf%lf%lf%lf", &y[i], &x[i], &f[i], &w[i]);
      scanf("%*[^\n] ");

      if (py > 8) {
        for (j = 4; j < py - 4; ++j)
          scanf("%lf", &mu[j]);
        scanf("%*[^\n] ");
      }
      if (px > 8) {
        for (j = 4; j < px - 4; ++j)
          scanf("%lf", &lamda[j]);
        scanf("%*[^\n] ");
      }
    }
    else {
      /* Allocate memory */
      if (!(lamda = NAG_ALLOC(px, double)) || !(mu = NAG_ALLOC(py, double))
             )
      {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      itemp = 0;
      for (i = 0; i < m; ++i)
        scanf("%lf%lf%lf%lf", &x[i], &y[i], &f[i], &w[i]);
      scanf("%*[^\n] ");

      if (px > 8) {
        for (j = 4; j < px - 4; ++j)
          scanf("%lf", &lamda[j]);
        scanf("%*[^\n] ");
      }
      if (py > 8) {
        for (j = 4; j < py - 4; ++j)
          scanf("%lf", &mu[j]);
        scanf("%*[^\n] ");
      }
    }

    printf("\nInterior %1.1s -knots\n", label + itemp);
    for (j = 4; j < px - 4; ++j)
      printf("%11.4f\n", lamda[j]);
    if (px == 8)
      printf("None\n");

    printf("\nInterior %1.1s -knots\n", label + (2 - itemp - 1));
    for (j = 4; j < py - 4; ++j)
      printf("%1s%11.4f\n", "", mu[j]);
    if (py == 8)
      printf("None\n");

    /* nag_2d_panel_sort (e02zac).
     * Sort two-dimensional data into panels for fitting bicubic
     * splines
     */
    nag_2d_panel_sort(px, py, lamda, mu, m, x, y, point, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_2d_panel_sort (e02zac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    /* Fit bicubic spline to data points */
    spline.nx = px;
    spline.ny = py;

    if (!(spline.c = NAG_ALLOC((spline.nx - 4) * (spline.ny - 4), double)) ||
        !(spline.lamda = NAG_ALLOC(spline.nx, double)) ||
        !(spline.mu = NAG_ALLOC(spline.ny, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    for (i = 0; i < spline.nx; i++)
      spline.lamda[i] = lamda[i];
    for (i = 0; i < spline.ny; i++)
      spline.mu[i] = mu[i];

    /* nag_2d_spline_fit_panel (e02dac).
     * Least squares surface fit, bicubic splines
     */
    nag_2d_spline_fit_panel(m, x, y, f, w, point, dl, eps, &sigma, &rank,
                            &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_2d_spline_fit_panel (e02dac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\nSum of squares of residual RHS%16.3e\n", sigma);
    printf("\nRank%5ld\n", rank);

    /* nag_2d_spline_eval (e02dec).
     * Evaluation of bicubic spline, at a set of points
     */
    nag_2d_spline_eval(m, x, y, ff, &spline, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_2d_spline_eval (e02dec).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    sum = 0.;
    if (itemp == 1)
      printf("\nx and y have been interchanged\n\n");

    /*Output data points, fitted values and residuals */
    printf("      X          Y          Data       Fit      Residual\n");
    for (i = 0; i < np; ++i) {
      iadres = i + m;
      while ((iadres = point[iadres] - 1) >= 0) {
        temp = ff[iadres] - f[iadres];

        printf("%11.4f%11.4f%11.4f%11.4f%12.3e\n", x[iadres],
               y[iadres], f[iadres], ff[iadres], temp);
        /* Computing 2nd power */
        d = temp * w[iadres];
        sum += d * d;
      }
    }

    printf("\nSum of squared residuals%16.3e\n", sum);
    printf("\nSpline coefficients\n");
    for (i = 0; i < px - 4; ++i) {
      for (j = 0; j < py - 4; ++j)
        printf("%11.4f", spline.c[i * (py - 4) + j]);
      printf("\n");
    }
  END:
    NAG_FREE(dl);
    NAG_FREE(f);
    NAG_FREE(ff);
    NAG_FREE(lamda);
    NAG_FREE(mu);
    NAG_FREE(w);
    NAG_FREE(x);
    NAG_FREE(y);
    NAG_FREE(point);
    NAG_FREE(spline.lamda);
    NAG_FREE(spline.mu);
    NAG_FREE(spline.c);
  }
  return exit_status;
}


Results matter. Trust NAG.

Privacy Policy | Trademarks