 |
 |
| Commit in sandbox/magfield/src on MAIN | |||
| spline.cpp | +987 | added 3075 | |
| spline.h | +177 | added 3075 | |
| spline1.cpp | +173 | added 3075 | |
| spline1.h | +74 | added 3075 | |
| +1411 | |||
core spline classes to be used for fitting the field map.
--- sandbox/magfield/src/spline.cpp (rev 0) +++ sandbox/magfield/src/spline.cpp 2014-04-03 15:54:52 UTC (rev 3075) @@ -0,0 +1,987 @@
+// Magic Software, Inc.
+// http://www.magic-software.com
+// Copyright (c) 2000, 2001. All Rights Reserved
+//
+// Source code from Magic Software is supplied under the terms of a license
+// agreement and may not be copied or disclosed except in accordance with the
+// terms of that agreement. The various license agreements may be found at
+// the Magic Software web site. This file is subject to the license
+//
+// FREE SOURCE CODE
+// http://www.magic-software.com/License/free.pdf
+
+#include <fstream>
+#include <math.h>
+#include <float.h>
+#include "spline.h"
+//===========================================================================
+const double mgcSpline::INVALID_DATA = FLT_MAX;
+
+// error handling
+int mgcSpline::verbose = 0;
+unsigned mgcSpline::error = 0;
+const unsigned mgcSpline::invalid_dimensions = 0x00000001;
+const unsigned mgcSpline::invalid_degree = 0x00000002;
+const unsigned mgcSpline::invalid_dim = 0x00000004;
+const unsigned mgcSpline::null_data = 0x00000008;
+const unsigned mgcSpline::allocation_failed = 0x00000010;
+const char* mgcSpline::message[5] = {
+ "invalid number of dimensions",
+ "degree must be positive",
+ "dimension must be at least DEGREE+1",
+ "input data pointer is null",
+ "failed to allocate memory"
+};
+//---------------------------------------------------------------------------
+void mgcSpline::
+InitNull ()
+{
+ dimensions = 0;
+ degree = 0;
+ Dp1 = 0;
+ Dp1_to_N = 0;
+ Dp1_to_2N = 0;
+ dim = 0;
+ data = 0;
+ dom_min = 0;
+ dom_max = 0;
+ grid_min = 0;
+ grid_max = 0;
+ base = 0;
+ old_base = 0;
+ cache = 0;
+ inter = 0;
+ poly = 0;
+ coeff = 0;
+ matrix = 0;
+ product = 0;
+ skip = 0;
+}
+//---------------------------------------------------------------------------
+void mgcSpline::
+FreeStorage ()
+{
+ if ( dim == 0 ) // nothing to free
+ return;
+
+ delete[] dim;
+ delete[] base;
+ delete[] old_base;
+ delete[] dom_min;
+ delete[] dom_max;
+ delete[] grid_min;
+ delete[] grid_max;
+ delete[] cache;
+ delete[] inter;
+ delete[] product;
+ delete[] skip;
+
+ int d;
+ for (d = 0; d < dimensions; d++)
+ delete[] poly[d];
+ delete[] poly;
+
+ for (d = 0; d <= degree; d++) {
+ delete[] matrix[d];
+ delete[] coeff[d];
+ }
+ delete[] matrix;
+ delete[] coeff;
+}
+//---------------------------------------------------------------------------
+void mgcSpline::
+Create ()
+{
+ int k, d;
+
+ // setup degree constants
+ Dp1 = degree+1;
+ Dp1_to_N = 1;
+ for (d = 0; d < dimensions; d++)
+ Dp1_to_N *= Dp1;
+ Dp1_to_2N = Dp1_to_N*Dp1_to_N;
+
+ // initialize base indices
+ base = new int[dimensions];
+ if ( base == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ old_base = new int[dimensions];
+ if ( old_base == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (d = 0; d < dimensions; d++)
+ old_base[d] = -1;
+
+ grid_min = new int[dimensions];
+ if ( grid_min == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ grid_max = new int[dimensions];
+ if ( grid_max == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (d = 0; d < dimensions; d++) {
+ grid_min[d] = -1;
+ grid_max[d] = -1;
+ }
+
+ // compute domain [min,max] for B-spline
+ dom_min = new double[dimensions];
+ if ( dom_min == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ dom_max = new double[dimensions];
+ if ( dom_max == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ for (d = 0; d < dimensions; d++) {
+ double tmpmin = 1.0f, tmpmax, dom_sup = double(dim[d]-degree+1);
+ double next = (tmpmin+dom_sup)/2;
+ do {
+ tmpmax = next;
+ next = (next+dom_sup)/2;
+ } while ( next < dom_sup );
+
+ dom_min[d] = tmpmin;
+ dom_max[d] = tmpmax;
+ }
+
+ // cache for optimizing compute_intermediate()
+ cache = new double[Dp1_to_N];
+ if ( cache == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ // storage for intermediate tensor product
+ inter = new double[Dp1_to_N];
+ if ( inter == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ // polynomial allocations
+ poly = new double*[dimensions];
+ if ( poly == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (d = 0; d < dimensions; d++) {
+ poly[d] = new double[Dp1];
+ if ( poly[d] == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ }
+
+ // coefficients for polynomial calculations
+ coeff = new double*[Dp1];
+ if ( coeff == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (int row = 0; row <= degree; row++) {
+ coeff[row] = new double[Dp1];
+ if ( coeff[row] == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (int col = row; col <= degree; col++) {
+ coeff[row][col] = 1;
+ for (d = 0; d <= row-1; d++)
+ coeff[row][col] *= col-d;
+ }
+ }
+
+ // generate spline blending matrix
+ matrix = BlendMatrix(degree);
+
+ // tensor product of m with itself N times
+ product = new double[Dp1_to_2N];
+ if ( product == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ skip = new int[Dp1_to_2N];
+ if ( skip == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ int* coord = new int[2*dimensions]; // for address decoding
+ if ( coord == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (k = 0; k < Dp1_to_2N; k++) {
+ int temp = k;
+ for (d = 0; d < 2*dimensions; d++) {
+ coord[d] = temp % Dp1;
+ temp /= Dp1;
+ }
+ product[k] = 1;
+ for (int s = 0; s < dimensions; s++)
+ product[k] *= matrix[coord[s]][coord[s+dimensions]];
+
+ skip[k] = 1;
+ }
+ delete[] coord;
+
+ // compute increments to skip zero values of mtensor
+ for (k = 0; k < Dp1_to_2N; ) {
+ for (d = k+1; d < Dp1_to_2N && product[d] == 0; d++)
+ skip[k]++;
+ k = d;
+ }
+
+ evaluate_callback = 0;
+}
+//---------------------------------------------------------------------------
+mgcSpline& mgcSpline::
+Renew (int _dimensions, int _degree, const int* _dim, double* _data)
+{
+ int d;
+
+ // blow away old data
+ FreeStorage();
+
+ // get number of dimensions
+ if ( (dimensions = _dimensions) <= 0 ) {
+ Report(invalid_dimensions);
+ return *this;
+ }
+
+ // get degree and check if valid
+ if ( (degree = _degree) <= 0 ) {
+ Report(invalid_degree);
+ return *this;
+ }
+
+ // get dimension bounds and check if valid
+ if ( _dim == 0 ) {
+ Report(invalid_dim);
+ return *this;
+ }
+
+ for (d = 0; d < dimensions; d++)
+ if ( _dim[d] <= degree+1 ) {
+ Report(invalid_dim);
+ return *this;
+ }
+
+ dim = new int[dimensions];
+ if ( dim == 0 ) {
+ Report(allocation_failed);
+ return *this;
+ }
+ for (d = 0; d < dimensions; d++)
+ dim[d] = _dim[d];
+
+ // get spline data
+ if ( (data = _data) == 0 ) {
+ Report(null_data);
+ return *this;
+ }
+
+ Create();
+
+ return *this;
+}
+//---------------------------------------------------------------------------
+mgcSpline::
+mgcSpline (int _dimensions, int _degree, const int* _dim, double* _data)
+{
+ InitNull();
+ Renew(_dimensions,_degree,_dim,_data);
+}
+//---------------------------------------------------------------------------
+void mgcSpline::
+SetPolynomial (int order, double diff, double* poly)
+{
+ double diffpower = 1;
+ for (int d = order; d <= degree; d++) {
+ poly[d] = coeff[order][d]*diffpower;
+ diffpower *= diff;
+ }
+}
+//---------------------------------------------------------------------------
+int mgcSpline::
+Choose (int n, int k)
+{
+ int i;
+
+ // computes combination "n choose k"
+
+ if ( n <= 1 || k >= n )
+ return 1;
+
+ int result = 1;
+ for (i = 0; i < k; i++)
+ result *= n-i;
+ for (i = 1; i <= k; i++)
+ result /= i;
+
+ return result;
+}
+//---------------------------------------------------------------------------
+double** mgcSpline::
+BlendMatrix (int deg)
+{
+ int degp1 = deg+1;
+ int row, col, i, j, d;
+
+ // allocate triple arrays
+ int*** a = new int**[degp1];
+ if ( a == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ int*** b = new int**[degp1];
+ if ( b == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (d = 0; d <= deg; d++) {
+ a[d] = new int*[degp1];
+ if ( a[d] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ b[d] = new int*[degp1];
+ if ( b[d] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (row = 0; row <= deg; row++) {
+ a[d][row] = new int[degp1];
+ if ( a[d][row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ b[d][row] = new int[degp1];
+ if ( b[d][row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (col = 0; col <= deg; col++) {
+ a[d][row][col] = 0;
+ b[d][row][col] = 0;
+ }
+ }
+ }
+
+ a[0][0][0] = 1;
+ b[0][0][0] = 1;
+
+ for (d = 1; d <= deg; d++) {
+ // compute a[]
+ for (row = 0; row <= d; row++) {
+ for (col = 0; col <= d; col++) {
+ a[d][row][col] = 0;
+ if ( col >= 1 ) {
+ a[d][row][col] += a[d-1][row][col-1];
+ if ( row >= 1 )
+ a[d][row][col] -= b[d-1][row-1][col-1];
+ }
+ if ( row >= 1 )
+ a[d][row][col] += (d+1)*b[d-1][row-1][col];
+ }
+ }
+
+ // compute b[]
+ for (row = 0; row <= d; row++) {
+ for (col = 0; col <= d; col++) {
+ b[d][row][col]= 0;
+ for (i = col; i <= d; i++)
+ if ( (i-col) % 2 )
+ b[d][row][col] -= Choose(i,col)*a[d][row][i];
+ else
+ b[d][row][col] += Choose(i,col)*a[d][row][i];
+ }
+ }
+ }
+
+ double** imat = new double*[degp1];
+ if ( imat == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (row = 0; row <= deg; row++) {
+ imat[row] = new double[degp1];
+ if ( imat[row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (col = 0; col <= deg; col++) {
+ imat[row][col]= 0;
+ for (i = col; i <= deg; i++) {
+ int prod = 1;
+ for (j = 1; j <= i-col; j++)
+ prod *= deg-row;
+ imat[row][col] += prod*Choose(i,col)*a[deg][deg-row][i];
+ }
+ }
+ }
+
+ double factorial = 1;
+ for (d = 1; d <= deg; d++)
+ factorial *= d;
+ double** mat = new double*[degp1];
+ if ( mat == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (row = 0; row <= deg; row++) {
+ mat[row] = new double[degp1];
+ if ( mat[row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (col = 0; col <= deg; col++)
+ mat[row][col] = imat[row][col]/factorial;
+ }
+
+ // deallocate triple arrays
+ for (d = 0; d <= deg; d++) {
+ for (row = 0; row <= deg; row++) {
+ delete[] b[d][row];
+ delete[] a[d][row];
+ }
+ delete[] b[d];
+ delete[] a[d];
+ }
+ delete[] b;
+ delete[] a;
+
+ // deallocate integer matrix
+ for (d = 0; d <= deg; d++)
+ delete[] imat[d];
+ delete[] imat;
+
+ return mat;
+}
+//---------------------------------------------------------------------------
+int mgcSpline::
+Number (unsigned single_error)
+{
+ int result;
+ for (result = -1; single_error; single_error >>= 1)
+ result++;
+ return result;
+}
+//---------------------------------------------------------------------------
+void mgcSpline::
+Report (unsigned single_error)
+{
+ if ( mgcSpline::verbose )
+ std::cout << "mgcSpline: " << message[Number(single_error)] << std::endl;
+ else
+ std::ofstream("spline.err",std::ios::out|std::ios::app)
+ << "mgcSpline: " << message[Number(single_error)] << std::endl;
+
+ error |= single_error;
+}
+//---------------------------------------------------------------------------
+void mgcSpline::
+Report (std::ostream& ostr)
+{
+ for (unsigned single_error = 1; single_error; single_error <<= 1)
+ if ( error & single_error )
+ ostr << "mgcSpline: " << message[Number(single_error)] << std::endl;
+
+ error = 0;
+}
+//===========================================================================
+const double mgcSplineD::INVALID_DATA = DBL_MAX;
+
+// error handling
+int mgcSplineD::verbose = 0;
+unsigned mgcSplineD::error = 0;
+const unsigned mgcSplineD::invalid_dimensions = 0x00000001;
+const unsigned mgcSplineD::invalid_degree = 0x00000002;
+const unsigned mgcSplineD::invalid_dim = 0x00000004;
+const unsigned mgcSplineD::null_data = 0x00000008;
+const unsigned mgcSplineD::allocation_failed = 0x00000010;
+const char* mgcSplineD::message[5] = {
+ "invalid number of dimensions",
+ "degree must be positive",
+ "dimension must be at least DEGREE+1",
+ "input data pointer is null",
+ "failed to allocate memory"
+};
+//---------------------------------------------------------------------------
+void mgcSplineD::
+InitNull ()
+{
+ dimensions = 0;
+ degree = 0;
+ Dp1 = 0;
+ Dp1_to_N = 0;
+ Dp1_to_2N = 0;
+ dim = 0;
+ data = 0;
+ dom_min = 0;
+ dom_max = 0;
+ grid_min = 0;
+ grid_max = 0;
+ base = 0;
+ old_base = 0;
+ cache = 0;
+ inter = 0;
+ poly = 0;
+ coeff = 0;
+ matrix = 0;
+ product = 0;
+ skip = 0;
+}
+//---------------------------------------------------------------------------
+void mgcSplineD::
+FreeStorage ()
+{
+ if ( dim == 0 ) // no storage to deallocate
+ return;
+
+ delete[] dim;
+ delete[] base;
+ delete[] old_base;
+ delete[] dom_min;
+ delete[] dom_max;
+ delete[] grid_min;
+ delete[] grid_max;
+ delete[] cache;
+ delete[] inter;
+ delete[] product;
+ delete[] skip;
+
+ int d;
+ for (d = 0; d < dimensions; d++)
+ delete[] poly[d];
+ delete[] poly;
+
+ for (d = 0; d <= degree; d++) {
+ delete[] matrix[d];
+ delete[] coeff[d];
+ }
+ delete[] matrix;
+ delete[] coeff;
+}
+//---------------------------------------------------------------------------
+void mgcSplineD::
+Create ()
+{
+ int k, d;
+
+ // setup degree constants
+ Dp1 = degree+1;
+ Dp1_to_N = 1;
+ for (d = 0; d < dimensions; d++)
+ Dp1_to_N *= Dp1;
+ Dp1_to_2N = Dp1_to_N*Dp1_to_N;
+
+ // initialize base indices
+ base = new int[dimensions];
+ if ( base == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ old_base = new int[dimensions];
+ if ( old_base == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (d = 0; d < dimensions; d++)
+ old_base[d] = -1;
+
+ grid_min = new int[dimensions];
+ if ( grid_min == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ grid_max = new int[dimensions];
+ if ( grid_max == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (d = 0; d < dimensions; d++) {
+ grid_min[d] = -1;
+ grid_max[d] = -1;
+ }
+
+ // compute domain [min,max] for B-spline
+ dom_min = new double[dimensions];
+ if ( dom_min == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ dom_max = new double[dimensions];
+ if ( dom_max == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ for (d = 0; d < dimensions; d++) {
+ double tmpmin = 1, tmpmax, dom_sup = dim[d]-degree+1;
+ double next = (tmpmin+dom_sup)/2;
+ do {
+ tmpmax = next;
+ next = (next+dom_sup)/2;
+ } while ( next < dom_sup );
+
+ dom_min[d] = tmpmin;
+ dom_max[d] = tmpmax;
+ }
+
+ // cache for optimizing compute_intermediate()
+ cache = new double[Dp1_to_N];
+ if ( cache == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ // storage for intermediate tensor product
+ inter = new double[Dp1_to_N];
+ if ( inter == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+
+ // polynomial allocations
+ poly = new double*[dimensions];
+ if ( poly == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (d = 0; d < dimensions; d++) {
+ poly[d] = new double[Dp1];
+ if ( poly[d] == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ }
+
+ // coefficients for polynomial calculations
+ coeff = new double*[Dp1];
+ if ( coeff == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (int row = 0; row <= degree; row++) {
+ coeff[row] = new double[Dp1];
+ if ( coeff[row] == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (int col = row; col <= degree; col++) {
+ coeff[row][col] = 1;
+ for (d = 0; d <= row-1; d++)
+ coeff[row][col] *= col-d;
+ }
+ }
+
+ // generate spline blending matrix
+ matrix = BlendMatrix(degree);
+
+ // tensor product of m with itself N times
+ product = new double[Dp1_to_2N];
+ if ( product == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ skip = new int[Dp1_to_2N];
+ if ( skip == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ int* coord = new int[2*dimensions]; // for address decoding
+ if ( coord == 0 ) {
+ Report(allocation_failed);
+ return;
+ }
+ for (k = 0; k < Dp1_to_2N; k++) {
+ int temp = k;
+ for (d = 0; d < 2*dimensions; d++) {
+ coord[d] = temp % Dp1;
+ temp /= Dp1;
+ }
+ product[k] = 1;
+ for (int s = 0; s < dimensions; s++)
+ product[k] *= matrix[coord[s]][coord[s+dimensions]];
+
+ skip[k] = 1;
+ }
+ delete[] coord;
+
+ // compute increments to skip zero values of mtensor
+ for (k = 0; k < Dp1_to_2N; ) {
+ for (d = k+1; d < Dp1_to_2N && product[d] == 0; d++)
+ skip[k]++;
+ k = d;
+ }
+
+ evaluate_callback = 0;
+}
+//---------------------------------------------------------------------------
+mgcSplineD& mgcSplineD::
+Renew (int _dimensions, int _degree, const int* _dim, double* _data)
+{
+ // blow away old data
+ FreeStorage();
+
+ // get number of dimensions
+ if ( (dimensions = _dimensions) <= 0 ) {
+ Report(invalid_dimensions);
+ return *this;
+ }
+
+ // get degree and check if valid
+ if ( (degree = _degree) <= 0 ) {
+ Report(invalid_degree);
+ return *this;
+ }
+
+ // get dimension bounds and check if valid
+ if ( _dim == 0 ) {
+ Report(invalid_dim);
+ return *this;
+ }
+
+ int d;
+ for (d = 0; d < dimensions; d++)
+ if ( _dim[d] <= degree+1 ) {
+ Report(invalid_dim);
+ return *this;
+ }
+
+ dim = new int[dimensions];
+ if ( dim == 0 ) {
+ Report(allocation_failed);
+ return *this;
+ }
+ for (d = 0; d < dimensions; d++)
+ dim[d] = _dim[d];
+
+ // get spline data
+ if ( (data = _data) == 0 ) {
+ Report(null_data);
+ return *this;
+ }
+
+ Create();
+
+ return *this;
+}
+//---------------------------------------------------------------------------
+mgcSplineD::
+mgcSplineD (int _dimensions, int _degree, const int* _dim, double* _data)
+{
+ InitNull();
+ Renew(_dimensions,_degree,_dim,_data);
+}
+//---------------------------------------------------------------------------
+void mgcSplineD::
+SetPolynomial (int order, double diff, double* poly)
+{
+ double diffpower = 1;
+ for (int d = order; d <= degree; d++) {
+ poly[d] = coeff[order][d]*diffpower;
+ diffpower *= diff;
+ }
+}
+//---------------------------------------------------------------------------
+int mgcSplineD::
+Choose (int n, int k)
+{
+ int i;
+
+ // computes combination "n choose k"
+
+ if ( n <= 1 || k >= n )
+ return 1;
+
+ int result = 1;
+ for (i = 0; i < k; i++)
+ result *= n-i;
+ for (i = 1; i <= k; i++)
+ result /= i;
+
+ return result;
+}
+//---------------------------------------------------------------------------
+double** mgcSplineD::
+BlendMatrix (int deg)
+{
+ int degp1 = deg+1;
+ int row, col, i, j, d;
+
+ // allocate triple arrays
+ int*** a = new int**[degp1];
+ if ( a == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ int*** b = new int**[degp1];
+ if ( b == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (d = 0; d <= deg; d++) {
+ a[d] = new int*[degp1];
+ if ( a[d] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ b[d] = new int*[degp1];
+ if ( b[d] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (row = 0; row <= deg; row++) {
+ a[d][row] = new int[degp1];
+ if ( a[d][row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ b[d][row] = new int[degp1];
+ if ( b[d][row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (col = 0; col <= deg; col++) {
+ a[d][row][col] = 0;
+ b[d][row][col] = 0;
+ }
+ }
+ }
+
+ a[0][0][0] = 1;
+ b[0][0][0] = 1;
+
+ for (d = 1; d <= deg; d++) {
+ // compute a[]
+ for (row = 0; row <= d; row++) {
+ for (col = 0; col <= d; col++) {
+ a[d][row][col] = 0;
+ if ( col >= 1 ) {
+ a[d][row][col] += a[d-1][row][col-1];
+ if ( row >= 1 )
+ a[d][row][col] -= b[d-1][row-1][col-1];
+ }
+ if ( row >= 1 )
+ a[d][row][col] += (d+1)*b[d-1][row-1][col];
+ }
+ }
+
+ // compute b[]
+ for (row = 0; row <= d; row++) {
+ for (col = 0; col <= d; col++) {
+ b[d][row][col]= 0;
+ for (i = col; i <= d; i++)
+ if ( (i-col) % 2 )
+ b[d][row][col] -= Choose(i,col)*a[d][row][i];
+ else
+ b[d][row][col] += Choose(i,col)*a[d][row][i];
+ }
+ }
+ }
+
+ double** imat = new double*[degp1];
+ if ( imat == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (row = 0; row <= deg; row++) {
+ imat[row] = new double[degp1];
+ if ( imat[row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (col = 0; col <= deg; col++) {
+ imat[row][col]= 0;
+ for (i = col; i <= deg; i++) {
+ int prod = 1;
+ for (j = 1; j <= i-col; j++)
+ prod *= deg-row;
+ imat[row][col] += prod*Choose(i,col)*a[deg][deg-row][i];
+ }
+ }
+ }
+
+ double factorial = 1;
+ for (d = 1; d <= deg; d++)
+ factorial *= d;
+ double** mat = new double*[degp1];
+ if ( mat == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (row = 0; row <= deg; row++) {
+ mat[row] = new double[degp1];
+ if ( mat[row] == 0 ) {
+ Report(allocation_failed);
+ return 0;
+ }
+ for (col = 0; col <= deg; col++)
+ mat[row][col] = imat[row][col]/factorial;
+ }
+
+ // deallocate triple arrays
+ for (d = 0; d <= deg; d++) {
+ for (row = 0; row <= deg; row++) {
+ delete[] b[d][row];
+ delete[] a[d][row];
+ }
+ delete[] b[d];
+ delete[] a[d];
+ }
+ delete[] b;
+ delete[] a;
+
+ // deallocate integer matrix
+ for (d = 0; d <= deg; d++)
+ delete[] imat[d];
+ delete[] imat;
+
+ return mat;
+}
+//---------------------------------------------------------------------------
+int mgcSplineD::
+Number (unsigned single_error)
+{
+ int result;
+ for (result = -1; single_error; single_error >>= 1)
+ result++;
+ return result;
+}
+//---------------------------------------------------------------------------
+void mgcSplineD::
+Report (unsigned single_error)
+{
+ if ( mgcSplineD::verbose )
+ std::cout << "mgcSplineD: " << message[Number(single_error)] << std::endl;
+ else
+ std::ofstream("spline.err",std::ios::out|std::ios::app)
+ << "mgcSplineD: " << message[Number(single_error)] << std::endl;
+
+ error |= single_error;
+}
+//---------------------------------------------------------------------------
+void mgcSplineD::
+Report (std::ostream& ostr)
+{
+ for (unsigned single_error = 1; single_error; single_error <<= 1)
+ if ( error & single_error )
+ ostr << "mgcSplineD: " << message[Number(single_error)] << std::endl;
+
+ error = 0;
+}
+//===========================================================================
+
--- sandbox/magfield/src/spline.h (rev 0) +++ sandbox/magfield/src/spline.h 2014-04-03 15:54:52 UTC (rev 3075) @@ -0,0 +1,177 @@
+// Magic Software, Inc.
+// http://www.magic-software.com
+// Copyright (c) 2000, 2001. All Rights Reserved
+//
+// Source code from Magic Software is supplied under the terms of a license
+// agreement and may not be copied or disclosed except in accordance with the
+// terms of that agreement. The various license agreements may be found at
+// the Magic Software web site. This file is subject to the license
+//
+// FREE SOURCE CODE
+// http://www.magic-software.com/License/free.pdf
+
+#ifndef SPLINE_H
+#define SPLINE_H
+
+#include <iostream>
+#include <math.h>
+
+class mgcSpline
+{
+public:
+ mgcSpline () { InitNull(); }
+ mgcSpline (int _dimensions, int _degree, const int* _dim, double* _data);
+ virtual ~mgcSpline () { FreeStorage(); }
+
+ int Dimensions () { return dimensions; }
+ int Degree () { return degree; }
+ const double** Matrix () { return (const double**)matrix; }
+ double DomainMin (int d = 0) { return dom_min[d]; }
+ double DomainMax (int d = 0) { return dom_max[d]; }
+ int GridMin (int d = 0) { return grid_min[d]; }
+ int GridMax (int d = 0) { return grid_max[d]; }
+ double (*evaluate_callback)(int);
+
+ static const double INVALID_DATA;
+
+ // spline evaluation for function interpolation (no derivatives)
+ virtual double operator() (double* x) = 0;
+
+ // spline evaluation, derivative counts given in dx[]
+ virtual double operator() (int* dx, double* x) = 0;
+
+ // reallocate spline with new parameters
+ mgcSpline& Renew ( int _dimensions, int _degree, const int* _dim,
+ double* _data );
+
+protected:
+ int dimensions; // N, number of dimensions
+ int degree; // D, degree of polynomial spline
+ int Dp1; // D+1
+ int Dp1_to_N; // power(D+1,N)
+ int Dp1_to_2N; // power(D+1,2N)
+ int* dim; // dimension sizes dim[0] through dim[N-1]
+ double* data; // N-dimensional array of data
+ double* dom_min; // minimum allowed value of spline input vector
+ double* dom_max; // maximum allowed value of spline input vector
+ int* grid_min; // minimum allowed value for current local grid
+ int* grid_max; // maximum allowed value for current local grid
+ int* base; // base indices for grid of local control points
+ int* old_base; // old base indices for grid of local control points
+ double* cache; // cache for subblock of data
+ double* inter; // intermediate product of data with blending matrix
+ double** poly; // poly[N][D+1] for storing polynomials and derivatives
+ double** coeff; // coefficients for polynomial construction
+ double** matrix; // (D+1)x(D+1) blending matrix
+ double* product; // outer tensor product of m with itself N times
+ int* skip; // for skipping zero values of mtensor
+
+ virtual void EvaluateUnknownData () = 0;
+ virtual void ComputeIntermediate () = 0;
+ void SetPolynomial (int order, double diff, double* poly);
+
+private:
+ void InitNull ();
+ void FreeStorage ();
+ void Create ();
+ static int Choose (int n, int k);
+ static double** BlendMatrix (int deg);
+
+// error handling
+public:
+ static int verbose;
+ static unsigned error;
+ static void Report (std::ostream& ostr);
+private:
+ static const unsigned invalid_dimensions;
+ static const unsigned invalid_degree;
+ static const unsigned invalid_dim;
+ static const unsigned null_data;
+ static const unsigned allocation_failed;
+ static const char* message[5];
+ static int Number (unsigned single_error);
+ static void Report (unsigned single_error);
+};
+
+
+
+class mgcSplineD
+{
+public:
+ mgcSplineD () { InitNull(); }
+ mgcSplineD (int _dimensions, int _degree, const int* _dim, double* _data);
+ virtual ~mgcSplineD () { FreeStorage(); }
+
+ int Dimensions () { return dimensions; }
+ int Degree () { return degree; }
+ const double** Matrix () { return (const double**)matrix; }
+ double DomainMin (int d = 0) { return dom_min[d]; }
+ double DomainMax (int d = 0) { return dom_max[d]; }
+ int GridMin (int d = 0) { return grid_min[d]; }
+ int GridMax (int d = 0) { return grid_max[d]; }
+ double (*evaluate_callback)(int);
+
+ static const double INVALID_DATA;
+
+ // spline evaluation for function interpolation (no derivatives)
+ virtual double operator() (double* x) = 0;
+
+ // spline evaluation, derivative counts given in dx[]
+ virtual double operator() (int* dx, double* x) = 0;
+
+ // reallocate spline with new parameters
+ mgcSplineD& Renew ( int _dimensions, int _degree, const int* _dim,
+ double* _data );
+
+protected:
+ int dimensions; // N, number of dimensions
+ int degree; // D, degree of polynomial spline
+ int Dp1; // D+1
+ int Dp1_to_N; // power(D+1,N)
+ int Dp1_to_2N; // power(D+1,2N)
+ int* dim; // dimension sizes dim[0] through dim[N-1]
+ double* data; // N-dimensional array of data
+ double* dom_min; // minimum allowed value of spline input vector
+ double* dom_max; // maximum allowed value of spline input vector
+ int* grid_min; // minimum allowed value for current local grid
+ int* grid_max; // maximum allowed value for current local grid
+ int* base; // base indices for grid of local control points
+ int* old_base; // old base indices for grid of local control points
+ double* cache; // cache for subblock of data
+ double* inter; // intermediate product of data with blending matrix
+ double** poly; // poly[N][D+1] for storing polynomials and derivatives
+ double** coeff; // coefficients for polynomial construction
+ double** matrix; // (D+1)x(D+1) blending matrix
+ double* product; // outer tensor product of m with itself N times
+ int* skip; // for skipping zero values of mtensor
+
+ virtual void EvaluateUnknownData () = 0;
+ virtual void ComputeIntermediate () = 0;
+ void SetPolynomial (int order, double diff, double* poly);
+
+private:
+ void InitNull ();
+ void FreeStorage ();
+ void Create ();
+ static int Choose (int n, int k);
+ static double** BlendMatrix (int deg);
+
+// error handling
+public:
+ static int verbose;
+ static unsigned error;
+ static void Report (std::ostream& ostr);
+private:
+ static const unsigned invalid_dimensions;
+ static const unsigned invalid_degree;
+ static const unsigned invalid_dim;
+ static const unsigned null_data;
+ static const unsigned allocation_failed;
+ static const char* message[5];
+ static int Number (unsigned single_error);
+ static void Report (unsigned single_error);
+};
+
+
+#endif
+
--- sandbox/magfield/src/spline1.cpp (rev 0) +++ sandbox/magfield/src/spline1.cpp 2014-04-03 15:54:52 UTC (rev 3075) @@ -0,0 +1,173 @@
+// Magic Software, Inc.
+// http://www.magic-software.com
+// Copyright (c) 2000, 2001. All Rights Reserved
+//
+// Source code from Magic Software is supplied under the terms of a license
+// agreement and may not be copied or disclosed except in accordance with the
+// terms of that agreement. The various license agreements may be found at
+// the Magic Software web site. This file is subject to the license
+//
+// FREE SOURCE CODE
+// http://www.magic-software.com/License/free.pdf
+
+#include "spline1.h"
+//===========================================================================
+const int mgcSpline1::N = 1;
+//---------------------------------------------------------------------------
+mgcSpline1::
+mgcSpline1 (int _degree, int _dim, double* _data) :
+ mgcSpline(N,_degree,&_dim,_data)
+{
+}
+//---------------------------------------------------------------------------
+void mgcSpline1::
+EvaluateUnknownData ()
+{
+ for (int k = grid_min[0]; k <= grid_max[0]; k++)
+ if ( data[k] == INVALID_DATA )
+ data[k] = evaluate_callback(k);
+}
+//---------------------------------------------------------------------------
+void mgcSpline1::
+ComputeIntermediate ()
+{
+ for (int k = 0; k <= degree; k++) {
+ inter[k] = 0;
+ for (int i = 0, j = base[0]-1; i <= degree; i++, j++)
+ inter[k] += data[j]*matrix[i][k];
+ }
+}
+//---------------------------------------------------------------------------
+double mgcSpline1::
+operator () (double x)
+{
+ base[0] = int(floor(x));
+ if ( old_base[0] != base[0] ) {
+ // switch to new local grid
+ old_base[0] = base[0];
+ grid_min[0] = base[0]-1;
+ grid_max[0] = grid_min[0]+degree;
+
+ // fill in missing grid data if necessary
+ if ( evaluate_callback )
+ EvaluateUnknownData();
+
+ ComputeIntermediate();
+ }
+
+ SetPolynomial(0,x-base[0],poly[0]);
+
+ double result = 0;
+ for (int k = 0; k <= degree; k++)
+ result += poly[0][k]*inter[k];
+ return result;
+}
+//---------------------------------------------------------------------------
+double mgcSpline1::
+operator () (int dx, double x)
+{
+ base[0] = int(floor(x));
+ if ( old_base[0] != base[0] ) {
+ // switch to new local grid
+ old_base[0] = base[0];
+ grid_min[0] = base[0]-1;
+ grid_max[0] = grid_min[0]+degree;
+
+ // fill in missing grid data if necessary
+ if ( evaluate_callback )
+ EvaluateUnknownData();
+
+ ComputeIntermediate();
+ }
+
+ SetPolynomial(dx,x-base[0],poly[0]);
+
+ double result = 0;
+ for (int k = dx; k <= degree; k++)
+ result += poly[0][k]*inter[k];
+ return result;
+}
+//===========================================================================
+const int mgcSplineD1::N = 1;
+//---------------------------------------------------------------------------
+mgcSplineD1::
+mgcSplineD1 (int _degree, int _dim, double* _data) :
+ mgcSplineD(N,_degree,&_dim,_data)
+{
+}
+//---------------------------------------------------------------------------
+void mgcSplineD1::
+EvaluateUnknownData ()
+{
+ for (int k = grid_min[0]; k <= grid_max[0]; k++)
+ if ( data[k] == INVALID_DATA )
+ data[k] = evaluate_callback(k);
+}
+//---------------------------------------------------------------------------
+void mgcSplineD1::
+ComputeIntermediate ()
+{
+ for (int k = 0; k <= degree; k++) {
+ inter[k] = 0;
+ for (int i = 0, j = base[0]-1; i <= degree; i++, j++)
+ inter[k] += data[j]*matrix[i][k];
+ }
+}
+//---------------------------------------------------------------------------
+double mgcSplineD1::
+operator () (double x)
+{
+ base[0] = int(floor(x));
+ if ( old_base[0] != base[0] ) {
+ // switch to new local grid
+ old_base[0] = base[0];
+ grid_min[0] = base[0]-1;
+ grid_max[0] = grid_min[0]+degree;
+
+ // fill in missing grid data if necessary
+ if ( evaluate_callback )
+ EvaluateUnknownData();
+
+ ComputeIntermediate();
+ }
+
+ SetPolynomial(0,x-base[0],poly[0]);
+
+ double result = 0;
+ for (int k = 0; k <= degree; k++)
+ result += poly[0][k]*inter[k];
+ return result;
+}
+//---------------------------------------------------------------------------
+double mgcSplineD1::
+operator () (int dx, double x)
+{
+ base[0] = int(floor(x));
+ if ( old_base[0] != base[0] ) {
+ // switch to new local grid
+ old_base[0] = base[0];
+ grid_min[0] = base[0]-1;
+ grid_max[0] = grid_min[0]+degree;
+
+ // fill in missing grid data if necessary
+ if ( evaluate_callback )
+ EvaluateUnknownData();
+
+ ComputeIntermediate();
+ }
+
+ SetPolynomial(dx,x-base[0],poly[0]);
+
+ double result = 0;
+ for (int k = dx; k <= degree; k++)
+ result += poly[0][k]*inter[k];
+ return result;
+}
+//===========================================================================
+
+
+#ifdef SPLINE1_TEST
+
+
+#endif
+
--- sandbox/magfield/src/spline1.h (rev 0) +++ sandbox/magfield/src/spline1.h 2014-04-03 15:54:52 UTC (rev 3075) @@ -0,0 +1,74 @@
+// Magic Software, Inc.
+// http://www.magic-software.com
+// Copyright (c) 2000, 2001. All Rights Reserved
+//
+// Source code from Magic Software is supplied under the terms of a license
+// agreement and may not be copied or disclosed except in accordance with the
+// terms of that agreement. The various license agreements may be found at
+// the Magic Software web site. This file is subject to the license
+//
+// FREE SOURCE CODE
+// http://www.magic-software.com/License/free.pdf
+
+#ifndef SPLINE1_H
+#define SPLINE1_H
+
+#include "spline.h"
+class mgcSpline1 : public mgcSpline
+{
+public:
+ mgcSpline1 () : mgcSpline() {;}
+ mgcSpline1 (int _degree, int _dim, double* _data);
+
+ // spline evaluation for function interpolation (no derivatives)
+ double operator() (double x);
+ double operator() (double* x) { return (*this)(*x); }
+
+ // spline evaluation, derivative count given in dx
+ double operator() (int dx, double x);
+ double operator() (int* dx, double* x) { return (*this)(*dx,*x); }
+
+ // reallocate spline with new parameters
+ mgcSpline1& Renew (int _degree, int _dim, double* _data)
+ {
+ mgcSpline::Renew(1,_degree,(const int*)&_dim,_data);
+ return *this;
+ }
+
+private:
+ static const int N;
+
+ void EvaluateUnknownData ();
+ void ComputeIntermediate ();
+};
+
+
+class mgcSplineD1 : public mgcSplineD
+{
+public:
+ mgcSplineD1 () : mgcSplineD() {;}
+ mgcSplineD1 (int _degree, int _dim, double* _data);
+
+ // spline evaluation for function interpolation (no derivatives)
+ double operator() (double x);
+ double operator() (double* x) { return (*this)(*x); }
+
+ // spline evaluation, derivative count given in dx
+ double operator() (int dx, double x);
+ double operator() (int* dx, double* x) { return (*this)(*dx,*x); }
+
+ // reallocate spline with new parameters
+ mgcSplineD1& Renew (int _degree, int _dim, double* _data)
+ {
+ mgcSplineD::Renew(1,_degree,(const int*)&_dim,_data);
+ return *this;
+ }
+
+private:
+ static const int N;
+
+ void EvaluateUnknownData ();
+ void ComputeIntermediate ();
+};
+#endif
+
Use REPLY-ALL to reply to list
To unsubscribe from the LCDET-SVN list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCDET-SVN&A=1