Skip to content

Commit

Permalink
Addressed Dylan's Comments from PR#283.
Browse files Browse the repository at this point in the history
  • Loading branch information
@ptranq
ptranq committed May 7, 2024
1 parent 706876b commit 841b47f
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 79 deletions.
65 changes: 32 additions & 33 deletions lib/algo/manifold_interp/SnapshotInterpolator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,13 @@

/**
* Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE
* PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE
* PIECEWISE CUBIC INTERPOLANTS", Fritsch and Butland (1984). as well as "MONOTONE
* PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980)
*
*/
#include "SnapshotInterpolator.h"

#include <cfloat>
#include <limits.h>
#include <cmath>
#include "linalg/Vector.h"
Expand All @@ -31,30 +32,24 @@ using namespace std;

namespace CAROM {

SnapshotInterpolator::SnapshotInterpolator()
{

}


std::vector<Vector*> SnapshotInterpolator::interpolate(std::vector<Vector*>
snapshot_ts,
std::vector<Vector*> snapshots,
std::vector<Vector*> output_ts)
void SnapshotInterpolator::interpolate(std::vector<Vector*>& snapshot_ts,
std::vector<Vector*>& snapshots,
std::vector<Vector*>& output_ts,
std::vector<Vector*>& output_snapshots)
{
CAROM_VERIFY(snapshot_ts.size() == snapshots.size());
CAROM_VERIFY(snapshot_ts.size() > 0);
CAROM_VERIFY(output_ts.size() > 0);
CAROM_VERIFY(output_snapshots.size() == 0);

int n_out = output_ts.size();
int n_snap = snapshots.size();
int n_dim = snapshots[0]->dim();

std::vector<Vector *> output_snapshots;
for(int i = 0; i < n_out; ++i)
{
Vector* temp_snapshot =new Vector(snapshots[0]->dim(),
snapshots[0]->distributed());
Vector* temp_snapshot = new Vector(snapshots[0]->dim(),
snapshots[0]->distributed());
output_snapshots.push_back(temp_snapshot);
}

Expand Down Expand Up @@ -113,7 +108,9 @@ std::vector<Vector*> SnapshotInterpolator::interpolate(std::vector<Vector*>
d_temp = 3*delta[n_snap-2];
}
d.push_back(d_temp);
while(counter < n_out && output_ts[counter]->getData()[0] <= t_in[n_snap-1])

while(counter < n_out
&& output_ts[counter]->getData()[0] - t_in[n_snap-1] <= FLT_EPSILON )
{
t = output_ts[counter]->getData()[0];
output_snapshots[counter]->getData()[i] = y_in[n_snap-2]*computeH1(t,
Expand All @@ -124,18 +121,20 @@ std::vector<Vector*> SnapshotInterpolator::interpolate(std::vector<Vector*>
counter += 1;
}
}
return output_snapshots;
}

std::vector<Vector*> SnapshotInterpolator::interpolate(std::vector<Vector*>
snapshot_ts,
std::vector<Vector*> snapshots,
int n_out,
std::vector<Vector*>* output_ts)
void SnapshotInterpolator::interpolate(std::vector<Vector*>&
snapshot_ts,
std::vector<Vector*>& snapshots,
int n_out,
std::vector<Vector*>& output_ts,
std::vector<Vector*>& output_snapshots)
{
CAROM_VERIFY(snapshot_ts.size() == snapshots.size());
CAROM_VERIFY(snapshot_ts.size() > 0);
CAROM_VERIFY(n_out > 0);
CAROM_VERIFY(output_ts.size() == 0 && output_snapshots.size() == 0);


int n_snap = snapshots.size();
int n_dim = snapshots[0]->dim();
Expand All @@ -145,21 +144,21 @@ std::vector<Vector*> SnapshotInterpolator::interpolate(std::vector<Vector*>
double t_min = snapshot_ts[0]->getData()[0];
double t_max = snapshot_ts[n_snap-1]->getData()[0];
double dt = (t_max-t_min)/(n_out-1);
output_ts->clear();

output_ts.clear();

for(int i = 0; i < n_out; ++i)
{
Vector* temp_t = new Vector(1, false);
temp_t->getData()[0] = t_min + i*dt;
output_ts->push_back(temp_t);
output_ts.push_back(temp_t);
}
std::cout << "Interpolating to " << n_out << " snapshots, from " << n_snap <<
std::endl;
return interpolate(snapshot_ts,snapshots,*output_ts);
interpolate(snapshot_ts,snapshots,output_ts, output_snapshots);
}

double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1,
const double SnapshotInterpolator::computeDerivative(double S1, double S2,
double h1,
double h2)
{
double d = 0.0;
Expand All @@ -175,41 +174,41 @@ double SnapshotInterpolator::computeDerivative(double S1, double S2, double h1,
return d;
}

double SnapshotInterpolator::computeH1(double x, double xl, double xr)
const double SnapshotInterpolator::computeH1(double x, double xl, double xr)
{
double h = xr - xl;
return computePhi((xr-x)/h);
}

double SnapshotInterpolator::computeH2(double x, double xl, double xr)
const double SnapshotInterpolator::computeH2(double x, double xl, double xr)
{
double h = xr - xl;
return computePhi((x-xl)/h);
}

double SnapshotInterpolator::computeH3(double x, double xl, double xr)
const double SnapshotInterpolator::computeH3(double x, double xl, double xr)
{
double h = xr-xl;
return -h*computePsi((xr-x)/h);
}

double SnapshotInterpolator::computeH4(double x, double xl, double xr)
const double SnapshotInterpolator::computeH4(double x, double xl, double xr)
{
double h = xr-xl;
return h*computePsi((x-xl)/h);
}

double SnapshotInterpolator::computePhi(double t)
const double SnapshotInterpolator::computePhi(double t)
{
return 3.*pow(t,2.) - 2*pow(t,3.);
}

double SnapshotInterpolator::computePsi(double t)
const double SnapshotInterpolator::computePsi(double t)
{
return pow(t,3.) - pow(t,2.);
}

int SnapshotInterpolator::sign(double a)
const int SnapshotInterpolator::sign(double a)
{
double TOL = 1e-15;
if(abs(a) < TOL)return 0;
Expand Down
67 changes: 47 additions & 20 deletions lib/algo/manifold_interp/SnapshotInterpolator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,39 +12,66 @@ class Vector;

/**
* Implements PCHIP algorithm. Based on "A METHOD FOR CONSTRUCTING LOCAL MONOTONE
* PIECEWISE CUBIC INTERPOLANTS", Fritchs and Butland (1984). as well as "MONOTONE
* PIECEWISE CUBIC INTERPOLANTS", Fritsch and Butland (1984). as well as "MONOTONE
* PIECEWISE CUBIC INTERPOLATION," Fritsch and Carlson (1980)
*
*/
class SnapshotInterpolator
{
public:
SnapshotInterpolator();
~SnapshotInterpolator();

std::vector<Vector*> interpolate(std::vector<Vector*> snapshot_ts,
std::vector<Vector*> snapshots,
std::vector<Vector*> output_ts);

std::vector<Vector*> interpolate(std::vector<Vector*> snapshot_ts,
std::vector<Vector*> snapshots,
int n_out,
std::vector<Vector*>* output_ts);

std::vector<Vector*> interpolate();
SnapshotInterpolator()
{}
~SnapshotInterpolator()
{}

/**
* @brief Compute new snapshots interpolated from snapshot_ts to
* output_ts.
*
* @param[in] snapshot_ts The parameter points.
* @param[in] snapshots The rotation matrices associated with
* each parameter point.
* @param[in] output_ts Requested times for interpolated
* snapshots
* @param[out] output_snapshots snapshots at output_ts interpolated
* from snapshot_ts
*/
void interpolate(std::vector<Vector*>& snapshot_ts,
std::vector<Vector*>& snapshots,
std::vector<Vector*>& output_ts,
std::vector<Vector*>&output_snapshots);

/**
* @brief Compute new snapshots interpolated from snapshot_ts to
* output_ts.
*
* @param[in] snapshot_ts The parameter points.
* @param[in] snapshots The rotation matrices associated with
* each parameter point.
* @param[in] n_out Number of output snapshots requested
* @param[out] output_ts std::vector of CAROM::Vectors that are
the times of the interpolated snapshots.
* @param[out] output_snapshots snapshots at output_ts interpolated
* from snapshot_ts
*/
void interpolate(std::vector<Vector*>& snapshot_ts,
std::vector<Vector*>& snapshots,
int n_out,
std::vector<Vector*>& output_ts,
std::vector<Vector*>& output_snapshots);

private:

double computeDerivative(double S1, double S2, double h1, double h2);
double computeH1(double x, double xl, double xr);
double computeH2(double x, double xl, double xr);
double computeH3(double x, double xl, double xr);
double computeH4(double x, double xl, double xr);
double computePhi(double t);
double computePsi(double t);
int sign(double a);
const double computeDerivative(double S1, double S2, double h1, double h2);
const double computeH1(double x, double xl, double xr);
const double computeH2(double x, double xl, double xr);
const double computeH3(double x, double xl, double xr);
const double computeH4(double x, double xl, double xr);
const double computePhi(double t);
const double computePsi(double t);
const int sign(double a);
};

}
Expand Down
66 changes: 40 additions & 26 deletions unit_tests/snapshot_interpolator_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#include "algo/manifold_interp/SnapshotInterpolator.h"
#include "linalg/Vector.h"
#include <cfloat>
#include <cmath>
#include <iostream>

Expand All @@ -27,6 +28,7 @@ int main(int argc, char *argv[])
{
int n_out = 25;
int n_snap = 11;
bool SUCCESS = true;
//input times
vector<double> t{ 0, 2, 3, 5, 6, 8, 9, 11, 12, 14, 15};
//Test function from original PCHIP paper
Expand Down Expand Up @@ -84,39 +86,51 @@ int main(int argc, char *argv[])

SnapshotInterpolator* interp = new SnapshotInterpolator();

std::cout << "Beginning base interpolation function" << std::endl;
out_snapshots = interp->interpolate(times,snapshots,out_times);
std::cout << "Finished interpolation " << std::endl;
/*
for(int i = 0; i < out_snapshots.size(); ++i)
{
std::cout << "Time " << i << " is " << out_times[i]->item(0) << std::endl;
std::cout << "Reference at " << i << " is (" << reference_snapshots[i]->item(0) <<
"," << reference_snapshots[i]->item(1) << ")" << std::endl;
std::cout << "Snapshot interpolation at " << i << " is (" << out_snapshots[i]->item(0) <<
"," << out_snapshots[i]->item(1) << ")" << std::endl;
}
*/

interp->interpolate(times,snapshots,out_times,out_snapshots);

for(int i = 0; i < out_snapshots.size(); ++i)
{
std::cout << "Error at time[" << i << "] = " << out_times[i]->item(
0) << " is (" <<
reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << ","
<< reference_snapshots[i]->item(1) - out_snapshots[i]->item(
1) << ") " << std::endl;
if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item(
0)) > 10.*FLT_EPSILON)
{
SUCCESS = false;
std::cout << "Test failed." << std::endl;
return SUCCESS;
}
else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item(
1)) > 10.*FLT_EPSILON)
{
SUCCESS = false;
std::cout << "Test failed." << std::endl;
return SUCCESS;
}
}


std::cout << "Beginning variant interpolation function" << std::endl;
out_snapshots = interp->interpolate(times,snapshots,26,&out_times);
std::cout << "Finished interpolation " << std::endl;

out_snapshots.clear();
out_times.clear();

interp->interpolate(times,snapshots,26,out_times,out_snapshots);

for(int i = 0; i < out_snapshots.size(); ++i)
{
std::cout << "Error at time[" << i << "] = " << out_times[i]->item(
0) << " is (" <<
reference_snapshots[i]->item(0) - out_snapshots[i]->item(0) << ","
<< reference_snapshots[i]->item(1) - out_snapshots[i]->item(
1) << ") " << std::endl;
if( abs(reference_snapshots[i]->item(0)-out_snapshots[i]->item(
0)) > 10.*FLT_EPSILON)
{
SUCCESS = false;
std::cout << "Test failed." << std::endl;
return SUCCESS;
}
else if(abs(reference_snapshots[i]->item(1)-out_snapshots[i]->item(
1)) > 10.*FLT_EPSILON)
{
SUCCESS = false;
std::cout << "Test failed." << std::endl;
return SUCCESS;
}
}
std::cout << "Test Succeeded" << std::endl;
return SUCCESS;
}

0 comments on commit 841b47f

Please sign in to comment.