Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: Add polynomial features and extra basis functions to scipy.interpolate.rbf #7157

Closed
wants to merge 1 commit into from

Conversation

evanlimanto
Copy link
Contributor

@evanlimanto evanlimanto commented Mar 11, 2017

See section 2.1 of this paper.

@evanlimanto
Copy link
Contributor Author

Surface plot of sample data with degree = 2:
figure_1

import numpy as np
from scipy.interpolate import Rbf
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm

fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
X, Y, Z = X[60:100], Y[60:100], Z[60:100]
ax.plot_surface(X, Y, Z, rstride=8, cstride=8, alpha=0.3, color='red')

x = np.linspace(np.min(X), np.max(X), 30)
y = np.linspace(np.min(Y), np.max(Y), 30)
xv, yv = np.meshgrid(x, y)
rbfi = Rbf(X, Y, Z, degree=2)
di = rbfi(xv.flatten(), yv.flatten())
ax.plot_surface(xv, yv, di.reshape(xv.shape[0], xv.shape[1]), rstride=8, cstride=8, alpha=0.5, cmap=cm.coolwarm)
plt.show()

@evanlimanto evanlimanto changed the title ENH: Add polynomials and basis functions to scipy.interpolate.rbf ENH: Add polynomials and extra basis functions to scipy.interpolate.rbf Mar 11, 2017
@evanlimanto evanlimanto changed the title ENH: Add polynomials and extra basis functions to scipy.interpolate.rbf ENH: Add polynomial features and extra basis functions to scipy.interpolate.rbf Mar 11, 2017
@ev-br
Copy link
Member

ev-br commented Mar 12, 2017

This now needs a rebase, as it crossed with gh-7161.

You can either do a rebase yourself, or grab the rebased branch: https://github.com/ev-br/scipy/commits/pr/7157

self.A = self._init_function(r) - np.eye(self.N)*self.smooth
self.nodes = linalg.solve(self.A, self.di)
self.A_aug = self._generate_augmented_matrix(self.xi.T)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should not be stored on the instance: the matrix is only used once and can be large. If really needs be, it can be made into a property, or just discarded.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed.

combs = combinations(n_dims, deg)
for i, c in enumerate(combs):
POLY[:, i] = xi[:, c].prod(axis=1)
return bmat([[self.A, POLY], [POLY.T, None]], format='csr')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Style nitpick: This function is a bit confusing: it receives xi as an argument, and is being called with self.xi; it also uses self.degree which it reads off the instance; it attaches self.n_features but returns self.A_aug. This could probably be combined into def _compute_nodes(self): which computes self.nodes and discards all local variables, e.g. A_aug.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed all the above.

n_points, n_dims = xi.shape
# closed-form formula for number of polynomial coefficients
deg = self.degree if self.degree else -1
self.n_features = int(comb(deg + n_dims, n_dims))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What are features?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

features are the generated polynomial coefficients. Updated the variable name and description.

@ev-br
Copy link
Member

ev-br commented Mar 12, 2017

It would be helpful to explain a bit what is going on here. Also references should be added to the docstring.

It would also be helpful to explain a connection to the suggestions in
#4523 (comment).

@ev-br ev-br added enhancement A new feature or improvement needs-work Items that are pending response from the author scipy.interpolate labels Mar 12, 2017
@ev-br
Copy link
Member

ev-br commented Mar 17, 2017

Oops, something went wrong on the git side. Don't merge master, please rebase on master instead, yhen force-push (the latter is key). Do you need help untangling this?

@evanlimanto
Copy link
Contributor Author

Yeah @ev-br trying to fix this now, but how do I rebase all the previous commits and keep my latest ones?

@ev-br
Copy link
Member

ev-br commented Mar 17, 2017

You squashed your work into this single commit, dbb346c
Then it's probably easiest to cherry-pick it into a new branch and force-push:

$ git co master
$ git co new-branch
$ git cherry-pick dbb346c56d43c9e949d0b24b52e735e61777cca8
$ git push origin new-branch:rbf-polynomials --force     # assuming your scipy fork is named "origin"

This will replace the the current contents of this PR by that single commit.

@evanlimanto
Copy link
Contributor Author

Done, thanks! Just couldn't figure out the new-branch:rbf-polynomials part.

@evanlimanto evanlimanto force-pushed the rbf-polynomials branch 3 times, most recently from c291427 to 5b9252f Compare March 23, 2017 06:44

If callable, then it must take 2 arguments (self, r). The epsilon
parameter will be available as self.epsilon. Other keyword
arguments passed in will be available as well.

degree : integer, optional
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This description needs to be clarified.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also "adjustable constant" reads a bit oddly.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Clarified.

epsilon : float, optional
Adjustable constant for gaussian or multiquadrics functions
- defaults to approximate average distance between nodes (which is
a good start).
m : float, optional
Adjustable constant for polyharmonic splines.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, what is m, how do you adjust it? Guidelines here would be helpful to the user.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see how this works now. Perhaps something like "A parameter used to construct polyharmonic splines. Large m yields ... Typical values lie between ... and ..."

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Updated. I specified that m refers to the degree of the polyharmonic spline.

m : float, optional
Adjustable constant for polyharmonic splines.
p : float, optional
Adjustable constant for bessel functions.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Updated.

def _compute_nodes(self):
"""
First, we generate all possible combinations of the input over all
dimensions for the polynomial features. Then, these are stacked to
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"polynomial features"?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changed to powers of input dimensions.

|C.T| 0 |
---------

where C is a matrix of size (n, num_coeffs) generated from the
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

More detail is needed to understand what happens here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Updated.

@pv pv removed the needs-work Items that are pending response from the author label Jul 30, 2017
@stefanv
Copy link
Member

stefanv commented Nov 2, 2017

@treverhines I just came across your RBF work at https://github.com/treverhines/RBF

I was wondering if you could take a look at this pull-request, and consider whether we could work together to include your existing work into SciPy?

Thanks!

@ev-br
Copy link
Member

ev-br commented Apr 29, 2022

@evanlimanto very sorry for this PR slipping under the radar. I don't think we will be merging this in the end, now that the RBFInterpolator is in. Your work is much appreciated regardless!

What is the status here, is there something from this PR to be ported to the RBFInterpolator interface? @treverhines @segasai @stefanv @evanlimanto

@stefanv
Copy link
Member

stefanv commented Apr 30, 2022

It would be good to ask @treverhines that question. (Woops, see you already did :) ).

@ev-br ev-br moved this from N-D scattered to In progress in scipy.interpolate Apr 30, 2022
@ev-br ev-br moved this from In progress to triaged, waiting-for-contributor in scipy.interpolate May 5, 2022
@ev-br
Copy link
Member

ev-br commented Jun 23, 2022

Ok, let's close this. If there's something to port to the current RBF framework, we can always reopen or cherry-pick. Thank you very much @evanlimanto , your work is much appreciated regardless!

@ev-br ev-br closed this Jun 23, 2022
scipy.interpolate automation moved this from triaged, waiting-for-contributor to Abandoned / rejected Jun 23, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement A new feature or improvement scipy.interpolate
Projects
Status: Abandoned / rejected / supeseded
scipy.interpolate
Abandoned / rejected / supeseded
Development

Successfully merging this pull request may close these issues.

None yet

4 participants