diff --git a/BSplineKokkos.cpp b/BSplineKokkos.cpp new file mode 100644 index 0000000..8ba3134 --- /dev/null +++ b/BSplineKokkos.cpp @@ -0,0 +1,53 @@ +#include "BSplineKokkos.h" +#include +/*template +BSplineKokkos::BSplineKokkos(int orderC, std::vector& ctrlPtsC, std::vector& knotsC, std::vector& weightsC) { + order = orderC; + //Allocate appropriate view space based on the number of control points, copy the data over to view + ctrlPts("ctrlPts", ctrlPtsC.size()); + for (int i = 0; i < ctrlPtsC.size(); i++) { + ctrlPts(i) = ctrlPtsC[i]; + } + //Same for knots and weights + knots("knots", knotsC.size()); + for (int i = 0; i < knotsC.size(); i++) { + knots(i) = knotsC[i]; + } + + weights("weights", weightsC.size()); + for (int i = 0; i < weightsC.size(); i++) { + weights(i) = weightsC[i]; + } + //Call the calculateDerivCoeff() to populate + //1st and 2nd deriavtive views + + //NOTE: UNCOMMENT THIS ONCE ALL IN FRONT ARE RESOLVED + //calculateDerivCoeff(); + +}*/ + + +/*template +void BSplineKokkos::calculateDerivCoeff() { + //Calculate first order derivative + //Allocate space for ctrlPts_1stD + ctrlPts_1stD("ctrlPts1Derivative", ctrlPts.extent(0)-1); + for (int i = 1; i < ctrlPts.extent(0); i++) { + double delta = double(order - 1)/(knots(i+order-1)-knots(i)); + ctrlPts_1stD(i) = ((ctrlPts(i) - ctrlPts(i-1)*delta)); + + } + + //Calculate second order derivative + //Allocate space for ctrlPts_2ndD + ctrlPts_2ndD("ctrlPts2Derivative", ctrlPts.extent(0)-2); + for (int i = 0; i < ctrlPts_1stD.extent(0); i++) { + double delta = double(order - 2) / (knots(i+order-1)-knots(i+1)); + ctrlPts_2ndD(i) = ((ctrlPts_1stD(i) - ctrlPts_1stD(i-1))*delta); + } + //TODO: find another way to verify the size of the second derivative view + +}*/ + +//Explicit instantiation of the templated class for Kokkos::serial +//template class BSplineKokkos; diff --git a/BSplineKokkos.h b/BSplineKokkos.h new file mode 100644 index 0000000..173c1a6 --- /dev/null +++ b/BSplineKokkos.h @@ -0,0 +1,432 @@ +#ifndef BSPLINEKOKKOS_H +#define BSPLINEKOKKOS_H + + +#include +#include +#include +#include + +//For the Kokkos version of the implementation, we will not be separating BSplineKokkos and BSplineKokkos2D, their respective functions will be in this file. + +template +class BSplineKokkos { +public: + using MemSpace = typename ExecutionSpace::memory_space; + + //Constructor for just 1 spline(Would this be used often?) + BSplineKokkos(int order_p, std::vector& ctrlPts_x, std::vector& ctrlPts_y, std::vector& knotsI) { + //In the original implementation, we had the x and y components separated and was storing redundant knots information + //The implementation now change to interleave the x and y coordinates + //We store just 1 set of knots rather than 2 + Kokkos::View orderV("Orders", 1); + auto host_orderV = Kokkos::create_mirror_view(orderV); + host_orderV(0) = order_p; + order = orderV; + Kokkos::deep_copy(order, host_orderV); + //Initialize control points and their offset + Kokkos::View ctrlPtsV("CtrlPoints", 2*ctrlPts_x.size()); + //Create the mirror on host so we can access the vector values + auto host_ctrlPtsV = Kokkos::create_mirror_view(ctrlPtsV); + + for (int i = 0; i < ctrlPts_x.size(); i++) { + host_ctrlPtsV(2*i) = ctrlPts_x[i]; + host_ctrlPtsV((2*i)+1) = ctrlPts_y[i]; + } + + ctrlPts = ctrlPtsV; + Kokkos::deep_copy(ctrlPts, host_ctrlPtsV); + Kokkos::View cPOffsetV("CtrlPointsOffset", 1); + auto host_cPOffsetV = Kokkos::create_mirror_view(cPOffsetV); + host_cPOffsetV(0) = 2*ctrlPts_x.size(); //End of ctrlPts + + cPOffset = cPOffsetV; + Kokkos::deep_copy(cPOffset, host_cPOffsetV); + + //Initialize knots and their offset + Kokkos::View knotsV("Knots", knotsI.size()); + Kokkos::View knotsOffsetV("KnotsOffset", 1); + auto host_knotsV = Kokkos::create_mirror_view(knotsV); + auto host_knotsOffsetV = Kokkos::create_mirror_view(knotsOffsetV); + host_knotsOffsetV(0) = knotsI.size(); //End of knots + + for (int i = 0; i < knotsI.size(); i++) { + host_knotsV(i) = knotsI[i]; + } + + knots = knotsV; + Kokkos::deep_copy(knots, host_knotsV); + + knotsOffset = knotsOffsetV; + Kokkos::deep_copy(knotsOffset, host_knotsOffsetV); + + //Derivative coefficients + calculateDerivCoeff(); + } + +//Create BSpline for a vector of BSplineKokkos objects + BSplineKokkos(std::vector& splines) { + //Preallocate space for Kokkos::View + //They are fixed in size and resize is expensive + int orderSize = 0; + int ctrlPSize = 0; + int knotsSize = 0; + int cPOffSize = 0; + int kOffSize = 0; + + //Input the offsets + //Track space needed to allocate for ctrlPts & knots + for (int i = 0; i < splines.size(); i++) { + cPOffSize += splines[i].getCPOffsetSize(); + kOffSize += splines[i].getKnotsOffsetSize(); + ctrlPSize += splines[i].getCtrlPtsSize(); + knotsSize += splines[i].getKnotsSize(); + orderSize += splines[i].getOrder().extent(0); + } + //Pre-Allocate space for views + Kokkos::View orderV ("Order", orderSize); + Kokkos::View ctrlPtsV ("CtrlPts", ctrlPSize); + Kokkos::View cPOffsetV ("CtrlPtsOffset", cPOffSize); + Kokkos::View knotsV("Knots", knotsSize); + Kokkos::View knotsOffsetV ("KnotsOffset", kOffSize); + auto host_orderV = Kokkos::create_mirror_view(orderV); + auto host_ctrlPtsV = Kokkos::create_mirror_view(ctrlPtsV); + auto host_cPOffsetV = Kokkos::create_mirror_view(cPOffsetV); + auto host_knotsV = Kokkos::create_mirror_view(knotsV); + auto host_knotsOffsetV = Kokkos::create_mirror_view(knotsOffsetV); + + int oidx = 0; + int cPidx = 0; + int cOidx = 0; + int kidx = 0; + int kOidx = 0; + + //Populate the views with data + int cpLast = 0; + int kLast = 0; + + //NOTE TO SELF: INDEXING MAY BE WRONG + for (int i = 0; i < splines.size(); i++) { + Kokkos::View intView = splines[i].getOrder(); + auto host_intView = Kokkos::create_mirror_view(intView); + Kokkos::deep_copy(host_intView, intView); + //Populate order + for (int j = 0; j < intView.extent(0); j++) { + host_orderV(oidx+j)= host_intView(j); + } + oidx += host_intView.extent(0); + //Populate ctrlpts offset + //Offset the offset by the last offset in the view + intView = splines[i].getCPOffset(); + host_intView = Kokkos::create_mirror_view(intView); + Kokkos::deep_copy(host_intView, intView); + for (int j = 0; j < host_intView.extent(0); j++) { + host_cPOffsetV(cOidx+j) = host_intView(j)+cpLast; + } + cOidx += host_intView.extent(0); + cpLast += host_intView(host_intView.extent(0)-1); + //Populate knots offset + //Offset the offset by the last offset in the view + + intView = splines[i].getKnotsOffset(); + host_intView = Kokkos::create_mirror_view(intView); + Kokkos::deep_copy(host_intView, intView); + for (int j = 0; j < host_intView.extent(0); j++) { + host_knotsOffsetV(kOidx+j) = host_intView(j)+kLast; + } + kOidx += host_intView.extent(0); + kLast += host_intView(host_intView.extent(0)-1); + + //Populate ctrlpts + Kokkos::View doubleView = splines[i].getCtrlPts(); + auto host_doubleView = Kokkos::create_mirror_view(doubleView); + Kokkos::deep_copy(host_doubleView, doubleView); + for (int j = 0; j < host_doubleView.extent(0); j++) { + host_ctrlPtsV(cPidx+j) = host_doubleView(j); + } + cPidx += host_doubleView.extent(0); + //Populate knots + doubleView = splines[i].getKnots(); + host_doubleView = Kokkos::create_mirror_view(doubleView); + Kokkos::deep_copy(host_doubleView, doubleView); + for (int j = 0; j < host_doubleView.extent(0); j++) { + host_knotsV(kidx+j) = host_doubleView(j); + } + kidx += host_doubleView.extent(0); + } + + //Update the member variables + order = orderV; + knots = knotsV; + ctrlPts = ctrlPtsV; + cPOffset = cPOffsetV; + knotsOffset = knotsOffsetV; + //Copy the modification made from the host data to the device view + Kokkos::deep_copy(order, host_orderV); + Kokkos::deep_copy(knots, host_knotsV); + Kokkos::deep_copy(ctrlPts, host_ctrlPtsV); + Kokkos::deep_copy(cPOffset, host_cPOffsetV); + Kokkos::deep_copy(knotsOffset, host_knotsOffsetV); + + //TO DO: Work on calculate derivative coeff + calculateDerivCoeff(); + + } + + //We perform De Boor recursion here + KOKKOS_FUNCTION void evalDeBoor (double x, int splineo, int lKnot, Kokkos::View result, const Kokkos::View order, const Kokkos::View knots, const Kokkos::View ctrlPts_1stD) const { + + int order_t = lKnot; + int leftPtX = 0; + int leftPtY = 1; + + while (x > knots(lKnot+1)) { + lKnot++; + leftPtX+=2; + leftPtY+=2; + } + + double ptsX[3]; + double ptsY[3]; + + int idx = 0; + for (int i = leftPtX; i < leftPtX+2*(order_t); i+=2) { + ptsX[idx] = ctrlPts_1stD(i); + ptsY[idx] = ctrlPts_1stD(i+1); + idx++; + } + + auto localKnots = Kokkos::subview(knots, Kokkos::pair(lKnot-order_t+2, lKnot+order_t)); + + //I am still working on the offset for the knots + for (int r = 1; r <= order_t; r++) { + for (int i = order_t-1; i >= r; i--) { + double alpha; + if (localKnots(i-1) == localKnots(i+order_t-r-1)) { + alpha = 0.; + } + else { + alpha = (x - localKnots(i-1))/(knots(i+order_t-r-1) - knots(i-1)); + } + ptsX[i] = (1. - alpha) * ptsX[i-1] + alpha * ptsX[i]; + ptsY[i] = (1. - alpha) * ptsY[i-1] + alpha * ptsY[i]; + } + } + + result(0) = ptsX[order_t-1]; + result(1) = ptsY[order_t-1]; + } + + std::vector eval1stDeriv(std::vector xVals, int splineo) const { + int lKnot; + Kokkos::deep_copy(lKnot, Kokkos::subview(order, splineo)); + lKnot--; + Kokkos::View res("result", 2); + Kokkos::parallel_for("parallel call for evalDeBoors", xVals.size(), KOKKOS_CLASS_LAMBDA(int i) { + evalDeBoor(xVals[i], splineo, lKnot, res, order, knots, ctrlPts_1stD); + }); + auto mvRes = Kokkos::create_mirror_view(res); + Kokkos::deep_copy(mvRes, res); + return {mvRes(0), mvRes(1)}; + } + + double eval2ndDeriv(double x, int splineo) const { + auto mv_order = Kokkos::create_mirror_view(order); + Kokkos::deep_copy(mv_order, order); + if (mv_order(splineo) == 2) { + return 0; + } + //Get the left knot, leftPtX, leftPtY + int lKnot = mv_order(splineo)-1; + int leftPtX = 0; + int leftPtY = 1; + + //Make sure that we did not go over to the next spline + auto mv_knotsOffset = Kokkos::create_mirror_view(knotsOffset); + Kokkos::deep_copy(mv_knotsOffset, knotsOffset); + + int bound = mv_knotsOffset(splineo); + + auto mv_knots = Kokkos::create_mirror_view(knots); + Kokkos::deep_copy(mv_knots, knots); + while (mv_knots(lKnot+1) < x) { + lKnot+=2; + leftPtX+=2; + leftPtY+=2; + if (lKnot == bound-1) { + break; + } + } + + //Populate pts and local knot views + int order_t = mv_order(splineo)-2; + auto mv_ctrlPts_2ndD = Kokkos::create_mirror_view(ctrlPts_2ndD); + int idx = 0; + Kokkos::deep_copy(mv_ctrlPts_2ndD, ctrlPts_2ndD); + + for (int i = 0; i < mv_ctrlPts_2ndD.extent(0); i++) { + std::cout << mv_ctrlPts_2ndD(i) << std::endl; + } + + Kokkos::View ptsX("ptsX", order_t); + Kokkos::View ptsY("ptsY", order_t); + auto mv_ptsX = Kokkos::create_mirror_view(ptsX); + auto mv_ptsY = Kokkos::create_mirror_view(ptsY); + + for (int i = leftPtX; i < leftPtX+order_t; i+=2) { + mv_ptsX(idx) = mv_ctrlPts_2ndD(i); + mv_ptsY(idx) = mv_ctrlPts_2ndD(i+1); + } + + //Start populating the local knots + idx = 0; + Kokkos::View localKnots("localKnots", 2*order_t-2); + auto mv_localKnots = Kokkos::create_mirror_view(localKnots); + for (int i = lKnot-order_t+2; i < lKnot+order_t; i++) { + mv_localKnots(idx) = mv_knots(i); + idx++; + } + + //Copy the data back to device + Kokkos::deep_copy(localKnots, mv_localKnots); + Kokkos::deep_copy(ptsX, mv_ptsX); + Kokkos::deep_copy(ptsY, mv_ptsY); + + Kokkos::parallel_for ("2nd derivative loop", order_t, KOKKOS_LAMBDA(int r){ + for (int i = order_t-1; i >= r+1; i--) { + double aLeft = localKnots(i-1); + double aRight = localKnots(i+order_t-(r+1)-1); + double alpha; + if (aLeft == aRight) { + alpha = 0; + } + else { + alpha = (x-aLeft)/(aRight-aLeft); + } + + ptsX(i) = (1. - alpha) * ptsX(i-1)+alpha*ptsX(i); + } + }); + + Kokkos::deep_copy(mv_ptsX, ptsX); + return mv_ptsX(order_t-1); + + } + +//Accessors + Kokkos::View getOrder() const {return order;} + Kokkos::View getCtrlPts() const {return ctrlPts;} + Kokkos::View getKnots() const {return knots;} + Kokkos::View getCPOffset() const {return cPOffset;} + Kokkos::View getKnotsOffset() const {return knotsOffset;} + Kokkos::View get1stD() const {return ctrlPts_1stD;} + Kokkos::View get2ndD() const {return ctrlPts_2ndD;} + + int getCPOffsetSize() const {return cPOffset.extent(0);} + int getKnotsOffsetSize() const {return knotsOffset.extent(0);} + int getCtrlPtsSize() const {return ctrlPts.extent(0);} + int getKnotsSize() const {return knots.extent(0);} + int getOrderSize() const {return order.extent(0);} + + double getCtrlPtCoor(int BspIdx, int cPIdx, char coor) const { + if (coor == 'x') { + //Find the correct offset + return ctrlPts(cPOffset(2*BspIdx)+cPIdx); + } + else { + return ctrlPts(cPOffset(2*BspIdx+1)+cPIdx); + } + } + double getKnotCoor(int BspIdx, int kIdx, char coor) const { + if (coor == 'x') { + return knots(knotsOffset(2*BspIdx)+kIdx); + } + else { + return knots(knotsOffset(2*BspIdx+1)+kIdx); + } + } + + void calculateDerivCoeff() { + //Moved here from the .cpp file since it uses + //Calculate first order derivative + //Allocate the space + + Kokkos::View ctrlPts_1stDV("ctrlPts1Derivative", ctrlPts.extent(0)-(2*cPOffset.extent(0))); + Kokkos::View cP1stDOffsetV("ctrlPts1DerivativeOffset", cPOffset.extent(0)); + Kokkos::View cP2ndDOffsetV("ctrlPts2DerivativeOffset", cPOffset.extent(0)); + + auto host_cP2ndDOffsetV = Kokkos::create_mirror_view(cP2ndDOffsetV); + auto host_cP1stDOffsetV = Kokkos::create_mirror_view(cP1stDOffsetV); + + //Adjust the offset to include be 1 & 2 less than the ctrlPtsOffset + auto host_cPOffset = Kokkos::create_mirror_view(cPOffset); + Kokkos::deep_copy(host_cPOffset, cPOffset); + if (host_cPOffset.extent(0) > 1) { + Kokkos::deep_copy(host_cPOffset, cPOffset); + for (int i = 1; i < cPOffset.extent(0); i++) { + host_cP1stDOffsetV(i) = host_cPOffset(i)-2; + host_cP2ndDOffsetV(i) = host_cPOffset(i)-4; + } + } + else { + host_cP1stDOffsetV(0) = host_cPOffset(0)-2; + host_cP2ndDOffsetV(0) = host_cPOffset(0)-4; + } + //We need to partition the x and y while we calculate the coefficient + int idx = 0; + int oidx = 0; + double deltaHold = 0; + int kidx = 2; + + auto host_knots = Kokkos::create_mirror_view(knots); + auto host_order = Kokkos::create_mirror_view(order); + auto host_ctrlPts = Kokkos::create_mirror_view(ctrlPts); + auto host_ctrlPts_1stDV = Kokkos::create_mirror_view(ctrlPts_1stDV); + + //Copy these value over + Kokkos::deep_copy(host_knots, knots); + Kokkos::deep_copy(host_order, order); + Kokkos::deep_copy(host_ctrlPts, ctrlPts); + + //We need to keep track of the indicies of the splines + //We ignore the first ctrl point, starting at idx = 2 + for (int i = 2; i < host_ctrlPts.extent(0); i++) { + //Even indices stores all the x coordinates + //Odd indices stores all the y coordinates + if (i == host_cPOffset(idx)) { + oidx++; + idx++; + } + if (i%2==0) { + double delta = double (host_order(oidx)-1)/(host_knots(kidx+host_order(oidx)-2)-host_knots(kidx-1)); + deltaHold = delta; + kidx++; + } + host_ctrlPts_1stDV(i-2) = (host_ctrlPts(i)-host_ctrlPts(i-2))*deltaHold; + } + + Kokkos::deep_copy(cP1stDOffsetV, host_cP1stDOffsetV); + Kokkos::deep_copy(ctrlPts_1stDV, host_ctrlPts_1stDV); + + ctrlPts_1stD = ctrlPts_1stDV; + cP1stDOffset = cP1stDOffsetV; + } + +private: + + Kokkos::View order; + + //CtrlPts, knots, weights and their offsets + Kokkos::View ctrlPts; + Kokkos::View cPOffset; + Kokkos::View knots; + Kokkos::View knotsOffset; + + //The 1st and 2nd derivatives + Kokkos::View ctrlPts_1stD; + Kokkos::View cP1stDOffset; + Kokkos::View ctrlPts_2ndD; + Kokkos::View cP2ndDOffset; + +}; +#endif diff --git a/BSplineKokkos2D.cpp b/BSplineKokkos2D.cpp new file mode 100644 index 0000000..d26ee24 --- /dev/null +++ b/BSplineKokkos2D.cpp @@ -0,0 +1 @@ +#include "BSplineKokkos2D.h" diff --git a/BSplineKokkos2D.h b/BSplineKokkos2D.h new file mode 100644 index 0000000..b3e4983 --- /dev/null +++ b/BSplineKokkos2D.h @@ -0,0 +1,372 @@ +#ifndef BSPLINEKOKKOS2D_H + +#define BSPLINEKOKKOS2D_H + +#include +#include +#include +#include + +template + +class BSplineKokkos2D { +public: + using MemSpace = typename ExecutionSpace::memory_space; + + BSplineKokkos2D(int order_p, std::vector& ctrlPts_x, std::vector& ctrlPts_y, std::vector& knotsI) { + Kokkos::View orderV("Orders", 1); + auto host_orderV = Kokkos::create_mirror_view(orderV); + host_orderV(0) = order_p; + order = orderV; + Kokkos::deep_copy(order, host_orderV); + + Kokkos::View ctrlPtsV("ctrlPts", ctrlPts_x.size()); + auto host_ctrlPtsV = Kokkos::create_mirror_view(ctrlPtsV); + for (int i = 0; i < ctrlPts_x.size(); i++) { + host_ctrlPtsV(i, 0) = ctrlPts_x[i]; + host_ctrlPtsV(i, 1) = ctrlPts_y[i]; + } + ctrlPts = ctrlPtsV; + Kokkos::deep_copy(ctrlPts, host_ctrlPtsV); + + Kokkos::View cpOffsetV("cpOffset", 1); + auto host_cpOffsetV = Kokkos::create_mirror_view(cpOffsetV); + host_cpOffsetV(0) = ctrlPts_x.size(); + host_cpOffsetV(0) = ctrlPts_x.size(); + + cpOffset = cpOffsetV; + Kokkos::deep_copy(cpOffset, host_cpOffsetV); + + Kokkos::View knotsV("knots", knotsI.size()); + auto host_knotsV = Kokkos::create_mirror_view(knotsV); + for (int i = 0; i < knotsI.size(); i++) { + host_knotsV(i) = knotsI[i]; + } + knots = knotsV; + Kokkos::deep_copy(knots, host_knotsV); + + Kokkos::View knotsOffsetV("knotsOffset", 1); + auto host_knotsOffsetV = Kokkos::create_mirror_view(knotsOffsetV); + host_knotsOffsetV(0) = knotsI.size(); + knotsOffset = knotsOffsetV; + Kokkos::deep_copy(knotsOffset, host_knotsOffsetV); + + //TO DO: Implement 1st&2nd Coef function + calculateDerivCoeff(); + } + + BSplineKokkos2D(std::vector& multiSplines) { + //Loop over all the splines to know how much space to allocate + int orderSize = 0; + int ctrlPtsSize = 0; + int cpOffsetSize = 0; + int knotSize = 0; + int knotsOffsetSize = 0; + + for (int i = 0; i < multiSplines.size(); i++) { + orderSize += multiSplines[i].getOrder().extent(0); + ctrlPtsSize += multiSplines[i].getCtrlPts().extent(0); + cpOffsetSize += multiSplines[i].getCPOffset().extent(0); + knotSize += multiSplines[i].getKnots().extent(0); + knotsOffsetSize += multiSplines[i].getKnotsOffset().extent(0); + } + + //Populate the views + Kokkos::View orderV ("order", orderSize); + Kokkos::View ctrlPtsV("ctrlPts", ctrlPtsSize); + Kokkos::View cPOffsetV("ctrlPtsOffset", cpOffsetSize); + Kokkos::View knotsV("knots", knotSize); + Kokkos::View knotsOffsetV("knotsOffset", knotsOffsetSize); + + //Create the mirror view on host to move the data over + auto mvOrderV = Kokkos::create_mirror_view(orderV); + auto mvCtrlPtsV = Kokkos::create_mirror_view(ctrlPtsV); + auto mvCPOffsetV = Kokkos::create_mirror_view(cPOffsetV); + auto mvKnotsV = Kokkos::create_mirror_view(knotsV); + auto mvKnotsOffsetV = Kokkos::create_mirror_view(knotsOffsetV); + //Copy the data over + int oidx = 0; + int cPidx = 0; + int cOidx = 0; + int kidx = 0; + int kOidx = 0; + for (int i = 0; i < multiSplines.size(); i++) { + //Populate orders + Kokkos::View intView = multiSplines[i].getOrder(); + auto mvIntView = Kokkos::create_mirror_view(intView); + Kokkos::deep_copy(mvIntView, intView); + for (int j = 0; j < intView.extent(0); j++) { + mvOrderV(oidx+j) = mvIntView(j); + } + oidx += mvIntView.extent(0); + + //Populate ctrlPts + Kokkos::View double2DView = multiSplines[i].getCtrlPts(); + auto mvDouble2DView = Kokkos::create_mirror_view(double2DView); + Kokkos::deep_copy(mvDouble2DView, double2DView); + for (int j = 0; j < mvDouble2DView.extent(0); j++) { + mvCtrlPtsV(cPidx+j, 0) = mvDouble2DView(j, 0); + mvCtrlPtsV(cPidx+j, 1) = mvDouble2DView(j, 1); + } + cPidx += mvDouble2DView.extent(0); + + //Populate ctrlPtsOffset + intView = multiSplines[i].getCPOffset(); + mvIntView = Kokkos::create_mirror_view(intView); + Kokkos::deep_copy(mvIntView, intView); + + for (int j = 0; j < intView.extent(0); j++) { + mvCPOffsetV(cOidx + j) = intView(j); + } + cOidx += intView.extent(0); + + //Populate knots + Kokkos::View doubleView = multiSplines[i].getKnots(); + auto mvDoubleView = Kokkos::create_mirror_view(doubleView); + Kokkos::deep_copy(mvDoubleView, doubleView); + + for (int j = 0; j < doubleView.extent(0); j++) { + mvKnotsV(kidx+j) = doubleView(j); + } + kidx += doubleView.extent(0); + + //Populate knots offset + intView = multiSplines[i].getKnotsOffset(); + mvIntView = Kokkos::create_mirror_view(intView); + Kokkos::deep_copy(mvIntView, intView); + + for (int j = 0; j < intView.extent(0); j++) { + mvKnotsOffsetV(kOidx+j) = mvIntView(j); + } + kOidx += mvIntView.extent(0); + } + //Copy the data on host to device + order = orderV; + ctrlPts = ctrlPtsV; + cpOffset = cPOffsetV; + knots = knotsV; + knotsOffset = knotsOffsetV; + + Kokkos::deep_copy(order, mvOrderV); + Kokkos::deep_copy(ctrlPts, mvCtrlPtsV); + Kokkos::deep_copy(cpOffset, mvCPOffsetV); + Kokkos::deep_copy(knots, mvKnotsV); + Kokkos::deep_copy(knotsOffset, mvKnotsOffsetV); + calculateDerivCoeff(); + } + + void calculateDerivCoeff() { + //Allocate the views we need + Kokkos::View ctrlPts1stDV("ctrlPts1stDeriv", ctrlPts.extent(0)-cpOffset.extent(0)); + Kokkos::View cP1stDOffsetV("cP1stDOffset", cpOffset.extent(0)); + Kokkos::View cP2ndDOffsetV("cP2ndDOffset", cpOffset.extent(0)); + + //Set up the offset views + auto mvCP1stDOffsetV = Kokkos::create_mirror_view(cP1stDOffsetV); + auto mvCP2ndDOffsetV = Kokkos::create_mirror_view(cP2ndDOffsetV); + + auto mvCPOffset = Kokkos::create_mirror_view(cpOffset); + + if (mvCPOffset.extent(0) > 1) { + Kokkos::deep_copy(mvCPOffset, cpOffset); + for (int i = 1; i < mvCPOffset.extent(0); i++) { + mvCP1stDOffsetV(i) = mvCPOffset(i)-1; + mvCP2ndDOffsetV(i) = mvCPOffset(i)-2; + } + } + else { + mvCP1stDOffsetV(0) = mvCPOffset(0)-1; + mvCP2ndDOffsetV(0) = mvCPOffset(0)-2; + } + + auto mvOrder = Kokkos::create_mirror_view(order); + auto mvKnots = Kokkos::create_mirror_view(knots); + auto mvCtrlPts = Kokkos::create_mirror_view(ctrlPts); + auto mvCtrlPts1stDV = Kokkos::create_mirror_view(ctrlPts1stDV); + Kokkos::deep_copy(mvOrder, order); + Kokkos::deep_copy(mvKnots, knots); + Kokkos::deep_copy(mvCtrlPts, ctrlPts); + + //Calculate 1st derivative coef + int offidx = 0; //Offset index + int oidx = 0; //Order index + for (int i = 1; i < mvCtrlPts1stDV.extent(0)+1; i++) { + if (i == mvCP1stDOffsetV(offidx)+1) { + oidx++; + offidx++; + continue; + } + //We need to check whether we are on the border for the next spline in our structure + double delta = double(mvOrder(oidx)-1) / (mvKnots(i+mvOrder(oidx)-1) - mvKnots(i)); + + mvCtrlPts1stDV(i-1, 0) = (mvCtrlPts(i, 0) - mvCtrlPts(i-1, 0))*delta; + mvCtrlPts1stDV(i-1, 1) = (mvCtrlPts(i, 1) - mvCtrlPts(i-1, 1))*delta; + } + + //Calculate 2nd derivative coef + Kokkos::View ctrlPts2ndDV("ctrlPts2ndDeriv", ctrlPts.extent(0)-(2*cpOffset.extent(0))); + auto mvCtrlPts2ndDV = Kokkos::create_mirror_view(ctrlPts2ndDV); + + offidx = 0; + oidx = 0; + for (int i = 1; i < ctrlPts1stDV.extent(0); i++) { + if (i == mvCP2ndDOffsetV(offidx)+2) { + oidx++; + offidx++; + continue; + } + double delta = double(mvOrder(oidx)-2) / (mvKnots(i+mvOrder(oidx)-1) - mvKnots(i+1)); + + mvCtrlPts2ndDV(i-1, 0) = (mvCtrlPts1stDV(i, 0) - mvCtrlPts1stDV(i-1, 0)) * delta; + mvCtrlPts2ndDV(i-1, 1) = (mvCtrlPts1stDV(i, 1) - mvCtrlPts1stDV(i-1, 1)) * delta; + } + Kokkos::deep_copy(cP1stDOffsetV, mvCP1stDOffsetV); + Kokkos::deep_copy(cP2ndDOffsetV, mvCP2ndDOffsetV); + Kokkos::deep_copy(ctrlPts1stDV, mvCtrlPts1stDV); + Kokkos::deep_copy(ctrlPts2ndDV, mvCtrlPts2ndDV); + + cp1stDOffset = cP1stDOffsetV; + cp2ndDOffset = cP2ndDOffsetV; + ctrlPts1stD = ctrlPts1stDV; + ctrlPts2ndD = ctrlPts2ndDV; + } + + KOKKOS_FUNCTION void evalDeBoor(double x, int splineo, int lKnot, Kokkos::View result, const Kokkos::View order, const Kokkos::View knots, const Kokkos::View ctrlPts_1stD) const { + //DeBoor's algorithm for BSpline 1st deriv calculation + int order_t = lKnot; + int leftPt = 0; + + while (x > knots(lKnot+1)) { + lKnot++; + leftPt++; + } + + //Allocate temporary points + double ptsX[3]; + double ptsY[3]; + + int idx = 0; + for (int i = leftPt; i < leftPt + order_t; i++) { + ptsX[idx] = ctrlPts_1stD(i, 0); + ptsY[idx] = ctrlPts_1stD(i, 1); + idx++; + } + + auto localKnots = Kokkos::subview(knots, Kokkos::pair(lKnot-order_t+2, lKnot+order_t)); + + //Calculation loop + for (int r = 1; r <= order_t; r++) { + for(int i = order_t-1; i >= r; i--) { + double alpha; + if (localKnots(i-1) - localKnots(i+order_t-r-1) > 1e-12) { + alpha = 0; + } + else { + alpha = (x - localKnots(i-1)) / (knots(i+order_t-r-1) - knots(i-1)); + } + ptsX[i] = (1. - alpha) * ptsX[i-1]+alpha * ptsX[i]; + ptsY[i] = (1. - alpha) * ptsY[i-1]+alpha * ptsY[i]; + } + + result(0) = ptsX[order_t-1]; + result(1) = ptsY[order_t-1]; + } + } + + Kokkos::View eval1stDeriv(Kokkos::View xVals, int splineo) const { + int lKnot; + Kokkos::deep_copy(lKnot, Kokkos::subview(order, splineo)); + lKnot--; + Kokkos::View res("result", 2); + Kokkos::parallel_for("parallel evalDeBoors", xVals.size(), KOKKOS_CLASS_LAMBDA(int i){ + evalDeBoor(xVals(i), splineo, lKnot, res, order, knots, ctrlPts1stD); + }); + return res; + } + + //Second derivative calculation + KOKKOS_FUNCTION void evalDeBoor2(double x, int splineo, int lKnot, Kokkos::View result, const Kokkos::View order, const Kokkos::View knots, const Kokkos::View ctrlPts_2ndD) const { + if (order(splineo) == 2) { + result(0) = 0; + result(1) = 0; + return; + } + + int order_t = lKnot-1; + int leftPt = 0; + + while (knots(lKnot+1) < x) { + lKnot++; + leftPt++; + } + + double ptsX[3]; + double ptsY[3]; + + int idx = 0; + for (int i = leftPt; i < leftPt+order_t; i++) { + ptsX[idx] = ctrlPts_2ndD(i, 0); + ptsY[idx] = ctrlPts_2ndD(i, 1); + idx++; + } + + auto localKnots = Kokkos::subview(knots, Kokkos::pair(lKnot-order_t+2, lKnot+order_t)); + + for (int r = 1; r <= order_t; r++) { + for (int i = order_t-1; i >= r; i--) { + double aLeft = localKnots(i-1); + double aRight = localKnots(i+order_t-r-1); + double alpha; + if (aLeft == aRight) { + alpha = 0.; + } + else { + alpha = (x - aLeft) / (aRight - aLeft); + } + ptsX[i] = (1. - alpha) * ptsX[i-1] + alpha * ptsX[i]; + ptsY[i] = (1. - alpha) * ptsY[i-1] + alpha * ptsY[i]; + } + } + + result(0) = ptsX[order_t-1]; + result(1) = ptsY[order_t-1]; + + } + + Kokkos::View eval2ndDeriv(Kokkos::View xVals, int splineo) const { + int lKnot; + Kokkos::deep_copy(lKnot, Kokkos::subview(order, splineo)); + lKnot--; + Kokkos::View res("result", 2); + Kokkos::parallel_for("parallel evalDeBoors 2", xVals.size(), KOKKOS_CLASS_LAMBDA(int i) { + evalDeBoor2(xVals(i), splineo, lKnot, res, order, knots, ctrlPts2ndD); + }); + return res; + } + + //Accessors + Kokkos::View getOrder() const {return order;} + Kokkos::View getCtrlPts() const {return ctrlPts;} + Kokkos::View getCPOffset() const {return cpOffset;} + Kokkos::View getKnots() const {return knots;} + Kokkos::View getKnotsOffset() const {return knotsOffset;} + Kokkos::View getCP1stD() const {return ctrlPts1stD;} + Kokkos::View getCP1stDOffset() const {return cp1stDOffset;} + Kokkos::View getCP2ndD() const {return ctrlPts2ndD;} + Kokkos::View getCP2ndDOffset() const {return cp2ndDOffset;} + + +private: + Kokkos::View order; + Kokkos::View ctrlPts; + Kokkos::View cpOffset; + Kokkos::View knots; + Kokkos::View knotsOffset; + Kokkos::View ctrlPts1stD; + Kokkos::View ctrlPts2ndD; + Kokkos::View cp1stDOffset; + Kokkos::View cp2ndDOffset; + +}; + +#endif diff --git a/CMakeLists.txt b/CMakeLists.txt index 1109d31..0479639 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,10 @@ project(simLandIceMeshGen VERSION 0.0.1 LANGUAGES CXX) enable_testing() list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake") +#Enabling Address Sanitizer (ASan) +#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fsanitize=address") +#set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -fsanitize=address") + option(USE_SIMMODSUITE "Enable SimModSuite support" ON) if(USE_SIMMODSUITE) @@ -26,6 +30,14 @@ target_include_directories(modelGen2d PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) target_include_directories(modelGen2d PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/Quadtree/include) target_link_libraries(modelGen2d PRIVATE splineInterp Omega_h::omega_h) +#NEED TO ADD THE FILES USED FOR BSpline_Kokkos HERE +add_library(splineKokkos BSplineKokkos.cpp) +target_link_libraries(splineKokkos PRIVATE Kokkos::kokkos) + +add_library(splineKokkos2D BSplineKokkos2D.cpp) +target_link_libraries(splineKokkos2D PRIVATE Kokkos::kokkos) + + set(testDataDir ${CMAKE_CURRENT_SOURCE_DIR}/contour_test_data) if(USE_SIMMODSUITE) @@ -169,3 +181,55 @@ add_test(NAME testDiscoverTopo_L11pts set_tests_properties(testDiscoverTopo_L11pts PROPERTIES PASS_REGULAR_EXPRESSION "number of model vertices: 3") +#THE FOLLOWING IS FOR BSPLINE_KOKKOS TESTS, MAY HAVE A LOT OF ERROR + +add_executable(testSplineKokkosBase testSplineKokkosBase.cpp) +target_link_libraries(testSplineKokkosBase PRIVATE splineInterp curveReader modelGen2d splineKokkos Kokkos::kokkos) +add_test(NAME testBSplineKokkosSerialStraightEdge + COMMAND ./testSplineKokkosBase ${testDataDir}/testCurveStrightLine.csv 9) +add_test(NAME testBSplineKokkosSerial_edge_2pts + COMMAND ./testSplineKokkosBase ${testDataDir}/testCurve2pts.csv 1) + +#TEST FOR BSPLINEKOKKOS 1ST DERIVATIVE +add_executable(testKokkos1stDerivative testKokkos1stDerivative.cpp) +target_link_libraries(testKokkos1stDerivative PRIVATE splineInterp curveReader modelGen2d splineKokkos Kokkos::kokkos) + +add_test(NAME testKokkos1stDerivativeStraightEdge + COMMAND ./testKokkos1stDerivative contour_test_data/testCurveStraightLine.csv 9) + +add_test(NAME testKokkos1stDerivativeEdge2Pts + COMMAND ./testKokkos1stDerivative contour_test_data/testCurve2pts.csv 1) + +add_test(NAME testKokkos1stDerivativeSingleEdge + COMMAND ./testKokkos1stDerivative ${testDataDir}/testCurveSplineFitting.csv 12.2685) + +#TEST FOR BSPLINEKOKKOS2D CONSTRUCTOR +add_executable(testSplineKokkos2DBase testSplineKokkos2DBase.cpp) +target_link_libraries(testSplineKokkos2DBase PRIVATE splineInterp curveReader modelGen2d splineKokkos2D Kokkos::kokkos) +add_test(NAME testSplineKokkos2DStraightEdge + COMMAND ./testSplineKokkos2DBase contour_test_data/testCurveStraightLine.csv 9) +add_test(NAME testSplineKokkos2D2Pts + COMMAND ./testSplineKokkos2DBase contour_test_data/testCurve2pts.csv 1) +add_test(NAME testSplineKokkos2DSingleEdge + COMMAND ./testSplineKokkos2DBase ${testDataDir}/testCurveSplineFitting.csv 12.2685) + +#TEST FOR BSPLINEKOKKOS2D 1ST DERIVATIVE +add_executable(testKokkos2D1stDeriv testKokkos2D1stDeriv.cpp) +target_link_libraries(testKokkos2D1stDeriv PRIVATE splineInterp curveReader modelGen2d splineKokkos2D Kokkos::kokkos) +add_test(NAME testDeriv2DStraightEdge + COMMAND ./testKokkos2D1stDeriv contour_test_data/testCurveStraightLine.csv 9) +add_test(NAME testDeriv2D2Pts + COMMAND ./testKokkos2D1stDeriv contour_test_data/testCurve2pts.csv 1) +add_test(NAME testDeriv2DSingleEdge + COMMAND ./testKokkos2D1stDeriv ${testDataDir}/testCurveSplineFitting.csv 12.2685) + +#TEST FOR BSPLINEKOKKOS2D 2ND DERIVATIVE +add_executable(testKokkos2D2ndDeriv testKokkos2D2ndDeriv.cpp) +target_link_libraries(testKokkos2D2ndDeriv PRIVATE splineInterp curveReader modelGen2d splineKokkos2D Kokkos::kokkos) +add_test(NAME test2NDDeriv2DStraightEdge + COMMAND ./testKokkos2D2ndDeriv contour_test_data/testCurveStraightLine.csv 9) +add_test(NAME test2NDDeriv2D2Pts + COMMAND ./testKokkos2D2ndDeriv contour_test_data/testCurve2pts.csv 1) +add_test(NAME test2NDDerive2DSingleEdge + COMMAND ./testKokkos2D2ndDeriv contour_test_data/testSplineFitting.csv 12.2685) + diff --git a/testKokkos1stDerivative.cpp b/testKokkos1stDerivative.cpp new file mode 100644 index 0000000..119014d --- /dev/null +++ b/testKokkos1stDerivative.cpp @@ -0,0 +1,134 @@ +#include +#include "BSplineKokkos.h" +#include "BSpline.h" +#include + +//These imports are from testSplineFitting.cc +#include "splineInterpolation.h" +#include "curveReader.h" +#include +#include +#include +#include +#include +#include + + +template +void printDView(Kokkos::View& dView) { + auto host_dView = Kokkos::create_mirror_view(dView); + Kokkos::deep_copy(host_dView, dView); + for (int i = 0; i < host_dView.extent(0); i++) { + std::cout << host_dView(i) << std::endl; + } +} + +template +void printIView(Kokkos::View& iView) { + auto host_iView = Kokkos::create_mirror_view(iView); + Kokkos::deep_copy(host_iView, iView); + for (int i = 0; i < host_iView.extent(0); i++) { + std::cout << host_iView(i) << std::endl; + } +} + +template +void printSpline(BSplineKokkos spline) { + std::cout << "Spline content" << std::endl; + Kokkos::View dView = spline.getCtrlPts(); + std::cout << "CtrlPts" << std::endl; + printDView(dView); + Kokkos::View iView = spline.getCPOffset(); + std::cout << "CtrlPts Offset" << std::endl; + printIView(iView); + dView = spline.getKnots(); + std::cout << "Knots" << std::endl; + printDView(dView); + iView = spline.getKnotsOffset(); + std::cout << "Knots Offset" << std::endl; + printIView(iView); + iView = spline.getOrder(); + std::cout << "Order" << std::endl; + printIView(iView); +} + +void printVector(const std::string& name, std::vector& vector) { + std::cout << name << ": " << vector.size() << std::endl; + for (int i = 0; i < vector.size(); i++) { + std::cout << vector[i] << std::endl; + } +} + +double EPSILON = 1e-12; + +int main(int argc, char* argv[]) { + int retVal = 0; + if (argc != 3) { + //Check for the arguments needed + std::cerr<< "Input arguments: " << std::endl; + std::cerr << "input csv need these columns: "; + std::cerr << "coordinate x, coordinate y, coordinate z,isOnCurve,angle,isMdlVtx" << std::endl; + return 1; + } + Kokkos::initialize(argc, argv); + { + //Check for what memory space we should use + //Sticking to just cuda space or host space for now + #ifdef KOKKOS_ENABLE_CUDA + #define MemSpace Kokkos::CudaSpace + #endif + #ifndef MemSpace + #define MemSpace Kokkos::HostSpace + #endif + using ExecutionSpace = MemSpace::execution_space; + + std::string inputCSV = argv[1]; + int extensionPos = inputCSV.rfind("."); + int slashPos = inputCSV.rfind("/"); + std::string fileNameNoExt = inputCSV.substr(slashPos+1, extensionPos); + double expectedCurveLength = std::stod(argv[2]); + auto curve = CurveReader::readCurveInfo(inputCSV); + + //Construct BSpline2d object + SplineInterp::BSpline2d serialBSP; + if (curve.x.size() == 2) { + serialBSP = SplineInterp::attach_piecewise_linear_curve(curve.x, curve.y); + } else { + serialBSP = SplineInterp::fitCubicSplineToPoints(curve.x, curve.y); + } + + std::vector ctrlPtsX, ctrlPtsY, knots, weight; + int order; + + //Get the info from serial spline, we will feed this to kokkos spline + serialBSP.x.getpara(order, ctrlPtsX, knots, weight); + serialBSP.y.getpara(order, ctrlPtsY, knots, weight); + + BSplineKokkos kokkosBSP(order, ctrlPtsX, ctrlPtsY, knots); + //printSpline(kokkosBSP); + std::vector evalAt = {0, 0.2, 0.41, 0.5, 0.66, 0.73, 0.75, 0.89, 0.94, 1}; + + for (int i = 0; i < 10; i++) { + double derivX = serialBSP.x.evalFirstDeriv(evalAt[i]); + double derivY = serialBSP.y.evalFirstDeriv(evalAt[i]); + std::vector kokkos1stDeriv = kokkosBSP.eval1stDeriv({evalAt[i]}, 0); + double kokkosDerivX = kokkos1stDeriv[0]; + double kokkosDerivY = kokkos1stDeriv[1]; + + double derivXDiff = std::fabs(derivX) - std::fabs(kokkosDerivX); + double derivYDiff = std::fabs(derivY) - std::fabs(kokkosDerivY); + + if (derivXDiff > EPSILON || derivYDiff > EPSILON) { + std::cout << "Test " << i+1 << " failed, eval at: " << evalAt[i] << std::endl; + std::cout << "EPSILON = " << EPSILON << std::endl; + std::cout << "SERIAL/KOKKOS DIFFERENCE: xcoor = " << derivXDiff << " ycoor = " << derivYDiff << std::endl; + std::cout << "SERIAL 1st Derivative: x = " << derivX << " y = " << derivY << std::endl; + std::cout << "KOKKOS 1st Derivative: x = " << kokkosDerivX << " y = " << kokkosDerivY << std::endl; + retVal = 1; + } + } + + } + Kokkos::finalize(); + return retVal; +} diff --git a/testKokkos2D1stDeriv.cpp b/testKokkos2D1stDeriv.cpp new file mode 100644 index 0000000..3fc17d5 --- /dev/null +++ b/testKokkos2D1stDeriv.cpp @@ -0,0 +1,88 @@ +#include +#include "BSplineKokkos2D.h" +#include "BSpline.h" +#include + +#include "splineInterpolation.h" +#include "curveReader.h" +#include +#include +#include +#include +#include +#include + +double EPSILON = 1e-12; + +int main(int argc, char* argv[]) { + int retVal; + if (argc != 3) { + std::cerr<< "Input arguments: " << std::endl; + std::cerr << "input csv need these columns: "; + std::cerr << "coordinate x, coordinate y, coordinate z,isOnCurve,angle,isMdlVtx" << std::endl; + return 1; + } + + Kokkos::initialize(argc, argv); + { + #ifdef KOKKOS_ENABLE_CUDA + #define MemSpace Kokkos::CudaSpace + #endif + #ifndef MemSpace + #define MemSpace Kokkos::HostSpace + #endif + + using ExecutionSpace = MemSpace::execution_space; + + std::string inputCSV = argv[1]; + int extensionPos = inputCSV.rfind("."); + int slashPos = inputCSV.rfind("/"); + std::string fileNameNoExt = inputCSV.substr(slashPos+1, extensionPos); + double expectedCurveLength = std::stod(argv[2]); + auto curve = CurveReader::readCurveInfo(inputCSV); + + //Construct BSpline2d object + SplineInterp::BSpline2d serialBSP; + if (curve.x.size() == 2) { + serialBSP = SplineInterp::attach_piecewise_linear_curve(curve.x, curve.y); + } else { + serialBSP = SplineInterp::fitCubicSplineToPoints(curve.x, curve.y); + } + + std::vector ctrlPtsX, ctrlPtsY, knots, weight; + int order; + + //Get the info from serial spline, we will feed this to kokkos spline + serialBSP.x.getpara(order, ctrlPtsX, knots, weight); + serialBSP.y.getpara(order, ctrlPtsY, knots, weight); + + BSplineKokkos2D kokkosBSP(order, ctrlPtsX, ctrlPtsY, knots); + + std::vector evalAt = {0, 0.2, 0.41, 0.5, 0.66, 0.73, 0.75, 0.89, 0.94, 1}; + for (int i = 0; i < 10; i++) { + double derivX = serialBSP.x.evalFirstDeriv(evalAt[i]); + double derivY = serialBSP.y.evalFirstDeriv(evalAt[i]); + + Kokkos::View res("deriv result", 2); + Kokkos::View xVals ("paraCoor", 1); + auto mvXVals = Kokkos::create_mirror_view(xVals); + mvXVals(0) = evalAt[i]; + Kokkos::deep_copy(xVals, mvXVals); + res = kokkosBSP.eval1stDeriv(xVals, 0); + auto mvRes = Kokkos::create_mirror_view(res); + Kokkos::deep_copy(mvRes, res); + + double xDiff = std::fabs(derivX) - std::fabs(mvRes(0)); + double yDiff = std::fabs(derivY) - std::fabs(mvRes(1)); + if (xDiff > EPSILON || yDiff > EPSILON) { + std::cout << "Test " << i+1 << " failed, eval at: " << evalAt[i] << std::endl; + std::cout << "Difference: x = " << xDiff << " y = " << yDiff << std::endl; + std::cout << "SERIAL 1st deriv: x = " << derivX << " y = " << derivY << std::endl; + std::cout << "KOKKOS 1st deriv: x = " << mvRes(0) << " y = " << mvRes(1) << std::endl; + retVal = 1; + } + } + } + Kokkos::finalize(); + return retVal; +} diff --git a/testKokkos2D2ndDeriv.cpp b/testKokkos2D2ndDeriv.cpp new file mode 100644 index 0000000..3ce2fa3 --- /dev/null +++ b/testKokkos2D2ndDeriv.cpp @@ -0,0 +1,87 @@ +#include +#include "BSplineKokkos2D.h" +#include "BSpline.h" +#include + +#include "splineInterpolation.h" +#include "curveReader.h" +#include +#include +#include +#include +#include +#include + +double EPSILON = 1e-12; + +int main(int argc, char* argv[]) { + int retVal = 0; + if (argc != 3) { + std::cerr<< "Input arguments: " << std::endl; + std::cerr << "input csv need these columns: "; + std::cerr << "coordinate x, coordinate y, coordinate z,isOnCurve,angle,isMdlVtx" << std::endl; + return 1; + } + + Kokkos::initialize(argc, argv); + { + #ifdef KOKKOS_ENABLE_CUDA + #define MemSpace Kokkos::CudaSpace + #endif + #ifndef MemSpace + #define MemSpace Kokkos::HostSpace + #endif + + using ExecutionSpace = MemSpace::execution_space; + + std::string inputCSV = argv[1]; + int extensionPos = inputCSV.rfind("."); + int slashPos = inputCSV.rfind("/"); + std::string fileNameNoExt = inputCSV.substr(slashPos+1, extensionPos); + double expectedCurveLength = std::stod(argv[2]); + auto curve = CurveReader::readCurveInfo(inputCSV); + + SplineInterp::BSpline2d serialBSP; + if (curve.x.size() == 2) { + serialBSP = SplineInterp::attach_piecewise_linear_curve(curve.x, curve.y); + } + else { + serialBSP = SplineInterp::fitCubicSplineToPoints(curve.x, curve.y); + } + std::vector ctrlPtsX, ctrlPtsY, knots, weight; + int order; + serialBSP.x.getpara(order, ctrlPtsX, knots, weight); + serialBSP.y.getpara(order, ctrlPtsY, knots, weight); + + BSplineKokkos2D kokkosBSP(order, ctrlPtsX, ctrlPtsY, knots); + + std::vector evalAt = {0, 0.2, 0.41, 0.5, 0.66, 0.73, 0.75, 0.89, 0.94, 1}; + for (int i = 0; i < 10; i++) { + double derivX = serialBSP.x.evalSecondDeriv(evalAt[i]); + double derivY = serialBSP.y.evalSecondDeriv(evalAt[i]); + + Kokkos::View res("Result", 2); + Kokkos::View xVals("paraCoor", 1); + auto mvXVals = Kokkos::create_mirror_view(xVals); + mvXVals(0) = evalAt[i]; + Kokkos::deep_copy(xVals, mvXVals); + + res = kokkosBSP.eval2ndDeriv(xVals, 0); + auto mvRes = Kokkos::create_mirror_view(res); + Kokkos::deep_copy(mvRes, res); + + double xDiff = std::fabs(derivX - mvRes(0)); + double yDiff = std::fabs(derivY - mvRes(1)); + + if (xDiff > EPSILON || yDiff > EPSILON) { + std::cout << "Test " << i+1 << " failed, eval at: " << evalAt[i] << std::endl; + std::cout << "Difference: x = " << xDiff << " y = " << derivY << std::endl; + std::cout << "SERIAL 2nd deriv: x = " << derivX << ", " << derivY << std::endl; + std::cout << "KOKKOS 2nd deriv: x = " << mvRes(0) << ", " << mvRes(1) << std::endl; + retVal = 1; + } + } + } + Kokkos::finalize(); + return retVal; +} diff --git a/testSplineKokkos2DBase.cpp b/testSplineKokkos2DBase.cpp new file mode 100644 index 0000000..11f5b97 --- /dev/null +++ b/testSplineKokkos2DBase.cpp @@ -0,0 +1,149 @@ +//Constructor test for new representation of Kokkos spline +//We will be comparing this against the serial version +#include +#include "BSplineKokkos2D.h" +#include "BSpline.h" +#include "splineInterpolation.h" +#include "curveReader.h" + +#include +#include +#include +#include +#include +#include +#include + +//For checking if the content of the splines are correct +// + +double EPSILON = 1e-12; + +int main(int argc, char* argv[]) { + //We check how many arguments are given + int retVal = 0; + if (argc != 3) { + std::cerr<< "Input arguments: " << std::endl; + std::cerr << "input csv need these columns: "; + std::cerr << "coordinate x, coordinate y, coordinate z,isOnCurve,angle,isMdlVtx" << std::endl; + return 1; + } + Kokkos::initialize(argc, argv); + { + #ifdef KOKKOS_ENABLE_CUDA + #define MemSpace Kokkos::CudaSpace + #endif + #ifndef MemSpace + #define MemSpace Kokkos::HostSpace + #endif + + using ExecutionSpace = MemSpace::execution_space; + + std::string inputCSV = argv[1]; + int extensionPos = inputCSV.rfind("."); + int slashPos = inputCSV.rfind("/"); + std::string fileNameNoExt = inputCSV.substr(slashPos+1, extensionPos); + double expectedCurveLength = std::stod(argv[2]); + auto curve = CurveReader::readCurveInfo(inputCSV); + + //Construct BSpline2d object + SplineInterp::BSpline2d serialBSP; + if (curve.x.size() == 2) { + serialBSP = SplineInterp::attach_piecewise_linear_curve(curve.x, curve.y); + } else { + serialBSP = SplineInterp::fitCubicSplineToPoints(curve.x, curve.y); + } + + std::vector ctrlPtsX, ctrlPtsY, knots, weight; + int order; + + serialBSP.x.getpara(order, ctrlPtsX, knots, weight); + serialBSP.y.getpara(order, ctrlPtsY, knots, weight); + + BSplineKokkos2D kokkosBSP(order, ctrlPtsX, ctrlPtsY, knots); + + auto intView = Kokkos::create_mirror_view(kokkosBSP.getOrder()); + Kokkos::deep_copy(intView, kokkosBSP.getOrder()); + //Testing the order initialization + double diff = std::fabs(order - intView(0)); + if (diff > EPSILON) { + std::cout << "Order difference : " << diff << std::endl; + std::cout << "Serial: " << order << " Kokkos: " << intView(0) << std::endl; + retVal = 1; + } + + auto double2DView = Kokkos::create_mirror_view(kokkosBSP.getCtrlPts()); + Kokkos::deep_copy(double2DView, kokkosBSP.getCtrlPts()); + + //Testing ctrlPts initialization + double xDiff, yDiff; + for (int i = 0; i < ctrlPtsX.size(); i++) { + xDiff = std::fabs(ctrlPtsX[i] - double2DView(i, 0)); + yDiff = std::fabs(ctrlPtsY[i] - double2DView(i, 1)); + if (xDiff > EPSILON || yDiff > EPSILON) { + std::cout << "CtrlPts difference: x = " << xDiff << " y = " << yDiff << std::endl; + std::cout << "Serial: " << ctrlPtsX[i] << ", " << ctrlPtsY[i] << std::endl; + std::cout << "Kokkos: " << double2DView(i, 0) << ", " << double2DView(i, 1) << std::endl; + retVal = 1; + } + } + + //Test knots initialization + auto doubleView = Kokkos::create_mirror_view(kokkosBSP.getKnots()); + Kokkos::deep_copy(doubleView, kokkosBSP.getKnots()); + for (int i = 0; i < knots.size(); i++) { + diff = std::fabs(knots[i] - doubleView(i)); + if (diff > EPSILON) { + std::cout << "Knots difference: " << diff << std::endl; + std::cout << "Serial: " << knots[i] << std::endl; + std::cout << "Kokkos: " << doubleView(i) << std::endl; + retVal = 1; + } + } + + //Test 1st deriv coef initialization + Kokkos::View coef ("getCPCoe2", kokkosBSP.getCP1stD().extent(0)); + coef = kokkosBSP.getCP1stD(); + auto mv_coef = Kokkos::create_mirror_view(coef); + + Kokkos::deep_copy(mv_coef, coef); + std::vector serialCoefX = serialBSP.x.getCtrlPts_1st(); + std::vector serialCoefY = serialBSP.y.getCtrlPts_1st(); + + for (int i = 0; i < serialCoefX.size(); i++) { + double diffX = std::fabs(serialCoefX[i] - mv_coef(i, 0)); + double diffY = std::fabs(serialCoefY[i] - mv_coef(i, 1)); + if (diffX > EPSILON || diffY > EPSILON) { + std::cout << "diffX: " << diffX << std::endl; + std::cout << "diffY: " << diffY << std::endl; + std::cout << "Serial 1st coef: " << serialCoefX[i] << ", " << serialCoefY[i] << std::endl; + std::cout << "Kokkos 1st coef: " << mv_coef(i, 0)<< ", " << mv_coef(i, 1) << std::endl; + retVal = 1; + } + } + + Kokkos::Viewcoef2 ("getCPCoe2", kokkosBSP.getCP2ndD().extent(0)); + coef2 = kokkosBSP.getCP2ndD(); + mv_coef = Kokkos::create_mirror_view(coef2); + + Kokkos::deep_copy(mv_coef, coef2); + + serialCoefX = serialBSP.x.getCtrlPts_2nd(); + serialCoefY = serialBSP.y.getCtrlPts_2nd(); + for (int i = 0; i < serialCoefX.size(); i++) { + double diffX = std::fabs(serialCoefX[i] - mv_coef(i, 0)); + double diffY = std::fabs(serialCoefY[i] - mv_coef(i, 1)); + if (diffX > EPSILON || diffY > EPSILON) { + std::cout << "diffX: " << diffX << std::endl; + std::cout << "diffY: " << diffY << std::endl; + std::cout << "Serial 2nd coef: " << serialCoefX[i] << ", " << serialCoefY[i] << std::endl; + std::cout << "Kokkos 2nd coef: " << mv_coef(i, 0)<< ", " << mv_coef(i, 1) << std::endl; + retVal = 1; + } + } + + } + Kokkos::finalize(); + return retVal; +} + diff --git a/testSplineKokkosBase.cpp b/testSplineKokkosBase.cpp new file mode 100644 index 0000000..37a78f6 --- /dev/null +++ b/testSplineKokkosBase.cpp @@ -0,0 +1,207 @@ +#include +#include "BSplineKokkos.h" +#include "BSpline.h" +#include + +//These imports are from testSplineFitting.cc +#include "splineInterpolation.h" +#include "curveReader.h" +#include +#include +#include +#include +#include +#include + +template +void printDView(Kokkos::View& dView) { + auto host_dView = Kokkos::create_mirror_view(dView); + Kokkos::deep_copy(host_dView, dView); + for (int i = 0; i < host_dView.extent(0); i++) { + std::cout << host_dView(i) << std::endl; + } +} + +template +void printIView(Kokkos::View& iView) { + auto host_iView = Kokkos::create_mirror_view(iView); + Kokkos::deep_copy(host_iView, iView); + for (int i = 0; i < host_iView.extent(0); i++) { + std::cout << host_iView(i) << std::endl; + } +} + +template +void printSpline(BSplineKokkos spline) { + std::cout << "Spline content" << std::endl; + Kokkos::View dView = spline.getCtrlPts(); + std::cout << "CtrlPts" << std::endl; + printDView(dView); + Kokkos::View iView = spline.getCPOffset(); + std::cout << "CtrlPts Offset" << std::endl; + printIView(iView); + dView = spline.getKnots(); + std::cout << "Knots" << std::endl; + printDView(dView); + iView = spline.getKnotsOffset(); + std::cout << "Knots Offset" << std::endl; + printIView(iView); + iView = spline.getOrder(); + std::cout << "Order" << std::endl; + printIView(iView); +} + +void printVector(const std::string& name, std::vector& vector) { + std::cout << name << ": " << vector.size() << std::endl; + for (int i = 0; i < vector.size(); i++) { + std::cout << vector[i] << std::endl; + } +} + +int main(int argc, char* argv[]) { + Kokkos::initialize(argc, argv); + { + //Check for what memory space we should use + //Sticking to just cuda space or host space for now + + #ifdef KOKKOS_ENABLE_CUDA + #define MemSpace Kokkos::CudaSpace + #endif + #ifndef MemSpace + #define MemSpace Kokkos::HostSpace + #endif + + using ExecutionSpace = MemSpace::execution_space; + //We will try to test with dummy data for BSplineKokkos creation + std::vector dummyCtrlPtsX{2.14, 3.98, 1.56, 9.10, 11.87}; + std::vector dummyCtrlPtsY{7.928, 12.18, 9.17, 6.67, 12.78}; + std::vector dummyKnotsX{5.99, 8.98, 17.71, 21.1, 2.154, 9.17, 6.32}; + int dOrder = 2; + //Initialize 1 Kokkos Spline + BSplineKokkos ex1(dOrder, dummyCtrlPtsX, dummyCtrlPtsY, dummyKnotsX); + + //Check if we actually initialized with correct values + std::cout << "ex1" << std::endl; + printSpline(ex1); + std::vector> splineList; + splineList.push_back(ex1); + + //Test constructor that takes in multiple BSplineKokkos objects + dummyCtrlPtsX.assign({1.20, 12.91, 7.11}); + dummyCtrlPtsY.assign({3.31, 7.89, 9.12}); + dummyKnotsX.assign({2.18, 5.56, 3.18, 16.15}); + dOrder = 3; + BSplineKokkos ex2(dOrder, dummyCtrlPtsX, dummyCtrlPtsY, dummyKnotsX); + splineList.push_back(ex2); + std::cout << "ex2" << std::endl; + printSpline(ex2); + + dummyCtrlPtsX.assign({18.29, 12.12, 5.92, 7.182, 1.16}); + dummyCtrlPtsY.assign({4.72, 5.78, 9.12, 1.29, 10.2}); + dummyKnotsX.assign({4.34, 29.3, 18.7, 15.14, 16.87, 9.99}); + dOrder = 5; + std::cout << "\t\tSTART OF MULTISPLINE CREATION" << std::endl; + BSplineKokkos ex3(dOrder, dummyCtrlPtsX, dummyCtrlPtsY, dummyKnotsX); + splineList.push_back(ex3); + std::cout << "ex3" << std::endl; + printSpline(ex3); + + BSplineKokkos fromList(splineList); + + std::cout << "Spline created from list of splines" << std::endl; + printSpline(fromList); + + std::cout << "\n\nChecking derivative calculation against serial version" << std::endl; + std::cout << "Number of commandline arguments: " << argc << std::endl; + for (int i = 0; i < argc; i++) { + std::cout << argv[i] << std::endl; + } + if (argc != 3) { + //Check for the arguments needed + std::cerr<<"Input arguments: " << std::endl; + std::cerr << "input csv need these columns: "; + std::cerr << "x,y,z,isOnCurve,angle,isMdlVtx" << std::endl; + return 1; + } + + //Set up before using curve reader + std::string inputCSV = argv[1]; + int extensionPos = inputCSV.rfind("."); + int slashPos = inputCSV.rfind("/"); + std::string fileNameNoExt = inputCSV.substr(slashPos+1, extensionPos); + double expectedCurveLength = std::stod(argv[2]); + auto curve = CurveReader::readCurveInfo(inputCSV); + + //Construct BSpline2d object + SplineInterp::BSpline2d serialBSP; + if (curve.x.size() == 2) { + serialBSP = SplineInterp::attach_piecewise_linear_curve(curve.x, curve.y); + } + else { + serialBSP = SplineInterp::fitCubicSplineToPoints(curve.x, curve.y); + } + + //Get the parameters and print the spline we initialized + std::vector ctrlPtsX, ctrlPtsY, knots, weight; + int order; + serialBSP.x.getpara(order, ctrlPtsX, knots, weight); + std::cout << "x component of the spline" << std::endl; + std::cout << "order: " << order << std::endl; + printVector("ctrlPtsX", ctrlPtsX); + printVector("knots", knots); + printVector("weight", weight); + serialBSP.y.getpara(order, ctrlPtsY, knots, weight); + std::cout << "y component of the spline" << std::endl; + std::cout << "order: " << order << std::endl; + printVector("ctrlPtsY", ctrlPtsY); + printVector("knots", knots); + printVector("weight", weight); + + //Now that we ensured that serial BSpline is created correctly we will now feed these data to our BSpline Kokkos Object + BSplineKokkos kokkosBSP(order, ctrlPtsX, ctrlPtsY, knots); + printSpline(kokkosBSP); + + //We can now test taking the first derivative + std::cout << "----Serial 1st derivative----" << std::endl; + for (int i = 0; i < ctrlPtsX.size(); i++) { + std::cout <<"Point: " << ctrlPtsX[i] << ", "; + std::cout << ctrlPtsY[i] << std::endl; + } + std::cout << "1st derivative, input x = 0: " << serialBSP.x.evalFirstDeriv(0) << ", " << serialBSP.y.evalFirstDeriv(0) << std::endl; + std::cout << "1st derivative, input x = 1: " << serialBSP.x.evalFirstDeriv(1) << ", " << serialBSP.y.evalFirstDeriv(1) << std::endl; + + std::cout << "----Kokkos 1st derivative----" << std::endl; + //Get and copy the control points here so we could print them + Kokkos::View dView = kokkosBSP.getCtrlPts(); + auto host_dView = Kokkos::create_mirror_view(dView); + Kokkos::deep_copy(host_dView, dView); + for (int i = 1; i < host_dView.extent(0); i+=2) { + std::cout <<"Point: " << host_dView(i-1) << ", " << host_dView(i) << std::endl; + } + std::cout << "1st derivative, input x = 0: " << kokkosBSP.eval1stDeriv({0}, 0)[0] << std::endl; + + std::cout << "1st derivative, input x = 1: " << kokkosBSP.eval1stDeriv({1}, 0)[0] << std::endl; + std::cout << "----Serial 2nd derivative----" << std::endl; + double x2ndD = serialBSP.x.evalSecondDeriv(0); + double y2ndD = serialBSP.y.evalSecondDeriv(0); + std::cout << "2nd derivative, input x = 0: " << x2ndD << ", " << y2ndD << std::endl; + + x2ndD = serialBSP.x.evalSecondDeriv(1); + y2ndD = serialBSP.y.evalSecondDeriv(1); + std::cout << "2nd derivative, input x = 1: " << x2ndD << ", " << y2ndD << std::endl; + + std::cout << "----Kokkos 2nd derivative----" << std::endl; + //Get and copy the control points here so we could print them + dView = kokkosBSP.getCtrlPts(); + host_dView = Kokkos::create_mirror_view(dView); + Kokkos::deep_copy(host_dView, dView); + for (int i = 1; i < host_dView.extent(0); i+=2) { + std::cout <<"Point: " << host_dView(i-1) << ", " << host_dView(i) << std::endl; + } + std::cout << "2nd derivative, input x = 0: " << kokkosBSP.eval2ndDeriv(0, 0) << std::endl; + + std::cout << "2nd derivative, input x = 1: " << kokkosBSP.eval2ndDeriv(1, 0) << std::endl; + + } + Kokkos::finalize(); +}