PPSolveFactory.cpp

Go to the documentation of this file.

/*
 *  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 <iostream>
#include <sstream>
#include <string>
#include <vector>
 
#include "Utilities/GSLCompat.h"
 
#include "Utility/Inform.h"  // ippl
 
#include "Utilities/GeneralOpalException.h"
 
#include "Fields/Interpolation/PPSolveFactory.h"
#include "Fields/Interpolation/PolynomialPatch.h"
#include "Fields/Interpolation/SolveFactory.h"
 
extern Inform* gmsg;
 
namespace interpolation {
 
    typedef PolynomialCoefficient Coeff;
 
    namespace {
        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;
        }
    }  // namespace
 
    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 == nullptr) {
            throw GeneralOpalException("PPSolveFactory::Solve", "nullptr 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 GeneralOpalException("PPSolveFactory::Solve", ss.str());
        }
        if (smoothingOrder < polyPatchOrder) {
            throw GeneralOpalException(
                    "PPSolveFactory::Solve",
                    "Polynomial order must be <= smoothing order in Solve");
        }
        polyMesh_m = points->dual();
        if (polyMesh_m == nullptr)
            throw GeneralOpalException("PPSolveFactory::Solve", "Dual of Mesh was nullptr");
        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> 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 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(const Mesh::Iterator& it) {
// calculate the derivatives near to polynomial at position indexed by it
#ifdef DEBUG
        bool verbose = false;
        if (verbose) {
            std::cerr << "PPSolveFactory::GetDerivs at position ";
            for (auto p : it.getPosition()) {
                std::cerr << p << " ";
            }
            std::cerr << std::endl;
        }
#endif
        // derivatives go in derivValues_m
        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 last = it.getMesh()->end() - 1;
        int deltaOrder      = smoothingOrder_m - polyPatchOrder_m;
        if (deltaOrder <= 0) {  // no smoothing needed
            return;
        }
        for (size_t i = 0; i < derivPoints_m.size(); ++i) {
            Mesh::Iterator nearest = it;
            bool inTheMesh         = true;
            for (int j = 0; j < posDim && inTheMesh; ++j) {
                nearest[j] += derivOrigins_m[i][j];
                inTheMesh = nearest[j] <= last[j];
            }
            if (!inTheMesh) {
                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);
                }
            }
#ifdef DEBUG
            if (verbose) {
                std::cerr << "Point ";
                for (auto p : derivPoints_m[i]) {
                    std::cerr << p << " ";
                }
                std::cerr << " values ";
                for (auto v : derivValues_m[i]) {
                    std::cerr << v << " ";
                }
                std::cerr << " derivIndices ";
                for (auto in : derivIndices_m[i]) {
                    std::cerr << in << " ";
                }
                std::cerr << std::endl;
            }
#endif
        }
    }
 
    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, nullptr);
        // get the list of points that are needed to make a given poly vector
        getPoints();
        getDerivPoints();
        SolveFactory solver(
                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.) {
                *gmsg << "    Done " << newPercentage << " %" << 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(GeneralOpalException(
                    "PPSolveFactory::getNearbyPointsSquares", "Point dimension must be > 0"));
        if (polyOrderLower > polyOrderUpper)
            throw(GeneralOpalException(
                    "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;
    }
}  // namespace interpolation