diff --git a/classic/5.0/src/AbsBeamline/SBend3D.cpp b/classic/5.0/src/AbsBeamline/SBend3D.cpp index e08c7811b21781974190bbc07dfc5bab63ca38e1..e8d3ca5c1f76d79a952e54ae3e96e613f7e2358d 100644 --- a/classic/5.0/src/AbsBeamline/SBend3D.cpp +++ b/classic/5.0/src/AbsBeamline/SBend3D.cpp @@ -32,13 +32,14 @@ SBend3D::SBend3D(const std::string &name) : Component(name), map_m(NULL), planarArcGeometry_m(1., 1.), fieldUnits_m(1.), lengthUnits_m(1.), - dummy() { + polyOrder_m(1), smoothOrder_m(1), dummy() { } SBend3D::SBend3D(const SBend3D &right) : Component(right), map_m(NULL), planarArcGeometry_m(right.planarArcGeometry_m), fieldUnits_m(right.fieldUnits_m), lengthUnits_m(right.lengthUnits_m), + polyOrder_m(right.polyOrder_m), smoothOrder_m(right.smoothOrder_m), dummy() { RefPartBunch_m = right.RefPartBunch_m; if (right.map_m != NULL) @@ -116,18 +117,18 @@ void SBend3D::setFieldMapFileName(std::string name) { map_m = NULL; } if (name != "") { - map_m = new SectorMagneticFieldMap - (name, "Dipole", lengthUnits_m, fieldUnits_m); + map_m = new SectorMagneticFieldMap( + name, + "Dipole", + lengthUnits_m, + fieldUnits_m, + polyOrder_m, + smoothOrder_m); double r_curv = (map_m->getPolarBoundingBoxMax()[0]+ map_m->getPolarBoundingBoxMin()[0])/2.; double delta_phi = map_m->getDeltaPhi(); planarArcGeometry_m.setElementLength(r_curv*delta_phi); planarArcGeometry_m.setCurvature(1./r_curv); - //std::cerr << "RCURV " << r_curv << " DPHI " << delta_phi << " DELTA X: " - // << planarArcGeometry_m.getTotalTransform().getVector()(0) << " Y: " - // << planarArcGeometry_m.getTotalTransform().getVector()(1) << " Z: " - // << planarArcGeometry_m.getTotalTransform().getVector()(2) << " phi: " - // << planarArcGeometry_m.getTotalTransform().getRotation().getAxis()(1) << std::endl; } } diff --git a/classic/5.0/src/AbsBeamline/SBend3D.h b/classic/5.0/src/AbsBeamline/SBend3D.h index 89b14f9d8c9c49a8ece9dd580706b7f89f34fe23..4983618b4e7647b118bdbbf7acc03523589724a0 100644 --- a/classic/5.0/src/AbsBeamline/SBend3D.h +++ b/classic/5.0/src/AbsBeamline/SBend3D.h @@ -156,6 +156,21 @@ class SBend3D : public Component { /** Set the scale factor */ void setFieldUnits(double fieldUnits) {fieldUnits_m = fieldUnits;} + /** Get the polynomial order */ + int getPolynomialOrder() const {return polyOrder_m;} + + /** Set the polynomial order */ + void setPolynomialOrder(int polyOrder) {polyOrder_m = polyOrder;} + + /** Get the smoothing factor */ + int getSmoothingOrder() const {return smoothOrder_m;} + + /** Set the smoothing factor */ + void setSmoothingOrder(int polyOrder) {smoothOrder_m = polyOrder;} + + /** Get the sector magnetic field map */ + SectorMagneticFieldMap* getSectorMagneticFieldMap() const {return map_m;} + private: SectorMagneticFieldMap* map_m; // Geometry @@ -163,6 +178,9 @@ class SBend3D : public Component { // Units used for field maps double fieldUnits_m; double lengthUnits_m; + // Polynomial interpolation + int polyOrder_m; + double smoothOrder_m; // Not implemented (I think this is something to do with error fields?) BMultipoleField dummy; }; diff --git a/classic/5.0/src/Fields/CMakeLists.txt b/classic/5.0/src/Fields/CMakeLists.txt index 6fdb211ba6c436570ef649f4058543bb3ca07bd7..f6f466344524ed25f15b672a240d8bebdb67cafb 100644 --- a/classic/5.0/src/Fields/CMakeLists.txt +++ b/classic/5.0/src/Fields/CMakeLists.txt @@ -37,11 +37,14 @@ set (_SRCS FM1DProfile2.cpp FMDummy.cpp SectorMagneticFieldMap.cpp + SectorField.cpp ) include_directories ( - SectorMagneticFieldMap + Interpolation ${CMAKE_CURRENT_SOURCE_DIR} ) +add_subdirectory (Interpolation) + add_cl_sources(${_SRCS}) diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/CMakeLists.txt b/classic/5.0/src/Fields/Interpolation/CMakeLists.txt similarity index 64% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/CMakeLists.txt rename to classic/5.0/src/Fields/Interpolation/CMakeLists.txt index 4e54721d1a5a5d8dbf3c61ac99a518393771437e..71c18eb45cc1b76cb679b28029e1b5d864d55346 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/CMakeLists.txt +++ b/classic/5.0/src/Fields/Interpolation/CMakeLists.txt @@ -2,14 +2,19 @@ set (_SRCS Interpolator3dGridTo1d.cpp Interpolator3dGridTo3d.cpp Mesh.cpp - SectorField.cpp + MMatrix.cpp + MVector.cpp + PolynomialPatch.cpp + PPSolveFactory.cpp + SolveFactory.cpp + SquarePolynomialVector.cpp ThreeDGrid.cpp TriLinearInterpolator.cpp VectorMap.cpp - ) +) include_directories ( ${CMAKE_CURRENT_SOURCE_DIR} - ) +) add_cl_sources(${_SRCS}) diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.cpp b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo1d.cpp similarity index 95% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.cpp rename to classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo1d.cpp index 494ab261bd30fa9e50f3a3dea479b1c00fc7fa9c..a988b7e8d4437d1c70f6dc95c86e60a56846ea7f 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.cpp +++ b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo1d.cpp @@ -25,8 +25,9 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h" +#include "Fields/Interpolation/Interpolator3dGridTo1d.h" +namespace interpolation { void Interpolator3dGridTo1d::deleteFunc(double*** func) { if (func == NULL) return; @@ -38,4 +39,4 @@ void Interpolator3dGridTo1d::deleteFunc(double*** func) { delete [] func; func = NULL; } - +} diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo1d.h similarity index 98% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h rename to classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo1d.h index b042459e4db1c7572a2be77b0dd5abf3fee46223..fc17c6e78ffae1a8db3880431ac5bef9281cfc88 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h +++ b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo1d.h @@ -28,8 +28,10 @@ #ifndef _CLASSIC_FIELDS_INTERPOLATOR3DGRIDTO1D_H_ #define _CLASSIC_FIELDS_INTERPOLATOR3DGRIDTO1D_H_ -#include "Fields/SectorMagneticFieldMap/VectorMap.h" -#include "Fields/SectorMagneticFieldMap/ThreeDGrid.h" +#include "Fields/Interpolation/VectorMap.h" +#include "Fields/Interpolation/ThreeDGrid.h" + +namespace interpolation { /** Interpolator3dGridTo1d is an abstraction for lookup on a 3D mesh to get a 1D * value. @@ -216,4 +218,5 @@ void Interpolator3dGridTo1d::clear() { coordinates_m->remove(this); } +} #endif // _CLASSIC_FIELDS_INTERPOLATOR3DGRIDTO1D_HH_ diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.cpp b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo3d.cpp similarity index 97% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.cpp rename to classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo3d.cpp index 31a04529982c46baf5dc94a9161f33ae1dbbdc16..3e4ee2d863498c8dbea3ceb9bf2e0f49044d41c5 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.cpp +++ b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo3d.cpp @@ -25,9 +25,11 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h" +#include "Fields/Interpolation/Interpolator3dGridTo3d.h" #include "Utilities/LogicalError.h" +namespace interpolation { + Interpolator3dGridTo3d::Interpolator3dGridTo3d (const Interpolator3dGridTo3d& rhs) { coordinates_m = new ThreeDGrid(*rhs.coordinates_m); @@ -77,3 +79,4 @@ void Interpolator3dGridTo3d::setAll(ThreeDGrid* grid, ); } } +} diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo3d.h similarity index 97% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h rename to classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo3d.h index 6f81461ee74b2418f84b7813e4000c3872a93fb0..8280b5d61a2199e1256d81fbcf56970c4fb073c9 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h +++ b/classic/5.0/src/Fields/Interpolation/Interpolator3dGridTo3d.h @@ -28,10 +28,12 @@ #ifndef _CLASSIC_FIELDS_INTERPOLATOR3DGRIDTO3D_HH_ #define _CLASSIC_FIELDS_INTERPOLATOR3DGRIDTO3D_HH_ -#include "Fields/SectorMagneticFieldMap/Mesh.h" -#include "Fields/SectorMagneticFieldMap/VectorMap.h" -#include "Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h" -#include "Fields/SectorMagneticFieldMap/TriLinearInterpolator.h" +#include "Fields/Interpolation/Mesh.h" +#include "Fields/Interpolation/VectorMap.h" +#include "Fields/Interpolation/Interpolator3dGridTo1d.h" +#include "Fields/Interpolation/TriLinearInterpolator.h" + +namespace interpolation { /** Interpolator3dGridTo3d interpolates from 3d grid to a 3d vector * @@ -219,4 +221,6 @@ void Interpolator3dGridTo3d::clear() { coordinates_m->remove(this); } +} + #endif diff --git a/classic/5.0/src/Fields/Interpolation/MMatrix.cpp b/classic/5.0/src/Fields/Interpolation/MMatrix.cpp new file mode 100644 index 0000000000000000000000000000000000000000..98b59684fd00ee499879ae0de5aa7a256e52a7f2 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/MMatrix.cpp @@ -0,0 +1,434 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "gsl/gsl_linalg.h" +#include "gsl/gsl_eigen.h" + +#include "Fields/Interpolation/MMatrix.h" + +///////////////////////////////// MMATRIX /////////////////////////////// + +//template class MMatrix<double>; + +/////////// CONSTRUCTORS, DESTRUCTORS +namespace interpolation { + +template <class Tmplt> +MMatrix<Tmplt>::MMatrix() : _matrix(NULL) +{} +template MMatrix<double> ::MMatrix(); +template MMatrix<m_complex>::MMatrix(); + + +template <> +void MMatrix<double>::delete_matrix() +{ + if(_matrix != NULL) gsl_matrix_free( (gsl_matrix*)_matrix ); + _matrix = NULL; +} + +template <> +void MMatrix<m_complex>::delete_matrix() +{ + if(_matrix != NULL) gsl_matrix_complex_free( (gsl_matrix_complex*)_matrix ); + _matrix = NULL; +} + + +template <> +MMatrix<double>& MMatrix<double>::operator= (const MMatrix<double>& mm) +{ + if (&mm == this) return *this; + delete_matrix(); + if(!mm._matrix) { _matrix = NULL; return *this; } + _matrix = gsl_matrix_alloc( mm.num_row(), mm.num_col() ); + gsl_matrix_memcpy((gsl_matrix*)_matrix, (const gsl_matrix*)mm._matrix); + return *this; +} + +template <> +MMatrix<m_complex>& MMatrix<m_complex>::operator= (const MMatrix<m_complex>& mm) +{ + if (&mm == this) return *this; + delete_matrix(); + if(!mm._matrix) { _matrix = NULL; return *this; } + _matrix = gsl_matrix_complex_alloc( mm.num_row(), mm.num_col() ); + gsl_matrix_complex_memcpy((gsl_matrix_complex*)_matrix, (const gsl_matrix_complex*)mm._matrix); + return *this; +} + +template <> +MMatrix<double>::MMatrix( const MMatrix<double>& mm ) : _matrix(NULL) +{ + if(mm._matrix) + { + _matrix = gsl_matrix_alloc(mm.num_row(), mm.num_col() ); + gsl_matrix_memcpy( (gsl_matrix*)_matrix, (gsl_matrix*)mm._matrix ); + } +} + +template <> +MMatrix<m_complex>::MMatrix( const MMatrix<m_complex>& mm ) : _matrix(NULL) +{ + if(mm._matrix) + { + _matrix = gsl_matrix_complex_alloc( mm.num_row(), mm.num_col() ); + gsl_matrix_complex_memcpy( (gsl_matrix_complex*)_matrix, (gsl_matrix_complex*)mm._matrix ); + } +} + +template <class Tmplt> +MMatrix<Tmplt>::MMatrix(size_t i, size_t j, Tmplt* data_beg ) : _matrix(NULL) +{ + build_matrix(i, j, data_beg); +} +template MMatrix<double> ::MMatrix(size_t i, size_t j, double* data_beg ); +template MMatrix<m_complex>::MMatrix(size_t i, size_t j, m_complex* data_beg ); + +template <class Tmplt> +MMatrix<Tmplt>::MMatrix(size_t i, size_t j, Tmplt value ) +{ + build_matrix(i, j); + for(size_t a=1; a<=i; a++) + for(size_t b=1; b<=j; b++) + operator()(a,b) = value; +} +template MMatrix<double> ::MMatrix(size_t i, size_t j, double value); +template MMatrix<m_complex>::MMatrix(size_t i, size_t j, m_complex value); + +template <class Tmplt> +MMatrix<Tmplt>::MMatrix(size_t i, size_t j ) +{ + build_matrix(i, j); +} +template MMatrix<double> ::MMatrix(size_t i, size_t j ); +template MMatrix<m_complex>::MMatrix(size_t i, size_t j ); + +template <class Tmplt> +MMatrix<Tmplt> MMatrix<Tmplt>::Diagonal(size_t i, Tmplt diag_value, Tmplt off_diag_value) +{ + MMatrix<Tmplt> mm(i,i); + for(size_t a=1; a<=i; a++) + { + for(size_t b=1; b< a; b++) mm(a, b) = off_diag_value; + for(size_t b=a+1; b<=i; b++) mm(a, b) = off_diag_value; + mm(a,a) = diag_value; + } + return mm; +} +template MMatrix<double> MMatrix<double> ::Diagonal(size_t i, double diag, double off_diag); +template MMatrix<m_complex> MMatrix<m_complex>::Diagonal(size_t i, m_complex diag, m_complex off_diag); + + +template <class Tmplt> +MMatrix<Tmplt>::~MMatrix() +{ + delete_matrix(); +} +template MMatrix<double>::~MMatrix(); +template MMatrix<m_complex>::~MMatrix(); + +template <> +void MMatrix<double>::build_matrix(size_t i, size_t j) +{ + _matrix = gsl_matrix_alloc(i, j); +} + +template <> +void MMatrix<m_complex>::build_matrix(size_t i, size_t j) +{ + _matrix = gsl_matrix_complex_alloc(i, j); +} + +template <> +void MMatrix<double>::build_matrix(size_t i, size_t j, double* data) +{ + _matrix = gsl_matrix_alloc(i, j); + for(size_t a=0; a<i; a++) + for(size_t b=0; b<j; b++) + operator()(a+1,b+1) = data[b*i + a]; +} + +template <> +void MMatrix<m_complex>::build_matrix(size_t i, size_t j, m_complex* data) +{ + _matrix = gsl_matrix_complex_alloc(i, j); + for(size_t a=0; a<i; a++) + for(size_t b=0; b<j; b++) + operator()(a+1,b+1) = data[b*i + a]; +} + + +////////////////// HIGH LEVEL FUNCTIONS + +template <> +m_complex MMatrix<m_complex>::determinant() const +{ + int signum = 1; + + if(num_row() != num_col()) throw(GeneralClassicException("MMatrix::determinant()", "Attempt to get determinant of non-square matrix")); + gsl_permutation * p = gsl_permutation_alloc (num_row()); + MMatrix<m_complex> copy(*this); + gsl_linalg_complex_LU_decomp ((gsl_matrix_complex*)copy._matrix, p, &signum); + gsl_permutation_free(p); + return gsl_linalg_complex_LU_det((gsl_matrix_complex*)copy._matrix, signum); +} + +template <> +double MMatrix<double>::determinant() const +{ + int signum = 1; + if(num_row() != num_col()) throw(GeneralClassicException("MMatrix::determinant()", "Attempt to get determinant of non-square matrix")); + gsl_permutation * p = gsl_permutation_alloc (num_row()); + MMatrix<double> copy(*this); + gsl_linalg_LU_decomp ((gsl_matrix*)copy._matrix, p, &signum); + double d = gsl_linalg_LU_det((gsl_matrix*)copy._matrix, signum); + gsl_permutation_free(p); + return d; +} + +template <class Tmplt> +MMatrix<Tmplt> MMatrix<Tmplt>::inverse() const +{ + MMatrix<Tmplt> copy(*this); + copy.invert(); + return copy; +} +template MMatrix<double> MMatrix<double> ::inverse() const; +template MMatrix<m_complex> MMatrix<m_complex>::inverse() const; + + +template <> +void MMatrix<m_complex>::invert() +{ + if(num_row() != num_col()) throw(GeneralClassicException("MMatrix::invert()", "Attempt to get inverse of non-square matrix")); + gsl_permutation* perm = gsl_permutation_alloc(num_row()); + MMatrix<m_complex> lu(*this); //hold LU decomposition + int s=0; + gsl_linalg_complex_LU_decomp ( (gsl_matrix_complex*)lu._matrix, perm, &s); + gsl_linalg_complex_LU_invert ( (gsl_matrix_complex*)lu._matrix, perm, (gsl_matrix_complex*)_matrix); + gsl_permutation_free( perm ); + for(size_t i=1; i<=num_row(); i++) + for(size_t j=1; j<=num_row(); j++) + if(operator()(i,j) != operator()(i,j)) throw(GeneralClassicException("MMatrix::invert()", "Failed to invert matrix - singular?")); +} + +template <> +void MMatrix<double>::invert() +{ + if(num_row() != num_col()) throw(GeneralClassicException("MMatrix::invert()", "Attempt to get inverse of non-square matrix")); + gsl_permutation* perm = gsl_permutation_alloc(num_row()); + MMatrix<double> lu(*this); //hold LU decomposition + int s=0; + gsl_linalg_LU_decomp ( (gsl_matrix*)lu._matrix, perm, &s); + gsl_linalg_LU_invert ( (gsl_matrix*)lu._matrix, perm, (gsl_matrix*)_matrix); + gsl_permutation_free( perm ); + for(size_t i=1; i<=num_row(); i++) + for(size_t j=1; j<=num_row(); j++) + if(operator()(i,j) != operator()(i,j)) throw(GeneralClassicException("MMatrix::invert()", "Failed to invert matrix - singular?")); +} + +template <> +MMatrix<double> MMatrix<double>::T() const +{ + MMatrix<double> in(num_col(), num_row()); + gsl_matrix_transpose_memcpy ((gsl_matrix*)in._matrix, (gsl_matrix*)_matrix); + return in; +} + +template <> +MMatrix<m_complex> MMatrix<m_complex>::T() const +{ + MMatrix<m_complex> in(num_col(), num_row()); + gsl_matrix_complex_transpose_memcpy ((gsl_matrix_complex*)in._matrix, (gsl_matrix_complex*)_matrix); + return in; +} + +template <class Tmplt> +MMatrix<Tmplt> MMatrix<Tmplt>::sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const +{ + MMatrix<Tmplt> sub_matrix(max_row-min_row+1, max_col-min_col+1); + for(size_t i=min_row; i<=max_row; i++) + for(size_t j=min_col; j<=max_col; j++) + sub_matrix(i-min_row+1, j-min_col+1) = operator()(i,j); + return sub_matrix; +} +template MMatrix<double> MMatrix<double> ::sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const; +template MMatrix<m_complex> MMatrix<m_complex>::sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const; + + +template <class Tmplt> +Tmplt MMatrix<Tmplt>::trace() const +{ + Tmplt t = operator()(1,1); + for(size_t i=2; i<=num_row() && i<=num_col(); i++) t+= operator()(i,i); + return t; +} +template double MMatrix<double>::trace() const; +template m_complex MMatrix<m_complex>::trace() const; + + +template <class Tmplt> +MVector<m_complex> MMatrix<Tmplt>::eigenvalues() const +{ + if(num_row() != num_col()) throw(GeneralClassicException("MMatrix<double>::eigenvalues", "Attempt to get eigenvectors of non-square matrix") ); + MMatrix<Tmplt> temp = *this; + MVector<m_complex> evals(num_row(), m_complex_build(2.,-1.)); + gsl_eigen_nonsymm_workspace * w = gsl_eigen_nonsymm_alloc(num_row()); + gsl_eigen_nonsymm_params(0, 1, w); + int ierr = gsl_eigen_nonsymm(get_matrix(temp), evals.get_vector(evals), w); + gsl_eigen_nonsymm_free(w); + if(ierr != 0) throw(GeneralClassicException("MMatrix<Tmplt>::eigenvalues", "Failed to calculate eigenvalue")); + return evals; +} +template MVector<m_complex> MMatrix<double> ::eigenvalues() const; + +template <class Tmplt> +std::pair<MVector<m_complex>, MMatrix<m_complex> > MMatrix<Tmplt>::eigenvectors() const +{ + if(num_row() != num_col()) throw(GeneralClassicException("MMatrix<double>::eigenvectors", "Attempt to get eigenvectors of non-square matrix") ); + MMatrix<Tmplt> temp = *this; + MVector<m_complex> evals(num_row()); + MMatrix<m_complex> evecs(num_row(), num_row()); + gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(num_row()); + int ierr = gsl_eigen_nonsymmv(get_matrix(temp), evals.get_vector(evals), get_matrix(evecs), w); + gsl_eigen_nonsymmv_free(w); + if(ierr != 0) throw(GeneralClassicException("MMatrix<Tmplt>::eigenvectors", "Failed to calculate eigenvalue")); + return std::pair<MVector<m_complex>, MMatrix<m_complex> > (evals, evecs) ; +} +template std::pair<MVector<m_complex>, MMatrix<m_complex> > MMatrix<double>::eigenvectors() const; + +///////////// OPERATORS + +MMatrix<double>& operator *= (MMatrix<double>& m1, MMatrix<double> m2) +{ + MMatrix<double> out(m1.num_row(), m2.num_col()); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., (gsl_matrix*)m1._matrix, (gsl_matrix*)m2._matrix, 0., (gsl_matrix*)out._matrix); + m1 = out; + return m1; +} + +MMatrix<m_complex> & operator *= (MMatrix<m_complex>& m1, MMatrix<m_complex> m2) +{ + MMatrix<m_complex> out(m1.num_row(), m2.num_col()); + gsl_blas_zgemm( CblasNoTrans, CblasNoTrans, m_complex_build(1.), (gsl_matrix_complex*)m1._matrix, + (gsl_matrix_complex*)m2._matrix, m_complex_build(0.), (gsl_matrix_complex*)out._matrix); + m1 = out; + return m1; +} + +template <class Tmplt> std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt> mat) +{ + out << mat.num_row() << " " << mat.num_col() << "\n"; + for(size_t i=1; i<=mat.num_row(); i++) + { + out << " "; + for(size_t j=1; j<mat.num_col(); j++) + out << mat(i,j) << " "; + out << mat(i,mat.num_col()) << "\n"; + } + return out; +} +template std::ostream& operator<<(std::ostream& out, MMatrix<double> mat); +template std::ostream& operator<<(std::ostream& out, MMatrix<m_complex> mat); + +template <class Tmplt> std::istream& operator>>(std::istream& in, MMatrix<Tmplt>& mat) +{ + size_t nr, nc; + in >> nr >> nc; + mat = MMatrix<Tmplt>(nr, nc); + for(size_t i=1; i<=nr; i++) + for(size_t j=1; j<=nc; j++) + in >> mat(i,j); + return in; +} +template std::istream& operator>>(std::istream& in, MMatrix<double>& mat); +template std::istream& operator>>(std::istream& in, MMatrix<m_complex>& mat); + + +///////////////// INTERFACES +const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m) +{ + if(m._matrix == NULL) throw(GeneralClassicException("MMatrix_to_gsl", "Attempt to reference uninitialised matrix")); + return (gsl_matrix*)m._matrix; +} + +const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<m_complex>& m) +{ + if(m._matrix == NULL) throw(GeneralClassicException("MMatrix_to_gsl", "Attempt to reference uninitialised matrix")); + return (gsl_matrix_complex*)m._matrix; +} + +MMatrix<double> re(MMatrix<m_complex> mc) +{ + MMatrix<double> md(mc.num_row(), mc.num_col()); + for(size_t i=1; i<=mc.num_row(); i++) + for(size_t j=1; j<=mc.num_col(); j++) + md(i,j) = re(mc(i,j)); + return md; +} + +MMatrix<double> im(MMatrix<m_complex> mc) +{ + MMatrix<double> md(mc.num_row(), mc.num_col()); + for(size_t i=1; i<=mc.num_row(); i++) + for(size_t j=1; j<=mc.num_col(); j++) + md(i,j) = im(mc(i,j)); + return md; +} + +MMatrix<m_complex> complex(MMatrix<double> real) +{ + MMatrix<m_complex> mc(real.num_row(), real.num_col()); + for(size_t i=1; i<=mc.num_row(); i++) + for(size_t j=1; j<=mc.num_col(); j++) + mc(i,j) = m_complex_build(real(i,j)); + return mc; +} + +MMatrix<m_complex> complex(MMatrix<double> real, MMatrix<double> imaginary) +{ + if(real.num_row() != imaginary.num_row() || real.num_col() != imaginary.num_col()) + throw(GeneralClassicException("MMatrix<m_complex>::complex", "Attempt to build complex matrix when real and imaginary matrix don't have the same size")); + MMatrix<m_complex> mc(real.num_row(), real.num_col()); + for(size_t i=1; i<=mc.num_row(); i++) + for(size_t j=1; j<=mc.num_col(); j++) + mc(i,j) = m_complex_build(real(i,j), imaginary(i,j)); + return mc; +} + +template <class Tmplt> +MVector<Tmplt> MMatrix<Tmplt>::get_mvector(size_t column) const +{ + MVector<Tmplt> temp(num_row()); + for(size_t i=1; i<=num_row(); i++) temp(i) = operator()(i,column); + return temp; +} +template MVector<double> MMatrix<double>::get_mvector(size_t column) const; +template MVector<m_complex> MMatrix<m_complex>::get_mvector(size_t column) const; + +} diff --git a/classic/5.0/src/Fields/Interpolation/MMatrix.h b/classic/5.0/src/Fields/Interpolation/MMatrix.h new file mode 100644 index 0000000000000000000000000000000000000000..305f244d37d9646cc4e069a8ba411d1489f4a283 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/MMatrix.h @@ -0,0 +1,374 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef _SRC_COMMON_CPP_MATHS_MMATRIX_HH_ +#define _SRC_COMMON_CPP_MATHS_MMATRIX_HH_ + +#include <iostream> +#include <vector> + +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_blas.h" + +#include "Utilities/GeneralClassicException.h" + +#include "Fields/Interpolation/MVector.h" + +namespace interpolation { + +///////////////// MMatrix ///////////////// + +/** @class MMatrix C++ wrapper for GSL matrix + * MMatrix class handles matrix algebra, maths operators and some + * higher level calculation like matrix inversion, eigenvector + * analysis etc + * + * Use template to define two types:\n + * (i) MMatrix<double> is a matrix of doubles\n + * (ii) MMatrix<m_complex> is a matrix of m_complex\n + * + * Maths operators and a few others are defined, but maths operators + * don't allow operations between types - i.e. you can't multiply + * a complex matrix by a double matrix. Instead use interface methods + * like real() and complex() to convert between types first + */ +template <class Tmplt> +class MMatrix +{ +public: + /** @brief default constructor makes an empty MMatrix of size (0,0) + */ + MMatrix(); + + /** @brief Copy constructor makes a deep copy of mv + */ + MMatrix(const MMatrix<Tmplt>& mv ); + + /** @brief Construct a matrix and fill with data from memory data_beginning + * + * @params nrows number of rows + * @params ncols number of columns + * @params data_beginning pointer to the start of a memory block of size + * nrows*ncols with data laid out <row> <row> <row>. Note MMatrix + * does not take ownership of memory in data_beginning. + */ + MMatrix(size_t nrows, size_t ncols, Tmplt* data_beginning ); + + /** @brief Construct a matrix and fill with identical data + * + * @params nrows number of rows + * @params ncols number of columns + * @params data variable to be copied into all items in the matrix + */ + MMatrix(size_t nrows, size_t ncols, Tmplt value); + + /** @brief Construct a matrix and fill all fields with 0 + * + * @params nrows number of rows + * @params ncols number of columns + */ + MMatrix(size_t nrows, size_t ncols); + + /** @brief Construct a square matrix filling on and off diagonal values + * + * @params i number of rows and number of columns + * @params diag_value fill values with row == column (i.e. on the diagonal) + * with this value + * @params off_diag_value fill values with row != column (i.e. off the + * diagonal) with this value + */ + static MMatrix<Tmplt> Diagonal(size_t i, Tmplt diag_value, Tmplt off_diag_value); + + /** @brief destructor + */ + ~MMatrix(); + + /** @brief returns number of rows in the matrix + */ + size_t num_row() const; + + /** @brief returns number of columns in the matrix + */ + size_t num_col() const; + + /** @brief returns sum of terms with row == column, even if matrix is not + * square + */ + Tmplt trace() const; + + /** @brief returns matrix determinant; throws an exception if matrix is not + * square + */ + Tmplt determinant() const; + + /** @brief returns matrix inverse leaving this matrix unchanged + */ + MMatrix<Tmplt> inverse() const; + + /** @brief turns this matrix into its inverse + */ + void invert(); + + /** @brief returns matrix transpose T (such that M(i,j) = T(j,i)) + */ + MMatrix<Tmplt> T() const; + + /** @brief returns a vector of eigenvalues. Throws an exception if either this + * matrix is not square or the eigenvalues could not be found (e.g. + * singular matrix or whatever). + */ + MVector<m_complex> eigenvalues() const; //return vector of eigenvalues + std::pair<MVector<m_complex>, MMatrix<m_complex> > + eigenvectors() const; //return pair of eigenvalues, eigenvectors (access using my_pair.first or my_pair.second) + + //return a submatrix, with data from min_row to max_row and min_row to max_row inclusive; indexing starts at 1 + MMatrix<Tmplt> sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const; + //create a column vector from the column^th column + MVector<Tmplt> get_mvector(size_t column) const; + //return a reference to the datum; indexing starts at 1, goes to num_row (inclusive) + const Tmplt& operator()(size_t row, size_t column) const; + Tmplt& operator()(size_t row, size_t column); + + MMatrix<Tmplt>& operator= (const MMatrix<Tmplt>& mm); + + //TODO - implement iterator + //class iterator + //{ + //} + + + friend MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m, m_complex c); + friend MMatrix<double>& operator *=(MMatrix<double>& m, double d); + friend MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m1, MMatrix<m_complex> m2); + friend MMatrix<double>& operator *=(MMatrix<double>& m1, MMatrix<double> m2); + friend MVector<m_complex> operator * (MMatrix<m_complex> m, MVector<m_complex> v); + friend MVector<double> operator * (MMatrix<double> m, MVector<double> v); + friend MMatrix<m_complex>& operator +=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2); + friend MMatrix<double>& operator +=(MMatrix<double>& m1, const MMatrix<double>& m2); + template <class Tmplt2> friend MMatrix<Tmplt2> operator + (MMatrix<Tmplt2> m1, const MMatrix<Tmplt2> m2); + + friend const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m); + friend const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<gsl_complex>& m); + + friend class MMatrix<double>; //To do the eigenvector problem, MMatrix<double> needs to see MMatrix<complex>'s _matrix + + +private: + //build the matrix with size i,j, data not initialised + void build_matrix(size_t i, size_t j); + //build the matrix with size i,j, data initialised to data in temp + void build_matrix(size_t i, size_t j, Tmplt* temp); + //delete the matrix and set it to null + void delete_matrix(); + //return the matrix overloaded to the correct type or throw + //if _matrix == NULL + static gsl_matrix* get_matrix(const MMatrix<double>& m); + static gsl_matrix_complex* get_matrix(const MMatrix<m_complex>& m); + //gsl object + void* _matrix; +}; + +//Operators do what you would expect - i.e. normal mathematical functions +MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m, m_complex c); +MMatrix<double>& operator *=(MMatrix<double>& m, double d); +MMatrix<m_complex>& operator *=(MMatrix<m_complex>& m1, MMatrix<m_complex> m2); //M1 *= M2 returns M1 = M1*M2 +MMatrix<double>& operator *=(MMatrix<double>& m1, MMatrix<double> m2); //M1 *= M2 returns M1 = M1*M2 +MVector<m_complex> operator * (MMatrix<m_complex> m, MVector<m_complex> v); +MVector<double> operator * (MMatrix<double> m, MVector<double> v); +MMatrix<m_complex>& operator +=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2); +MMatrix<double>& operator +=(MMatrix<double>& m1, const MMatrix<double>& m2); + +template <class Tmplt> MMatrix<Tmplt>& operator -=(MMatrix<double>& m1, MMatrix<double> m2); + +template <class Tmplt> MMatrix<Tmplt> operator*(MMatrix<Tmplt>, MMatrix<Tmplt>); +template <class Tmplt> MMatrix<Tmplt> operator*(MMatrix<Tmplt>, MVector<Tmplt>); +template <class Tmplt> MMatrix<Tmplt> operator*(MMatrix<Tmplt>, Tmplt); +template <class Tmplt> MMatrix<Tmplt> operator/(MMatrix<Tmplt>, Tmplt); + +template <class Tmplt> MMatrix<Tmplt> operator+(MMatrix<Tmplt>, MMatrix<Tmplt>); +template <class Tmplt> MMatrix<Tmplt> operator-(MMatrix<Tmplt>, MMatrix<Tmplt>); +template <class Tmplt> MMatrix<Tmplt> operator-(MMatrix<Tmplt>); //unary minus returns matrix*(-1) + +template <class Tmplt> bool operator==(const MMatrix<Tmplt>& c1, const MMatrix<Tmplt>& c2); +template <class Tmplt> bool operator!=(const MMatrix<Tmplt>& c1, const MMatrix<Tmplt>& c2); + +//format is +// num_row num_col +// m11 m12 m13 ... +// m21 m22 m23 ... +// ... +template <class Tmplt> std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt> mat); +template <class Tmplt> std::istream& operator>>(std::istream& in, MMatrix<Tmplt>& mat); + +//return matrix of doubles filled with either real or imaginary part of m +MMatrix<double> re(MMatrix<m_complex> m); +MMatrix<double> im(MMatrix<m_complex> m); +//return matrix of m_complex filled with real and imaginary parts +MMatrix<m_complex> complex(MMatrix<double> real); +MMatrix<m_complex> complex(MMatrix<double> real, MMatrix<double> imaginary); + +//return pointer to gsl_matrix objects that store matrix data in m +const gsl_matrix* MMatrix_to_gsl(const MMatrix<double>& m); +const gsl_matrix_complex* MMatrix_to_gsl(const MMatrix<gsl_complex>& m); + +//////////////////////////// MMatrix declaration end /////////////// + + + +//////////////////////////// MMatrix Inlined Functions ////////////// + +template <> +inline size_t MMatrix<double>::num_row() const +{ + if(_matrix) return ((gsl_matrix*)_matrix)->size1; + return 0; +} +template <> +inline size_t MMatrix<m_complex>::num_row() const +{ + if(_matrix) return ((gsl_matrix_complex*)_matrix)->size1; + return 0; +} +template <> +inline size_t MMatrix<double>::num_col() const +{ + if(_matrix) return ((gsl_matrix*)_matrix)->size2; + return 0; +} + +template <> +inline size_t MMatrix<m_complex>::num_col() const +{ + if(_matrix) return ((gsl_matrix_complex*)_matrix)->size2; + return 0; +} + +MMatrix<m_complex> inline & operator *=(MMatrix<m_complex>& m, m_complex c) +{gsl_matrix_complex_scale((gsl_matrix_complex*)m._matrix, c); return m;} + +MMatrix<double> inline & operator *=(MMatrix<double>& m, double d) +{gsl_matrix_scale((gsl_matrix*)m._matrix, d); return m;} + +MVector<m_complex> inline operator * (MMatrix<m_complex> m, MVector<m_complex> v) +{ + MVector<m_complex> v0(m.num_row()); + gsl_blas_zgemv(CblasNoTrans, m_complex_build(1.), (gsl_matrix_complex*)m._matrix, (gsl_vector_complex*)v._vector, m_complex_build(0.), (gsl_vector_complex*)v0._vector); + return v0; +} + +MVector<double> inline operator * (MMatrix<double> m, MVector<double> v) +{ + MVector<double> v0(m.num_row()); + gsl_blas_dgemv(CblasNoTrans, 1., (gsl_matrix*)m._matrix, (gsl_vector*)v._vector, 0., (gsl_vector*)v0._vector); + return v0; +} + +MMatrix<m_complex> inline & operator +=(MMatrix<m_complex>& m1, const MMatrix<m_complex>& m2) +{gsl_matrix_complex_add((gsl_matrix_complex*)m1._matrix, (gsl_matrix_complex*)m2._matrix); return m1;} + +MMatrix<double> inline & operator +=(MMatrix<double>& m1, const MMatrix<double>& m2) +{gsl_matrix_add((gsl_matrix*)m1._matrix, (gsl_matrix*)m2._matrix); return m1;} + +template <class Tmplt> MMatrix<Tmplt> inline operator - (MMatrix<Tmplt> m1) +{for(size_t i=1; i<=m1.num_row(); i++) for(size_t j=1; j<=m1.num_col(); j++) m1(i,j) = -1.*m1(i,j); return m1; } + +template <class Tmplt> MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2) +{return m1*=m2;} + +template <class Tmplt> MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m, MVector<Tmplt> v) +{return m*=v;} + +template <class Tmplt> MMatrix<Tmplt> inline operator*(MMatrix<Tmplt> m, Tmplt t) +{return m*=t;} + +template <class Tmplt> MMatrix<Tmplt> inline & operator/=(MMatrix<Tmplt>& m, Tmplt t) +{return m*=1./t;} +template <class Tmplt> MMatrix<Tmplt> inline operator/(MMatrix<Tmplt> m, Tmplt t) +{return m/=t;} + +template <class Tmplt> MMatrix<Tmplt> inline operator+(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2) +{ MMatrix<Tmplt> m3 = m1+=m2; return m3; } + +template <class Tmplt> MMatrix<Tmplt> inline & operator-=(MMatrix<Tmplt>& m1, MMatrix<Tmplt> m2) +{return m1 += -m2;} +template <class Tmplt> MMatrix<Tmplt> inline operator-(MMatrix<Tmplt> m1, MMatrix<Tmplt> m2) +{return m1 -= m2;} + + + +template <> +const double inline & MMatrix<double>::operator()(size_t i, size_t j) const +{ return *gsl_matrix_ptr((gsl_matrix*)_matrix, i-1, j-1); } +template <> +const m_complex inline & MMatrix<m_complex>::operator()(size_t i, size_t j) const +{ return *gsl_matrix_complex_ptr((gsl_matrix_complex*)_matrix, i-1, j-1); } + +template <> +double inline & MMatrix<double>::operator()(size_t i, size_t j) +{ return *gsl_matrix_ptr((gsl_matrix*)_matrix, i-1, j-1); } +template <> +m_complex inline & MMatrix<m_complex>::operator()(size_t i, size_t j) +{ return *gsl_matrix_complex_ptr((gsl_matrix_complex*)_matrix, i-1, j-1); } + +template <class Tmplt> +bool operator==(const MMatrix<Tmplt>& lhs, const MMatrix<Tmplt>& rhs) +{ + if( lhs.num_row() != rhs.num_row() || lhs.num_col() != rhs.num_col() ) return false; + for(size_t i=1; i<=lhs.num_row(); i++) + for(size_t j=1; j<=lhs.num_col(); j++) + if(lhs(i,j) != rhs(i,j)) return false; + return true; +} + +template <class Tmplt> +bool inline operator!=(const MMatrix<Tmplt>& lhs, const MMatrix<Tmplt>& rhs) {return ! (lhs == rhs);} + +//iostream +template <class Tmplt> std::ostream& operator<<(std::ostream& out, MMatrix<Tmplt>); +//template <class Tmplt> std::istream& operator>>(std::istream& in, MMatrix<Tmplt>); + +template <class Tmplt> +gsl_matrix inline* MMatrix<Tmplt>::get_matrix(const MMatrix<double>& m) +{ + if(m._matrix == NULL) + throw(GeneralClassicException("MMatrix::get_matrix", "Attempt to access uninitialised matrix")); + return (gsl_matrix*)m._matrix; +} + +template <class Tmplt> +gsl_matrix_complex inline* MMatrix<Tmplt>::get_matrix(const MMatrix<m_complex>& m) +{ + if(m._matrix == NULL) + throw(GeneralClassicException("MMatrix::get_matrix", "Attempt to access uninitialised matrix")); + return (gsl_matrix_complex*)m._matrix; + +} + +////////////////////////// MMatrix inlined functions end ////////////////////// +} + +#endif diff --git a/classic/5.0/src/Fields/Interpolation/MVector.cpp b/classic/5.0/src/Fields/Interpolation/MVector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7ecd6e1d28aaf13e216553a4baa818d8fe607058 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/MVector.cpp @@ -0,0 +1,228 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "Fields/Interpolation/MVector.h" +#include "Fields/Interpolation/MMatrix.h" + +namespace interpolation { +///////////////////////// m_complex //////////////////////////// + +std::ostream& operator<<(std::ostream& out, m_complex c) +{ + out << re(c) << " r " << im(c) << " i"; + return out; +} + +std::istream& operator>>(std::istream& in, m_complex& c) +{ + std::string dummy; + in >> re(c) >> dummy >> im(c) >> dummy; + return in; +} + + +///////////////////////// MVector ////////////////////////////// + + +template <typename Tmplt> +MVector<Tmplt>::MVector( size_t i ) : _vector(NULL) +{ + build_vector(i); +} +template MVector<double> ::MVector(size_t i); +template MVector<m_complex>::MVector(size_t i); + +template <typename Tmplt> +MVector<Tmplt>::MVector( const MVector<Tmplt>& mv) : _vector(NULL) +{ *this = mv; } +template MVector<double> ::MVector(const MVector<double>&); +template MVector<m_complex>::MVector(const MVector<m_complex>&); + + +template <typename Tmplt> +MVector<Tmplt>::MVector( size_t size, Tmplt value ) : _vector(NULL) +{ + build_vector(size); + for(size_t i=0; i<size; i++) operator()(i+1) = value; +} +template MVector<double> ::MVector(size_t i, double value); +template MVector<m_complex>::MVector(size_t i, m_complex value); + + +template <> +void MVector<double>::build_vector ( size_t size ) +{ + if(_vector != NULL) gsl_vector_free((gsl_vector*)_vector); + _vector = gsl_vector_alloc(size); +} + +template <> +void MVector<m_complex>::build_vector( size_t size ) +{ + if(_vector != NULL) gsl_vector_complex_free((gsl_vector_complex*)_vector); + _vector = gsl_vector_complex_alloc(size); +} + +template <class Tmplt> +void MVector<Tmplt>::build_vector ( const Tmplt* data_begin, const Tmplt* data_end ) +{ + build_vector(data_end - data_begin); + for(size_t i=0; i<num_row(); i++) operator()(i+1) = data_begin[i]; +} +template void MVector<double>::build_vector( const double* data_begin, const double* data_end ); +template void MVector<m_complex>::build_vector( const m_complex* data_begin, const m_complex* data_end ); + + +template <class Tmplt> +size_t MVector<Tmplt>::num_row() const +{ if(_vector != NULL) return ((gsl_vector*)_vector)->size; else return 0;} +template size_t MVector<double> ::num_row() const; +template size_t MVector<m_complex>::num_row() const; + +template <> +const double& MVector<double>::operator()(const size_t i) const +{return *gsl_vector_ptr( (gsl_vector*)_vector, i-1);} +template <> +const m_complex& MVector<m_complex>::operator()(const size_t i) const +{return *gsl_vector_complex_ptr((gsl_vector_complex*)_vector, i-1);} +template <> +double& MVector<double>::operator()(const size_t i) +{return *gsl_vector_ptr( (gsl_vector*)_vector, i-1);} +template <> +m_complex& MVector<m_complex>::operator()(const size_t i) +{return *gsl_vector_complex_ptr((gsl_vector_complex*)_vector, i-1);} + + +template <class Tmplt> +MMatrix<Tmplt> MVector<Tmplt>::T() const +{ + MMatrix<Tmplt> mat = MMatrix<Tmplt>(1,num_row()); + for(size_t i=1; i<=num_row(); i++) + mat(1, i) = this->operator()(i); + return mat; +} +template MMatrix<double> MVector<double> ::T() const; +template MMatrix<m_complex> MVector<m_complex>::T() const; + +template <class Tmplt> bool operator==(const MVector<Tmplt>& c1, const MVector<Tmplt>& c2) +{ + if(c1.num_row() != c2.num_row()) + return false; + for(size_t i=0; i<c1.num_row(); i++) + if(c1(i+1) != c2(i+1)) return false; + return true; +} +template bool operator==(const MVector<double>& c1, const MVector<double>& c2); +template bool operator==(const MVector<m_complex>& c1, const MVector<m_complex>& c2); + +template <> +MVector<double>& MVector<double>::operator= (const MVector<double>& mv) +{ + if (&mv == this) return *this; + delete_vector(); + if(!mv._vector) { _vector = NULL; return *this; } + _vector = gsl_vector_alloc( mv.num_row() ); + gsl_vector_memcpy((gsl_vector*)_vector, (const gsl_vector*)mv._vector); + return *this; +} + +template <> +MVector<m_complex>& MVector<m_complex>::operator= (const MVector<m_complex>& mv) +{ + if (&mv == this) return *this; + delete_vector(); + if(!mv._vector) { _vector = NULL; return *this; } + _vector = gsl_vector_complex_alloc( mv.num_row() ); + gsl_vector_complex_memcpy((gsl_vector_complex*)_vector, (const gsl_vector_complex*)mv._vector); + return *this; +} + +template <class Tmplt> std::ostream& operator<<(std::ostream& out, MVector<Tmplt> v) +{ + out << v.num_row() << "\n"; + for(size_t i=0; i<v.num_row(); i++) out << " " << v(i+1) << "\n"; + return out; +} +template std::ostream& operator<<(std::ostream& out, MVector<double> v); +template std::ostream& operator<<(std::ostream& out, MVector<m_complex> v); + +template <class Tmplt> std::istream& operator>>(std::istream& in, MVector<Tmplt>& v) +{ + size_t n; + in >> n; + v = MVector<Tmplt>(n); + for(size_t i=1; i<=v.num_row(); i++) in >> v(i); + return in; +} +template std::istream& operator>>(std::istream& out, MVector<double>& v); +template std::istream& operator>>(std::istream& out, MVector<m_complex>& v); + +const gsl_vector* MVector_to_gsl(const MVector<double>& vd) +{return vd.get_vector(vd);} +const gsl_vector_complex* MVector_to_gsl(const MVector<gsl_complex>& vc) +{return vc.get_vector(vc);} + +template <class Tmplt> +MVector<Tmplt> MVector<Tmplt>::sub(size_t n1, size_t n2) const +{ + MVector<Tmplt> temp(n2-n1+1); + for(size_t i=n1; i<=n2; i++) temp(i-n1+1) = operator()(i); + return temp; +} +template MVector<double> MVector<double> ::sub(size_t n1, size_t n2) const; +template MVector<m_complex> MVector<m_complex>::sub(size_t n1, size_t n2) const; + +MVector<m_complex> complex(MVector<double> real) +{ + MVector<m_complex> c(real.num_row()); + for(size_t i=1; i<=real.num_row(); i++) c(i) = m_complex_build(real(i)); + return c; +} + +MVector<m_complex> complex(MVector<double> real, MVector<double> im) +{ + MVector<m_complex> c(real.num_row()); + for(size_t i=1; i<=real.num_row(); i++) c(i) = m_complex_build(real(i), im(i)); + return c; +} + +MVector<double> re (MVector<m_complex> c) +{ + MVector<double> d(c.num_row()); + for(size_t i=1; i<=c.num_row(); i++) d(i) = re( c(i) ); + return d; +} + +MVector<double> im (MVector<m_complex> c) +{ + MVector<double> d(c.num_row()); + for(size_t i=1; i<=c.num_row(); i++) d(i) = im( c(i) ); + return d; +} + +} + diff --git a/classic/5.0/src/Fields/Interpolation/MVector.h b/classic/5.0/src/Fields/Interpolation/MVector.h new file mode 100644 index 0000000000000000000000000000000000000000..88829aa6c42d581718239d4eeb49b01fb4efb786 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/MVector.h @@ -0,0 +1,273 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef MVector_h +#define MVector_h + +#include <iostream> +#include <vector> + +#include "gsl/gsl_complex_math.h" +#include "gsl/gsl_vector.h" +#include "gsl/gsl_vector_complex_double.h" + +#include "Utilities/GeneralClassicException.h" + + +namespace interpolation { + +/** Wrapper for GSL vector and gsl_complex + * + * MVector class handles maths operators and a few other operators + * + * Use template to define two types: + * (i) MVector<double> is a vector of doubles + * (ii) MVector<m_complex> is a vector of m_complex + * + * Maths operators and a few others are defined, but maths operators + * don't allow operations between types - i.e. you can't multiply + * a complex matrix by a double matrix. Instead use interface methods + * like real() and complex() to convert between types first + */ + + +/////////////// Complex Start ////////////////////// + +typedef gsl_complex m_complex; //typedef because I guess at some point I may want to do something else + +inline m_complex m_complex_build(double r, double i) { m_complex c = {{r,i}}; return c;} +inline m_complex m_complex_build(double r) { m_complex c = {{r,0.}}; return c;} +//Just overload the standard operators +inline m_complex operator *(m_complex c, double d) {return gsl_complex_mul_real(c,d);} +inline m_complex operator *(double d, m_complex c) {return gsl_complex_mul_real(c,d);} +inline m_complex operator *(m_complex c1, m_complex c2) {return gsl_complex_mul(c1,c2);} + +inline m_complex operator /(m_complex c, double d) {return gsl_complex_div_real(c,d);} +inline m_complex operator /(m_complex c1, m_complex c2) {return gsl_complex_div(c1,c2);} +inline m_complex operator /(double d, m_complex c) {m_complex c1 = m_complex_build(d); return gsl_complex_div(c1,c);} + +inline m_complex operator +(m_complex c, double d) {return gsl_complex_add_real(c, d);} +inline m_complex operator +(double d, m_complex c) {return gsl_complex_add_real(c, d);} +inline m_complex operator +(m_complex c1, m_complex c2) {return gsl_complex_add (c1,c2);} + +inline m_complex operator -(m_complex c) {c.dat[0] = -c.dat[0]; c.dat[1] = -c.dat[1]; return c;} +inline m_complex operator -(m_complex c, double d) {return gsl_complex_sub_real(c, d);} +inline m_complex operator -(double d, m_complex c) {return -gsl_complex_sub_real(c, d);} +inline m_complex operator -(m_complex c1, m_complex c2) {return gsl_complex_sub (c1,c2);} + +//suboptimal *= ... should do the substitution in place +inline m_complex& operator *=(m_complex& c, double d) {c = gsl_complex_mul_real(c,d); return c;} +inline m_complex& operator *=(m_complex& c1, m_complex c2) {c1 = gsl_complex_mul(c1,c2); return c1;} + +inline m_complex& operator /=(m_complex& c, double d) {c = gsl_complex_div_real(c,d); return c;} +inline m_complex& operator /=(m_complex& c1, m_complex c2) {c1 = gsl_complex_div(c1,c2); return c1;} + +inline m_complex& operator +=(m_complex& c, double d) {c = gsl_complex_add_real(c, d); return c;} +inline m_complex& operator +=(m_complex& c1, m_complex c2) {c1 = gsl_complex_add (c1,c2); return c1;} + +inline m_complex& operator -=(m_complex& c, double d) {c = gsl_complex_sub_real(c, d); return c;} +inline m_complex& operator -=(m_complex& c1, m_complex c2) {c1 = gsl_complex_sub (c1,c2); return c1;} +//comparators +inline bool operator ==(m_complex c1, m_complex c2) {return c1.dat[0] == c2.dat[0] && c1.dat[1] == c2.dat[1];} +inline bool operator !=(m_complex c1, m_complex c2) {return !(c1==c2);} + +//real and imaginary +inline double& re(m_complex& c) {return c.dat[0];} +inline double& im(m_complex& c) {return c.dat[1];} +inline const double& re(const m_complex& c) {return c.dat[0];} +inline const double& im(const m_complex& c) {return c.dat[1];} + +//complex conjugate +inline m_complex conj(const m_complex& c) {return m_complex_build(re(c), -im(c));} + +std::ostream& operator<<(std::ostream& out, m_complex c); +std::istream& operator>>(std::istream& in, m_complex& c); + + +//////////////// Complex End /////////////////// + +template <class Tmplt> +class MMatrix; //forward declaration because MVector needs to see MMatrix for T() method (and others) + +//////////////// MVector Start /////////////////// + +template <class Tmplt> +class MVector +{ +public: + MVector() : _vector(NULL) {;} + MVector( const MVector<Tmplt>& mv); + MVector( const Tmplt* ta_beg, const Tmplt* ta_end ) : _vector(NULL) { build_vector(ta_beg, ta_end); } //copy from data and put it in the vector + MVector( std::vector<Tmplt> tv) : _vector(NULL) { build_vector(&tv[0], &tv.back()+1); } + MVector( size_t i ); + MVector( size_t i, Tmplt value ); + template <class Tmplt2> MVector(MVector<Tmplt2>); + ~MVector(); + + //return number of elements + size_t num_row() const; + //return reference to ith element; indexing starts from 1, runs to num_row (inclusive); not bound checked + Tmplt& operator()(size_t i); + const Tmplt& operator()(size_t i) const; + + MVector<Tmplt> sub(size_t n1, size_t n2) const; //return a vector running from n1 to n2 (inclusive) + MMatrix<Tmplt> T() const; //return a matrix that is the transpose of the vector + MVector<Tmplt>& operator= (const MVector<Tmplt>& mv); //set *this to be equal to mv + + friend MVector<m_complex>& operator *=(MVector<m_complex>& v, m_complex c); + friend MVector<double>& operator *=(MVector<double>& v, double d); + + friend MVector<m_complex>& operator +=(MVector<m_complex>& v1, MVector<m_complex> v2); + friend MVector<double>& operator +=(MVector<double>& v1, MVector<double> v2); + + friend MVector<m_complex> operator * (MMatrix<m_complex> m, MVector<m_complex> v); + friend MVector<double> operator * (MMatrix<double> m, MVector<double> v); + + friend class MMatrix<Tmplt>; + friend class MMatrix<double>; + +// friend gsl_vector* MVectorToGSL(MVector<double>& ); +// friend gsl_vector_complex* MVectorToGSL(MVector<gsl_complex>&); + friend const gsl_vector* MVector_to_gsl(const MVector<double>& ); + friend const gsl_vector_complex* MVector_to_gsl(const MVector<gsl_complex>&); + +private: + void build_vector ( size_t size ); //copy from data and put it in the vector + void build_vector ( const Tmplt* data_start, const Tmplt* data_end ); //copy from data and put it in the vector + + void delete_vector( ); + + static gsl_vector* get_vector(const MVector<double>& m); + static gsl_vector_complex* get_vector(const MVector<m_complex>& m); + + + void* _vector; +}; + +//Operators do what you would expect - i.e. normal mathematical functions +template <class Tmplt> bool operator ==(const MVector<Tmplt>& v1, const MVector<Tmplt>& v2); +template <class Tmplt> bool operator !=(const MVector<Tmplt>& v1, const MVector<Tmplt>& v2); + +MVector<m_complex>& operator *=(MVector<m_complex>& v, m_complex c); +MVector<double>& operator *=(MVector<double>& v, double d); +template <class Tmplt> MVector<Tmplt> operator *(MVector<Tmplt> v, Tmplt t); +template <class Tmplt> MVector<Tmplt> operator *(Tmplt t, MVector<Tmplt> v); + +template <class Tmplt> MVector<Tmplt>& operator /=(MVector<Tmplt>& v, Tmplt t); +template <class Tmplt> MVector<Tmplt> operator / (MVector<Tmplt> v, Tmplt t); +template <class Tmplt> MVector<Tmplt> operator / (Tmplt t, MVector<Tmplt> v); + +MVector<m_complex>& operator +=(MVector<m_complex>& v1, MVector<m_complex> v2); +MVector<double>& operator +=(MVector<double>& v1, MVector<double> v2); +template <class Tmplt> MVector<Tmplt> operator + (MVector<Tmplt> v1, MVector<Tmplt> v2); + +template <class Tmplt> MVector<Tmplt>& operator -=(MVector<Tmplt>& v1, MVector<Tmplt> v2); +template <class Tmplt> MVector<Tmplt> operator - (MVector<Tmplt> v1, MVector<Tmplt> v2); +template <class Tmplt> MVector<Tmplt> operator - (const MVector<Tmplt>& v); + +//format is +//num_row +//v1 +//v2 +//... +template <class Tmplt> std::ostream& operator<<(std::ostream& out, MVector<Tmplt> v); +template <class Tmplt> std::istream& operator>>(std::istream& out, MVector<Tmplt>& v); + +//Interfaces between MVector<m_complex> and MVector<double> +MVector<m_complex> complex(MVector<double> real); +MVector<m_complex> complex(MVector<double> real, MVector<double> imaginary); +//return vector of doubles filled with either real or imaginary part of mv +MVector<double> re (MVector<m_complex> mv); +MVector<double> im (MVector<m_complex> mv); + +//Interface to gsl +const gsl_vector* MVector_to_gsl(const MVector<double>& vd); +const gsl_vector_complex* MVector_to_gsl(const MVector<gsl_complex>& vc); + +///////////////// MVector End ///////////////// Nb: some inlined functions below... + +//////////////////////////// MVector Inlined Functions ////////////// + +template <class Tmplt> +inline MVector<Tmplt>::~MVector() { delete_vector();} +template <> +inline void MVector<double>::delete_vector() { if(_vector != NULL) gsl_vector_free( (gsl_vector*)_vector); } +template <> +inline void MVector<m_complex>::delete_vector() { if(_vector != NULL) gsl_vector_complex_free( (gsl_vector_complex*)_vector); } + +MVector<m_complex> inline & operator *=(MVector<m_complex>& v, gsl_complex c) +{gsl_vector_complex_scale((gsl_vector_complex*)v._vector, c); return v;} +MVector<double> inline & operator *=(MVector<double>& v, double d) +{gsl_vector_scale((gsl_vector*)v._vector, d); return v;} + +template <class Tmplt> MVector<Tmplt> inline operator *(MVector<Tmplt> v, Tmplt d) +{return v*=d;} +template <class Tmplt> MVector<Tmplt> inline operator *(Tmplt d, MVector<Tmplt> v) +{return v*d;} + +template <class Tmplt> MVector<Tmplt> inline & operator /=(MVector<Tmplt>& v, Tmplt d) +{return v*=1./d;} +template <class Tmplt> MVector<Tmplt> inline operator / (MVector<Tmplt> v, Tmplt d) +{return v/=d;} + +MVector<m_complex> inline & operator +=(MVector<m_complex>& v1, MVector<m_complex> v2) +{gsl_vector_complex_add((gsl_vector_complex*)v1._vector, (gsl_vector_complex*)v2._vector); return v1;} +MVector<double> inline & operator +=(MVector<double>& v1, MVector<double> v2) +{gsl_vector_add((gsl_vector*)v1._vector, (gsl_vector*)v2._vector); return v1;} +template <class Tmplt> MVector<Tmplt> inline operator +(MVector<Tmplt> v1, MVector<Tmplt> v2) +{return v1+=v2;} + +template <class Tmplt> MVector<Tmplt> inline operator - (const MVector<Tmplt>& v) +{MVector<Tmplt> vo(v); for(size_t i=1; i<=vo.num_row(); i++) vo(i) = -vo(i); return vo;} +template <class Tmplt> MVector<Tmplt> inline & operator -=(MVector<Tmplt>& v1, MVector<Tmplt> v2) +{return v1+= -v2;} +template <class Tmplt> MVector<Tmplt> inline operator - (MVector<Tmplt> v1, MVector<Tmplt> v2) +{return v1-=v2;} + +template <class Tmplt> bool inline operator!=(const MVector<Tmplt>& c1, const MVector<Tmplt>& c2) +{return !(c1==c2);} + +template <class Tmplt> +gsl_vector inline* MVector<Tmplt>::get_vector(const MVector<double>& m) +{ + if(m._vector == NULL) + throw(GeneralClassicException("MVector::get_vector", "Attempt to access uninitialised matrix")); + return (gsl_vector*)m._vector; +} + +template <class Tmplt> +gsl_vector_complex inline* MVector<Tmplt>::get_vector(const MVector<m_complex>& m) +{ + if(m._vector == NULL) + throw(GeneralClassicException("MVector::get_vector", "Attempt to access uninitialised vector")); + return (gsl_vector_complex*)m._vector; + +} + +} + +#endif diff --git a/classic/5.0/src/Fields/Interpolation/Makefile b/classic/5.0/src/Fields/Interpolation/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..ece7233be9ee9c26f4d8b65e0bbf4821068aa133 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/Makefile @@ -0,0 +1,156 @@ +# CMAKE generated file: DO NOT EDIT! +# Generated by "Unix Makefiles" Generator, CMake Version 2.8 + +# Default target executed when no arguments are given to make. +default_target: all +.PHONY : default_target + +#============================================================================= +# Special targets provided by cmake. + +# Disable implicit rules so canonical targets will work. +.SUFFIXES: + +# Remove some rules from gmake that .SUFFIXES does not remove. +SUFFIXES = + +.SUFFIXES: .hpux_make_needs_suffix_list + +# Suppress display of executed commands. +$(VERBOSE).SILENT: + +# A target that is always out of date. +cmake_force: +.PHONY : cmake_force + +#============================================================================= +# Set environment variables for the build. + +# The shell in which to execute make rules. +SHELL = /bin/sh + +# The CMake executable. +CMAKE_COMMAND = /usr/bin/cmake + +# The command to remove a file. +RM = /usr/bin/cmake -E remove -f + +# Escaping for special characters. +EQUALS = = + +# The top-level source directory on which CMake was run. +CMAKE_SOURCE_DIR = /home/cr67/OPAL/source/opal_git + +# The top-level build directory on which CMake was run. +CMAKE_BINARY_DIR = /home/cr67/OPAL/source/opal_git + +#============================================================================= +# Targets provided globally by CMake. + +# Special rule for the target edit_cache +edit_cache: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running interactive CMake command-line interface..." + /usr/bin/cmake -i . +.PHONY : edit_cache + +# Special rule for the target edit_cache +edit_cache/fast: edit_cache +.PHONY : edit_cache/fast + +# Special rule for the target install +install: preinstall + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Install the project..." + /usr/bin/cmake -P cmake_install.cmake +.PHONY : install + +# Special rule for the target install +install/fast: preinstall/fast + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Install the project..." + /usr/bin/cmake -P cmake_install.cmake +.PHONY : install/fast + +# Special rule for the target install/local +install/local: preinstall + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing only the local directory..." + /usr/bin/cmake -DCMAKE_INSTALL_LOCAL_ONLY=1 -P cmake_install.cmake +.PHONY : install/local + +# Special rule for the target install/local +install/local/fast: install/local +.PHONY : install/local/fast + +# Special rule for the target list_install_components +list_install_components: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Available install components are: \"Unspecified\"" +.PHONY : list_install_components + +# Special rule for the target list_install_components +list_install_components/fast: list_install_components +.PHONY : list_install_components/fast + +# Special rule for the target rebuild_cache +rebuild_cache: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..." + /usr/bin/cmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) +.PHONY : rebuild_cache + +# Special rule for the target rebuild_cache +rebuild_cache/fast: rebuild_cache +.PHONY : rebuild_cache/fast + +# The main all target +all: cmake_check_build_system + cd /home/cr67/OPAL/source/opal_git && $(CMAKE_COMMAND) -E cmake_progress_start /home/cr67/OPAL/source/opal_git/CMakeFiles /home/cr67/OPAL/source/opal_git/classic/5.0/src/Fields/Interpolation/CMakeFiles/progress.marks + cd /home/cr67/OPAL/source/opal_git && $(MAKE) -f CMakeFiles/Makefile2 classic/5.0/src/Fields/Interpolation/all + $(CMAKE_COMMAND) -E cmake_progress_start /home/cr67/OPAL/source/opal_git/CMakeFiles 0 +.PHONY : all + +# The main clean target +clean: + cd /home/cr67/OPAL/source/opal_git && $(MAKE) -f CMakeFiles/Makefile2 classic/5.0/src/Fields/Interpolation/clean +.PHONY : clean + +# The main clean target +clean/fast: clean +.PHONY : clean/fast + +# Prepare targets for installation. +preinstall: all + cd /home/cr67/OPAL/source/opal_git && $(MAKE) -f CMakeFiles/Makefile2 classic/5.0/src/Fields/Interpolation/preinstall +.PHONY : preinstall + +# Prepare targets for installation. +preinstall/fast: + cd /home/cr67/OPAL/source/opal_git && $(MAKE) -f CMakeFiles/Makefile2 classic/5.0/src/Fields/Interpolation/preinstall +.PHONY : preinstall/fast + +# clear depends +depend: + cd /home/cr67/OPAL/source/opal_git && $(CMAKE_COMMAND) -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 1 +.PHONY : depend + +# Help Target +help: + @echo "The following are some of the valid targets for this Makefile:" + @echo "... all (the default if no target is provided)" + @echo "... clean" + @echo "... depend" + @echo "... edit_cache" + @echo "... install" + @echo "... install/local" + @echo "... list_install_components" + @echo "... rebuild_cache" +.PHONY : help + + + +#============================================================================= +# Special targets to cleanup operation of make. + +# Special rule to run CMake to check the build system integrity. +# No rule that depends on this can have commands that come from listfiles +# because they might be regenerated. +cmake_check_build_system: + cd /home/cr67/OPAL/source/opal_git && $(CMAKE_COMMAND) -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 0 +.PHONY : cmake_check_build_system + diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh-inl.icc b/classic/5.0/src/Fields/Interpolation/Mesh-inl.icc similarity index 99% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh-inl.icc rename to classic/5.0/src/Fields/Interpolation/Mesh-inl.icc index 676821be270c605e386d46cc150d9eec9f45d106..9915abc3111d0340def8c143614da4ad2dff061d 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh-inl.icc +++ b/classic/5.0/src/Fields/Interpolation/Mesh-inl.icc @@ -27,6 +27,7 @@ #include <vector> +namespace interpolation { // inline declarations for Mesh class - keep them in a separate file to keep the // header file clean. @@ -179,4 +180,5 @@ bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) { return false; return true; } +} diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.cpp b/classic/5.0/src/Fields/Interpolation/Mesh.cpp similarity index 96% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.cpp rename to classic/5.0/src/Fields/Interpolation/Mesh.cpp index 58d72f6aeb59e21fb8cc6e3721d0b0c88047c074..6b0a43da283d7818059291c3f60fa11727baa242 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.cpp +++ b/classic/5.0/src/Fields/Interpolation/Mesh.cpp @@ -25,10 +25,11 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/Mesh.h" +#include "Fields/Interpolation/Mesh.h" #include <iomanip> +namespace interpolation { std::ostream& operator<<(std::ostream& out, const Mesh::Iterator& it) { out << std::setw(5) << it.toInteger() << " ** "; for (unsigned int i = 0; i < it.getState().size(); i++) @@ -39,3 +40,5 @@ std::ostream& operator<<(std::ostream& out, const Mesh::Iterator& it) { << it.getPosition()[i] << " "; return out; } +} + diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.h b/classic/5.0/src/Fields/Interpolation/Mesh.h similarity index 99% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.h rename to classic/5.0/src/Fields/Interpolation/Mesh.h index c4a78d5923bbcb971fa49656033aa7d7ae2a971c..bba87cb326e7d06ac23bfe188979b465ceec3a32 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/Mesh.h +++ b/classic/5.0/src/Fields/Interpolation/Mesh.h @@ -31,6 +31,8 @@ #ifndef _CLASSIC_FIELDS_MESH_H_ #define _CLASSIC_FIELDS_MESH_H_ +namespace interpolation { + class VectorMap; /** Mesh Base class for meshing routines @@ -63,7 +65,7 @@ class Mesh { * * Dual is a polyhedron that has centre of each face as a point on the mesh */ - virtual Mesh* dual() = 0; + virtual Mesh* dual() const = 0; /** Destructor - does nothing */ inline virtual ~Mesh(); @@ -338,7 +340,9 @@ inline bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); */ std::ostream& operator<<(std::ostream& out, const Mesh::Iterator& it); -#include "Fields/SectorMagneticFieldMap/Mesh-inl.icc" +} + +#include "Fields/Interpolation/Mesh-inl.icc" #endif diff --git a/classic/5.0/src/Fields/Interpolation/PPSolveFactory.cpp b/classic/5.0/src/Fields/Interpolation/PPSolveFactory.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e4ecccf100872489946f012e14d85d47fcd57702 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PPSolveFactory.cpp @@ -0,0 +1,369 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include <vector> +#include <iostream> +#include <string> +#include <sstream> + +#include <gsl/gsl_sf_pow_int.h> +#include <gsl/gsl_sf_gamma.h> + +#include "Utilities/GeneralClassicException.h" + +#include "Fields/Interpolation/SolveFactory.h" +#include "Fields/Interpolation/PolynomialPatch.h" +#include "Fields/Interpolation/PPSolveFactory.h" + +namespace interpolation { + +typedef PolynomialCoefficient Coeff; + +PPSolveFactory::PPSolveFactory(Mesh* points, + std::vector<std::vector<double> > values, + int polyPatchOrder, + int smoothingOrder) + : polyPatchOrder_m(polyPatchOrder), + smoothingOrder_m(smoothingOrder), + polyDim_m(0), + polyMesh_m(), + points_m(points), + values_m(values), + polynomials_m(), + thisPoints_m(), + thisValues_m(), + derivPoints_m(), + derivValues_m(), + derivIndices_m(), + edgePoints_m(), + smoothingPoints_m() { + if (points == NULL) { + throw GeneralClassicException( + "PPSolveFactory::Solve", + "NULL Mesh input to Solve"); + } + if (points->end().toInteger() != int(values.size())) { + std::stringstream ss; + ss << "Mismatch between Mesh and values in Solve. Mesh size was " + << points->end().toInteger() << " but field values size was " + << values.size(); + throw GeneralClassicException( + "PPSolveFactory::Solve", + ss.str()); + } + if (smoothingOrder < polyPatchOrder) { + throw GeneralClassicException( + "PPSolveFactory::Solve", + "Polynomial order must be <= smoothing order in Solve"); + } + polyMesh_m = points->dual(); + if (polyMesh_m == NULL) + throw GeneralClassicException( + "PPSolveFactory::Solve", + "Dual of Mesh was NULL"); + polyDim_m = values[0].size(); + int posDim = points->getPositionDimension(); + + edgePoints_m = std::vector< std::vector< std::vector<int> > >(posDim+1); + for (int i = 1; i <= posDim; ++i) + edgePoints_m[i] = getNearbyPointsSquares + (i, 0, smoothingOrder-polyPatchOrder); + smoothingPoints_m = getNearbyPointsSquares + (posDim, polyPatchOrder_m-1, polyPatchOrder_m); +} + +std::vector<double> getDeltaPos(Mesh* mesh) { + Mesh::Iterator it = mesh->begin(); + std::vector<double> pos1 = it.getPosition(); + for (size_t i = 0; i < it.getState().size(); ++i) + it[i]++; + std::vector<double> pos2 = it.getPosition(); + for (size_t i = 0; i < pos2.size(); ++i) + pos2[i] -= pos1[i]; + return pos2; +} + +std::vector<double> PPSolveFactory::outOfBoundsPosition( + Mesh::Iterator outOfBoundsIt) { + const Mesh* mesh = outOfBoundsIt.getMesh(); + std::vector<int> state = outOfBoundsIt.getState(); + std::vector<int> delta = std::vector<int>(state.size(), 2); + Mesh::Iterator end = mesh->end()-1; + std::vector<double> pos1 = mesh->begin().getPosition(); + std::vector<double> pos2 = Mesh::Iterator(delta, mesh).getPosition(); + for (size_t i = 0; i < state.size(); ++i) + if (state[i] < 1) { + state[i] = 1; + } else if (state[i] > end[i]) { + state[i] = end[i]; + } + std::vector<double> position = Mesh::Iterator(state, mesh).getPosition(); + state = outOfBoundsIt.getState(); + for (size_t i = 0; i < state.size(); ++i) + if (state[i] < 1) { + position[i] = (pos2[i] - pos1[i])*(state[i]-1) + position[i]; + } else if (state[i] > end[i]) { + position[i] = (pos2[i] - pos1[i])*(state[i]-end[i]) + position[i]; + } + return position; +} + +void PPSolveFactory::getPoints() { + int posDim = points_m->getPositionDimension(); + std::vector< std::vector<int> > dataPoints = + getNearbyPointsSquares(posDim, -1, polyPatchOrder_m); + thisPoints_m = std::vector< std::vector<double> >(dataPoints.size()); + std::vector<double> deltaPos = getDeltaPos(points_m); + // now iterate over the index_by_power_list elements - offset; each of + // these makes a point in the polynomial fit + for (size_t i = 0; i < dataPoints.size(); ++i) { + thisPoints_m[i] = std::vector<double>(posDim); + for (int j = 0; j < posDim; ++j) + thisPoints_m[i][j] = (0.5-dataPoints[i][j])*deltaPos[j]; + } +} + +void PPSolveFactory::getValues(Mesh::Iterator it) { + thisValues_m = std::vector< std::vector<double> >(); + int posDim = it.getState().size(); + std::vector< std::vector<int> > dataPoints = + getNearbyPointsSquares(posDim, -1, polyPatchOrder_m); + size_t dataPointsSize = dataPoints.size(); + Mesh::Iterator end = points_m->end()-1; + // now iterate over the indexByPowerList elements - offset; each of + // these makes a point in the polynomial fit + for (size_t i = 0; i < dataPointsSize; ++i) { + std::vector<int> itState = it.getState(); + for (int j = 0; j < posDim; ++j) + itState[j]++; + Mesh::Iterator itCurrent = Mesh::Iterator(itState, points_m); + bool outOfBounds = false; // element is off the edge of the mesh + for (int j = 0; j < posDim; ++j) { + itCurrent[j] -= dataPoints[i][j]; + outOfBounds = outOfBounds || + itCurrent[j] < 1 || + itCurrent[j] > end[j]; + } + if (outOfBounds) { // if off the edge, then just constrain to zero + thisValues_m.push_back(values_m[it.toInteger()]); + } else { // else fit using values + thisValues_m.push_back(values_m[itCurrent.toInteger()]); + } + } +} + +void PPSolveFactory::getDerivPoints() { + derivPoints_m = std::vector< std::vector<double> >(); + derivIndices_m = std::vector< std::vector<int> >(); + derivOrigins_m = std::vector< std::vector<int> >(); + derivPolyVec_m = std::vector<MVector<double> >(); + int posDim = points_m->getPositionDimension(); + // get the outer layer of points + int deltaOrder = smoothingOrder_m - polyPatchOrder_m; + if (deltaOrder <= 0) + return; + std::vector<double> deltaPos = getDeltaPos(points_m); + // smoothingPoints_m is the last layer of points used for polynomial fit + // walk around smoothingPoints_m finding the points used for smoothing + for (size_t i = 0; i < smoothingPoints_m.size(); ++i) { + // make a list of the axes that are on the edge of the space + std::vector<int> equalAxes; + for (size_t j = 0; j < smoothingPoints_m[i].size(); ++j) + if (smoothingPoints_m[i][j] == polyPatchOrder_m) + equalAxes.push_back(j); + std::vector<std::vector<int> > edgePoints = + edgePoints_m[equalAxes.size()]; + // note the first point, 0,0, is ignored + for (size_t j = 0; j < edgePoints.size(); ++j) { + derivIndices_m.push_back(std::vector<int>(posDim)); + derivOrigins_m.push_back(smoothingPoints_m[i]); + derivPoints_m.push_back(std::vector<double>(posDim)); + + for (size_t k = 0; k < smoothingPoints_m[i].size(); ++k) { + derivPoints_m.back()[k] = (smoothingPoints_m[i][k]-0.5)*deltaPos[k]; + } + for (size_t k = 0; k < edgePoints[j].size(); ++k) { + derivIndices_m.back()[equalAxes[k]] = edgePoints[j][k]; + } + int index = 0; + // potential bug ... derivIndexByPower_m is never used, is it + // needed? + while (true) { + bool isEqual = true; + std::vector<int> polyIndex = + SquarePolynomialVector::IndexByPower(index, posDim); + for (int k = 0; k < posDim && isEqual; ++k) { + isEqual = polyIndex[k] == derivIndices_m.back()[k]; + } + if (isEqual) { + derivIndexByPower_m.push_back(index); + break; + } + ++index; + } + } + } + + int nPolyCoeffs = SquarePolynomialVector::NumberOfPolynomialCoefficients + (posDim, smoothingOrder_m); + for (size_t i = 0; i < derivPoints_m.size(); ++i) { + derivPolyVec_m.push_back(MVector<double>(nPolyCoeffs, 1.)); + // Product[x^{p_j-q_j}*p_j!/(p_j-q_j)!] + for (int j = 0; j < nPolyCoeffs; ++j) { + std::vector<int> pVec = + SquarePolynomialVector::IndexByPower(j, posDim); + std::vector<int> qVec = derivIndices_m[i]; + for (int l = 0; l < posDim; ++l) { + int pow = pVec[l]-qVec[l]; + if (pow < 0) { + derivPolyVec_m.back()(j+1) = 0.; + break; + } + derivPolyVec_m.back()(j+1) *= gsl_sf_fact(pVec[l])/gsl_sf_fact(pow)*gsl_sf_pow_int(-deltaPos[l]/2., pow); + } + } + } +} + +void PPSolveFactory::getDerivs(Mesh::Iterator it) { + derivValues_m = std::vector< std::vector<double> >(derivPoints_m.size()); + + std::vector<double> itPos = it.getPosition(); + int posDim = it.getState().size(); + // get the outer layer of points + Mesh::Iterator end = it.getMesh()->end()-1; + int deltaOrder = smoothingOrder_m - polyPatchOrder_m; + if (deltaOrder <= 0) + return; + for (size_t i = 0; i < derivPoints_m.size(); ++i) { + std::vector<double> point = std::vector<double>(posDim, 0.); + Mesh::Iterator nearest = it; + for (int j = 0; j < posDim; ++j) { + point[j] = itPos[j]-derivPoints_m[i][j]; + } + for (int j = 0; j < posDim; ++j) { + nearest[j] += derivOrigins_m[i][j]; + if (nearest[j] > end[j]) { + nearest = it; + break; + } + } + if (polynomials_m[nearest.toInteger()] == NULL) { + derivValues_m[i] = std::vector<double>(polyDim_m, 0.); + } else { + MMatrix<double> coeffs = + polynomials_m[nearest.toInteger()]->GetCoefficientsAsMatrix(); + MVector<double> values = coeffs*derivPolyVec_m[i]; + derivValues_m[i] = std::vector<double>(polyDim_m); + for(int j = 0; j < posDim; ++j) { + derivValues_m[i][j] = values(j+1); + } + } + } +} + +PolynomialPatch* PPSolveFactory::solve() { + Mesh::Iterator begin = points_m->begin(); + Mesh::Iterator end = points_m->end()-1; + int meshSize = polyMesh_m->end().toInteger(); + polynomials_m = std::vector<SquarePolynomialVector*>(meshSize, NULL); + // get the list of points that are needed to make a given poly vector + getPoints(); + getDerivPoints(); + SolveFactory solver(polyPatchOrder_m, + smoothingOrder_m, + polyMesh_m->getPositionDimension(), + polyDim_m, + thisPoints_m, + derivPoints_m, + derivIndices_m); + int total = (polyMesh_m->end()-1).toInteger(); + double oldPercentage = 0.; + for (Mesh::Iterator it = polyMesh_m->end()-1; + it >= polyMesh_m->begin(); + --it) { + double newPercentage = (total-it.toInteger())/double(total)*100.; + if (newPercentage - oldPercentage > 10.) { + std::cout << " Done " << newPercentage << " %" << std::endl; + oldPercentage = newPercentage; + } + // find the set of values and derivatives for this point in the mesh + getValues(it); + getDerivs(it); + // The polynomial is found using simultaneous equation solve + polynomials_m[it.toInteger()] = solver.PolynomialSolve + (thisValues_m, derivValues_m); + } + return new PolynomialPatch(polyMesh_m, points_m, polynomials_m); +} + +void PPSolveFactory::nearbyPointsRecursive( + std::vector<int> check, + size_t checkIndex, + size_t polyPower, + std::vector<std::vector<int> >& nearbyPoints) { + check[checkIndex] = polyPower; + nearbyPoints.push_back(check); + if (checkIndex+1 == check.size()) + return; + for (int i = 1; i < int(polyPower); ++i) + nearbyPointsRecursive(check, checkIndex+1, i, nearbyPoints); + for (int i = 0; i <= int(polyPower); ++i) { + check[checkIndex] = i; + nearbyPointsRecursive(check, checkIndex+1, polyPower, nearbyPoints); + } +} + +std::vector<std::vector<int> > PPSolveFactory::getNearbyPointsSquares( + int pointDim, + int polyOrderLower, + int polyOrderUpper) { + if (pointDim < 1) + throw(GeneralClassicException("PPSolveFactory::getNearbyPointsSquares", + "Point dimension must be > 0")); + if (polyOrderLower > polyOrderUpper) + throw(GeneralClassicException("PPSolveFactory::getNearbyPointsSquares", + "Polynomial lower bound must be <= polynomial upper bound")); + // polyOrder -1 (or less) means no terms + if (polyOrderUpper == polyOrderLower || polyOrderUpper < 0) + return std::vector<std::vector<int> >(); + // polyOrder 0 means constant term + std::vector<std::vector<int> > nearbyPoints(1, std::vector<int>(pointDim, 0)); + // polyOrder > 0 means corresponding poly terms (1 is linear, 2 is quadratic...) + for (size_t thisPolyOrder = 0; + thisPolyOrder < size_t(polyOrderUpper); + ++thisPolyOrder) + nearbyPointsRecursive(nearbyPoints[0], 0, thisPolyOrder+1, nearbyPoints); + if (polyOrderLower < 0) + return nearbyPoints; // no terms in lower bound + int nPointsLower = 1; + for (int dim = 0; dim < pointDim; ++dim) + nPointsLower *= polyOrderLower+1; + nearbyPoints.erase(nearbyPoints.begin(), nearbyPoints.begin()+nPointsLower); + return nearbyPoints; +} +} diff --git a/classic/5.0/src/Fields/Interpolation/PPSolveFactory.h b/classic/5.0/src/Fields/Interpolation/PPSolveFactory.h new file mode 100644 index 0000000000000000000000000000000000000000..3d1c5c9cb1a8259e57ecccc05f13789cf99c482e --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PPSolveFactory.h @@ -0,0 +1,160 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include <vector> + +#include "Fields/Interpolation/Mesh.h" +#include "Fields/Interpolation/PolynomialCoefficient.h" +#include "Fields/Interpolation/SquarePolynomialVector.h" + +namespace interpolation { + +class PolynomialPatch; + +/** \class[PPSolveFactory] solves the system of linear equations to interpolate + * from a grid of points using higher order polynomials, creating a + * PolynomialPatch object. The order of the polynomial is user-defined. + * + * The PPSolveFactory can be used to match nearby points for smoothing or for + * straight fitting. + * + * PPSolveFactory generates polynomials on a grid that is at the centre point + * of each of the mesh squares; PPSolveFactory picks nearby points in a + * square around the PPSolveFactory, picked by getNearbyPointsSquares. Two sets + * of points are chosen; the nearest points are used for fitting; smoothing is + * performed by choosing the points on the edge of the grid, and smoothing to + * derivatives on these points. + * + * Because we work in a square grid, we have matching polynomials and + * derivatives; so e.g. 2D points are (0, 0), (1, 0), (0, 1), (1, 1) so + * corresponding polynomial coefficients that we solve for are + * x^0 y^0, x^1 y^0, x^0 y^1, x^1 y^1; note that hence we define "first order" + * to be all terms with no single x/y vector power > 1; i.e. "first order" in + * x OR first order in y, etc. We allow products so long as the maximum power + * in any one dimension <= 1 + * + * The PPSolveFactory sits on top of SolveFactory; PPSolveFactory has the job + * of picking values and derivatives for fitting, SolveFactory then does the + * actual solve for each individual polynomial. + */ + +class PPSolveFactory { + public: + /** Constructor + * + * \param[points] Set of points on which values are stored. Must be a + * rectangular grid. PPSolveFactory takes ownership of + * points (will delete on exit). + * \param[values] Set of values to which we fit. Must be one value per mesh + * point and each value must have the same size. + * \param[polyPatchOrder] The order of the fitted part of the polynomial. + * \param[smoothingOrder] The total order of the fitted and the smoothed + * part of the polynomial. Smoothing order should always be + * >= polyPatchOrder. So for example if polyPatchOrder + * is 2 and smoothing order is 2, we get a 2nd order + * function with no derivative matching. If polyPatchOrder + * is 2 and smoothing order is 3, we get a 3rd order + * function with derivative matching in first order at two + * mesh points from the centre. + */ + PPSolveFactory(Mesh* points, + std::vector<std::vector<double> > values, + int polyPatchOrder, + int smoothingOrder); + + /** Destructor deletes points_ */ + ~PPSolveFactory() {} + + /** Solve the system of equations to generate a PolynomialPatch object. + * + * Returns a PolynomialPatch object, caller owns the returned memory. + */ + PolynomialPatch* solve(); + + /** Get nearby points in a square pattern + * + * Get nearby points at least poly_order_lower distance away and less than + * poly_order_upper distance away; e.g. getNearbyPointsSquares(2, 4, 6) + * will get nearby points on a 2D grid at least 4 grid points and at most + * 6 grid points, like: + * + * x x x o o x ... + * x x x o o x + * x x x o o x + * o o o o o x + * o o o o o x + * x x x x x x + * ... + */ + static std::vector<std::vector<int> > getNearbyPointsSquares( + int pointDim, + int polyOrderLower, + int polyOrderUpper); + + private: + void getPoints(); + void getValues(Mesh::Iterator it); + void getDerivPoints(); + void getDerivs(Mesh::Iterator it); + + // nothing calls this method but I don't quite field brave enought to remove + // it... + std::vector<double> outOfBoundsPosition(Mesh::Iterator outOfBoundsIt); + static void nearbyPointsRecursive( + std::vector<int> check, + size_t checkIndex, + size_t polyPower, + std::vector<std::vector<int> >& nearbyPoints); + + PolynomialCoefficient getDeltaIterator(Mesh::Iterator it1, + Mesh::Iterator it2, int valueIndex); + + int polyPatchOrder_m; + int smoothingOrder_m; + int polyDim_m; + int valueDim_m; + Mesh* polyMesh_m; + Mesh* points_m; + std::vector<std::vector<double> > values_m; + std::vector<SquarePolynomialVector*> polynomials_m; + + std::vector< std::vector<double> > thisPoints_m; + std::vector< std::vector<double> > thisValues_m; + std::vector< std::vector<double> > derivPoints_m; + std::vector< std::vector<double> > derivValues_m; + std::vector< std::vector<int> > derivOrigins_m; + std::vector< std::vector<int> > derivIndices_m; + std::vector< MVector<double> > derivPolyVec_m; + std::vector<int> derivIndexByPower_m; + + std::vector<std::vector<std::vector<int> > > edgePoints_m; + std::vector< std::vector<int> > smoothingPoints_m; + +}; + +} + diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialCoefficient.cpp b/classic/5.0/src/Fields/Interpolation/PolynomialCoefficient.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c2d824b31b8b14e0d1bc0801c7fd1218d9e1eed3 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialCoefficient.cpp @@ -0,0 +1,43 @@ +namespace interpolation { + +std::ostream& operator << (std::ostream& out, + const PolynomialCoefficient& coeff) { + std::vector<int> in_var = coeff.InVariables(); + for (size_t i = 0; i < coeff.InVariables().size(); ++i) + out << in_var[i] << " "; + out << "| " << coeff.OutVariable(); + out << " | " << coeff.Coefficient(); + return out; +} + +void PolynomialCoefficient::SpaceTransform(std::vector<int> space_in, + std::vector<int> space_out) { + std::map<int, int> mapping; //probably optimise this + for (unsigned int i = 0; i < space_out.size(); i++) + for (unsigned int j = 0; j < space_in.size(); j++) + if (space_out[i] == space_in[j]) + mapping[j] = i; //mapping[space_in_index] returns space_out_index + + std::vector<int> in_variables(_inVarByVec.size()); + for (unsigned int con=0; con<in_variables.size(); con++) + { + if (mapping.find(in_variables[con]) != mapping.end()) + in_variables[con] = mapping[in_variables[con]]; + else + throw(GeneralClassicException( + "PolynomialVector::PolynomialCoefficient::SpaceTransform", + "Input variable not found in space transform" + )); + } + + if(mapping.find(_outVar) != mapping.end()) + _outVar = mapping[_outVar]; + else + throw(GeneralClassicException( + "PolynomialVector::PolynomialCoefficient::SpaceTransform", + "Output variable not found in space transform" + )); + _inVarByVec = in_variables; +} +} + diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialCoefficient.h b/classic/5.0/src/Fields/Interpolation/PolynomialCoefficient.h new file mode 100644 index 0000000000000000000000000000000000000000..6c99f55774ab98e7365658e33b9229d4a4b23dfb --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialCoefficient.h @@ -0,0 +1,106 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef PolynomialCoefficient_h +#define PolynomialCoefficient_h + +namespace interpolation { + +/** \class PolynomialCoefficient represents a coefficient in a multi-dimensional + * polynomial + * + * PolynomialCoefficient has three member data + * \param[inVarByVec_m] is the x power to which the coefficient pertains, + * indexed by vector e.g. \f$x_1^4 x_2^3 =\f$ {1,1,1,1,2,2,2} + * \param[outVar_m] is the output y variable to which the coefficient pertains + * \param[coefficient] is the value of the coefficient + */ +class PolynomialCoefficient +{ + public: + /** Construct the coefficient + * \param[inVariablesByVector] x power indexed like e.g. + * \param[outVariable] y index that the coefficient pertains to + * \param[coefficient] value of the coefficient + */ + PolynomialCoefficient(std::vector<int> inVariablesByVector, + int outVariable, + double coefficient) + : _inVarByVec(inVariablesByVector), _outVar(outVariable), + _coefficient(coefficient) { + } + + + /** Copy constructor */ + PolynomialCoefficient(const PolynomialCoefficient& pc) + : _inVarByVec(pc._inVarByVec), _outVar(pc._outVar), + _coefficient(pc._coefficient) { + } + + /** Return the vector of input variables, indexed by vector */ + std::vector<int> InVariables() const {return _inVarByVec;} + + /** Return the output variable */ + int OutVariable() const {return _outVar;} + + /** Return the coefficient value */ + double Coefficient() const {return _coefficient;} + + /** Set the vector of input variables, indexed by vector */ + std::vector<int> InVariables(std::vector<int> inVar ) {_inVarByVec = inVar; return _inVarByVec;} + + /** Set the output variable */ + int OutVariable(int outVar) {_outVar = outVar; return _outVar;} + /** Set the coefficient */ + double Coefficient(double coeff ) {_coefficient = coeff; return _coefficient;} + + /** Transform coefficient from subspace space_in to subspace space_out, both + * subspaces of some larger space + * + * \param{spaceIn} Describes the input subspace + * \param{spaceOut} Describes the output subspace + * Throws an exception if the coefficient is not in the input or output + * subspace + * So for coeff({1,2},0,1.1), coeff.space_transform({0,2,3,5}, {4,7,1,2,3,0}) + * would return coeff({3,4},5,1.1) + */ + void SpaceTransform(std::vector<int> spaceIn, + std::vector<int> spaceOut); + /** Return true if var is in inVariables */ + bool HasInVariable(int var) const {for(unsigned int i=0; i<_inVarByVec.size(); i++) if(_inVarByVec[i] == var) return true; return false;} + + private: + std::vector<int> _inVarByVec; + int _outVar; + double _coefficient; +}; + +std::ostream& operator << (std::ostream&, const PolynomialCoefficient&); + +} + +#endif diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialPatch.cpp b/classic/5.0/src/Fields/Interpolation/PolynomialPatch.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5600e6c56f937a3f249fb3f756f4d66cfa2f02a2 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialPatch.cpp @@ -0,0 +1,135 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "Utilities/GeneralClassicException.h" + +#include "Fields/Interpolation/SolveFactory.h" +#include "Fields/Interpolation/PolynomialPatch.h" + +namespace interpolation { +PolynomialPatch::PolynomialPatch(Mesh* grid_points, + Mesh* validity_region, + std::vector<SquarePolynomialVector*> polynomials) { + grid_points_ = grid_points; + validity_region_ = validity_region; + points_ = polynomials; + // validate input data + if (grid_points_ == NULL) + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "PolynomialPatch grid_points_ was NULL"); + if (validity_region_ == NULL) + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "PolynomialPatch validity_region_ was NULL"); + if (points_.size() == 0) { + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "Could not make PolynomialPatch with 0 length polynomials vector"); + } + Mesh::Iterator end = grid_points_->end(); + if (int(points_.size()) != end.toInteger()) { + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "Could not make PolynomialPatch with bad length polynomials vector" + ); + } + point_dimension_ = grid_points_->getPositionDimension(); + if (validity_region_->getPositionDimension() != point_dimension_) { + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "PolynomialPatch validity_region_ has bad dimension"); + } + value_dimension_ = points_[0]->ValueDimension(); + for (size_t i = 0; i < points_.size(); ++i) { + if (points_[i] == NULL) + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "PolynomialPatch points_ element was NULL"); + if (points_[i]->PointDimension() != point_dimension_) { + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "Polynomial with mismatched PointDimension in PolynomialPatch"); + } + if (points_[i]->ValueDimension() != value_dimension_) + throw GeneralClassicException( + "PolynomialPatch::PolynomialPatch", + "Polynomial with mismatched ValueDimension in PolynomialPatch"); + } +} + +VectorMap* PolynomialPatch::clone() const { + Mesh* new_mesh = NULL; + if (grid_points_ != NULL) { + new_mesh = grid_points_->clone(); + } + Mesh* new_validity = NULL; + if (validity_region_ != NULL) { + new_validity = validity_region_->clone(); + } + std::vector<SquarePolynomialVector*> new_points(points_.size()); + for (size_t i = 0; i < points_.size(); ++i) { + if (points_[i] == NULL) { + new_points[i] = NULL; + } else { + new_points[i] = new SquarePolynomialVector(*points_[i]); + } + } + return new PolynomialPatch(new_mesh, new_validity, new_points); +} + +PolynomialPatch::PolynomialPatch() + : grid_points_(NULL), validity_region_(NULL), points_(), point_dimension_(0), + value_dimension_(0) { +} + +PolynomialPatch::~PolynomialPatch() { + if (grid_points_ != NULL) + delete grid_points_; + if (validity_region_ != NULL) + delete validity_region_; + for (size_t i = 0; i < points_.size(); ++i) + delete points_[i]; +} + +void PolynomialPatch::function(const double* point, double* value) const { + Mesh::Iterator nearest = grid_points_->getNearest(point); + std::vector<double> point_temp(point_dimension_); + std::vector<double> nearest_pos = nearest.getPosition(); + for (size_t i = 0; i < point_dimension_; ++i) + point_temp[i] = point[i] - nearest_pos[i]; + int points_index = nearest.toInteger(); + points_[points_index]->F(&point_temp[0], value); +} + +SquarePolynomialVector* PolynomialPatch::getPolynomialVector(const double* point) const { + Mesh::Iterator nearest = grid_points_->getNearest(point); + int points_index = nearest.toInteger(); + return points_[points_index]; +} +} + diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialPatch.h b/classic/5.0/src/Fields/Interpolation/PolynomialPatch.h new file mode 100644 index 0000000000000000000000000000000000000000..484845c98664f59069fc4dfa10501a26511ca1a7 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialPatch.h @@ -0,0 +1,138 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "Fields/Interpolation/SquarePolynomialVector.h" +#include "Fields/Interpolation/VectorMap.h" +#include "Fields/Interpolation/Mesh.h" + +#ifndef PolynomialPatch_hh +#define PolynomialPatch_hh + +namespace interpolation { + +/** \class[PolynomialPatch] patches together many SquarePolynomialVectors to + * make a multidimensional polynomial spline. + * + * SquarePolynomialVectors are distributed on a mesh. + * \param[grid_points_] The grid points on which patches are defined. + * SquarePolynomialVector owns this memory (deletes on destructor). + * \param[validity_region_] Defines whether a region is valid or not. This is + * defined as a mesh for compatibility with the + * SectorMagneticFieldMap (but possibly this is the wrong way to do + * it - better as a lookup, like bool is_valid(point), which could be + * inherited from the VectorMap class) + * SquarePolynomialVector owns this memory (deletes on destructor). + * \param[points_] The PolynomialVector data + * SquarePolynomialVector owns this memory (deletes on destructor). + * \param[point_dimension_] Dimension of the ordinate + * \param[value_dimension_] Dimension of the abscissa + */ + + +class PolynomialPatch : public VectorMap { + public: + /** Construct the PolynomialPatch + * + * \param[grid_points_] the mesh over which the polynomials are + * distributed. PolynomialPatch takes ownership of this memory. + * \param[validity_region_] the mesh over which the grid points are valid; + * the validity region is defined by min/max values of the grid. + * PolynomialPatch takes ownership of this memory. + * \param[polynomials_] the set of polynomials which make up the patch. + * PolynomialPatch takes ownership of this memory. + */ + PolynomialPatch(Mesh* grid_points_, + Mesh* validity_region, + std::vector<SquarePolynomialVector*> polynomials_); + + /** Default constructor leaves validity_region_ and grid_points_ as NULL + * + * Note it is not safe to call some member functions e.g. function(...) on + * default constructed PolynomialPatch. + */ + PolynomialPatch(); + + /** Destructor delets validity_region_, grid_points_ and points_ */ + ~PolynomialPatch(); + + /** Deep copy the PolynomialPatch + * + * This is the polymorphic copy constructor. Always returns a + * PolynomialPatch. + */ + VectorMap* clone() const; //copy function + + /** Get the value at a given point + * + * \param[point] array of length point_dimension_. Caller owns this memory. + * Not bound checked - it wants to be fast. + * \param[value] array of length value_dimension_. Data gets overwritten by + * PolynomialPatch. Caller owns this memory. Not bound checked. + */ + virtual void function(const double* point, double* value) const; + + /** Get the point dimension (length of the ordinate) */ + inline unsigned int getPointDimension() const {return point_dimension_;} + + /** Get the value dimension (length of the abscissa) */ + inline unsigned int getValueDimension() const {return value_dimension_;} + + /** Get the nearest SquarePolynomialVector to point + * + * \param[point] array of length point_dimension_. Caller owns this memory. + * Not bound checked. + */ + SquarePolynomialVector* getPolynomialVector(const double* point) const; + + /** Get the validity region mesh + * + * Note that the PolynomialVectors are distributed on the dual of this mesh + */ + Mesh* getMesh() const {return validity_region_;} + + /** Get all polynomials in the PolynomialPatch + * + * This is a borrowed reference - PolynomialPatch still owns this memory. + */ + std::vector<SquarePolynomialVector*> getPolynomials() const {return points_;} + + private: + Mesh* validity_region_; + Mesh* grid_points_; + std::vector<SquarePolynomialVector*> points_; + unsigned int point_dimension_; + unsigned int value_dimension_; + + PolynomialPatch(const PolynomialPatch& rhs); + PolynomialPatch& operator=(const PolynomialPatch& rhs); +}; + +} + +#endif // PolynomialPatch_hh + + diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialVector.cpp b/classic/5.0/src/Fields/Interpolation/PolynomialVector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..56f06a99660e419a87dddf9249c750510789879b --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialVector.cpp @@ -0,0 +1,589 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + + + +#include <iomanip> +#include <sstream> +#include <math.h> + +#include "gsl/gsl_sf_gamma.h" + +#include "Utilities/GeneralClassicException.h" + +#include "Fields/Interpolation/MMatrix.h" +#include "Fields/Interpolation/MVector.h" + +#include "Fields/Interpolation/PolynomialVector.h" + + +namespace interpolation { + +////////// POLYNOMIALVECTOR START /////////// + +bool PolynomialVector::_printHeaders=true; + +PolynomialVector::PolynomialVector() + : _pointDim(0), _polyKeyByPower(), + _polyKeyByVector(), _polyCoeffs(), + _polyVector() { +} + +PolynomialVector::PolynomialVector (const PolynomialVector& pv) + : _pointDim(pv._pointDim), _polyKeyByPower(pv._polyKeyByPower), + _polyKeyByVector(pv._polyKeyByVector), + _polyCoeffs(pv._polyCoeffs.num_row(), pv._polyCoeffs.num_col(), 0.), + _polyVector(pv._polyVector) { + SetCoefficients(pv._polyCoeffs); +} + + +PolynomialVector::PolynomialVector(int numberOfInputVariables, MMatrix<double> polynomialCoefficients) + : _pointDim(numberOfInputVariables), _polyKeyByPower(), _polyKeyByVector(), _polyCoeffs(polynomialCoefficients) +{ + SetCoefficients(numberOfInputVariables, polynomialCoefficients); +} + +PolynomialVector::PolynomialVector(std::vector<PolynomialCoefficient> coefficients) + : _pointDim(0), _polyKeyByPower(), _polyKeyByVector(), _polyCoeffs() +{ + SetCoefficients(coefficients); +} + +void PolynomialVector::SetCoefficients(int pointDim, MMatrix<double> coeff) +{ + int nPCoeff = coeff.num_col(); + _pointDim = pointDim; + _polyKeyByPower = std::vector< std::vector<int> >(); + _polyKeyByVector = std::vector< std::vector<int> >(); + _polyVector = std::vector<PolynomialCoefficient>(); + + for(int i=0; i<nPCoeff; i++) + _polyKeyByPower.push_back(IndexByPower (i, pointDim)); + for(int i=0; i<nPCoeff; i++) + _polyKeyByVector.push_back(IndexByVector(i, pointDim)); + + for(size_t i=0; i<_polyCoeffs.num_row(); i++) + for(int j=0; j<nPCoeff; j++) + _polyVector.push_back( PolynomialVector::PolynomialCoefficient(_polyKeyByVector[j], i, _polyCoeffs(i+1, j+1) ) ); +} + +void PolynomialVector::SetCoefficients(std::vector<PolynomialCoefficient> coeff) +{ + _pointDim = 0; + _polyKeyByPower = std::vector< std::vector<int> >(); + _polyKeyByVector = std::vector< std::vector<int> >(); + _polyVector = coeff; + + int maxPolyOrder = 0; + int valueDim = 0; + for(unsigned int i=0; i<coeff.size(); i++) { + int polyOrder = coeff[i].InVariables().size(); + for(unsigned int j=0; j<coeff[i].InVariables().size(); j++) + if(coeff[i].InVariables()[j] > _pointDim) _pointDim = coeff[i].InVariables()[j]; + if(coeff[i].OutVariable() > valueDim) valueDim = coeff[i].OutVariable(); + if(polyOrder > maxPolyOrder) maxPolyOrder = polyOrder; + } + _pointDim++; //PointDim indexes from 0 + _polyCoeffs = MMatrix<double>(valueDim+1, NumberOfPolynomialCoefficients(_pointDim, maxPolyOrder+1), 0.); + + for(size_t i=0; i<_polyCoeffs.num_col(); i++) + _polyKeyByPower.push_back(IndexByPower (i, _pointDim)); + for(size_t i=0; i<_polyCoeffs.num_col(); i++) + { + _polyKeyByVector.push_back(IndexByVector(i, _pointDim)); + for(unsigned int j=0; j<coeff.size(); j++) + if(coeff[j].InVariables() == _polyKeyByVector.back()) + { + int dim = coeff[j].OutVariable(); + _polyCoeffs(dim+1,i+1) = coeff[j].Coefficient(); + } + } +} + +void PolynomialVector::SetCoefficients(MMatrix<double> coeff) +{ + for(size_t i=0; i<coeff.num_row() && i<_polyCoeffs.num_row(); i++) + for(size_t j=0; j<coeff.num_col() && j<_polyCoeffs.num_col(); j++) + _polyCoeffs(i+1,j+1) = coeff(i+1,j+1); +} + +PolynomialVector* PolynomialVector::Recentred(double * point) const +{ + throw(GeneralClassicException("PolynomialVector::Recentred", "Recentred not implemented")); +} + +void PolynomialVector::F(const double* point, double* value) const +{ + MVector<double> pointV(PointDimension(), 1), valueV(ValueDimension(), 1); + for(unsigned int i=0; i<PointDimension(); i++) pointV(i+1) = point[i]; + F(pointV, valueV); + for(unsigned int i=0; i<ValueDimension(); i++) value[i] = valueV(i+1); +} + +void PolynomialVector::F(const MVector<double>& point, MVector<double>& value) const +{ + MVector<double> polyVector(_polyKeyByVector.size(), 1); + MakePolyVector(point, polyVector); + MVector<double> answer = _polyCoeffs*polyVector; + for(unsigned int i=0; i<ValueDimension(); i++) value(i+1) = answer(i+1); +} + + +MVector<double>& PolynomialVector::MakePolyVector(const MVector<double>& point, MVector<double>& polyVector) const +{ + for(unsigned int i=0; i<_polyKeyByVector.size(); i++) + for(unsigned int j=0; j<_polyKeyByVector[i].size(); j++) + polyVector(i+1) *= point( _polyKeyByVector[i][j]+1 ); + return polyVector; +} + +double* PolynomialVector::MakePolyVector(const double* point, double* polyVector) const +{ + for(unsigned int i=0; i<_polyKeyByVector.size(); i++) + { + polyVector[i] = 1.; + for(unsigned int j=0; j<_polyKeyByVector[i].size(); j++) + polyVector[i] *= point[ _polyKeyByVector[i][j] ]; + } + return polyVector; +} + +double* PolynomialVector::MakeDerivVector(const double* positions, const int* deriv_indices, double* deriv_vec) const { + size_t pkey_size = _polyKeyByPower.size(); + for (size_t i = 0; i < pkey_size; ++i) { + deriv_vec[i] = 1.; + for (int j = 0; j < _pointDim; ++j) { + int power = _polyKeyByPower[i][j] - deriv_indices[j]; + if (power < 0) { + deriv_vec[i] = 0.; + } else { + // x^(n-m) + for (int k = 0; k < power; ++k) + deriv_vec[i] *= positions[j]; + } + // n*(n-1)*(n-2)*...*(n-m) + for (int k = _polyKeyByPower[i][j]; k > 0 && k > power; --k) { + deriv_vec[i] *= k; + } + } + } + return deriv_vec; +} + +//Turn int index into a std::vector<int> 'a' of length 'd' with values x_1^a_1 x_2^a_2 ... x_d^a_d +// e.g. x_1^4 x_2^3 = {4,3,0,0} +std::vector<int> PolynomialVector::IndexByPower(int index, int nInputVariables) +{ + std::vector<int> powerIndex(nInputVariables, 0); + std::vector<int> vectorIndex = IndexByVector(index, nInputVariables); + for(unsigned int i=0; i<vectorIndex.size(); i++) + powerIndex[vectorIndex[i]]++; + return powerIndex; +} + +//Turn int index into an index for element of polynomial +// e.g. x_1^4 x_2^3 = {1,1,1,1,2,2,2} +std::vector<int> PolynomialVector::IndexByVector(int index, int nInputVariables) +{ + if(index == 0) return std::vector<int>(); + std::vector<int> indices(1,-1); + nInputVariables--; + for(int i=0; i<index; i++) + { + if (indices.front() == nInputVariables) { indices = std::vector<int>(indices.size()+1, 0); indices.back()--;} + if (indices.back() == nInputVariables) + { + int j=indices.size()-1; + while(indices[j] == nInputVariables) j--; + for(int k=indices.size()-1; k>=j; k--) indices[k] = indices[j]+1; + } + else indices.back()++; + } + return indices; +} + +unsigned int PolynomialVector::NumberOfPolynomialCoefficients(int pointDimension, int order) +{ + int n=0; + if(order<=0) return 0; + for(int i=1; i<order; i++) + n += gsl_sf_choose(pointDimension+i-1, i); + return n+1; +} + +std::ostream& operator<<(std::ostream& out, const PolynomialVector& pv) +{ + if(PolynomialVector::_printHeaders) pv.PrintHeader(out, '.', ' ', 14, true); + out << '\n' << pv.GetCoefficientsAsMatrix(); + return out; +} + +std::ostream& operator << (std::ostream& out, const PolynomialVector::PolynomialCoefficient& coeff) { + std::vector<int> in_var = coeff.InVariables(); + for (size_t i = 0; i < coeff.InVariables().size(); ++i) + out << in_var[i] << " "; + out << "| " << coeff.OutVariable(); + out << " | " << coeff.Coefficient(); + return out; +} + +template <class Container > +void PolynomialVector::PrintContainer(std::ostream& out, const Container& container, char T_separator, char str_separator, int length, bool pad_at_start) +{ +// class Container::iterator it; + std::stringstream strstr1(""); + std::stringstream strstr2(""); + class Container::const_iterator it1 = container.begin(); + class Container::const_iterator it2 = it1; + while(it1!=container.end()) + { + it2++; + if(it1 != container.end() && it2 != container.end()) + strstr1 << (*it1) << T_separator; + else strstr1 << (*it1); + it1++; + } + + if(int(strstr1.str().size()) > length && length > 0) strstr2 << strstr1.str().substr(1, length); + else strstr2 << strstr1.str(); + + if(!pad_at_start) out << strstr2.str(); + for(int i=0; i<length - int(strstr2.str().size()); i++) out << str_separator; + if(pad_at_start) out << strstr2.str(); +} + +void PolynomialVector::PrintHeader(std::ostream& out, char int_separator, char str_separator, int length, bool pad_at_start) const +{ + if(_polyKeyByPower.size() > 0) PrintContainer< std::vector<int> >(out, _polyKeyByPower[0], int_separator, str_separator, length-1, pad_at_start); + for(unsigned int i=1; i<_polyKeyByPower.size(); i++) + PrintContainer<std::vector<int> >(out, _polyKeyByPower[i], int_separator, str_separator, length, pad_at_start); +} + +MVector<double> PolynomialVector::Means(std::vector<std::vector<double> > values, std::vector<double> weights) +{ + if(values.size() < 1) throw(GeneralClassicException("PolynomialVector::Means", "No input values")); + if(values[0].size() < 1) throw(GeneralClassicException("PolynomialVector::Means", "Dimension < 1")); + if(weights.size() != values.size()) weights = std::vector<double>(values.size(), 1.); + int npoints = values.size(); + int dim = values[0].size(); + MVector<double> means(dim,0); + double totalWeight = 0; + for(int x=0; x<npoints; x++) totalWeight += weights[x]; + for(int x=0; x<npoints; x++) + for(int i=0; i<dim; i++) + { + means(i+1) += values[x][i]*weights[x]/totalWeight; + } + + return means; +} + +MMatrix<double> PolynomialVector::Covariances(std::vector<std::vector<double> > values, std::vector<double> weights, MVector<double> means) +{ + size_t npoints = values.size(); + if(npoints < 1) throw(GeneralClassicException("PolynomialVector::Covariances()", "No input values")); + size_t dim = values[0].size(); + if(dim < 1) throw(GeneralClassicException("PolynomialVector::Covariances()", "Dimension < 1")); + if(means.num_row() != dim) means = Means(values, weights); + if(weights.size() != values.size()) weights = std::vector<double>(values.size(), 1.); + + MMatrix<double> cov(dim,dim,0.); + + double totalWeight = 0; + for(size_t x=0; x<npoints; x++) + totalWeight += weights[x]; + + for(size_t x=0; x<npoints; x++) + for(size_t i=0; i<dim; i++) + for(size_t j=0; j<dim; j++) //make a sym matrix for speed + cov(i+1,j+1) += values[x][i]*values[x][j]*weights[x]/totalWeight; + + for(size_t i=0; i<dim; i++) + for(size_t j=0; j<dim; j++) + cov(i+1,j+1) -= means(i+1)*means(j+1); + + return cov; +} + +void PolynomialVector::FDeriv(const double* point, const int* derivative_by_power, double* value) const { + MVector<double> polyVector(_polyKeyByVector.size(), 1); + MakeDerivVector(point, derivative_by_power, &polyVector(1)); + MVector<double> answer = _polyCoeffs*polyVector; + for(unsigned int i=0; i<ValueDimension(); i++) + value[i] = answer(i+1); +} + + + +PolynomialVector* PolynomialVector::PolynomialLeastSquaresFit(const std::vector< std::vector<double> >& points, const std::vector< std::vector<double> >& values, + int polynomialOrder, const std::vector<double>& weights) +{ + //Algorithm: We have F2 = sum_i ( f_k f_l) where f are polynomial terms; FY = sum_i (f_) + + int pointDim = points[0].size(); + int valueDim = values[0].size(); + int nPoints = points.size(); + int nCoeffs = NumberOfPolynomialCoefficients(pointDim, polynomialOrder); + + MMatrix<double> Fy(nCoeffs, valueDim, 0.); + MMatrix<double> F2(nCoeffs, nCoeffs, 0.); + + MMatrix<double> dummy(valueDim, nCoeffs, 0.); + for(int i=0; i<valueDim; i++) for(int j=0; j<nCoeffs; j++) dummy(i+1, j+1) = 1; + PolynomialVector* temp = new PolynomialVector(pointDim, dummy); + + std::vector<double> tempFx(nCoeffs,0); //tempfx = x_i^j + std::vector<double> wt (nPoints,1); + if(weights.size() > 0) wt = weights; + + //sum over points, values + for(int i=0; i<nPoints; i++) + { + temp->MakePolyVector(&points[i][0], &tempFx[0]); + for(int k=0; k<nCoeffs; k++) + for(int j=0; j<nCoeffs; j++) + F2(j+1, k+1) += tempFx[k]*tempFx[j]*wt[i]; + for(int j=0; j<nCoeffs; j++) + for(int k=0; k<valueDim; k++) + Fy(j+1,k+1) += values[i][k]*tempFx[j]*wt[i]; + } + + try { + F2.invert(); + } catch(GeneralClassicException exc) { + delete temp; + throw(GeneralClassicException( + "PolynomialVector::PolynomialLeastSquaresFit", + "Could not find least squares fit for data")); + } + MMatrix<double> A = F2 * Fy; + delete temp; + temp = new PolynomialVector(pointDim, A.T()); + return temp; +} + +PolynomialVector* PolynomialVector::ConstrainedPolynomialLeastSquaresFit( + const std::vector< std::vector<double> >& points, + const std::vector< std::vector<double> >& values, + int polynomialOrder, + std::vector< PolynomialVector::PolynomialCoefficient > coeffs, + const std::vector<double>& weights) { + //Algorithm: we want g(x) = old_f(x) + new_f(x), + //where old_f has polynomial terms in poly, new_f has all the rest + //Use value - f(point) as input and do LLS forcing old_f(x) polynomial terms to 0 + //Nb: in this constrained case we have to go through the output vector and treat each like a 1D vector + int pointDim = points[0].size(); + int valueDim = values[0].size(); + int nPoints = points.size(); + if(nPoints<1) throw(GeneralClassicException("PolynomialVector::ConstrainedPolynomialLeastSquaresFit(...)", "No data for LLS Fit")); + if(int(values.size())!=nPoints) throw(GeneralClassicException("PolynomialVector::ConstrainedPolynomialLeastSquaresFit(...)", "Misaligned array in LLS Fit")); + int nCoeffsNew = NumberOfPolynomialCoefficients(pointDim, polynomialOrder); + for (size_t i = 0; i < coeffs.size(); ++i) { + std::vector<int> in_var = coeffs[i].InVariables(); + if (coeffs[i].OutVariable() >= valueDim) + throw(GeneralClassicException( + "PolynomialVector::ConstrainedPolynomialLeastSquaresFit(...)", + "Coefficient output variable out of range")); + if ( int(in_var.size()) >= polynomialOrder) + throw(GeneralClassicException( + "PolynomialVector::ConstrainedPolynomialLeastSquaresFit(...)", + "Coefficient input variables of too high order")); + for (size_t j = 0; j < coeffs[i].InVariables().size(); ++j) { + if (in_var[j] >= pointDim) + throw(GeneralClassicException( + "PolynomialVector::ConstrainedPolynomialLeastSquaresFit(...)", + "Coefficient input variable of too high dimension")); + } + } + + MMatrix<double> A(valueDim, nCoeffsNew, 0.); + + std::vector<double> wt (nPoints, 1); + if(weights.size() > 0) wt = weights; + + PolynomialVector newPolyVector1D(pointDim, MMatrix<double>(1, nCoeffsNew)); //guaranteed to be of right order + std::vector< std::vector<double> > tempFx(nPoints, std::vector<double>(nCoeffsNew) ); + for(int i=0; i<nPoints; i++) + newPolyVector1D.MakePolyVector(&points[i][0], &tempFx[i][0]); + + for(int dim=0; dim<valueDim ; dim++) + { + std::vector<PolynomialVector::PolynomialCoefficient> oldCoeffs1D; + for(unsigned int i=0; i<coeffs.size(); i++) + if(coeffs[i].OutVariable() == dim && int(coeffs[i].InVariables().size()) < polynomialOrder) + oldCoeffs1D.push_back(coeffs[i]); + PolynomialVector oldPolyVector1D(oldCoeffs1D); + int nCoeffsOld = oldCoeffs1D.size(); + int deltaCoeff = nCoeffsNew - nCoeffsOld; + if(deltaCoeff <= 0) break; + + + std::vector<int> needsCalculation; //index of variables that need calculation + for(int i=0; i<nCoeffsNew; i++) + { + bool exists = true; + for(unsigned int j=0; j<oldCoeffs1D.size(); j++) + if(IndexByVector(i, pointDim) == oldCoeffs1D[j].InVariables()) exists = false; + if(exists) needsCalculation.push_back(i); + } + + MVector<double> Fy(deltaCoeff, 0); + MMatrix<double> F2(deltaCoeff, deltaCoeff, 0.); + + for(int i=0; i<nPoints && needsCalculation.size()>0; i++) //optimisation note - this algorithm spends all its time in this loop + { + std::vector<double> oldValue(valueDim); + oldPolyVector1D.F(&points[i][0], &oldValue[0]); + for(int k=0; k<deltaCoeff; k++) + { + Fy(k+1) += (values[i][dim] - oldValue[dim])*tempFx[i][needsCalculation[k]]*wt[i]; //ynew - yold + for(int j=0; j<deltaCoeff; j++) + F2(j+1,k+1) += tempFx[i][needsCalculation[k]]*tempFx[i][needsCalculation[j]]*wt[i]; + } + } + try{ F2.invert(); } + catch(GeneralClassicException exc) { + throw(GeneralClassicException("PolynomialVector::ConstrainedPolynomialLeastSquaresFit", "Could not find least squares fit for data")); } + MVector<double> AVec = F2 * Fy; + for(int i=0; i<deltaCoeff; i++) A(dim+1, needsCalculation[i]+1) = AVec(i+1); + } + + PolynomialVector tempPoly(coeffs); + for(unsigned int i=0; i<coeffs.size(); i++) + for(unsigned int j=0; j<tempPoly._polyKeyByVector.size(); j++) + if(tempPoly._polyKeyByVector[j] == coeffs[i].InVariables()) + { + A(coeffs[i].OutVariable()+1,j+1) = coeffs[i].Coefficient(); + } + + return new PolynomialVector(pointDim, A); +} + + +PolynomialVector* PolynomialVector::Chi2ConstrainedLeastSquaresFit + (const std::vector< std::vector<double> >& xin, const std::vector< std::vector<double> >& xout, + int polynomialOrder, std::vector< PolynomialVector::PolynomialCoefficient > coeffs, + double chi2Start, double discardStep, double* chi2End, double chi2Limit, + std::vector<double> weights, bool firstIsMean) +{ + int nParticles = xin.size(); + int dimensionOut = xout[0].size(); + if(int(weights.size()) != nParticles) weights = std::vector<double>(xin.size(), 1.); + std::vector<double> weights_in = std::vector<double>(weights); + std::vector<double> amps(nParticles); + MVector<double> means(dimensionOut,0); + if(!firstIsMean) + means = Means (xout, weights); + else for(int i=0; i<dimensionOut; i++) means(i+1) = xout[0][i]; + + MMatrix<double> covInv = Covariances(xout, weights, means); + try{ covInv.invert(); } + catch(GeneralClassicException exc) {throw(GeneralClassicException("PolynomialVector::Chi2ConstrainedLeastSquaresFit", "Failed to find least squares fit for data"));} + + double totalWeight = 0.; + double discard = chi2Start; + double chi2 = chi2Limit*2.; + int nCoefficients = PolynomialVector::NumberOfPolynomialCoefficients(xin[0].size(), polynomialOrder)*xout[0].size(); + int nGood = nParticles; + + if(discard <= 0.) for(int i=0; i<nParticles; i++) //if chi2Start ill-defined, just use maximum chi2 in sample + { + MVector<double> xoutVec (dimensionOut, 0); + for(int j=0; j<dimensionOut; j++) + xoutVec(j+1) = xout[i][j] - means(j+1); + amps[i] = (xoutVec.T() * covInv * xoutVec)(1); + if(amps[i] > discard) discard = amps[i]; + totalWeight += weights[i]; + } + discard /= discardStep; //Set discard to the largest amplitude in the data + + PolynomialVector* pvec = new PolynomialVector(std::vector<PolynomialVector::PolynomialCoefficient>() ); + while( nGood >= nCoefficients && chi2 > chi2Limit) + { + chi2 = 0.; + nGood = nParticles; + if(pvec != NULL) delete pvec; + + //optimisation note - algorithm spends all its time doing the CPLS Fit + try { pvec = PolynomialVector::ConstrainedPolynomialLeastSquaresFit(xin, xout, polynomialOrder, coeffs, weights); } + catch(GeneralClassicException exc) { pvec = NULL; chi2 = chi2Limit*2.;} + for(int i=0; i<nParticles && pvec!=NULL; i++) + { + if(fabs(weights[i]) > 1.e-9) + { + std::vector<double> pout(dimensionOut, 0.); + pvec->F(&xin[i][0], &pout[0]); + MVector<double> residualVec(dimensionOut,0); + for(int j=0; j<dimensionOut; j++) residualVec(j+1) = pout[j] - xout[i][j]; + chi2 += weights[i]*weights[i]/totalWeight/totalWeight * (residualVec.T() * covInv * residualVec)(1); //watch weights handling here + } + } + for(int i=0; i<nParticles; i++) + if(amps[i] > discard) {totalWeight -= weights[i]; weights[i] = 0.; nGood--;} + discard /= discardStep; + std::cout << "ConstrainedPolynomialLeastSquaresFit - chi2: " << chi2 << " cut: " << discard << " cut_step: " << discardStep << " weight: " << totalWeight + << " pvec: " << pvec << std::endl; + } + if(chi2 > chi2Limit || pvec == NULL) + { + std::stringstream err; + err << "PolynomialVector::Chi2ConstrainedLeastSquaresFit failed to converge. Best fit had <chi2>=" + << chi2 << " compared to limit " << chi2Limit << std::endl; + throw(GeneralClassicException("PolynomialVector::Chi2ConstrainedLeastSquaresFit(...)", err.str())); + } + if(chi2End != NULL) *chi2End = discard; //save discard at end for future use + return pvec; +} + +void PolynomialVector::PolynomialCoefficient::SpaceTransform(std::vector<int> space_in, std::vector<int> space_out) +{ + std::map<int, int> mapping; //probably optimise this + for(unsigned int i=0; i<space_out.size(); i++) + for(unsigned int j=0; j<space_in.size(); j++) + if(space_out[i] == space_in[j]) mapping[j] = i; //mapping[space_in_index] returns space_out_index + + std::vector<int> in_variables(_inVarByVec.size()); + for(unsigned int con=0; con<in_variables.size(); con++) + { + if(mapping.find(in_variables[con]) != mapping.end()) in_variables[con] = mapping[in_variables[con]]; + else throw(GeneralClassicException("PolynomialVector::PolynomialCoefficient::SpaceTransform", "Input variable not found in space transform")); + } + + if(mapping.find(_outVar) != mapping.end()) _outVar = mapping[_outVar]; + else throw(GeneralClassicException( + "PolynomialVector::PolynomialCoefficient::SpaceTransform", + "Output variable not found in space transform" +)); + + _inVarByVec = in_variables; +} +////////// POLYNOMIALVECTOR END //////////// + +} + + diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialVector.h b/classic/5.0/src/Fields/Interpolation/PolynomialVector.h new file mode 100644 index 0000000000000000000000000000000000000000..ebb6a2f5bf22e74edbc893ac1755ecb64e9903ff --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialVector.h @@ -0,0 +1,290 @@ +// MAUS WARNING: THIS IS LEGACY CODE. +// Copyright 2010-2011 Chris Rogers +// +// This file is a part of G4MICE +// +// G4MICE is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// G4MICE is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with G4MICE in the doc folder. If not, see +// <http://www.gnu.org/licenses/>. + +#ifndef PolynomialVector_hh +#define PolynomialVector_hh 1 + +#include "MMatrix.h" +#include "MVector.h" + +#include <map> + +/// PolynomialVector, an arbitrary order polynomial vector class. + +/// PolynomialVector describes a vector of multivariate polynomials \f$y_i = a_0 + Sum (a_j x^j)\f$ +/// i.e. maps a vector \f$\vec{x}\f$ onto a vector \f$\vec{y}\f$ with \f$\vec{y} = a_0 + sum( a_{j_1j_2...j_n} x_1^{j1} x_2^{j2} ... x_n^{jn} )\f$. +/// Also algorithms to map a single index J into \f$j_1...j_n\f$.\n +/// \n +/// PolynomialVector represents the polynomial coefficients as a matrix of doubles so that\n +/// \f$ \vec{y} = \mathbf{A} \vec{w} \f$\n +/// where \f$\vec{w}\f$ is a vector with elements given by \n +/// \f$ w_i = x_1^{i_1} x_2^{i_2} ... x_n^{i_n} \f$ \n +/// So the index \f$ i \f$ is actually itself a vector. The vectorisation of \f$ i \f$ is handled by IndexByPower and IndexByVector methods. +/// IndexByPower gives an index like:\n +/// \f$ w_i = x_1^{i_1} x_2^{i_2} ... x_n^{i_n} \f$ \n +/// While IndexByVector gives an index like:\n +/// \f$ w_i = x_{i_1}x_{i_2} \ldots x_{i_n} \f$ \n +/// \n +/// Polynomial Vector includes several functions to do least squares fits. +/// Here we find a polynomial of arbitrary dimension and order from a set of points by the method of least squares. +/// Note that the problem must be overspecified, so the number of points must be >= number of polynomial coefficients. +/// A few interesting variations on this theme... have a look!\n +///\n + +class PolynomialVector +{ +public: + + class PolynomialCoefficient; + /// Default constructor (leaves everything empty, call SetCoefficients to set up properly) + PolynomialVector(); + /// Copy constructor + PolynomialVector (const PolynomialVector& pv); + /// Construct polynomial vector passing polynomial coefficients as a matrix. + PolynomialVector (int pointDimension, MMatrix<double> polynomialCoefficients); + /// Construct polynomial vector passing polynomial coefficients as a list of PolynomialCoefficient objects. + + /// Any coefficients missing from the vector are set to 0. Maximum order of the polynomial is given by the + /// maximum order of the coefficients in the vector. + PolynomialVector (std::vector<PolynomialCoefficient> coefficients); + /// Destructor - no memory allocated so doesn't do anything + ~PolynomialVector() {;} + + /// Reinitialise the polynomial vector with new point (x) dimension and coefficients + void SetCoefficients(int pointDim, MMatrix<double> coeff); + /// Reinitialise the polynomial vector with new point (x) dimension and coefficients. + + /// Any coefficients missing from the vector are set to 0. Maximum order of the polynomial is given by the + /// maximum order of the coefficients in the vector. + void SetCoefficients(std::vector<PolynomialCoefficient> coeff); + /// Set the coefficients. + + /// This method can't be used to change point dimension or + /// value dimension - any range mismatch will be ignored. + void SetCoefficients(MMatrix<double> coeff); + /// Return the coefficients as a matrix of doubles. + MMatrix<double> GetCoefficientsAsMatrix() const {return _polyCoeffs;} + /// Return the coefficients as a vector of PolynomialCoefficients, including 0 values. + std::vector<PolynomialCoefficient> GetCoefficientsAsVector() const {return _polyVector;} + + /// Fill value with \f$ y_i \f$ at some set of \f$ x_i \f$ (point). + + /// point should be array of length PointDimension(). + /// value should be array of length ValueDimension(). + void F (const double* point, double* value) const; + + /** Calculate derivative of the field at a given point + */ + void FDeriv(const double* point, const int* derivative_by_power, double* value) const; + + + /// Fill value with \f$ y_i \f$ at some set of \f$ x_i \f$ (point). + + /// point should be vector of length PointDimension(). + /// value should be vector of length ValueDimension(). + /// Note that there is no bound checking here. + void F (const MVector<double>& point, MVector<double>& value) const; + /// Length of the input point (x) vector. + unsigned int PointDimension() const {return _pointDim;} + /// Length of the output value (y) vector. + unsigned int ValueDimension() const {return _polyCoeffs.num_row();} + /// Index of highest power - 0 is const, 1 is linear, 2 is quadratic etc... + unsigned int PolynomialOrder() const {if(_polyKeyByVector.size()>0) return _polyKeyByVector.back().size()+1; else return 0;} + /// Polymorphic copy constructor. + + /// This is a special copy constructor for inheritance structures + PolynomialVector* Clone() const {return new PolynomialVector(*this);} + /// Return a copy, centred on point. + PolynomialVector* Recentred(double* point) const; + /// Return operator $\f R = P(Q) \f$ - NOT IMPLEMENTED + PolynomialVector Chain(PolynomialVector Q) const; + /// Return inverse of the polynomial truncated at order n (in general an infinite series that does not converge, so beware!) - NOT IMPLEMENTED + PolynomialVector Inverse(int max_order) const; + + /// Make a vector like \f$(c, x, x^2, x^3...)\f$. + /// polyVector should be of size NumberOfPolynomialCoefficients(). + /// could be static but faster as member function (use lookup table for _polyKey). + MVector<double>& MakePolyVector(const MVector<double>& point, MVector<double>& polyVector) const; + /// Make a vector like \f$(c, x, x^2, x^3...)\f$. + /// PolyVector should be of size NumberOfPolynomialCoefficients(). + /// Could be static but faster as member function (use lookup table for _polyKey). + double* MakePolyVector(const double* point, double* polyVector) const; + + /** Make a vector like \f$d^\vec{p}(c, x, x^2, x^3...)/d\vec{x}^\vec{p}\f$. + * - positions: array of size PointDimension() containing the position at + * which the vector should be calculated, \f$\vec{x}\f$ + * - deriv_indices: array in the "index by power" format defining the + * derivative index \vec{p}. + * - deriv_vec: array that will hold the return value; should be + * initialised to size at least NumberOfPolynomialCoefficients() + */ + double* MakeDerivVector(const double* positions, const int* deriv_indices, double* deriv_vec) const; + + /// Transforms from a 1d index of polynomial coefficients to an nd index. + /// This is slow - you should use it to build a lookup table. + /// For polynomial term \f$x_1^i x_2^j ... x_d^n\f$ index like [i][j] ... [n] + /// e.g. \f$x_1^4 x_2^3 = x_1^4 x_2^3 x_3^0 x_4^0\f$ = {4,3,0,0}. + static std::vector<int> IndexByPower (int index, int nInputVariables); + /// Transforms from a 1d index of polynomial coefficients to an nd index. + /// This is slow - you should use it to build a lookup table. + /// For polynomial term \f$x_i x_j...x_n\f$ index like [i][j] ... [n] + /// e.g. \f$x_1^4 x_2^3\f$ = {1,1,1,1,2,2,2} + static std::vector<int> IndexByVector(int index, int nInputVariables); + /// Returns the number of coefficients required for an arbitrary dimension, order polynomial + /// e.g. \f$a_0 + a_1 x + a_2 y + a_3 z + a_4 x^2 + a_5 xy + a_6 y^2 + a_7 xz + a_8 yz + a_9 z^2\f$ + /// => NumberOfPolynomialCoefficients(3,2) = 9 (indexing starts from 0). + static unsigned int NumberOfPolynomialCoefficients(int pointDimension, int order); + unsigned int NumberOfPolynomialCoefficients() {return NumberOfPolynomialCoefficients(_pointDim, PolynomialOrder());} + + /// Write out the PolynomialVector (effectively just polyCoeffs). + friend std::ostream& operator<<(std::ostream&, const PolynomialVector&); + + /// Control the formatting of the polynomial vector. If PrintHeaders is true, then every time I write + /// a PolynomialVector it will write the header also (default). + static void PrintHeaders(bool willPrintHeaders) {_printHeaders = willPrintHeaders;} + /// Write out the header (PolynomialByPower index for each element). + void PrintHeader(std::ostream& out, char int_separator='.', char str_separator=' ', int length=14, bool pad_at_start=true) const; + + /// Return a PolynomialVector found from a set of points by method of least squares. + + /// Finds the polynomial pol that minimises sum_i ( w_i (\vec{y_i} - Pol(\vec{x_i}))^2 + /// y_i = values[i], x_i = points[i], w_i = weightFunction(x_i) OR weights[i], polynomialOrder is the order of pol + /// Nb if w_i not defined, I assume that w_i = 1 i.e. unweighted + /// points and values should be the same length; points[i] should all be the same length; values[i] should all be the same length + /// weight function should be a vector map that has PointDimension of points and returns a weight of dimension 1 + static PolynomialVector* PolynomialLeastSquaresFit(const std::vector< std::vector<double> >& points, const std::vector< std::vector<double> >& values, + int polynomialOrder, const std::vector<double>& weights=std::vector<double>()); + /// Find a polynomial least squares fit given that I know certain coefficients already, (stored in coeffs) + /// (e.g. I know the polynomial to 2nd order, I want to extend to 3rd order using a least squares) + static PolynomialVector* ConstrainedPolynomialLeastSquaresFit(const std::vector< std::vector<double> >& points, const std::vector< std::vector<double> >& values, + int polynomialOrder, std::vector< PolynomialCoefficient > coeffs, + const std::vector<double>& weights=std::vector<double>()); + /// Find a polynomial least squares fit given that I know certain coefficients already + /// After finding the fit, calculate "chi2" for each particle. Compare x_out with polynomial fit of x_in, x_pol, + /// using delta_chi2 = \vec(x_out - x_pol).T() V^{-1} \vec(x_out - x_pol). If sum(delta_chi2)/totalWeight > chi2Limit, + /// calculate chi2_out = \vec(x_out).T() V^{-1} \vec(x_out) and discard particles with chi2_out > chi2Start*(discardStep^n) + /// repeat until sum(delta_chi2)/totalWeight <= chi2Limit, each time incrementing n by 1 + /// if chi2Start < 0, start with maximum chi2; if chi2End != NULL, put the final value in chi2End + static PolynomialVector* Chi2ConstrainedLeastSquaresFit + (const std::vector< std::vector<double> >& xin, const std::vector< std::vector<double> >& xout, + int polynomialOrder, std::vector< PolynomialVector::PolynomialCoefficient > coeffs, + double chi2Start, double discardStep, double* chi2End, double chi2Limit, + std::vector<double> weights, bool firstIsMean=false); +/* + static PolynomialVector* PolynomialSolve( + int polynomialOrder, + const std::vector< std::vector<double> >& positions, + const std::vector< std::vector<double> >& values, + const std::vector< std::vector<double> > &deriv_positions, + const std::vector< std::vector<double> >& deriv_values, + const std::vector< std::vector<int> >& deriv_indices); +*/ + /// Now some utility functions + /// Should probably go in a "utility function" namespace, that does not currently exist. + + /// Return the mean of some set of values, using some set of weights + static MVector<double> Means (std::vector<std::vector<double> > values, std::vector<double> weights); + /// Return the covariance matrix of some set of values, using some set of weights and some mean + /// If length of weights is not equal to length of values use weights = 1. + /// If length of mean is not equal to length of first value recalculate mean + static MMatrix<double> Covariances(std::vector<std::vector<double> > values, std::vector<double> weights, MVector<double> mean); + /// Return chi2 of residuals of some set of output variables compared with the polynomial vector of some set of input variables + /// Here Chi2 Avg is defined as Sum( v^T M^{-1} v ) where v is the vector of residuals and M is covariance matrix of v + double GetAvgChi2OfDifference(std::vector< std::vector<double> > in, std::vector< std::vector<double> > out); + /// Should probably go in a "utility function" namespace, that does not currently exist. + /// Print a sequence: with some string that separates elements of index and some character that pads the total vector + /// length is the total length of the output, set to < 1 to ignore (means no padding) + /// pad_at_start determines whether to pad at start (i.e. right align) or end (i.e. left align) + template < class Container > + static void PrintContainer(std::ostream& out, const Container& container, char T_separator, char str_separator, int length, bool pad_at_start); + /// Return true if a and b are identical; a must have an iterator object dereferenced by *it + /// must have increment prefix operator defined and with != operator defined by the object pointed to by it + template <class TEMP_CLASS, class TEMP_ITER> + static bool IterableEquality(const TEMP_CLASS& a, const TEMP_CLASS& b, TEMP_ITER a_begin, TEMP_ITER a_end, TEMP_ITER b_begin, TEMP_ITER b_end); + + /// PolynomialCoefficient sub-class - useful, could be extended if needed (e.g. become an iterator etc) + class PolynomialCoefficient + { + public: + /// inVariablesByVector is x power indexed like e.g. \f$x_1^4 x_2^3 =\f$ {1,1,1,1,2,2,2} + PolynomialCoefficient(std::vector<int> inVariablesByVector, int outVariable, double coefficient) + : _inVarByVec(inVariablesByVector), _outVar(outVariable), _coefficient(coefficient) {} + PolynomialCoefficient(const PolynomialCoefficient& pc) + : _inVarByVec(pc._inVarByVec), _outVar(pc._outVar), _coefficient(pc._coefficient) {} + /// Accessors + /// Pass a parameter to *set* + std::vector<int> InVariables() const {return _inVarByVec;} + int OutVariable() const {return _outVar;} + double Coefficient() const {return _coefficient;} + + std::vector<int> InVariables(std::vector<int> inVar ) {_inVarByVec = inVar; return _inVarByVec;} + int OutVariable(int outVar) {_outVar = outVar; return _outVar;} + double Coefficient(double coeff ) {_coefficient = coeff; return _coefficient;} + + /// transform coefficient from subspace space_in to subspace space_out, both subspaces of some larger space + /// if any of coeff variables is not in space_out OR not in space_in, leave this coeff untouched and Squeal + /// so for coeff({1,2},0,1.1), coeff.space_transform({0,2,3,5}, {4,7,1,2,3,0}) would return coeff({3,4},5,1.1) + void SpaceTransform(std::vector<int> space_in, std::vector<int> space_out); + /// if var is in inVariables return true + bool HasInVariable(int var) const {for(unsigned int i=0; i<_inVarByVec.size(); i++) if(_inVarByVec[i] == var) return true; return false;} + + + + private: + std::vector<int> _inVarByVec; + int _outVar; + double _coefficient; + }; + + +private: + int _pointDim; + std::vector< std::vector<int> > _polyKeyByPower; //std::vector<int>[i_1] = j_1 + std::vector< std::vector<int> > _polyKeyByVector; //std::vector<int>[j_1] = i_1 + MMatrix<double> _polyCoeffs; + std::vector<PolynomialCoefficient> _polyVector; + static bool _printHeaders; +}; + +std::ostream& operator<<(std::ostream&, const PolynomialVector&); +std::ostream& operator << (std::ostream&, const PolynomialVector::PolynomialCoefficient&); + +//template function + +template <class TEMP_CLASS, class TEMP_ITER> +bool PolynomialVector::IterableEquality(const TEMP_CLASS& a, const TEMP_CLASS& b, TEMP_ITER a_begin, TEMP_ITER a_end, TEMP_ITER b_begin, TEMP_ITER b_end) +{ + TEMP_ITER a_it=a_begin; + TEMP_ITER b_it=b_begin; + while(a_it!=a_end && b_it!=b_end) { + if( *a_it != *b_it ) return false; + ++a_it; + ++b_it; + } + if ( a_it != a_end || b_it != b_end ) return false; + return true; +} + + + +#endif + + + + diff --git a/classic/5.0/src/Fields/Interpolation/PolynomialVectorFactory.h b/classic/5.0/src/Fields/Interpolation/PolynomialVectorFactory.h new file mode 100644 index 0000000000000000000000000000000000000000..e526234d26994862aebe93e5c42732cbad105a4e --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/PolynomialVectorFactory.h @@ -0,0 +1,37 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +namespace interpolation { + +class PolynomialVectorFactory { + public: + PolynomialVectorFactory() {} + virtual ~PolynomialVectorFactory() {} +}; + +} + diff --git a/classic/5.0/src/Fields/Interpolation/SolveFactory.cpp b/classic/5.0/src/Fields/Interpolation/SolveFactory.cpp new file mode 100644 index 0000000000000000000000000000000000000000..85249619ca9fc0f60d6edda02d99c28da226ef9e --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/SolveFactory.cpp @@ -0,0 +1,165 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include <gsl/gsl_sf_pow_int.h> + +#include "Utilities/GeneralClassicException.h" + +#include "Fields/Interpolation/PPSolveFactory.h" +#include "Fields/Interpolation/SolveFactory.h" + +namespace interpolation { + +SolveFactory::SolveFactory(int polynomial_order, + int smoothing_order, + int point_dim, + int value_dim, + std::vector< std::vector<double> > positions, + std::vector< std::vector<double> > deriv_positions, + std::vector< std::vector<int> >& deriv_indices) + : polynomial_order_(polynomial_order), smoothing_order_(smoothing_order) { + n_poly_coeffs_ = SquarePolynomialVector::NumberOfPolynomialCoefficients(point_dim, smoothing_order); + square_points_ = PPSolveFactory::getNearbyPointsSquares(point_dim, -1, smoothing_order); + square_deriv_nearby_points_ = PPSolveFactory::getNearbyPointsSquares(point_dim, -1, smoothing_order); + BuildHInvMatrix(positions, deriv_positions, deriv_indices); + MMatrix<double> A_temp(value_dim, n_poly_coeffs_, 0.); + square_temp_ = SquarePolynomialVector(point_dim, A_temp); +} + +void SolveFactory::BuildHInvMatrix( + std::vector< std::vector<double> > positions, + std::vector< std::vector<double> > deriv_positions, + std::vector< std::vector<int> >& deriv_indices) { + int nCoeffs = positions.size(); + h_inv_ = MMatrix<double>(n_poly_coeffs_, n_poly_coeffs_, 0.); + for (int i = 0; i < nCoeffs; ++i) { + std::vector<double> poly_vec = MakeSquareVector(positions[i]); + for (int j = 0; j < n_poly_coeffs_; ++j) { + h_inv_(i+1, j+1) = poly_vec[j]; + } + } + for (int i = 0; i < deriv_positions.size(); ++i) { + std::vector<double> deriv_vec = MakeSquareDerivVector(deriv_positions[i], + deriv_indices[i]); + for (int j = 0; j < n_poly_coeffs_; ++j) { + h_inv_(i+1+nCoeffs, j+1) = deriv_vec[j]; + } + } + h_inv_.invert(); +} + +std::vector<double> SolveFactory::MakeSquareVector(std::vector<double> x) { + std::vector<double> square_vector(square_points_.size(), 1.); + for (size_t i = 0; i < square_points_.size(); ++i) { + std::vector<int>& point = square_points_[i]; + for (size_t j = 0; j < point.size(); ++j) + square_vector[i] *= gsl_sf_pow_int(x[j], point[j]); + } + return square_vector; +} + +std::vector<double> SolveFactory::MakeSquareDerivVector(std::vector<double> positions, std::vector<int> deriv_indices) { + std::vector<double> deriv_vec(square_deriv_nearby_points_.size(), 1.); + int square_deriv_nearby_points_size = square_deriv_nearby_points_.size(); + int dim = square_deriv_nearby_points_[0].size(); + for (int i = 0; i < square_deriv_nearby_points_size; ++i) { + for (int j = 0; j < dim; ++j) { + int power = square_deriv_nearby_points_[i][j] - deriv_indices[j]; // p_j - q_j + if (power < 0) { + deriv_vec[i] = 0.; + break; + } else { + // x^(p_j-q_j) + deriv_vec[i] *= gsl_sf_pow_int(positions[j], power); + } + // p_j*(p_j-1)*(p_j-2)*...*(p_j-q_j) + for (int k = square_deriv_nearby_points_[i][j]; k > power; --k) { + deriv_vec[i] *= k; + } + } + } + return deriv_vec; +} + +SquarePolynomialVector* SolveFactory::PolynomialSolve( + const std::vector< std::vector<double> >& values, + const std::vector< std::vector<double> >& deriv_values) { + // Algorithm: + // define G_i = vector of F_i and d^pF/dF^p with values taken from coeffs + // and derivs respectively + // define H_ij = vector of f_i and d^pf/df^p) + // Then use G = H A => A = H^-1 G + // where A is vector of polynomial coefficients + // First set up G_i from coeffs and derivs; then calculate H_ij; + // then do the matrix inversion + // It is an error to include the same polynomial index more than once in + // coeffs or derivs or both; any polynomial indices that are missing will be + // assigned 0 value; polynomial order will be taken as the maximum + // polynomial order from coeffs and derivs. + // PointDimension and ValueDimension will be taken from coeffs and derivs; + // it is an error if these do not all have the same dimensions. + + // OPTIMISATION - if we are doing this many times and only changing values, + // can reuse H + int nCoeffs = values.size(); + int nDerivs = deriv_values.size(); + if (values.size()+deriv_values.size() != size_t(n_poly_coeffs_)) { + throw GeneralClassicException( + "SolveFactory::PolynomialSolve", + "Values and derivatives over or under constrained" + ); + } + int pointDim = square_temp_.PointDimension(); + int valueDim = 0; + if (values.size() != 0) { + valueDim = values[0].size(); + } else if (deriv_values.size() != 0) { + valueDim = deriv_values[0].size(); + } + + MMatrix<double> A(valueDim, n_poly_coeffs_, 0.); + for (size_t y_index = 0; y_index < values[0].size(); ++y_index) { + MVector<double> G(n_poly_coeffs_, 0.); + // First fill the values + for (int i = 0; i < nCoeffs && i < n_poly_coeffs_; ++i) { + G(i+1) = values[i][y_index]; + } + // Now fill the derivatives + for (int i = 0; i < nDerivs; ++i) { + G(i+nCoeffs+1) = deriv_values[i][y_index]; + } + MVector<double> A_vec = h_inv_*G; + for (int j = 0; j < n_poly_coeffs_; ++j) { + A(y_index+1, j+1) = A_vec(j+1); + } + } + SquarePolynomialVector* temp = new SquarePolynomialVector(square_temp_); + temp->SetCoefficients(A); + return temp; +} +} + diff --git a/classic/5.0/src/Fields/Interpolation/SolveFactory.h b/classic/5.0/src/Fields/Interpolation/SolveFactory.h new file mode 100644 index 0000000000000000000000000000000000000000..2b1a09f240295e7a8f0bc5c7994058472aaee1d3 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/SolveFactory.h @@ -0,0 +1,133 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "Fields/Interpolation/MMatrix.h" +#include "Fields/Interpolation/SquarePolynomialVector.h" + +#ifndef SolveFactory_hh +#define SolveFactory_hh + +namespace interpolation { + +/** \class[SolveFactory] is a factory class for solving a set of linear + * equations to generate a polynomial based on nearby points. + * + * See also \class[PPSolveFactory] for more details + * + */ + +class SolveFactory { + public: + /** Construct a new SolveFactory. + * + * \param[polynomial_order] Order of the polynomial part of the fit + * \param[smoothing_order] Order of the smoothing part of the fit; should + * always be >= polynomial_order + * \param[point_dim] Dimension of the points (i.e. space of x), for + * y = f(x) + * \param[value_dim] Dimension of the values (i.e. space of y), for + * y = f(x) + * \param[positions] Position of the values; should be a vector of length + * polynomial_order^point_dim. Each element should be a + * vector of length point_dim + * \param[deriv_positions] Position of the derivatives; should be a vector + * of length + * smoothing_oder^point_dim - polynomial_order^point_dim. + * Each element should be a vector of length point_dim. + * \param[deriv_indices] Index the derivatives. Should be a vector with + * same length as deriv_positions. Each element should + * be a vector of length point_dim. Each element tells us + * what derivative will be defined, indexing by the order + * of the derivative. + * + * The constructor takes the set of positions; then performs a matrix + * inversion. Each call to PolynomialSolve uses the same matrix inversion + * in solution of the linear equations (hence when we do this lots of times + * as in PolynomialPatch, we get a lot faster). The matrix inversion can + * fail for badly formed sets of positions/deriv_positions; caveat emptor! + */ + SolveFactory(int polynomial_order, + int smoothing_order, + int point_dim, + int value_dim, + std::vector< std::vector<double> > positions, + std::vector< std::vector<double> > deriv_positions, + std::vector< std::vector<int> >& deriv_indices); + + /** Destructor; does nothing */ + ~SolveFactory() {} + + /** Solve to get a SquarePolynomialVector + * + * \param[values] Set of y values for the solve, giving function values at + * the positions defined in the constructor + * \param[deriv_values] Set of y derivatives for the solve, giving q^th + * derivatives, q defined by deriv_indices, at positions + * defined by deriv_positions (in constructor) + */ + SquarePolynomialVector* PolynomialSolve( + const std::vector< std::vector<double> >& values, + const std::vector< std::vector<double> >& deriv_values); + + /** Convert a position vector to a set of polynomial products + * + * \param[position] a position vector \f$ x_{b} \f$ + * + * Return value is a set of polynomial products \f$ x_{b}^{p} \f$ + */ + std::vector<double> MakeSquareVector(std::vector<double> position); + + /** Convert a position vector to a derivative of a set of polynomial products + * + * \param[position] a position vector \f$ x_{b} \f$ + * \param[derivIndices] indexes the derivative, \f$\vec{q} \f$ + * + * Return value is a set of derivatives of polynomial products + * \f$ p_j!/(p_j-q_j)! x_{b_j}^{p_j-q_j} \f$ + */ + std::vector<double> MakeSquareDerivVector(std::vector<double> position, std::vector<int> derivIndices); + + private: + void BuildHInvMatrix(std::vector< std::vector<double> > positions, + std::vector< std::vector<double> > deriv_positions, + std::vector< std::vector<int> >& deriv_indices); + + int polynomial_order_; + int smoothing_order_; + int n_poly_coeffs_; + std::vector< std::vector<int> > square_points_; + std::vector< std::vector<int> > square_deriv_nearby_points_; + + SquarePolynomialVector square_temp_; + MMatrix<double> h_inv_; + +}; +} + +#endif // SolveFactory_hh + + diff --git a/classic/5.0/src/Fields/Interpolation/SquarePolynomialVector.cpp b/classic/5.0/src/Fields/Interpolation/SquarePolynomialVector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f72862d9727c7b98374a17b870515a23dbca7410 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/SquarePolynomialVector.cpp @@ -0,0 +1,250 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include <iomanip> +#include <sstream> +#include <math.h> + +#include "gsl/gsl_sf_gamma.h" + +#include "Utilities/GeneralClassicException.h" + +#include "Fields/Interpolation/MMatrix.h" +#include "Fields/Interpolation/MVector.h" +#include "Fields/Interpolation/SquarePolynomialVector.h" + +namespace interpolation { + +bool SquarePolynomialVector::_printHeaders=true; +std::vector< std::vector< std::vector<int> > > SquarePolynomialVector::_polyKeyByPower; +std::vector< std::vector< std::vector<int> > > SquarePolynomialVector::_polyKeyByVector; + +SquarePolynomialVector::SquarePolynomialVector() + : _pointDim(0), _polyCoeffs() { +} + +SquarePolynomialVector::SquarePolynomialVector (const SquarePolynomialVector& pv) + : _pointDim(pv._pointDim), + _polyCoeffs(pv._polyCoeffs.num_row(), pv._polyCoeffs.num_col(), 0.) { + SetCoefficients(pv._polyCoeffs); +} + + +SquarePolynomialVector::SquarePolynomialVector(int numberOfInputVariables, MMatrix<double> polynomialCoefficients) + : _pointDim(numberOfInputVariables), _polyCoeffs(polynomialCoefficients) { + SetCoefficients(numberOfInputVariables, polynomialCoefficients); +} + +SquarePolynomialVector::SquarePolynomialVector(std::vector<PolynomialCoefficient> coefficients) + : _pointDim(0), _polyCoeffs() { + SetCoefficients(coefficients); +} + +void SquarePolynomialVector::SetCoefficients(int pointDim, MMatrix<double> coeff) { + _pointDim = pointDim; + _polyCoeffs = coeff; + + if (int(_polyKeyByVector.size()) < pointDim || _polyKeyByVector[pointDim-1].size() < coeff.num_col()) + IndexByVector(coeff.num_col(), pointDim); // sets _polyKeyByVector and _polyKeyByPower +} + +void SquarePolynomialVector::SetCoefficients(std::vector<PolynomialCoefficient> coeff) { + _pointDim = 0; + int maxPolyOrder = 0; + int valueDim = 0; + for(unsigned int i=0; i<coeff.size(); i++) { + int polyOrder = coeff[i].InVariables().size(); + for(unsigned int j=0; j<coeff[i].InVariables().size(); j++) + if(coeff[i].InVariables()[j] > _pointDim) _pointDim = coeff[i].InVariables()[j]; + if(coeff[i].OutVariable() > valueDim) valueDim = coeff[i].OutVariable(); + if(polyOrder > maxPolyOrder) maxPolyOrder = polyOrder; + } + _pointDim++; //PointDim indexes from 0 + _polyCoeffs = MMatrix<double>(valueDim+1, NumberOfPolynomialCoefficients(_pointDim, maxPolyOrder+1), 0.); + if (int(_polyKeyByVector.size()) < _pointDim || _polyKeyByVector[_pointDim-1].size() < _polyCoeffs.num_col()) + IndexByVector(_polyCoeffs.num_col(), _pointDim); // sets _polyKeyByVector and _polyKeyByPower + + for(size_t i=0; i<_polyCoeffs.num_col(); i++) { + for(unsigned int j=0; j<coeff.size(); j++) + if(coeff[j].InVariables() == _polyKeyByVector[_pointDim-1].back()) { + int dim = coeff[j].OutVariable(); + _polyCoeffs(dim+1,i+1) = coeff[j].Coefficient(); + } + } +} + +void SquarePolynomialVector::SetCoefficients(MMatrix<double> coeff) +{ + for(size_t i=0; i<coeff.num_row() && i<_polyCoeffs.num_row(); i++) + for(size_t j=0; j<coeff.num_col() && j<_polyCoeffs.num_col(); j++) + _polyCoeffs(i+1,j+1) = coeff(i+1,j+1); +} + +void SquarePolynomialVector::F(const double* point, double* value) const +{ + MVector<double> pointV(PointDimension(), 1), valueV(ValueDimension(), 1); + for(unsigned int i=0; i<PointDimension(); i++) pointV(i+1) = point[i]; + F(pointV, valueV); + for(unsigned int i=0; i<ValueDimension(); i++) value[i] = valueV(i+1); +} + +void SquarePolynomialVector::F(const MVector<double>& point, MVector<double>& value) const +{ + MVector<double> polyVector(_polyCoeffs.num_col(), 1); + MakePolyVector(point, polyVector); + MVector<double> answer = _polyCoeffs*polyVector; + for(unsigned int i=0; i<ValueDimension(); i++) value(i+1) = answer(i+1); +} + + +MVector<double>& SquarePolynomialVector::MakePolyVector(const MVector<double>& point, MVector<double>& polyVector) const +{ + for(unsigned int i=0; i<_polyCoeffs.num_col(); i++) { + polyVector(i+1) = 1.; + for(unsigned int j=0; j<_polyKeyByVector[_pointDim-1][i].size(); j++) + polyVector(i+1) *= point( _polyKeyByVector[_pointDim-1][i][j]+1 ); + } + return polyVector; +} + +double* SquarePolynomialVector::MakePolyVector(const double* point, double* polyVector) const +{ + for(unsigned int i=0; i<_polyCoeffs.num_col(); i++) + { + polyVector[i] = 1.; + for(unsigned int j=0; j<_polyKeyByVector[i].size(); j++) + polyVector[i] *= point[_polyKeyByVector[_pointDim-1][i][j] ]; + } + return polyVector; +} + + +void SquarePolynomialVector::IndexByPowerRecursive(std::vector<int> check, size_t check_index, size_t poly_power, std::vector<std::vector<int> >& nearby_points) { + check[check_index] = poly_power; + nearby_points.push_back(check); + if (check_index+1 == check.size()) + return; + for (int i = 1; i < int(poly_power); ++i) + IndexByPowerRecursive(check, check_index+1, i, nearby_points); + for (int i = 0; i <= int(poly_power); ++i) { + check[check_index] = i; + IndexByPowerRecursive(check, check_index+1, poly_power, nearby_points); + } +} + +std::vector<int> SquarePolynomialVector::IndexByPower(int index, int point_dim) { + if (point_dim < 1) + throw(GeneralClassicException( + "PPSolveFactory::GetNearbyPointsSquares", + "Point dimension must be > 0" + )); + while (int(_polyKeyByPower.size()) < point_dim) + _polyKeyByPower.push_back(std::vector< std::vector<int> >()); + if (index < int(_polyKeyByPower[point_dim-1].size())) + return _polyKeyByPower[point_dim-1][index]; + // poly_order 0 means constant term + std::vector<std::vector<int> > nearby_points(1, std::vector<int>(point_dim, 0)); + // poly_order > 0 means corresponding poly terms (1 is linear, 2 is quadratic...) + int this_poly_order = 0; + while (int(nearby_points.size()) < index+1) { + IndexByPowerRecursive(nearby_points[0], 0, this_poly_order+1, nearby_points); + this_poly_order += 1; + } + _polyKeyByPower[point_dim-1] = nearby_points; + return _polyKeyByPower[point_dim-1][index]; +} + +//Turn int index into an index for element of polynomial +// e.g. x_1^4 x_2^3 = {1,1,1,1,2,2,2} +std::vector<int> SquarePolynomialVector::IndexByVector(int index, int point_dim) { + while (int(_polyKeyByVector.size()) < point_dim) + _polyKeyByVector.push_back(std::vector< std::vector<int> >()); + // if _polyKeyByVector is big enough, just return the value + if (index < int(_polyKeyByVector[point_dim-1].size())) + return _polyKeyByVector[point_dim-1][index]; + // make sure _polyKeyByPower is big enough + std::vector<int> index_by_power = IndexByPower(index, point_dim); + // update _polyKeyByVector with values from _polyKeyByPower + for (size_t i = _polyKeyByVector[point_dim-1].size(); i < _polyKeyByPower[point_dim-1].size(); ++i) { + _polyKeyByVector[point_dim-1].push_back(std::vector<int>()); + for (size_t j = 0; j < _polyKeyByPower[point_dim-1][i].size(); ++j) + for (int k = 0; k < _polyKeyByPower[point_dim-1][i][j]; ++k) + _polyKeyByVector[point_dim-1][i].push_back(j); + } + return _polyKeyByVector[point_dim-1][index]; +} + +unsigned int SquarePolynomialVector::NumberOfPolynomialCoefficients(int pointDimension, int order) { + int number = 1; + for (int i = 0; i < pointDimension; ++i) + number *= order+1; + return number; +} + +std::ostream& operator<<(std::ostream& out, const SquarePolynomialVector& spv) +{ + if(SquarePolynomialVector::_printHeaders) spv.PrintHeader(out, '.', ' ', 14, true); + out << '\n' << spv.GetCoefficientsAsMatrix(); + return out; +} + +void SquarePolynomialVector::PrintHeader(std::ostream& out, char int_separator, char str_separator, int length, bool pad_at_start) const +{ + if(_polyKeyByPower[_pointDim-1].size() > 0) PrintContainer< std::vector<int> >(out, _polyKeyByPower[_pointDim-1][0], int_separator, str_separator, length-1, pad_at_start); + for(unsigned int i=1; i<_polyCoeffs.num_col(); ++i) + PrintContainer<std::vector<int> >(out, _polyKeyByPower[_pointDim-1][i], int_separator, str_separator, length, pad_at_start); +} + + +template <class Container > +void SquarePolynomialVector::PrintContainer(std::ostream& out, const Container& container, char T_separator, char str_separator, int length, bool pad_at_start) +{ +// class Container::iterator it; + std::stringstream strstr1(""); + std::stringstream strstr2(""); + class Container::const_iterator it1 = container.begin(); + class Container::const_iterator it2 = it1; + while(it1!=container.end()) + { + it2++; + if(it1 != container.end() && it2 != container.end()) + strstr1 << (*it1) << T_separator; + else strstr1 << (*it1); + it1++; + } + + if(int(strstr1.str().size()) > length && length > 0) strstr2 << strstr1.str().substr(1, length); + else strstr2 << strstr1.str(); + + if(!pad_at_start) out << strstr2.str(); + for(int i=0; i<length - int(strstr2.str().size()); i++) out << str_separator; + if(pad_at_start) out << strstr2.str(); +} + +} + + diff --git a/classic/5.0/src/Fields/Interpolation/SquarePolynomialVector.h b/classic/5.0/src/Fields/Interpolation/SquarePolynomialVector.h new file mode 100644 index 0000000000000000000000000000000000000000..8abf6518b2188d95d1ea416f99645a6ce16a770d --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/SquarePolynomialVector.h @@ -0,0 +1,246 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef SquarePolynomialVector_hh +#define SquarePolynomialVector_hh 1 + +#include "Fields/Interpolation/MMatrix.h" +#include "Fields/Interpolation/MVector.h" +#include "Fields/Interpolation/PolynomialCoefficient.h" + +#include <map> + +namespace interpolation { + +/** SquarePolynomialVector, an arbitrary order polynomial vector class. + * + * SquarePolynomialVector describes a vector of multivariate polynomials\n + * \f$y_i = a_0 + Sum (a_j x^j)\f$ + * i.e. maps a vector \f$\vec{x}\f$ onto a vector \f$\vec{y}\f$ with\n + * \f$\vec{y} = a_0 + sum( a_{j_1j_2...j_n} x_1^{j_1} x_2^{j_2} ... x_n^{j_n})\f$. + * Also algorithms to map a single index J into \f$j_1...j_n\f$.\n + * \n + * PolynomialVector represents the polynomial coefficients as a matrix of + * doubles so that\n + * \f$ \vec{y} = \mathbf{A} \vec{w} \f$\n + * where \f$\vec{w}\f$ is a vector with elements given by \n + * \f$ w_i = x_1^{i_1} x_2^{i_2} ... x_n^{i_n} \f$ \n + * So the index \f$ i \f$ is actually itself a vector. The vectorisation of + * \f$ i \f$ is handled by IndexByPower and IndexByVector methods. + * IndexByPower gives an index like:\n + * \f$ w_i = x_1^{i_1} x_2^{i_2} ... x_n^{i_n} \f$ \n + * While IndexByVector gives an index like:\n + * \f$ w_i = x_{i_1}x_{i_2} \ldots x_{i_n} \f$ \n + * \n + * + * Nb: it is a SquarePolynomialVector because coefficients include all + * polynomial coefficients with \f$ i_ j <= n $\f; coefficients sit within an + * n-dimensional square. The distinction should be made with a PolynomialVector + * where coefficients include all polynomial coefficients with + * \f$ \Sigma(i_j) <= n. + */ + +class SquarePolynomialVector { + public: + /** Default constructor + * + * Leaves everything empty, call SetCoefficients to set up properly. + */ + SquarePolynomialVector(); + + /** Copy constructor + */ + SquarePolynomialVector (const SquarePolynomialVector& pv); + + /** Construct polynomial vector passing polynomial coefficients as a matrix. + */ + SquarePolynomialVector (int pointDimension, + MMatrix<double> polynomialCoefficients); + + /** Construct polynomial vector passing polynomial coefficients as a list of + * PolynomialCoefficient objects. + * + * Any coefficients missing from the vector are set to 0. Maximum order of + * the polynomial is given by the maximum order of the coefficients in the + * vector. + */ + SquarePolynomialVector (std::vector<PolynomialCoefficient> coefficients); + + /** Destructor - no memory allocated so doesn't do anything */ + ~SquarePolynomialVector() {;} + + /** Reinitialise the polynomial vector with new point (x) dimension and + * coefficients + */ + void SetCoefficients(int pointDim, MMatrix<double> coeff); + + /** Reinitialise the polynomial vector with new point (x) dimension and + * coefficients. + * + * Any coefficients missing from the vector are set to 0. Maximum order of + * the polynomial is given by the maximum order of the coefficients in the + * vector. + */ + void SetCoefficients(std::vector<PolynomialCoefficient> coeff); + + /** Set the coefficients. + * + * This method can't be used to change point dimension or value dimension - + * any range mismatch will be ignored. + */ + void SetCoefficients(MMatrix<double> coeff); + + /** Return the coefficients as a matrix of doubles. */ + MMatrix<double> GetCoefficientsAsMatrix() const + {return _polyCoeffs;} + + /** Fill value with \f$ y_i \f$ at some set of \f$ x_i \f$ (point). + * + * point should be array of length PointDimension(). + * value should be array of length ValueDimension(). + */ + void F(const double* point, double* value) const; + + /** Fill value with \f$ y_i \f$ at some set of \f$ x_i \f$ (point). + * + * point should be vector of length PointDimension(). + * value should be vector of length ValueDimension(). + * Note that there is no bound checking here. + */ + void F(const MVector<double>& point, MVector<double>& value) const; + + /** Length of the input point (x) vector. */ + unsigned int PointDimension() const {return _pointDim;} + + /** Length of the output value (y) vector. */ + unsigned int ValueDimension() const {return _polyCoeffs.num_row();} + + /** Index of highest power - 0 is const, 1 is linear, 2 is quadratic etc... + */ + unsigned int PolynomialOrder() const; + + /** Polymorphic copy constructor. + * + * This is a special copy constructor for inheritance structures + */ + SquarePolynomialVector* Clone() const + {return new SquarePolynomialVector(*this);} + + /** Make a vector like \f$(c, x, x^2, x^3...)\f$. + * polyVector should be of size NumberOfPolynomialCoefficients(). + * could be static but faster as member function (use lookup table for + * _polyKey). + */ + MVector<double>& MakePolyVector(const MVector<double>& point, + MVector<double>& polyVector) const; + + /** Make a vector like \f$(c, x, x^2, x^3...)\f$. + * + * PolyVector should be of size NumberOfPolynomialCoefficients(). + * Could be static but faster as member function (use lookup table for + * _polyKey). + */ + double* MakePolyVector(const double* point, double* polyVector) const; + + /** Transforms from a 1d index of polynomial coefficients to an nd index. + * + * This is slow - you should use it to build a lookup table. + * For polynomial term \f$x_1^i x_2^j ... x_d^n\f$ index like [i][j] ... [n] + * e.g. \f$x_1^4 x_2^3 = x_1^4 x_2^3 x_3^0 x_4^0\f$ = {4,3,0,0}. + */ + static std::vector<int> IndexByPower (int index, int nInputVariables); + + /** Transforms from a 1d index of polynomial coefficients to an nd index. + * + * This is slow - you should use it to build a lookup table. + * For polynomial term \f$x_i x_j...x_n\f$ index like [i][j] ... [n] + * e.g. \f$x_1^4 x_2^3\f$ = {1,1,1,1,2,2,2} + */ + static std::vector<int> IndexByVector(int index, int nInputVariables); + + /** Returns the number of coefficients required for an arbitrary dimension, + * order polynomial + * e.g. \f$a_0 + a_1 x + a_2 y + a_3 z + a_4 x^2 + a_5 xy + a_6 y^2 + a_7 xz + a_8 yz + a_9 z^2\f$ + * => NumberOfPolynomialCoefficients(3,2) = 9 (indexing starts from 0). + */ + static unsigned int NumberOfPolynomialCoefficients(int pointDimension, int order); + + /** Write out the PolynomialVector (effectively just polyCoeffs). */ + friend std::ostream& operator<<(std::ostream&, const SquarePolynomialVector&); + + /** Control the formatting of the polynomial vector. + * + * If PrintHeaders is true, then every time I write a PolynomialVector it + * will write the header also (default). + */ + static void PrintHeaders(bool willPrintHeaders) {_printHeaders = willPrintHeaders;} + + /** Write out the header (PolynomialByPower index for each element). */ + void PrintHeader(std::ostream& out, + char int_separator='.', + char str_separator=' ', + int length=14, + bool pad_at_start=true) const; + + /** Print a sequence: with some string that separates elements of index and + * some character that pads + * + * The total vector length is the total length of the output, set to < 1 to + * ignore (means no padding) + * pad_at_start determines whether to pad at start (i.e. right align) or + * end (i.e. left align) + */ + template < class Container > + static void PrintContainer(std::ostream& out, + const Container& container, + char T_separator, + char str_separator, + int length, + bool pad_at_start); + +private: + static void IndexByPowerRecursive( + std::vector<int> check, + size_t check_index, + size_t poly_power, + std::vector<std::vector<int> >& nearby_points); + + int _pointDim; + MMatrix<double> _polyCoeffs; + static std::vector< std::vector< std::vector<int> > > _polyKeyByPower; + static std::vector< std::vector< std::vector<int> > > _polyKeyByVector; + static bool _printHeaders; +}; + +std::ostream& operator<<(std::ostream&, const SquarePolynomialVector&); + +} +#endif // SquarePolynomialVector_hh + + + + diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.cpp b/classic/5.0/src/Fields/Interpolation/ThreeDGrid.cpp similarity index 83% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.cpp rename to classic/5.0/src/Fields/Interpolation/ThreeDGrid.cpp index 0fce740caa55387ad966ff31f967334523aeefcf..e3895cbd8803bc3e1e827ca76fbbd65e09b7025d 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.cpp +++ b/classic/5.0/src/Fields/Interpolation/ThreeDGrid.cpp @@ -25,12 +25,14 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/ThreeDGrid.h" +#include "Fields/Interpolation/ThreeDGrid.h" #include "Utilities/LogicalError.h" +namespace interpolation { + ThreeDGrid::ThreeDGrid() - : x_m(2, 0), y_m(2, 0), z_m(2, 0), xSize_m(0), ySize_m(0), zSize_m(0), maps_m(0), - constantSpacing_m(false) { + : x_m(2, 0), y_m(2, 0), z_m(2, 0), xSize_m(0), ySize_m(0), zSize_m(0), + maps_m(0), constantSpacing_m(false) { setConstantSpacing(); x_m[1] = y_m[1] = z_m[1] = 1.; } @@ -70,8 +72,8 @@ ThreeDGrid::ThreeDGrid(double dX, double dY, double dZ, int numberOfXCoords, int numberOfYCoords, int numberOfZCoords) : x_m(numberOfXCoords), y_m(numberOfYCoords), z_m(numberOfZCoords), - xSize_m(numberOfXCoords), ySize_m(numberOfYCoords), zSize_m(numberOfZCoords), - maps_m(), constantSpacing_m(true) { + xSize_m(numberOfXCoords), ySize_m(numberOfYCoords), + zSize_m(numberOfZCoords), maps_m(), constantSpacing_m(true) { for (int i = 0; i < numberOfXCoords; i++) x_m[i] = minX+i*dX; for (int j = 0; j < numberOfYCoords; j++) @@ -96,12 +98,28 @@ Mesh::Iterator ThreeDGrid::end() const { return Mesh::Iterator(end, this); } -void ThreeDGrid::getPosition(const Mesh::Iterator& it, double * position) const { +void ThreeDGrid::getPosition(const Mesh::Iterator& it, + double * position) const { position[0] = x(it.state_m[0]); position[1] = y(it.state_m[1]); position[2] = z(it.state_m[2]); } +Mesh* ThreeDGrid::dual() const { + std::vector<double> new_x(x_m.size()-1); + std::vector<double> new_y(y_m.size()-1); + std::vector<double> new_z(z_m.size()-1); + for (size_t i = 0; i < x_m.size()-1; ++i) { + new_x[i] = (x_m[i]+x_m[i+1])/2.; + } + for (size_t i = 0; i < y_m.size()-1; ++i) { + new_y[i] = (y_m[i]+y_m[i+1])/2.; + } + for (size_t i = 0; i < z_m.size()-1; ++i) { + new_z[i] = (z_m[i]+z_m[i+1])/2.; + } + return new ThreeDGrid(new_x, new_y, new_z); +} Mesh::Iterator& ThreeDGrid::addEquals (Mesh::Iterator& lhs, int difference) const { @@ -147,6 +165,29 @@ Mesh::Iterator& ThreeDGrid::subEquals return lhs; } +void ThreeDGrid::vectorLowerBound(std::vector<double> vec, + double x, + int& index) { + if (x < vec[0]) { + index = -1; + return; + } + if (x >= vec.back()) { + index = vec.size()-1; + return; + } + size_t xLower = 0; + size_t xUpper = vec.size()-1; + while (xUpper - xLower > 1) { + index = (xUpper+xLower)/2; + if (x >= vec[index]) { + xLower = index; + } else { + xUpper = index; + } + } + index = xLower; +} Mesh::Iterator& ThreeDGrid::addEquals (Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const { @@ -190,7 +231,8 @@ bool ThreeDGrid::isGreater (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const { if (lhs.state_m[0] > rhs.state_m[0]) return true; - else if (lhs.state_m[0] == rhs.state_m[0] && lhs.state_m[1] > rhs.state_m[1]) + else if (lhs.state_m[0] == rhs.state_m[0] && + lhs.state_m[1] > rhs.state_m[1]) return true; else if (lhs.state_m[0] == rhs.state_m[0] && lhs.state_m[1] == rhs.state_m[1] && @@ -202,7 +244,7 @@ bool ThreeDGrid::isGreater // remove *map if it exists; delete this if there are no more VectorMaps void ThreeDGrid::remove(VectorMap* map) { std::vector<VectorMap*>::iterator it = - std::find(maps_m.begin(), maps_m.end(), map); + std::find(maps_m.begin(), maps_m.end(), map); if (it < maps_m.end()) { maps_m.erase(it); } @@ -214,7 +256,7 @@ void ThreeDGrid::remove(VectorMap* map) { // add *map if it has not already been added void ThreeDGrid::add(VectorMap* map) { std::vector<VectorMap*>::iterator it = - std::find(maps_m.begin(), maps_m.end(), map); + std::find(maps_m.begin(), maps_m.end(), map); if (it == maps_m.end()) { maps_m.push_back(map); } @@ -273,3 +315,5 @@ Mesh::Iterator ThreeDGrid::getNearest(const double* position) const { index[2] = zSize_m; return Mesh::Iterator(index, this); } +} + diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.h b/classic/5.0/src/Fields/Interpolation/ThreeDGrid.h similarity index 89% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.h rename to classic/5.0/src/Fields/Interpolation/ThreeDGrid.h index 91f06947167ace1c5e0a605e2ea0de38f6a718eb..fa047b22deb59dfe6ab224200f70d18abd21be65 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/ThreeDGrid.h +++ b/classic/5.0/src/Fields/Interpolation/ThreeDGrid.h @@ -33,7 +33,9 @@ #include <algorithm> #include <vector> -#include "Fields/SectorMagneticFieldMap/Mesh.h" +#include "Fields/Interpolation/Mesh.h" + +namespace interpolation { /** ThreeDGrid holds grid information for a rectangular grid used in e.g. * interpolation @@ -61,7 +63,7 @@ class ThreeDGrid : public Mesh { Mesh * clone() {return new ThreeDGrid(*this);} /** Not implemented (returns NULL) */ - Mesh * dual () {return NULL;} + Mesh * dual () const; /** Default constructor */ ThreeDGrid(); @@ -291,10 +293,22 @@ class ThreeDGrid : public Mesh { /** Return the point in the mesh nearest to position */ Mesh::Iterator getNearest(const double* position) const; + /** Custom LowerBound routine like std::lower_bound + * + * Make a binary search to find the element in vec with vec[i] <= x. + * \param vec STL vector object (could make this any sorted iterator...) + * \param x object for which to search. + * \param index vectorLowerBound sets index to the position of the element. If + * x < vec[0], vectorLowerBound fills with -1. + */ + static void vectorLowerBound(std::vector<double> vec, double x, int& index); + protected: // Change position - virtual Mesh::Iterator& addEquals(Mesh::Iterator& lhs, int difference) const; - virtual Mesh::Iterator& subEquals(Mesh::Iterator& lhs, int difference) const; + virtual Mesh::Iterator& addEquals(Mesh::Iterator& lhs, + int difference) const; + virtual Mesh::Iterator& subEquals(Mesh::Iterator& lhs, + int difference) const; virtual Mesh::Iterator& addEquals (Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const; virtual Mesh::Iterator& subEquals @@ -328,12 +342,18 @@ class ThreeDGrid : public Mesh { friend Mesh::Iterator& operator+= (Mesh::Iterator& lhs, const Mesh::Iterator& rhs); - friend bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); - friend bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); - friend bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); - friend bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); - friend bool operator< (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); - friend bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs); + friend bool operator==(const Mesh::Iterator& lhs, + const Mesh::Iterator& rhs); + friend bool operator!=(const Mesh::Iterator& lhs, + const Mesh::Iterator& rhs); + friend bool operator>=(const Mesh::Iterator& lhs, + const Mesh::Iterator& rhs); + friend bool operator<=(const Mesh::Iterator& lhs, + const Mesh::Iterator& rhs); + friend bool operator< (const Mesh::Iterator& lhs, + const Mesh::Iterator& rhs); + friend bool operator> (const Mesh::Iterator& lhs, + const Mesh::Iterator& rhs); }; double* ThreeDGrid::newXArray() { @@ -357,27 +377,6 @@ double* ThreeDGrid::newZArray() { return z; } -void ThreeDGrid::xLowerBound(const double& x, int& xIndex) const { - if (constantSpacing_m) - xIndex = static_cast<int>(floor((x - x_m[0])/(x_m[1]-x_m[0]) )); - else - xIndex = std::lower_bound(x_m.begin(), x_m.end(), x)-x_m.begin()-1; -} - -void ThreeDGrid::yLowerBound(const double& y, int& yIndex) const { - if (constantSpacing_m) - yIndex = static_cast<int>(floor((y - y_m[0])/(y_m[1]-y_m[0]))); - else - yIndex = std::lower_bound(y_m.begin(), y_m.end(), y)-y_m.begin()-1; -} - -void ThreeDGrid::zLowerBound(const double& z, int& zIndex) const { - if (constantSpacing_m) - zIndex = static_cast<int>(floor((z - z_m[0])/(z_m[1]-z_m[0]))); - else - zIndex = std::lower_bound(z_m.begin(), z_m.end(), z)-z_m.begin()-1; -} - void ThreeDGrid::lowerBound(const double& x, int& xIndex, const double& y, int& yIndex, const double& z, int& zIndex) const { @@ -387,9 +386,9 @@ void ThreeDGrid::lowerBound(const double& x, int& xIndex, } void ThreeDGrid::lowerBound(const double& x, - const double& y, - const double& z, - Mesh::Iterator& it) const { + const double& y, + const double& z, + Mesh::Iterator& it) const { xLowerBound(x, it[0]); yLowerBound(y, it[1]); zLowerBound(z, it[2]); @@ -398,6 +397,29 @@ void ThreeDGrid::lowerBound(const double& x, it[2]++; } +void ThreeDGrid::xLowerBound(const double& x, int& xIndex) const { + if (constantSpacing_m) + xIndex = static_cast<int>(floor((x - x_m[0])/(x_m[1]-x_m[0]) )); + else { + vectorLowerBound(x_m, x, xIndex); + } +} + + +void ThreeDGrid::yLowerBound(const double& y, int& yIndex) const { + if (constantSpacing_m) + yIndex = static_cast<int>(floor((y - y_m[0])/(y_m[1]-y_m[0]))); + else + vectorLowerBound(y_m, y, yIndex); +} + +void ThreeDGrid::zLowerBound(const double& z, int& zIndex) const { + if (constantSpacing_m) + zIndex = static_cast<int>(floor((z - z_m[0])/(z_m[1]-z_m[0]))); + else + vectorLowerBound(z_m, z, zIndex); +} + double ThreeDGrid::minX() const { return x_m[0]; } @@ -451,6 +473,6 @@ void ThreeDGrid::setConstantSpacing(bool spacing) { bool ThreeDGrid::getConstantSpacing() const { return constantSpacing_m; } - +} #endif diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.cpp b/classic/5.0/src/Fields/Interpolation/TriLinearInterpolator.cpp similarity index 97% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.cpp rename to classic/5.0/src/Fields/Interpolation/TriLinearInterpolator.cpp index 0ebf9ac170f2b04ceb7c785b06f296f93f875ef9..bf897f54387d2f9dd6f39d9c7e390b26a4a202ae 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.cpp +++ b/classic/5.0/src/Fields/Interpolation/TriLinearInterpolator.cpp @@ -25,7 +25,9 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/TriLinearInterpolator.h" +#include "Fields/Interpolation/TriLinearInterpolator.h" + +namespace interpolation { TriLinearInterpolator::TriLinearInterpolator(const TriLinearInterpolator& lhs) { coordinates_m = new ThreeDGrid(*lhs.coordinates_m); @@ -75,4 +77,5 @@ void TriLinearInterpolator::function (coordinates_m->z(k+2) - coordinates_m->z(k+1))* (Point[2] - coordinates_m->z(k+1)) + f_xy[0]; } +} diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.h b/classic/5.0/src/Fields/Interpolation/TriLinearInterpolator.h similarity index 97% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.h rename to classic/5.0/src/Fields/Interpolation/TriLinearInterpolator.h index b7e5c3c15e1b0a6de5841f4330755d40bd12ff96..cea8eb2821a1bb5d760f43ea6acb6ef859416a63 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/TriLinearInterpolator.h +++ b/classic/5.0/src/Fields/Interpolation/TriLinearInterpolator.h @@ -28,7 +28,8 @@ #ifndef _SRC_LEGACY_INTERFACE_INTERPOLATION_TRILINEARINTERPOLATOR_HH_ #define _SRC_LEGACY_INTERFACE_INTERPOLATION_TRILINEARINTERPOLATOR_HH_ -#include "Fields/SectorMagneticFieldMap/Interpolator3dGridTo1d.h" +#include "Fields/Interpolation/Interpolator3dGridTo1d.h" +namespace interpolation { /** TriLinearInterpolator performs a linear interpolation in x then y then z * @@ -89,7 +90,7 @@ TriLinearInterpolator::~TriLinearInterpolator() { TriLinearInterpolator* TriLinearInterpolator::clone() const { return new TriLinearInterpolator(*this); } - +} #endif diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.cpp b/classic/5.0/src/Fields/Interpolation/VectorMap.cpp similarity index 96% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.cpp rename to classic/5.0/src/Fields/Interpolation/VectorMap.cpp index c03983f2e7ffa8c043431273027d77ce8c6a1617..f7676f7900dd52606ba0ddf040ae26ecc33e4cce 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.cpp +++ b/classic/5.0/src/Fields/Interpolation/VectorMap.cpp @@ -25,5 +25,5 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/VectorMap.h" +#include "Fields/Interpolation/VectorMap.h" diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.h b/classic/5.0/src/Fields/Interpolation/VectorMap.h similarity index 97% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.h rename to classic/5.0/src/Fields/Interpolation/VectorMap.h index 9c286d1861591d56609b2544c3e579812ce641f2..f93f7388dd25090e371836c48d7e829334cde2cc 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/VectorMap.h +++ b/classic/5.0/src/Fields/Interpolation/VectorMap.h @@ -30,7 +30,9 @@ #include <vector> -#include "Fields/SectorMagneticFieldMap/ThreeDGrid.h" +#include "Fields/Interpolation/ThreeDGrid.h" + +namespace interpolation { /** VectorMap is an abstract class that defines mapping from one vector to * another. @@ -86,7 +88,7 @@ class VectorMap { virtual ~VectorMap() {;} /** Return the mesh used by the vector map or NULL if no mesh */ - virtual ThreeDGrid* getMesh() const {return NULL;} + virtual Mesh* getMesh() const {return NULL;} private: }; @@ -113,5 +115,5 @@ void VectorMap::functionAppend function(&point_vec[i][0], &value_vec[i][0]); } } - +} #endif // _CLASSIC_FIELDS_VECTORMAP_HH_ diff --git a/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.aux b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.aux new file mode 100644 index 0000000000000000000000000000000000000000..8428308c1faa71561f930384a388cbfccac2dc8f --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.aux @@ -0,0 +1,6 @@ +\relax +\@writefile{toc}{\contentsline {section}{\numberline {1}Arbitrary Order, Arbitrary Dimensional Polynomial Solve With Smoothing}{1}} +\@writefile{toc}{\contentsline {subsection}{\numberline {1.1}Generalised Indexing and Notation}{1}} +\newlabel{eq:quadratic}{{1}{1}} +\@writefile{toc}{\contentsline {subsection}{\numberline {1.2}Polynomial Solve with no Smoothing}{2}} +\@writefile{toc}{\contentsline {subsection}{\numberline {1.3}Polynomial Solve with Smoothing}{2}} diff --git a/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.log b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.log new file mode 100644 index 0000000000000000000000000000000000000000..c203c5ea84bf6941128b032857870c27d35d0a57 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.log @@ -0,0 +1,173 @@ +This is pdfTeX, Version 3.1415926-2.5-1.40.14 (TeX Live 2013/Debian) (format=pdflatex 2014.9.17) 1 MAY 2015 17:54 +entering extended mode + restricted \write18 enabled. + %&-line parsing enabled. +**polynomial_solve_with_smoothing.tex +(./polynomial_solve_with_smoothing.tex +LaTeX2e <2011/06/27> +Babel <3.9h> and hyphenation patterns for 2 languages loaded. +(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls +Document Class: article 2007/10/19 v1.4h Standard LaTeX document class +(/usr/share/texlive/texmf-dist/tex/latex/base/size10.clo +File: size10.clo 2007/10/19 v1.4h Standard LaTeX file (size option) +) +\c@part=\count79 +\c@section=\count80 +\c@subsection=\count81 +\c@subsubsection=\count82 +\c@paragraph=\count83 +\c@subparagraph=\count84 +\c@figure=\count85 +\c@table=\count86 +\abovecaptionskip=\skip41 +\belowcaptionskip=\skip42 +\bibindent=\dimen102 +) +(/usr/share/texlive/texmf-dist/tex/latex/graphics/graphicx.sty +Package: graphicx 1999/02/16 v1.0f Enhanced LaTeX Graphics (DPC,SPQR) + +(/usr/share/texlive/texmf-dist/tex/latex/graphics/keyval.sty +Package: keyval 1999/03/16 v1.13 key=value parser (DPC) +\KV@toks@=\toks14 +) +(/usr/share/texlive/texmf-dist/tex/latex/graphics/graphics.sty +Package: graphics 2009/02/05 v1.0o Standard LaTeX Graphics (DPC,SPQR) + +(/usr/share/texlive/texmf-dist/tex/latex/graphics/trig.sty +Package: trig 1999/03/16 v1.09 sin cos tan (DPC) +) +(/usr/share/texlive/texmf-dist/tex/latex/latexconfig/graphics.cfg +File: graphics.cfg 2010/04/23 v1.9 graphics configuration of TeX Live +) +Package graphics Info: Driver file: pdftex.def on input line 91. + +(/usr/share/texlive/texmf-dist/tex/latex/pdftex-def/pdftex.def +File: pdftex.def 2011/05/27 v0.06d Graphics/color for pdfTeX + +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/infwarerr.sty +Package: infwarerr 2010/04/08 v1.3 Providing info/warning/error messages (HO) +) +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ltxcmds.sty +Package: ltxcmds 2011/11/09 v1.22 LaTeX kernel commands for general use (HO) +) +\Gread@gobject=\count87 +)) +\Gin@req@height=\dimen103 +\Gin@req@width=\dimen104 +) +(./polynomial_solve_with_smoothing.aux) +\openout1 = `polynomial_solve_with_smoothing.aux'. + +LaTeX Font Info: Checking defaults for OML/cmm/m/it on input line 4. +LaTeX Font Info: ... okay on input line 4. +LaTeX Font Info: Checking defaults for T1/cmr/m/n on input line 4. +LaTeX Font Info: ... okay on input line 4. +LaTeX Font Info: Checking defaults for OT1/cmr/m/n on input line 4. +LaTeX Font Info: ... okay on input line 4. +LaTeX Font Info: Checking defaults for OMS/cmsy/m/n on input line 4. +LaTeX Font Info: ... okay on input line 4. +LaTeX Font Info: Checking defaults for OMX/cmex/m/n on input line 4. +LaTeX Font Info: ... okay on input line 4. +LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 4. +LaTeX Font Info: ... okay on input line 4. + +(/usr/share/texlive/texmf-dist/tex/context/base/supp-pdf.mkii +[Loading MPS to PDF converter (version 2006.09.02).] +\scratchcounter=\count88 +\scratchdimen=\dimen105 +\scratchbox=\box26 +\nofMPsegments=\count89 +\nofMParguments=\count90 +\everyMPshowfont=\toks15 +\MPscratchCnt=\count91 +\MPscratchDim=\dimen106 +\MPnumerator=\count92 +\makeMPintoPDFobject=\count93 +\everyMPtoPDFconversion=\toks16 +) (/usr/share/texlive/texmf-dist/tex/generic/oberdiek/pdftexcmds.sty +Package: pdftexcmds 2011/11/29 v0.20 Utility functions of pdfTeX for LuaTeX (HO +) + +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifluatex.sty +Package: ifluatex 2010/03/01 v1.3 Provides the ifluatex switch (HO) +Package ifluatex Info: LuaTeX not detected. +) +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/ifpdf.sty +Package: ifpdf 2011/01/30 v2.3 Provides the ifpdf switch (HO) +Package ifpdf Info: pdfTeX in PDF mode is detected. +) +Package pdftexcmds Info: LuaTeX not detected. +Package pdftexcmds Info: \pdf@primitive is available. +Package pdftexcmds Info: \pdf@ifprimitive is available. +Package pdftexcmds Info: \pdfdraftmode found. +) +(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/epstopdf-base.sty +Package: epstopdf-base 2010/02/09 v2.5 Base part for package epstopdf + +(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/grfext.sty +Package: grfext 2010/08/19 v1.1 Manage graphics extensions (HO) + +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/kvdefinekeys.sty +Package: kvdefinekeys 2011/04/07 v1.3 Define keys (HO) +)) +(/usr/share/texlive/texmf-dist/tex/latex/oberdiek/kvoptions.sty +Package: kvoptions 2011/06/30 v3.11 Key value format for package options (HO) + +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/kvsetkeys.sty +Package: kvsetkeys 2012/04/25 v1.16 Key value parser (HO) + +(/usr/share/texlive/texmf-dist/tex/generic/oberdiek/etexcmds.sty +Package: etexcmds 2011/02/16 v1.5 Avoid name clashes with e-TeX commands (HO) +Package etexcmds Info: Could not find \expanded. +(etexcmds) That can mean that you are not using pdfTeX 1.50 or +(etexcmds) that some package has redefined \expanded. +(etexcmds) In the latter case, load this package earlier. +))) +Package grfext Info: Graphics extension search list: +(grfext) [.png,.pdf,.jpg,.mps,.jpeg,.jbig2,.jb2,.PNG,.PDF,.JPG,.JPE +G,.JBIG2,.JB2,.eps] +(grfext) \AppendGraphicsExtensions on input line 452. + +(/usr/share/texlive/texmf-dist/tex/latex/latexconfig/epstopdf-sys.cfg +File: epstopdf-sys.cfg 2010/07/13 v1.3 Configuration of (r)epstopdf for TeX Liv +e +)) +Overfull \hbox (1.308pt too wide) in paragraph at lines 6--6 +[]\OT1/cmr/bx/n/14.4 Arbitrary Or-der, Ar-bi-trary Di-men-sional Poly- + [] + +LaTeX Font Info: External font `cmex10' loaded for size +(Font) <7> on input line 19. +LaTeX Font Info: External font `cmex10' loaded for size +(Font) <5> on input line 19. +[1 + +{/var/lib/texmf/fonts/map/pdftex/updmap/pdftex.map}] [2] +(./polynomial_solve_with_smoothing.aux) ) +Here is how much of TeX's memory you used: + 1473 strings out of 495028 + 21242 string characters out of 6181497 + 77159 words of memory out of 5000000 + 4699 multiletter control sequences out of 15000+600000 + 6123 words of font info for 22 fonts, out of 8000000 for 9000 + 14 hyphenation exceptions out of 8191 + 37i,5n,23p,257b,115s stack positions out of 5000i,500n,10000p,200000b,80000s +</usr/share/texlive/texmf-dist/fonts/t +ype1/public/amsfonts/cm/cmbx10.pfb></usr/share/texlive/texmf-dist/fonts/type1/p +ublic/amsfonts/cm/cmbx12.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/ +amsfonts/cm/cmex10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfon +ts/cm/cmmi10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/ +cmmi5.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmmi7.p +fb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmr10.pfb></us +r/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmr5.pfb></usr/share/ +texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmr7.pfb></usr/share/texlive/ +texmf-dist/fonts/type1/public/amsfonts/cm/cmsy10.pfb></usr/share/texlive/texmf- +dist/fonts/type1/public/amsfonts/cm/cmsy7.pfb></usr/share/texlive/texmf-dist/fo +nts/type1/public/amsfonts/cm/cmti10.pfb> +Output written on polynomial_solve_with_smoothing.pdf (2 pages, 124220 bytes). +PDF statistics: + 59 PDF objects out of 1000 (max. 8388607) + 42 compressed objects within 1 object stream + 0 named destinations out of 1000 (max. 500000) + 1 words of extra memory for PDF output out of 10000 (max. 10000000) + diff --git a/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.pdf b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.pdf new file mode 100644 index 0000000000000000000000000000000000000000..28c5ad17a31dcc2c69ef5bf4749685aa97c29351 Binary files /dev/null and b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.pdf differ diff --git a/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.tex b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.tex new file mode 100644 index 0000000000000000000000000000000000000000..86ab7f461ae5226d8512c4ac796ff568a95a5d87 --- /dev/null +++ b/classic/5.0/src/Fields/Interpolation/polynomial_solve_with_smoothing.tex @@ -0,0 +1,119 @@ +\documentclass{article} +\usepackage{graphicx} + +\begin{document} + +\section{Arbitrary Order, Arbitrary Dimensional Polynomial Solve With Smoothing} + +\emph{Chris Rogers, STFC Rutherford Appleton Laboratory, 2015} + +Below I outline the mathematical foundation for higher order polynomial solving. + +\subsection{Generalised Indexing and Notation} + +The polynomial solve is a quite straightforward consequence of a simultaneous +equation solve. The only really tricky thing is one of notation. For a +polynomial in higher dimensional space, we have to be quite careful how things +are indexed. For example, consider the generalised quadratic polynomial in two +dimensions +\begin{equation} +\label{eq:quadratic} +y = a_{00} + a_{10} x_0 + a_{01} x_1 + a_{11} x_0 x_1 + a_{20} x_0^2 + a_{02} x_1^2 +\end{equation} + +The polynomial coefficients $a_{jk}$ have been indexed by the power on the +corresponding products of $x_i$. Explicitly, + +\begin{equation} +y = a_{00} x_0^0 x_1^0 + a_{10} x_0^1 x_1^0 + a_{01} x_0^0 x_1^1 + a_{11} x_0^1 x_1^1 + a_{20} x_0^2 x_1^0 + a_{02} x_0^0 x_1^2 +\end{equation} + +which is identical to equation \ref{eq:quadratic}. In the code, this is termed +\emph{index by power}. Occasionally we also use the concept of \emph{index by +vector}. In this indexing scheme, we index the $i$ in the product of $x_i$, so +for example the quadratic equation above becomes + +\begin{equation} +y = b+ b_{0} x_0 + b_{1} x_1 + b_{01} x_0 x_1 + b_{00} x_0 x_0 + b_{11} x_1 x_1 +\end{equation} + +Conventionally, a quadratic equation in higher dimension is one in which the sum +of the powers $<=$ 2, i.e. the sum of \emph{index by power} $<=$ 2 and the +length of \emph{index by vector} $<=$ 2. + +For reasons discussed below, in this note a polynomial of order $i$ will be one +in which no term in \emph{index by power} is more than $i$, or equivalently no +integer in \emph{index by vector} is repeated more than $i$ times. + +The coefficients $a_{00}, a_{01}, \ldots$ can be written in shorthand as +$a_{\vec{p}}$ where $\vec{p}$ is a vector of integers as discussed above. + +This notation can also be used to extend to derivatives. Consider the derivative +\begin{equation} +\frac{\partial^ny_i}{\partial x_0^{q_0} \partial x_1^{q_1} \ldots} +\end{equation} +which can be written in shorthand as +\begin{equation} +y^{(\vec{q})} +\end{equation} + +\subsection{Polynomial Solve with no Smoothing} +Consider a polynomial in variables $\vec{x} = x_0, x_1, \ldots x_n$ such that +\begin{equation} +\vec{y}(x) = \sum_{\vec{p}} a_{\vec{p}} \prod^{j=n}_{j=1} x_j^{p_j} +\end{equation} +If we have a set of known values $\vec{y}_{\vec{b}}(\vec{x}_{\vec{b}}$ at some +known positions $\vec{x}_{\vec{b}}$ e.g. on a rectangular grid; then we can +solve for $a_{\vec{p}}$ in the usual way as a linear system of simultaneous +equations. + +Then +\begin{equation} +\vec{y}_{\vec{b}}(\vec{x}_{\vec{b}}) = \sum_{\vec{p}} a_{\vec{p}} \prod^{j=n}_{j=1} x_{b_j}^{p_j} +\end{equation} +where the measured values to be fitted, $\vec{y}_{\vec{b}}$ are known, the +positions at which those values were measured, $\vec{x}$, are known but the +polynomial coefficients $a_{\vec{p}}$ are not known. + +This can be written as a system of linear equations like +\begin{equation} +\vec{G} = \mathbf{H} \vec{A} +\end{equation} + +where $\vec{G}$ is the vector of $\vec{y}_{\vec{b}}$, $\vec{a}$ is the vector of +$a_{\vec{p}}$ and $\mathbf{H}$ is the matrix of $\prod x_{b_j}^{p_j}$. + +In order for the polynomial fit to be successful, $\mathbf{H}$ must be invertible, +so the points $\vec{x}_{\vec{b}}$ need to be chosen carefully. Otherwise no +presumption has been made on the dimension or order of the fit. + +\subsection{Polynomial Solve with Smoothing} +Derivatives of the polynomial may be known, in which case they can be used in +addition to measured points. + +Consider the measured derivatives +$\vec{y}_{\vec{b}}^{(\vec{q})}(\vec{x}_{\vec{b}})$. Then +\begin{equation} +y^{\vec{q}}_{\vec{b}} = \sum_{\vec{p}} a_{\vec{p}} + \prod^j \frac{p_j!}{(p_j-q_j)!} \vec{x}_{\vec{b}}^{(p_j-q_j)} +\end{equation} +The simultaneous equation can then be reformulated as +\begin{equation} +\vec{G'} = \mathbf{H'} \vec{A} +\end{equation} +where $\vec{G'}$ is the vector of $y^{\vec{p}}$ and $y^{(\vec{q})}$; $\mathbf{H'}$ +is the matrix of $\prod x_{b_j}^{p_j}$ and $\prod p_j!/(p_j-q_j)! x_{b_j}^{p_j-q_j}$ + +In order for the polynomial fit to be successful, $\mathbf{H'}$ must be +invertible. + +\subsection{Issues} +In this note, no attempt is made to demonstrate that a given set of points +result in an invertible matrix. This is found empirically. So there may be a set +of pathological cases that break the fitting. + +Boundary conditions can be an issue. At the moment, boundary condition is taken +as derivative is zero. + +\end{document} + diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.cpp b/classic/5.0/src/Fields/SectorField.cpp similarity index 99% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.cpp rename to classic/5.0/src/Fields/SectorField.cpp index 8f0c20c10594f72b8ca606f281650d9871ce68f9..3bd2f5198d655c39c923a1420267ce087bc30ffe 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.cpp +++ b/classic/5.0/src/Fields/SectorField.cpp @@ -25,7 +25,7 @@ * POSSIBILITY OF SUCH DAMAGE. */ -#include "Fields/SectorMagneticFieldMap/SectorField.h" +#include "Fields/SectorField.h" #include <math.h> diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.h b/classic/5.0/src/Fields/SectorField.h similarity index 100% rename from classic/5.0/src/Fields/SectorMagneticFieldMap/SectorField.h rename to classic/5.0/src/Fields/SectorField.h diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap.cpp b/classic/5.0/src/Fields/SectorMagneticFieldMap.cpp index 987ac7f00fa7e8ceb6f8ca7662a524aa7e49e744..dfb1201f484079ba87969b9bc581faf32d09783b 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap.cpp +++ b/classic/5.0/src/Fields/SectorMagneticFieldMap.cpp @@ -25,10 +25,17 @@ * POSSIBILITY OF SUCH DAMAGE. */ + #include "Fields/SectorMagneticFieldMap.h" -#include "Fields/SectorMagneticFieldMap/Mesh.h" -#include "Fields/SectorMagneticFieldMap/ThreeDGrid.h" -#include "Fields/SectorMagneticFieldMap/Interpolator3dGridTo3d.h" +// Grid on which field values are stored +#include "Fields/Interpolation/Mesh.h" +#include "Fields/Interpolation/ThreeDGrid.h" +// Linear interpolation routines +#include "Fields/Interpolation/VectorMap.h" +#include "Fields/Interpolation/Interpolator3dGridTo3d.h" +// Higher order interpolation routines +#include "Fields/Interpolation/PolynomialPatch.h" +#include "Fields/Interpolation/PPSolveFactory.h" #include "Utilities/LogicalError.h" @@ -40,18 +47,23 @@ #include <fstream> #include <string> +using namespace interpolation; + // allow a fairly generous phi tolerance - we don't care about phi much and // calculation can be flaky due to ascii truncation of double and conversions // from polar to Cartesian -const double SectorMagneticFieldMap::fractionalBBPhiTolerance_m(1e-2); +const double SectorMagneticFieldMap::fractionalBBPhiTolerance_m(1e-12); std::map<std::string, SectorMagneticFieldMap*> SectorMagneticFieldMap::_fields; SectorMagneticFieldMap::SectorMagneticFieldMap(std::string file_name, std::string symmetry, double length_units, - double field_units) + double field_units, + int polynomial_order, + int smoothing_order) : SectorField(file_name), interpolator_m(NULL), symmetry_m(dipole), - units_m(6, 1.), filename_m(file_name), phiOffset_m(0.) { + units_m(6, 1.), filename_m(file_name), phiOffset_m(0.), + poly_order_m(polynomial_order), smoothing_order_m(smoothing_order) { units_m[0] *= length_units; units_m[1] *= length_units; units_m[2] *= length_units; @@ -64,7 +76,8 @@ SectorMagneticFieldMap::SectorMagneticFieldMap(std::string file_name, _fields[file_name] = new SectorMagneticFieldMap(*this); } else { SectorMagneticFieldMap* tgt = _fields[file_name]; - if (symmetry_m != tgt->symmetry_m || units_m != tgt->units_m || + if (symmetry_m != tgt->symmetry_m || + units_m != tgt->units_m || filename_m != tgt->filename_m) { throw(LogicalError( "SectorMagneticFieldMap::SectorMagneticFieldMap", @@ -72,46 +85,64 @@ SectorMagneticFieldMap::SectorMagneticFieldMap(std::string file_name, std::string("file but different settings") )); } - setInterpolator(tgt->interpolator_m->clone()); + setInterpolator(tgt->interpolator_m); } - phiOffset_m = -interpolator_m->getMesh()->minZ(); + ThreeDGrid* grid = dynamic_cast<ThreeDGrid*>(interpolator_m->getMesh()); + phiOffset_m = -grid->minZ(); } SectorMagneticFieldMap::SectorMagneticFieldMap - (const SectorMagneticFieldMap& field) + (const SectorMagneticFieldMap& field) : SectorField(field), interpolator_m(NULL), symmetry_m(field.symmetry_m), units_m(field.units_m), filename_m(field.filename_m), phiOffset_m(field.phiOffset_m) { - Interpolator3dGridTo3d* interpolator = NULL; + VectorMap* interpolator = NULL; if (field.interpolator_m != NULL) { - interpolator = field.interpolator_m->clone(); + interpolator = field.interpolator_m; } setInterpolator(interpolator); } SectorMagneticFieldMap::~SectorMagneticFieldMap() { - delete interpolator_m; + // delete interpolator_m; We reuse the interpolators... hack! should use smart pointer } -Interpolator3dGridTo3d* SectorMagneticFieldMap::getInterpolator() { +VectorMap* SectorMagneticFieldMap::getInterpolator() { return interpolator_m; } -void SectorMagneticFieldMap::setInterpolator - (Interpolator3dGridTo3d* interpolator) { +void SectorMagneticFieldMap::setInterpolator(VectorMap* interpolator) { if (interpolator_m != NULL) { delete interpolator_m; } interpolator_m = interpolator; if (interpolator_m != NULL) { - ThreeDGrid* grid = interpolator_m->getMesh(); + ThreeDGrid* fp_grid = dynamic_cast<ThreeDGrid*>(interpolator_m->getMesh()); + + if (interpolator_m->getPointDimension() != 3) + throw(LogicalError( + "SectorMagneticFieldMap::setInterpolator", + "Attempt to load interpolator with PointDimension != 3") + ); + if (interpolator_m->getValueDimension() != 3) + throw(LogicalError( + "SectorMagneticFieldMap::setInterpolator", + "Attempt to load interpolator with ValueDimension != 3" + )); + ThreeDGrid* grid = dynamic_cast<ThreeDGrid*>(interpolator_m->getMesh()); + if (grid == NULL) + throw(LogicalError( + "SectorMagneticFieldMap::setInterpolator", + "Attempt to load interpolator with grid not ThreeDGrid" + )); double phiTol = (grid->maxZ()-grid->minZ())*fractionalBBPhiTolerance_m; SectorField::setPolarBoundingBox - (grid->minX(), grid->minY(), grid->minZ(), + (grid->minX(), grid->minY()-grid->maxY(), grid->minZ(), grid->maxX(), grid->maxY(), grid->maxZ(), - 0., 0., phiTol); + 0., 0., 0.); } - delete interpolator_m->clone(); // why? + Mesh* mesh = interpolator_m->getMesh(); + print(std::cerr); } std::string SectorMagneticFieldMap::getSymmetry() const { @@ -123,7 +154,8 @@ void SectorMagneticFieldMap::setSymmetry(std::string name) { } void SectorMagneticFieldMap::readMap() { - setInterpolator(IO::readMap(filename_m, units_m, symmetry_m)); + setInterpolator(IO::readMap + (filename_m, units_m, symmetry_m, poly_order_m, smoothing_order_m)); } void SectorMagneticFieldMap::freeMap() { @@ -180,32 +212,44 @@ bool SectorMagneticFieldMap::getFieldstrength // correct sense but this should be implemented in the end by Ring (i.e. // should be able to do off-midplane transformations) double radius = (getPolarBoundingBoxMin()[0]+getPolarBoundingBoxMax()[0])/2; + double midplane = (getPolarBoundingBoxMin()[1]+getPolarBoundingBoxMax()[1])/2; double R_temp[3] = {R_c(0)+radius, R_c(1), R_c(2)}; double B_temp[3] = {0., 0., 0.}; SectorField::convertToPolar(R_temp); + if (R_temp[1] < midplane) { + R_temp[1] = midplane + (midplane - R_temp[1]); + } // interpolator has phi in 0. to dphi R_temp[2] -= phiOffset_m; + // bool symmetryWasApplied = applySymmetry(R_temp); if (!isInBoundingBox(R_temp)) { - // std::cerr << "SectorMagneticFieldMap r: " << radius << " R_c " << R_c << " R_temp " - // << R_temp[0] << " " << R_temp[1] << " " << R_temp[2] - // << std::endl; return true; } interpolator_m->function(R_temp, B_temp); // and here we transform back // we want phi in 0. to dphi - R_temp[2] += phiOffset_m; - SectorField::convertToCartesian(R_temp, B_temp); - B_c(0) = B_temp[0]; - B_c(1) = B_temp[1]; - B_c(2) = B_temp[2]; - // B_c *= 10.; - // std::cerr << "\nSMFM::getFieldStrength r: " << radius << " R_c: " << R_c - // << " R_p: " << R_temp[0] << " " << R_temp[1] << " " << R_temp[2] - // << " B: " << B_c[0] << " " << B_c[1] << " " << B_c[2] << std::endl; + if (R_temp[1] > midplane) { + B_temp[0] *= +1; // reflect Bx + B_temp[2] *= +1; // reflect Bz + } + if (R_temp[1] < midplane) { + B_temp[0] *= -1; // reflect Bx + B_temp[2] *= -1; // reflect Bz + } + B_c(0) = B_temp[0]*cos(phiOffset_m)-B_temp[2]*sin(phiOffset_m); // x + B_c(1) = B_temp[1]; // axial + B_c(2) = B_temp[0]*sin(phiOffset_m)+B_temp[2]*cos(phiOffset_m); // z return false; } +bool SectorMagneticFieldMap::applySymmetry(double* R_temp) const { + double ymin = SectorField::getPolarBoundingBoxMin()[1]; + if (symmetry_m == dipole && R_temp[1] <= ymin) { + R_temp[1] = 2*ymin-R_temp[1]; + return true; + } + return false; +} void SectorMagneticFieldMap::clearFieldCache() { for (std::map<std::string, SectorMagneticFieldMap*>::iterator it = @@ -219,8 +263,9 @@ void SectorMagneticFieldMap::getInfo(Inform* msg) { std::vector<double> bbmin = SectorField::getPolarBoundingBoxMin(); std::vector<double> bbmax = SectorField::getPolarBoundingBoxMax(); (*msg) << Filename_m << " (3D Sector Magnetostatic); " - << "zini= " << bbmin[2] << " m; zfinal= " << bbmax[2] << " m;" - << "rini= " << bbmin[0] << " m; rfinal= " << bbmax[0] << " m;" + << "zini= " << bbmin[2] << " mm; zfinal= " << bbmax[2] << " mm;" + << "phiini= " << bbmin[1] << " rad; phifinal= " << bbmax[1] << " rad;" + << "rini= " << bbmin[0] << " mm; rfinal= " << bbmax[0] << " mm;" << endl; } @@ -229,6 +274,7 @@ void SectorMagneticFieldMap::print(std::ostream& out) { std::vector<double> bbmax = SectorField::getPolarBoundingBoxMax(); out << Filename_m << " (3D Sector Magnetostatic); " << "zini= " << bbmin[2] << " m; zfinal= " << bbmax[2] << " m;" + << "phiini= " << bbmin[1] << " rad; phifinal= " << bbmax[1] << " rad;" << "rini= " << bbmin[0] << " m; rfinal= " << bbmax[0] << " m;" << std::endl; } @@ -243,25 +289,39 @@ void SectorMagneticFieldMap::swap() {} double SectorMagneticFieldMap::getDeltaPhi() const { - ThreeDGrid* grid = interpolator_m->getMesh(); + ThreeDGrid* grid = reinterpret_cast<ThreeDGrid*>(interpolator_m->getMesh()); return grid->maxZ() - grid->minZ(); } const double SectorMagneticFieldMap::IO::floatTolerance_m = 1e-3; const int SectorMagneticFieldMap::IO::sortOrder_m[3] = {0, 1, 2}; -Interpolator3dGridTo3d* SectorMagneticFieldMap::IO::readMap - (std::string file_name, std::vector<double> units, - SectorMagneticFieldMap::symmetry sym) { +VectorMap* SectorMagneticFieldMap::IO::readMap( + std::string file_name, + std::vector<double> units, + SectorMagneticFieldMap::symmetry sym, + int polynomial_order, + int smoothing_order) { try { - INFOMSG("Opening sector field map " << file_name << endl); + INFOMSG("Opening sector field map " << file_name + << " fit order " << polynomial_order + << " smoothing order " << smoothing_order << endl); // get raw data std::vector< std::vector<double> > field_points = readLines (file_name, units); // build grid ThreeDGrid* grid = generateGrid(field_points, sym); // build interpolator (convert grid to useful units) - return getInterpolator(field_points, grid, sym); + if (polynomial_order == 1 && smoothing_order == 1) { + return getInterpolator(field_points, grid, sym); + } else { + return getInterpolatorPolyPatch( + field_points, + grid, + sym, + polynomial_order, + smoothing_order); + } } catch(std::exception& exc) { throw(LogicalError( "SectorMagneticFieldMap::IO::ReadMap", @@ -271,14 +331,42 @@ Interpolator3dGridTo3d* SectorMagneticFieldMap::IO::readMap return NULL; } -Interpolator3dGridTo3d* SectorMagneticFieldMap::IO::getInterpolator +VectorMap* SectorMagneticFieldMap::IO::getInterpolatorPolyPatch( + std::vector< std::vector<double> > field_points, + ThreeDGrid* grid, + SectorMagneticFieldMap::symmetry sym, + int polynomial_order, + int smoothing_order) { + // too lazy to write code to handle this case - not available to user anyway + if (sym != SectorMagneticFieldMap::dipole) { + throw(LogicalError( + "SectorMagneticFieldMap::IO::ReadMap", + "Failed to recognise symmetry type" + )); + } + std::vector< std::vector<double> > data(field_points.size(), + std::vector<double>(3)); + for (size_t i = 0; i < field_points.size(); ++i) { + data[i][0] = field_points[i][3]; + data[i][1] = field_points[i][4]; + data[i][2] = field_points[i][5]; + } + // symmetry is dipole + try { + INFOMSG("Calculating polynomials..." << endl); + PPSolveFactory solver(grid, data, polynomial_order, smoothing_order); + PolynomialPatch* patch = solver.solve(); + INFOMSG(" ... done" << endl); + return patch; + } catch (GeneralClassicException& exc) { + throw exc; + } +} + +VectorMap* SectorMagneticFieldMap::IO::getInterpolator (const std::vector< std::vector<double> > field_points, ThreeDGrid* grid, SectorMagneticFieldMap::symmetry sym) { - int y_start = 0; - if (sym == SectorMagneticFieldMap::dipole) { - y_start = (grid->ySize()+1)/2-1; - } // build field arrays double *** Bx, *** By, *** Bz; int index = 0; @@ -293,7 +381,7 @@ Interpolator3dGridTo3d* SectorMagneticFieldMap::IO::getInterpolator Bx[i][j] = new double[grid->zSize()]; By[i][j] = new double[grid->zSize()]; Bz[i][j] = new double[grid->zSize()]; - for (int k = 0; k < grid->zSize() && j >= y_start; ++k) { + for (int k = 0; k < grid->zSize(); ++k) { Bx[i][j][k] = field_points[index][3]; By[i][j][k] = field_points[index][4]; Bz[i][j][k] = field_points[index][5]; @@ -301,22 +389,8 @@ Interpolator3dGridTo3d* SectorMagneticFieldMap::IO::getInterpolator } } } - - // extend field array downwards if field is dipole - if (sym == SectorMagneticFieldMap::dipole) { - for (int i = 0; i < grid->xSize(); ++i) { - for (int j = 0; j < y_start; ++j) { - for (int k = 0; k < grid->zSize(); ++k) { - Bx[i][j][k] = -Bx[i][grid->ySize()-j-1][k]; - By[i][j][k] = By[i][grid->ySize()-j-1][k]; - Bz[i][j][k] = -Bz[i][grid->ySize()-j-1][k]; - ++index; - } - } - } - } - - return new Interpolator3dGridTo3d(grid, Bx, By, Bz); + Interpolator3dGridTo3d* interpolator = new Interpolator3dGridTo3d(grid, Bx, By, Bz); + return dynamic_cast<VectorMap*>(interpolator); } bool SectorMagneticFieldMap::IO::comparator(std::vector<double> field_item1, @@ -390,7 +464,7 @@ ThreeDGrid* SectorMagneticFieldMap::IO::generateGrid (const std::vector< std::vector<double> > field_points, SectorMagneticFieldMap::symmetry sym) { std::vector<double> r_grid(1, field_points[0][0]); - std::vector<double> y_grid(1, field_points[0][1]); + std::vector<double> y_grid(1, field_points[0][1]), y_grid_neg; std::vector<double> phi_grid(1, field_points[0][2]); for (size_t i = 0; i < field_points.size(); ++i) { if (floatGreaterEqual(field_points[i][0], r_grid.back())) { @@ -403,15 +477,8 @@ ThreeDGrid* SectorMagneticFieldMap::IO::generateGrid phi_grid.push_back(field_points[i][2]); } } - + // reflect about y if symmetry is dipole - int y_end = static_cast<int>(y_grid.size()); - if (sym == SectorMagneticFieldMap::dipole) { - for (int i = 0; i < y_end-1; ++i) { - y_grid.insert(y_grid.begin(), 2.*y_grid[i]-y_grid[2*i+1]); - } - y_end = y_grid.size(); - } INFOMSG("Grid size (r, y, phi) = (" << r_grid.size() << ", " << y_grid.size() << ", " << phi_grid.size() << ")" << endl); @@ -423,6 +490,7 @@ ThreeDGrid* SectorMagneticFieldMap::IO::generateGrid << ")" << endl); ThreeDGrid* grid = new ThreeDGrid(r_grid, y_grid, phi_grid); + grid->setConstantSpacing(true); return grid; } diff --git a/classic/5.0/src/Fields/SectorMagneticFieldMap.h b/classic/5.0/src/Fields/SectorMagneticFieldMap.h index 7bb59fff29d10c87a43a385a511469beacc23212..be8b4961b7305ba776883161163796227694aa4f 100644 --- a/classic/5.0/src/Fields/SectorMagneticFieldMap.h +++ b/classic/5.0/src/Fields/SectorMagneticFieldMap.h @@ -33,10 +33,12 @@ #include <map> #include <iostream> -#include "Fields/SectorMagneticFieldMap/SectorField.h" +#include "Fields/SectorField.h" -class Interpolator3dGridTo3d; -class ThreeDGrid; +namespace interpolation { + class VectorMap; + class ThreeDGrid; +} /** \class SectorMagneticFieldMap handles field map grids with sector geometry * @@ -77,9 +79,15 @@ class SectorMagneticFieldMap : public SectorField { * field map - either "Dipole" or "None" * \param length_units multiplier for lengths * \param field_units multiplier for fields + * \param polynomial_order order of polynomial fit to the mesh points + * \param smoothing_order order of smoothing (should be >= polynomial_order) */ - SectorMagneticFieldMap(std::string file_name, std::string symmetry, - double length_units, double field_units); + SectorMagneticFieldMap(std::string file_name, + std::string symmetry, + double length_units, + double field_units, + int polynomial_order, + int smoothing_order); /** Copy constructor - makes a deep copy of the field map */ SectorMagneticFieldMap(const SectorMagneticFieldMap& field); @@ -125,14 +133,14 @@ class SectorMagneticFieldMap : public SectorField { * * Note SectorMagneticFieldMap still owns this memory. */ - Interpolator3dGridTo3d* getInterpolator(); + interpolation::VectorMap* getInterpolator(); /** Set the interpolator * * \param interpolator set the field map interpolator. Note * SectorMagneticFieldMap now owns the memory pointed to by interpolator */ - void setInterpolator(Interpolator3dGridTo3d* interpolator); + void setInterpolator(interpolation::VectorMap* interpolator); /** Get the field map file name */ @@ -188,23 +196,28 @@ class SectorMagneticFieldMap : public SectorField { /** Get change in azimuthal angle between entrance and exit */ double getDeltaPhi() const; - void test_f(); - friend class FieldMap; private: enum symmetry {none, dipole}; + /** Reflect R_temp about y if below bbmin + * \returns true if symmetry transformation was applied + */ + bool applySymmetry(double* R_temp) const; + static symmetry StringToSymmetry(std::string name); static std::string SymmetryToString(symmetry sym); void Rotate(double* value, double angle); - Interpolator3dGridTo3d* interpolator_m; + interpolation::VectorMap* interpolator_m; symmetry symmetry_m; std::vector<double> units_m; const std::string filename_m; double phiOffset_m; + int poly_order_m; + int smoothing_order_m; static const double fractionalBBPhiTolerance_m; static std::map<std::string, SectorMagneticFieldMap*> _fields; @@ -232,9 +245,14 @@ class SectorMagneticFieldMap::IO { * \param units units of the file - should be 6-vector * \param sym symmetry of the file - either "none" (full field map) or * "dipole" (field map is reflected about y = 0) + * \param poly_order order of the polynomial fit + * \param smoothing_order order of the polynomial smoothing */ - static Interpolator3dGridTo3d* readMap(std::string file_name, - std::vector<double> units, SectorMagneticFieldMap::symmetry sym); + static interpolation::VectorMap* readMap(std::string file_name, + std::vector<double> units, + SectorMagneticFieldMap::symmetry sym, + int poly_order, + int smoothing_order); private: static const double floatTolerance_m; @@ -245,16 +263,23 @@ class SectorMagneticFieldMap::IO { (std::string file_name, std::vector<double> units); // generate a grid based on the input map file - static ThreeDGrid* generateGrid - (const std::vector< std::vector<double> > field_points, + static interpolation::ThreeDGrid* generateGrid( + const std::vector< std::vector<double> > field_points, SectorMagneticFieldMap::symmetry sym); // get the interpolator based on field points and grid information - static Interpolator3dGridTo3d* getInterpolator - (const std::vector< std::vector<double> > field_points, - ThreeDGrid* grid, + static interpolation::VectorMap* getInterpolator( + const std::vector< std::vector<double> > field_points, + interpolation::ThreeDGrid* grid, SectorMagneticFieldMap::symmetry sym); + static interpolation::VectorMap* getInterpolatorPolyPatch( + const std::vector< std::vector<double> > field_points, + interpolation::ThreeDGrid* grid, + SectorMagneticFieldMap::symmetry sym, + int poly_order, + int smoothing_order); + // Compare two floats are same within tolerance static bool floatGreaterEqual(double in1, double in2); diff --git a/src/unit_tests/classic_src/AbsBeamline/SBend3DTest.cpp b/src/unit_tests/classic_src/AbsBeamline/SBend3DTest.cpp index 764622a9c9e07164ef13d6aa77ed675b647670b4..f879de17a6c9e1ecaffbb57bb27ec66e6e7ebc13 100644 --- a/src/unit_tests/classic_src/AbsBeamline/SBend3DTest.cpp +++ b/src/unit_tests/classic_src/AbsBeamline/SBend3DTest.cpp @@ -35,74 +35,97 @@ #include "Structure/BoundaryGeometry.h" #include "AbsBeamline/SBend3D.h" -void writeFieldMap() { - std::ofstream out("SBend3D_map.field"); - out << " 422280 422280 422280 1\n" - << " 1 X [LENGU]\n" - << " 2 Y [LENGU]\n" - << " 3 Z [LENGU]\n" - << " 4 BX [FLUXU]\n" - << " 5 BY [FLUXU]\n" - << " 6 BZ [FLUXU]\n" - << " 0\n" - << " 194.014700000 0.00000000000 80.3635200000 0.682759323466E-07 -5.37524925776 0.282807068056E-07\n" - << " 210.000000000 0.00000000000 0.00000000000 0.00000000000 5038.98826666 0.00000000000 \n" - << " 194.014700000 0.00000000000 -80.3635200000 0.682759055339E-07 -5.37524929545 -0.282807663112E-07\n" - << " 194.014700000 0.250000000000 80.3635200000 -0.250876890041E-01 -5.37468265486 -0.105656867773E-01\n" - << " 210.000000000 0.250000000000 0.00000000000 81.6064485044 5043.86454663 -0.162646432266E-02\n" - << " 194.014700000 0.250000000000 -80.3635200000 -0.252106804794E-01 -5.37468268550 0.102687594900E-01\n" - << " 194.014700000 0.500000000000 80.3635200000 -0.501213811240E-01 -5.37048637738 -0.211087378232E-01\n" - << " 210.000000000 0.500000000000 0.00000000000 163.479461107 5051.16103059 -0.319620022674E-02\n" - << " 194.014700000 0.500000000000 -80.3635200000 -0.503671731133E-01 -5.37048640042 0.205153440729E-01\n" - << " 194.938580000 0.00000000000 80.7462000000 -0.109193207295E-06 -5.47695895319 -0.452292335858E-07\n" - << " 211.000000000 0.00000000000 0.00000000000 0.00000000000 5346.80707376 0.00000000000 \n" - << " 194.938580000 0.00000000000 -80.7462000000 -0.109194150249E-06 -5.47695923342 0.452297706086E-07\n" - << " 194.938580000 0.250000000000 80.7462000000 -0.228949123764E-01 -5.47615607268 -0.965935126817E-02\n" - << " 211.000000000 0.250000000000 0.00000000000 71.6525328925 5351.82535885 -0.156196844700E-02\n" - << " 194.938580000 0.250000000000 -80.7462000000 -0.230215180683E-01 -5.47615608020 0.935369738637E-02\n" - << " 194.938580000 0.500000000000 80.7462000000 -0.457451523318E-01 -5.47168044969 -0.193017688931E-01\n" - << " 211.000000000 0.500000000000 0.00000000000 143.398642339 5359.42974721 -0.306846139443E-02\n" - << " 194.938580000 0.500000000000 -80.7462000000 -0.459981937280E-01 -5.47168045716 0.186908719298E-01\n" - << " 195.862460000 0.00000000000 81.1288900000 0.00000000000 -5.56838923159 0.00000000000 \n" - << " 212.000000000 0.00000000000 0.00000000000 0.00000000000 5612.41675566 0.00000000000 \n" - << " 195.862460000 0.00000000000 -81.1288900000 0.00000000000 -5.56838923767 0.00000000000 \n" - << " 195.862460000 0.250000000000 81.1288900000 -0.209034621315E-01 -5.56743204962 -0.884197597858E-02\n" - << " 212.000000000 0.250000000000 0.00000000000 61.0257302103 5617.50584022 -0.152992556601E-02\n" - << " 195.862460000 0.250000000000 -81.1288900000 -0.210330844408E-01 -5.56743205033 0.852904004824E-02\n" - << " 195.862460000 0.500000000000 81.1288900000 -0.416276950059E-01 -5.56278428219 -0.176093837788E-01\n" - << " 212.000000000 0.500000000000 0.00000000000 121.985754491 5625.18212699 -0.300403876153E-02\n" - << " 195.862460000 0.500000000000 -81.1288900000 -0.418867765797E-01 -5.56278428318 0.169839055390E-01" - << std::endl; - - if (!out.good()) - throw std::runtime_error("could not write field map 'SBend3D_map.field'"); - - out.close(); -} +#include "Utilities/LogicalError.h" -void removeFieldMap() { - if (std::remove("SBend3D_map.field") != 0) - throw std::runtime_error("could not delete field map 'SBend3D_map.field'"); +class LoadFieldMap { + public: + LoadFieldMap(std::string name, int polynomial_order, int smoothing_order) + : fname_m(name), sbend3d_m(NULL) { + writeFieldMap(); + getFieldMap(polynomial_order, smoothing_order); + } -} + ~LoadFieldMap() { + delete sbend3d_m; + removeFieldMap(); + } -SBend3D* loadFieldMap() { - try { - SBend3D* sbend3D = new SBend3D("name"); - sbend3D->setLengthUnits(10.); // cm - sbend3D->setFieldUnits(1.e-4); // ? - sbend3D->setFieldMapFileName("SBend3D_map.field"); - // sbend3D->setFieldMapFileName("src/unit_tests/classic_src/AbsBeamline/SBend3D_map.field"); - return sbend3D; - } catch (std::exception& exc) { - return NULL; + void getFieldMap(int polynomial_order, int smoothing_order) { + try { + sbend3d_m = new SBend3D("name"); + sbend3d_m->setLengthUnits(10.); // cm + sbend3d_m->setFieldUnits(1.e-4); // ? + sbend3d_m->setLengthUnits(10.); // cm + sbend3d_m->setPolynomialOrder(polynomial_order); // ? + sbend3d_m->setSmoothingOrder(smoothing_order); // ? + sbend3d_m->setFieldMapFileName(fname_m); + } catch (std::exception& exc) { + std::cerr << exc.what() << std::endl; + sbend3d_m = NULL; + throw std::runtime_error("Failed to get field map"); + } } -} -TEST(SBend3DTest, SBend3DGeometryTest) { - writeFieldMap(); - SBend3D* field = loadFieldMap(); + void writeFieldMap() { + std::ofstream out(fname_m); + out << " 422280 422280 422280 1\n" + << " 1 X [LENGU]\n" + << " 2 Y [LENGU]\n" + << " 3 Z [LENGU]\n" + << " 4 BX [FLUXU]\n" + << " 5 BY [FLUXU]\n" + << " 6 BZ [FLUXU]\n" + << " 0\n" + << " 194.014700000 0.00000000000 80.3635200000 0.682759323466E-07 -5.37524925776 0.282807068056E-07\n" + << " 210.000000000 0.00000000000 0.00000000000 0.00000000000 5038.98826666 0.00000000000 \n" + << " 194.014700000 0.00000000000 -80.3635200000 0.682759055339E-07 -5.37524929545 -0.282807663112E-07\n" + << " 194.014700000 0.250000000000 80.3635200000 -0.250876890041E-01 -5.37468265486 -0.105656867773E-01\n" + << " 210.000000000 0.250000000000 0.00000000000 81.6064485044 5043.86454663 -0.162646432266E-02\n" + << " 194.014700000 0.250000000000 -80.3635200000 -0.252106804794E-01 -5.37468268550 0.102687594900E-01\n" + << " 194.014700000 0.500000000000 80.3635200000 -0.501213811240E-01 -5.37048637738 -0.211087378232E-01\n" + << " 210.000000000 0.500000000000 0.00000000000 163.479461107 5051.16103059 -0.319620022674E-02\n" + << " 194.014700000 0.500000000000 -80.3635200000 -0.503671731133E-01 -5.37048640042 0.205153440729E-01\n" + << " 194.938580000 0.00000000000 80.7462000000 -0.109193207295E-06 -5.47695895319 -0.452292335858E-07\n" + << " 211.000000000 0.00000000000 0.00000000000 0.00000000000 5346.80707376 0.00000000000 \n" + << " 194.938580000 0.00000000000 -80.7462000000 -0.109194150249E-06 -5.47695923342 0.452297706086E-07\n" + << " 194.938580000 0.250000000000 80.7462000000 -0.228949123764E-01 -5.47615607268 -0.965935126817E-02\n" + << " 211.000000000 0.250000000000 0.00000000000 71.6525328925 5351.82535885 -0.156196844700E-02\n" + << " 194.938580000 0.250000000000 -80.7462000000 -0.230215180683E-01 -5.47615608020 0.935369738637E-02\n" + << " 194.938580000 0.500000000000 80.7462000000 -0.457451523318E-01 -5.47168044969 -0.193017688931E-01\n" + << " 211.000000000 0.500000000000 0.00000000000 143.398642339 5359.42974721 -0.306846139443E-02\n" + << " 194.938580000 0.500000000000 -80.7462000000 -0.459981937280E-01 -5.47168045716 0.186908719298E-01\n" + << " 195.862460000 0.00000000000 81.1288900000 0.00000000000 -5.56838923159 0.00000000000 \n" + << " 212.000000000 0.00000000000 0.00000000000 0.00000000000 5612.41675566 0.00000000000 \n" + << " 195.862460000 0.00000000000 -81.1288900000 0.00000000000 -5.56838923767 0.00000000000 \n" + << " 195.862460000 0.250000000000 81.1288900000 -0.209034621315E-01 -5.56743204962 -0.884197597858E-02\n" + << " 212.000000000 0.250000000000 0.00000000000 61.0257302103 5617.50584022 -0.152992556601E-02\n" + << " 195.862460000 0.250000000000 -81.1288900000 -0.210330844408E-01 -5.56743205033 0.852904004824E-02\n" + << " 195.862460000 0.500000000000 81.1288900000 -0.416276950059E-01 -5.56278428219 -0.176093837788E-01\n" + << " 212.000000000 0.500000000000 0.00000000000 121.985754491 5625.18212699 -0.300403876153E-02\n" + << " 195.862460000 0.500000000000 -81.1288900000 -0.418867765797E-01 -5.56278428318 0.169839055390E-01" + << std::endl; + + if (!out.good()) + throw std::runtime_error("could not write field map 'SBend3D_map.field'"); + + out.close(); + } + + void removeFieldMap() { + if (std::remove(fname_m.c_str()) != 0) + std::cout << "Failed to delete field map " << fname_m + << " - it's probably okay" << std::endl; + } + + + SBend3D* sbend3d_m; + std::string fname_m; +}; + +TEST(SBend3DTest, SBend3DGeometryTest) { + LoadFieldMap fieldLoader("field1", 1, 1); + SBend3D* field = fieldLoader.sbend3d_m; if (field == NULL) return; // skip the test Vector_t B, E, centroid; @@ -133,13 +156,13 @@ TEST(SBend3DTest, SBend3DGeometryTest) { << std::endl; } } - removeFieldMap(); } void testField(double r, double y, double phi, double bx, double by, double bz, double tol) { - SBend3D* field = loadFieldMap(); + LoadFieldMap fieldLoader("field2", 1, 1); + SBend3D* field = fieldLoader.sbend3d_m; if (field == NULL) { std::cerr << "SKIPPING SBEND3D TEST - FAILED TO LOAD FIELD" << std::endl; return; // skip the test @@ -166,14 +189,43 @@ void testField(double r, double y, double phi, } TEST(SBend3DTest, SBend3DFieldTest) { - writeFieldMap(); - // 211.000000000 0.00000000000 0.00000000000 0.00000000000 0.00000000000 testField(0., 0., Physics::pi/8., 0., 5346.80707376*1e-4, 0., 1e-6); - // 211.000000000 0.250000000000 0.00000000000 71.6525328925 5351.82535885 -0.156196844700E-02 testField(0., 2.5, Physics::pi/8., 71.65253*1e-4, 5351.8253*1e-4, -0.1561E-02*1e-4, 1e-7); +} - removeFieldMap(); +TEST(SBend3DTest, SBend3DPolyPatchTest) { + // make the poly order > 1; this puts SBend3D in PolynomialPatch mode; + // check the field map loads correctly including e.g. dipole symmetry + try { + LoadFieldMap fieldLoader1("field3", 2, 2); + LoadFieldMap fieldLoader2("field3", 2, 2); + LoadFieldMap fieldLoader3("field4", 1, 1); + SBend3D* field1 = fieldLoader1.sbend3d_m; + SBend3D* field2 = fieldLoader2.sbend3d_m; + SBend3D* field3 = fieldLoader3.sbend3d_m; + const double dphi = Physics::pi/8; + double radius = 2110.; + for (double y = -5.; y < 5.1; y += 2.5) + for (double r = 2100.; r < 2120.1; r += 5.) + for (double phi =-dphi/2.; phi < dphi*3.1; phi += dphi/2.) { + Vector_t B, E, centroid; + Vector_t pos(r*cos(phi)-radius, y, r*sin(phi)); + std::cerr << "pos: " << pos << " r: " << r << " phi: " << phi/Physics::pi*180. << std::endl; + field1->apply(pos, centroid, 0, E, B); + std::cerr << " field1: " << B << std::endl; + field2->apply(pos, centroid, 0, E, B); + std::cerr << " field2: " << B << std::endl; + field3->apply(pos, centroid, 0, E, B); + std::cerr << " field3: " << B << std::endl; + } + } catch (LogicalError& err) { + std::cerr << err.what() << std::endl; + } catch (GeneralClassicException& err) { + std::cerr << err.what() << std::endl; + } catch (std::exception& err) { + std::cerr << err.what() << std::endl; + } } diff --git a/src/unit_tests/classic_src/Fields/Interpolation/PPSolveFactoryTest.cpp b/src/unit_tests/classic_src/Fields/Interpolation/PPSolveFactoryTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..11aad0de1ed6ba871b842ee47f499a97b937e975 --- /dev/null +++ b/src/unit_tests/classic_src/Fields/Interpolation/PPSolveFactoryTest.cpp @@ -0,0 +1,379 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include <sstream> + +#include "gtest/gtest.h" + +#define PolynomialPatchTest_MakePlots +#ifdef PolynomialPatchTest_MakePlots +#include "TCanvas.h" +#include "TH2D.h" +#include "TGraph.h" +#endif + +#include "Fields/Interpolation/SquarePolynomialVector.h" +#include "Fields/Interpolation/PolynomialPatch.h" +#include "Fields/Interpolation/PolynomialCoefficient.h" +#include "Fields/Interpolation/PPSolveFactory.h" + +using namespace interpolation; + +void test_points(int dim, int lower, int upper, std::vector< std::vector<int> > pts) { + int upper_size = 1; + int lower_size = 1; + for (int i = 0; i < dim; ++i) + upper_size *= upper+1; + for (int i = 0; i < dim; ++i) + lower_size *= lower+1; + if (lower < 0) + lower_size = 0; + if (upper < 0) + upper_size = 0; + // size should be difference in area of the squares + EXPECT_EQ(pts.size(), upper_size - lower_size); + for (size_t i = 0; i < pts.size(); ++i) { + // each pts element should have length dim + EXPECT_EQ(pts[i].size(), dim); + // each pts element should have indices with lower < size <= upper + bool in_bounds = true; + for (int j = 0; j < dim; ++j) { + in_bounds = in_bounds || (pts[i][j] > lower && pts[i][j] <= upper); + } + EXPECT_TRUE(in_bounds); + // each pts element should be unique + for (size_t j = 0; j < pts.size(); ++j) { + if (j == i) + continue; + bool equal = true; + for (int k = 0; k < dim; ++k) + equal &= pts[i][k] == pts[j][k]; + EXPECT_FALSE(equal); + } + } +} + +TEST(PPSolveFactoryTest, TestNearbyPointsSquares) { + for (int upper = 0; upper < 5; ++upper) { + for (int lower = 0; lower < upper; ++lower) { + for (int dim = 1; dim < 5; ++dim) { + std::vector< std::vector<int> > pts = + PPSolveFactory::getNearbyPointsSquares(dim, lower, upper); + test_points(dim, lower, upper, pts); + } + } + } +} + +class PPSolveFactoryTestFixture : public ::testing::Test { + protected: + std::vector<std::vector<double> > values; + SquarePolynomialVector ref; + ThreeDGrid* grid; + MMatrix<double> refCoeffs; + int np; + + PPSolveFactoryTestFixture() { + // we make a reference poly vector + std::vector<double> data(2*27); // 2^3 + for (size_t i = 0; i < data.size(); ++i) + data[i] = i/10.; + refCoeffs = MMatrix<double>(2, 27, &data[0]); + ref = SquarePolynomialVector(3, refCoeffs); + // we get some data + np = 6; + // choose the grid so that a midpoint falls at 0/0/0; we want to compare it + // later on... + grid = new ThreeDGrid(1., 2., 3., -2.5, -5.0, -7.5, np, np, np); + for (Mesh::Iterator it = grid->begin(); it < grid->end(); ++it) { + std::vector<double> value(2); + ref.F(&it.getPosition()[0], &value[0]); + values.push_back(value); + } + } +}; + +// check linear fit exactly reproduces data on grid points +TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialLinear) { + PPSolveFactory fac1(grid->clone(), + values, + 1, + 1); + PolynomialPatch* patch1 = fac1.solve(); + for (Mesh::Iterator it = grid->begin(); it < grid->end(); ++it) { + std::vector<double> value(2); + patch1->function(&it.getPosition()[0], &value[0]); + for (size_t i = 0; i < 2; ++i) + EXPECT_NEAR(values[it.toInteger()][i], value[i], 1e-6); + } + delete patch1; +} + +// check quadratic fit exactly reproduces data on and off grid points (data +// comes from a quadratic polynomial as source) +TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadratic) { + PPSolveFactory fac2(grid->clone(), + values, + 2, + 2); + PolynomialPatch* patch2 = fac2.solve(); + // first check that the polynomial at 0, 0, 0 is the same as ref + std::vector<double> zero(3, 0.); + SquarePolynomialVector* test = patch2->getPolynomialVector(&zero[0]); + EXPECT_EQ(patch2->getValueDimension(), 2); + MMatrix<double> testCoeffs = test->GetCoefficientsAsMatrix(); + std::cout << "Ref" << std::endl; + std::cout << refCoeffs << std::endl; + std::cout << "Test" << std::endl; + std::cout << testCoeffs << std::endl; + ASSERT_EQ(testCoeffs.num_row(), refCoeffs.num_row()); + ASSERT_EQ(testCoeffs.num_col(), refCoeffs.num_col()); + for (size_t i = 0; i < testCoeffs.num_row(); ++i) + for (size_t j = 0; j < testCoeffs.num_row(); ++j) { + EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs(i+1, j+1), 1e-6) + << "col " << i << " row " << j; + EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs(i+1, j+1), 1e-6) + << "col " << i << " row " << j; + } + + // now check that the values in each polynomial are correct + ThreeDGrid testGrid(1./4., 2./4., 3./4., -1., -2., -3., np*4-1, np*4-1, np*4-1); + for (Mesh::Iterator it = testGrid.begin(); it < testGrid.end(); ++it) { + std::vector<double> refValue(2); + std::vector<double> testValue(2); + ref.F(&it.getPosition()[0], &refValue[0]); + patch2->function(&it.getPosition()[0], &testValue[0]); + for (size_t i = 0; i < 2; ++i) + EXPECT_NEAR(refValue[i], testValue[i], 1e-6) << std::endl << it; + } +} + +// check smoothed quadratic fit exactly reproduces data on and off grid points +// except near to the boundary +TEST_F(PPSolveFactoryTestFixture, TestSolvePolynomialQuadraticSmoothed) { + ASSERT_TRUE(false); + PPSolveFactory fac2(grid->clone(), + values, + 1, + 2); + PolynomialPatch* patch2 = fac2.solve(); + // first check that the polynomial at 0, 0, 0 is the same as ref + std::vector<double> zero(3, 0.); + SquarePolynomialVector* test = patch2->getPolynomialVector(&zero[0]); + MMatrix<double> testCoeffs = test->GetCoefficientsAsMatrix(); + std::cout << "Ref" << std::endl; + std::cout << refCoeffs << std::endl; + std::cout << "Test" << std::endl; + std::cout << testCoeffs << std::endl; + ASSERT_EQ(testCoeffs.num_row(), refCoeffs.num_row()); + ASSERT_EQ(testCoeffs.num_col(), refCoeffs.num_col()); + for (size_t i = 0; i < testCoeffs.num_row(); ++i) + for (size_t j = 0; j < testCoeffs.num_row(); ++j) { + EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs(i+1, j+1), 1e-6) + << "col " << i << " row " << j; + EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs(i+1, j+1), 1e-6) + << "col " << i << " row " << j; + } + + // now check that the values in each polynomial are correct + ThreeDGrid testGrid(1./4., 2./4., 3./4., -1., -2., -3., np*4-1, np*4-1, np*4-1); + for (Mesh::Iterator it = testGrid.begin(); it < testGrid.end(); ++it) { + std::vector<double> refValue(2); + std::vector<double> testValue(2); + ref.F(&it.getPosition()[0], &refValue[0]); + patch2->function(&it.getPosition()[0], &testValue[0]); + for (size_t i = 0; i < 2; ++i) + EXPECT_NEAR(refValue[i], testValue[i], 1e-6) << std::endl << it; + } +} + + +std::vector<double> get_value(std::vector<double> x, int np) { + static const double twopi = asin(1.)*4; + std::vector<double> a_value(3); + a_value[0] = sin(twopi*x[0]/np)*sin(twopi*x[1]/np)*sin(twopi*x[2]/np); + a_value[1] = cos(twopi*x[0]/np)*cos(twopi*x[1]/np)*cos(twopi*x[2]/np); + a_value[2] = 1.; + return a_value; +} + +#ifdef PolynomialPatchTest_MakePlots +void plot(int n_points, std::vector<double> start, std::vector<double> end, PolynomialPatch* patch, int n_grid_points, std::string title) { + double delta_sq = 0; + for (size_t i = 0; i < start.size(); ++i) + delta_sq += (start[i]-end[i])*(start[i]-end[i]); + std::vector<TGraph*> graphs(3*patch->getValueDimension(), NULL); + for (size_t j = 0; j < graphs.size(); ++j) { + graphs[j] = new TGraph(n_points); + } + std::vector<double> errors(3, 0.); + for (int i = 0; i < n_points; ++i) { + std::vector<double> point = start; + for (size_t j = 0; j < point.size(); ++j) { + point[j] += (end[j]-start[j])*i/n_points; + } + std::vector<double> fit_value(patch->getValueDimension(), 0.); + patch->function(&point[0], &fit_value[0]); + for (size_t j = 0; j < patch->getValueDimension(); ++j) { + graphs[j]->SetPoint(i, point[0], fit_value[j]); + } + std::vector<double> true_value = get_value(point, n_grid_points); + for (size_t j = 0; j < patch->getValueDimension(); ++j) { + graphs[j+patch->getValueDimension()]->SetPoint(i, point[0], true_value[j]); + graphs[j+2*patch->getValueDimension()]->SetPoint(i, point[0], fabs(fit_value[j]-true_value[j])); + errors[j] += fabs(fit_value[j]-true_value[j]); + } + } + graphs[0]->SetTitle("Fit Bx"); + graphs[1]->SetTitle("Fit By"); + graphs[2]->SetTitle("Fit Bz"); + graphs[3]->SetTitle("True Bx"); + graphs[4]->SetTitle("True By"); + graphs[5]->SetTitle("True Bz"); + + TCanvas* c1 = new TCanvas("canvas", "canvas", 1024, 768); + c1->Draw(); + TH2D* h1 = new TH2D("h1", (title+";pos [AU]; Fit [AU]").c_str(), 10000, start[0]-(end[0]-start[0])/10., end[0]+(end[0]-start[0])/10., 1000, -5., 5.); + h1->SetStats(false); + h1->Draw(); + for (size_t j = 0; j < 2*patch->getValueDimension(); ++j) { + graphs[j]->SetFillColor(0); + graphs[j]->SetLineColor(j+1); + graphs[j]->SetMarkerColor(j+1); + graphs[j]->Draw("l"); + } + c1->BuildLegend(); + static int test_index = 0; + std::stringstream test_name; + test_name << "test_" << test_index; + c1->SetGridx(); + c1->Print((test_name.str()+".png").c_str()); + c1->Print((test_name.str()+".root").c_str()); + + TCanvas* c2 = new TCanvas("deltas canvas", "deltas canvas", 1024, 768); + c2->Draw(); + TH2D* h2 = new TH2D("h2", ("Residuals "+title+";pos [AU]; Fit [AU]").c_str(), 10000, start[0]-(end[0]-start[0])/10., end[0]+(end[0]-start[0])/10., 1000, 1e-6, 5.); + h2->SetStats(false); + h2->Draw(); + + std::vector<std::stringstream> strings(errors.size()); + for (size_t i = 0; i < errors.size(); ++i) { + strings[i] << " <#epsilon> " << errors[i]/n_points; + } + + graphs[6]->SetTitle(("Delta Bx "+strings[0].str()).c_str()); + graphs[7]->SetTitle(("Delta By "+strings[1].str()).c_str()); + graphs[8]->SetTitle(("Delta Bz "+strings[2].str()).c_str()); + for (size_t j = 2*patch->getValueDimension(); j < 3*patch->getValueDimension(); ++j) { + graphs[j]->SetFillColor(0); + graphs[j]->SetLineColor(j+1); + graphs[j]->SetMarkerColor(j+1); + graphs[j]->Draw("l"); + } + c2->BuildLegend(); + c2->SetGridx(); + c2->SetLogy(); + c2->Print(("residuals_"+test_name.str()+".png").c_str()); + c2->Print(("residuals_"+test_name.str()+".root").c_str()); + test_index++; +} +#else // PolynomialPatchTest_MakePlots +void plot(int n_points, std::vector<double> start, std::vector<double> end, PolynomialPatch* patch, int n_grid_points, std::string title) { +} +#endif // PolynomialPatchTest_MakePlots + + +TEST(PPSolveFactoryTest, TestThreeDSolveSinCos) { + int np = 11; + ThreeDGrid grid(1., 1., 1., 1., 2., 3., np, np, np); + // test grid starts a few points in to avoid boundary effects + ThreeDGrid fineGrid(1./4., 1./4., 1./4., 3., 4., 5., (np-2)*4, (np-2)*4, (np-2)*4); + std::vector<double> start = grid.begin().getPosition(); + std::vector<double> end = (grid.end()-1).getPosition(); + + std::vector<std::vector<double> > values; + for (Mesh::Iterator it = grid.begin(); it < grid.end(); ++it) { + values.push_back(get_value(it.getPosition(), np)); + } + std::vector<PolynomialPatch*> ppVec; + for (int smooth_order = 1; smooth_order < 4; ++smooth_order) { + for (int pp_order = smooth_order; pp_order <= smooth_order; ++pp_order) { + std::cout << "Building pp of order " << pp_order << " smooth " << smooth_order << std::endl; + PPSolveFactory fac(grid.clone(), + values, + pp_order, + smooth_order); + ppVec.push_back(fac.solve()); + std::stringstream ss; + ss << "smooth: " << smooth_order << " poly: " << pp_order; + plot(100, start, end, ppVec.back(), np, ss.str()); + EXPECT_EQ(ppVec.back()->getValueDimension(), 3); + } + } + std::cout << "Testing" << std::endl; + // check we get exactly right answer on mesh points + for (Mesh::Iterator it = grid.begin(); it < grid.end(); ++it) { + for (size_t ppi = 0; ppi < ppVec.size(); ppi++) { + std::vector<double> testValue(3, 0.); + ppVec[ppi]->function(&it.getPosition()[0], &testValue[0]); + for (size_t i = 0; i < 3; ++i) { + EXPECT_NEAR(get_value(it.getPosition(), np)[i], testValue[i], 1e-6); + } + } + } + // check we get nearly right answer between mesh points + std::vector<std::vector<double> > sumError(ppVec.size(), std::vector<double>(3, 0.)); + for (Mesh::Iterator it = fineGrid.begin(); it < fineGrid.end(); ++it) { + std::vector<double> prevTestValue(3, 0.); + std::vector<double> refValue = get_value(it.getPosition(), np); + ppVec[0]->function(&it.getPosition()[0], &prevTestValue[0]); + for (size_t ppi = 0; ppi < ppVec.size(); ppi++) { + std::vector<double> testValue(3, 0.); + ppVec[ppi]->function(&it.getPosition()[0], &testValue[0]); + for (size_t i = 0; i < 3; ++i) { + sumError[ppi][i] += fabs(testValue[i]-refValue[i]); + } + prevTestValue = testValue; + } + } + + for (size_t i = 0; i < sumError.size(); ++i) { + std::cout << i << ": "; + for (size_t j = 0; j < sumError[i].size(); ++j) { + std::cout << sumError[i][j] << " "; + // higher order interpolations should fit better than lower order + // except a little bit of floating point error can come in for e.g. + // Bz where we are fitting a straight line (pathological case)... + if (i > 0) { + EXPECT_LT(sumError[i][j], sumError[i-1][j]+1e-6); + } + } + std::cout << std::endl; + } +} + + diff --git a/src/unit_tests/classic_src/Fields/Interpolation/PolynomialPatchTest.cpp b/src/unit_tests/classic_src/Fields/Interpolation/PolynomialPatchTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..50a81c8e6a8c10612fbac300821bee2d6dc8ddf2 --- /dev/null +++ b/src/unit_tests/classic_src/Fields/Interpolation/PolynomialPatchTest.cpp @@ -0,0 +1,63 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "gtest/gtest.h" + +#include "Fields/Interpolation/PolynomialPatch.h" + +using namespace interpolation; + +TEST(PolynomialPatchTest, TestPolynomialPatch) { + ThreeDGrid grid(1, 2, 3, -1, -2, -3, 4, 3, 2); + // we make a reference poly vector + std::vector<double> data(54); + for (size_t i = 0; i < data.size(); ++i) + data[i] = i/10.; + MMatrix<double> refCoeffs(2, 27, &data[0]); + SquarePolynomialVector ref(3, refCoeffs); + // copy it into the grid + std::vector<SquarePolynomialVector*> poly; + for (size_t i = 0; i < grid.end().toInteger(); ++i) + poly.push_back(new SquarePolynomialVector(ref)); + PolynomialPatch patch(grid.clone(), grid.clone(), poly); + ThreeDGrid testGrid(1/4., 2/4., 3/4., -1, -2, -3, 4*4, 3*4, 2*4); + for (Mesh::Iterator it = testGrid.begin(); it < testGrid.begin()+1; ++it) { + std::vector<double> testValue(2); + patch.function(&it.getPosition()[0], &testValue[0]); + Mesh::Iterator nearest = grid.getNearest(&it.getPosition()[0]); + std::vector<double> localPosition(3); + for (size_t i = 0; i < 3; ++i) { + localPosition[i] = nearest.getPosition()[i] - it.getPosition()[i]; + } + std::vector<double> refValue(2); + ref.F(&localPosition[0], &refValue[0]); + for (size_t i = 0; i < 2; ++i) + EXPECT_NEAR(testValue[i], refValue[i], 1e-6); + } +} + + diff --git a/src/unit_tests/classic_src/Fields/Interpolation/SolveFactoryTest.cpp b/src/unit_tests/classic_src/Fields/Interpolation/SolveFactoryTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..2e42495ec422c3e7c0014dc05705a5051fbd9e00 --- /dev/null +++ b/src/unit_tests/classic_src/Fields/Interpolation/SolveFactoryTest.cpp @@ -0,0 +1,70 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "gtest/gtest.h" + +#include "Fields/Interpolation/SolveFactory.h" +#include "Fields/Interpolation/SquarePolynomialVector.h" +#include "Fields/Interpolation/MMatrix.h" + +using namespace interpolation; + +TEST(SolveFactoryTest, TestSolveNoDerivs) { + // we make a reference poly vector + std::vector<double> data(27); + for (size_t i = 0; i < data.size(); ++i) + data[i] = i; + MMatrix<double> refCoeffs(3, 9, &data[0]); + SquarePolynomialVector ref(2, refCoeffs); + // Make a set of points + std::vector< std::vector<double> > positions; + std::vector< std::vector<double> > values; + std::vector<double> pos(2); + std::vector<double> val(3); + for (pos[0] = 0.; pos[0] < 2.5; pos[0] += 1.) + for (pos[1] = 0.; pos[1] < 2.5; pos[1] += 1.) { + ref.F(&pos[0], &val[0]); + positions.push_back(pos); + values.push_back(val); + } + // now try to solve for the test poly vector; should be that the test and + // ref are identical + std::vector<std::vector<double> > deriv_pos; + std::vector< std::vector<int> > deriv_index; + SolveFactory fac(2, 2, 2, 3, positions, deriv_pos, deriv_index); + SquarePolynomialVector* vec = fac.PolynomialSolve(values, + deriv_pos); + MMatrix<double> testCoeffs = vec->GetCoefficientsAsMatrix(); + ASSERT_EQ(testCoeffs.num_row(), 3); + ASSERT_EQ(testCoeffs.num_col(), 9); + for (size_t i = 0; i < 3; ++i) + for (size_t j = 0; j < 3; ++j) + EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs(i+1, j+1), 1e-6); +} + + + diff --git a/src/unit_tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp b/src/unit_tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..74788a52270de2861d250408b015e0cb930af67f --- /dev/null +++ b/src/unit_tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp @@ -0,0 +1,106 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "gtest/gtest.h" + +#include "Fields/Interpolation/SquarePolynomialVector.h" +#include "Fields/Interpolation/PolynomialCoefficient.h" + +using namespace interpolation; + +TEST(SquarePolynomialVectorTest, TestConstructorDestructor) { + SquarePolynomialVector test1; + EXPECT_EQ(test1.PointDimension(), 0); + EXPECT_EQ(test1.ValueDimension(), 0); + EXPECT_EQ(test1.GetCoefficientsAsMatrix().num_col(), 0); + EXPECT_EQ(test1.GetCoefficientsAsMatrix().num_row(), 0); + + std::vector<double> data(18); + for (size_t i = 0; i < data.size(); ++i) + data[i] = i; + MMatrix<double> refCoeffs(2, 9, &data[0]); // 2x9 -> c, x, y, xy, xx, xxy, xxyy, xyy, yy, + SquarePolynomialVector test2(2, refCoeffs); + MMatrix<double> testCoeffs = test2.GetCoefficientsAsMatrix(); + ASSERT_EQ(testCoeffs.num_row(), 2); + ASSERT_EQ(testCoeffs.num_col(), 9); + for (size_t i = 0; i < testCoeffs.num_row(); ++i) { + for (size_t j = 0; j < testCoeffs.num_col(); ++j) { + EXPECT_EQ(testCoeffs(i+1, j+1), refCoeffs(i+1, j+1)); + } + } +} + +TEST(SquarePolynomialVectorTest, TestMakePolyVector) { + std::vector<double> data(18, 0.); + MMatrix<double> refCoeffs(2, 9, &data[0]); + SquarePolynomialVector ref(2, refCoeffs); + MVector<double> point(2, 0.); + MVector<double> polyVector(9, -99); + ref.MakePolyVector(point, polyVector); + EXPECT_EQ(polyVector(1), 1.); + for (size_t i = 1; i < 9; ++i) + EXPECT_EQ(polyVector(i+1), 0.); + point(1) = 3; + point(2) = 2; + // c x0 x1 x0x1 x0x0 x0x0x1 x1x1 x1x1x0, x1x1x0x0 + double refVector[] = {1., 3., 2., 6., 9., 18., 4., 12., 36.}; + ref.MakePolyVector(point, polyVector); + for (size_t i = 0; i < 9; ++i) + EXPECT_EQ(polyVector(i+1), refVector[i]); +} + +TEST(SquarePolynomialVectorTest, TestF) { + std::vector<double> data(18); + for (size_t i = 0; i < data.size(); ++i) + data[i] = i; + MMatrix<double> refCoeffs(2, 9, &data[0]); + SquarePolynomialVector ref(2, refCoeffs); + std::vector<double> point1(2, 0.); + std::vector<double> value(2); + ref.F(&point1[0], &value[0]); + EXPECT_EQ(value[0], refCoeffs(1, 1)); + EXPECT_EQ(value[1], refCoeffs(2, 1)); + MVector<double> point2(2); + point2(1) = 2.; + point2(2) = 3.; + MVector<double> polyVector(9, -99); + ref.MakePolyVector(point2, polyVector); + MVector<double> refValue = refCoeffs*polyVector; + MVector<double> testValue(2, -99); + ref.F(point2, testValue); + EXPECT_EQ(testValue(1), refValue(1)); // sum (0, ..., 8) + EXPECT_EQ(testValue(2), refValue(2)); // sum (9, ..., 17) + std::vector<double> point3(2); + point3[0] = 2.; + point3[1] = 3.; + ref.F(&point3[0], &value[0]); + EXPECT_EQ(value[0], refValue(1)); // sum (0, ..., 8) + EXPECT_EQ(value[1], refValue(2)); // sum (9, ..., 17) +} + + + diff --git a/src/unit_tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp b/src/unit_tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f7af14d28d107beb0197fe932f74a749d16146e3 --- /dev/null +++ b/src/unit_tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp @@ -0,0 +1,84 @@ +/* + * Copyright (c) 2015, Chris Rogers + * All rights reserved. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * 3. Neither the name of STFC nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "gtest/gtest.h" +#include "Fields/Interpolation/ThreeDGrid.h" + +using interpolation::ThreeDGrid; + +TEST(ThreeDGridTest, LowerBoundTest) { + std::vector<double> xVar(7), yVar(8), zVar(9); + for (size_t i = 0; i < xVar.size(); ++i) { + xVar[i] = 4.+1.*i; + } + for (size_t i = 0; i < yVar.size(); ++i) { + yVar[i] = -5.+2.*i; + } + for (size_t i = 0; i < zVar.size(); ++i) { + zVar[i] = 6.+3.*i*i; + } + zVar[5] = (zVar[4]+zVar[6])*0.51; // force not constant spacing + ThreeDGrid gridConst(1., 2., 3., // dx/dy/dz + 4., 5., 6., // min x/y/z + 7, 8, 9); // n x/y/z coordinates in the grid + ThreeDGrid gridVar(xVar, yVar, zVar); + gridVar.setConstantSpacing(false); + + ThreeDGrid* gridArray[] = {&gridConst, &gridVar}; + for (size_t i = 0; i < 2; ++i) { + ThreeDGrid* grid = gridArray[i]; + int index = -99; + for (size_t j = 0; j < grid->xSize(); ++j) { + grid->xLowerBound(grid->x(j+1)-1e-9, index); + EXPECT_EQ(index, int(j)-1) << "grid" << i << " " << j; + grid->xLowerBound(grid->x(j+1)+0e-9, index); + EXPECT_EQ(index, int(j)) << "grid" << i << " " << j; + grid->xLowerBound(grid->x(j+1)+1e-9, index); + EXPECT_EQ(index, int(j)) << "grid" << i << " " << j; + } + + index = -99; + for (size_t j = 0; j < grid->ySize(); ++j) { + grid->yLowerBound(grid->y(j+1)-1e-9, index); + EXPECT_EQ(index, int(j)-1) << "grid" << i << " " << j; + grid->yLowerBound(grid->y(j+1)+0e-9, index); + EXPECT_EQ(index, int(j)) << "grid" << i << " " << j; + grid->yLowerBound(grid->y(j+1)+1e-9, index); + EXPECT_EQ(index, int(j)) << "grid" << i << " " << j; + } + + index = -99; + for (size_t j = 0; j < grid->zSize(); ++j) { + grid->zLowerBound(grid->z(j+1)-1e-9, index); + EXPECT_EQ(index, int(j)-1) << "grid" << i << " " << j; + grid->zLowerBound(grid->z(j+1)+0e-9, index); + EXPECT_EQ(index, int(j)) << "grid" << i << " " << j; + grid->zLowerBound(grid->z(j+1)+1e-9, index); + EXPECT_EQ(index, int(j)) << "grid" << i << " " << j; + } + } +} diff --git a/src/unit_tests/classic_src/Utilities/RingSectionTest.cpp b/src/unit_tests/classic_src/Utilities/RingSectionTest.cpp index 5b5c6639057bbb610271347ad9e4aaea54ab17b8..fbc02587c1e3f37942f0013665e8e9631987814f 100644 --- a/src/unit_tests/classic_src/Utilities/RingSectionTest.cpp +++ b/src/unit_tests/classic_src/Utilities/RingSectionTest.cpp @@ -216,3 +216,10 @@ TEST(RingSectionTest, TestDoesOverlap) { EXPECT_FALSE(ors3.doesOverlap(f1, f1)); EXPECT_FALSE(ors3.doesOverlap(f4, f4)); } + +TEST(RingSectionTest, TestToDoFails) { + EXPECT_TRUE(false) << "SectorMagneticFieldMap: fixed bug that I rotated from Cartesian B field to Polar B field when I want Cartesian B field; moved Symmetry operation (reflection in Bphi/Bz) to field map lookup, it should be done at build time; added some print statements.\n" + << "ParallelCyclotronTracker: added the WRite3DFieldMap routines that should be removed; would it be possible to make a routine to write the field map?\n" + << "RingSection: changed the transformation routines in getFieldSTregnth; add unit test to follow\n" + << std::endl; +}