Interpolation Schemes
Interpolators goven how the model parameters are interpolated across the field of view.
The Interp base class
- class piff.Interp[source]
The base class for interpolating a set of data vectors across the field of view.
In general, the interpolator is agnostic as to the meaning of the parameter vectors. These parameter vectors are passed as simple numpy arrays. They are imbued meaning by a Model instance. Thus, the same interpolators may be used with many different Model types.
The principal ways that interpolators will differ are:
Which properties of the star are used for their interpolation.
What functional form (or algorithm) is used to interpolate between measurements.
Whether the interpolator assumes each sample has a non-degenerate parameter fit, vs getting a differential quadratic form for chisq from each sample.
The answer to #3 is given in a boolean property degenerate_points.
This is essentially an abstract base class intended to define the methods that should be implemented by any derived class.
- getProperties(star)[source]
Extract the appropriate properties to use as the independent variables for the interpolation.
The base class implementation returns the field position (u,v) as a 1d numpy array.
- Parameters:
star – A Star instance from which to extract the properties to use.
- Returns:
A numpy vector of these properties.
- initialize(stars, logger=None)[source]
Initialize both the interpolator to some state prefatory to any solve iterations and initialize the stars for use with this interpolator.
The nature of the initialization is specific to the derived classes.
The base class implentation calls interpolateList, which will set the stars to have the right type object in its star.fit.params attribute.
- Parameters:
stars – A list of Star instances to use to initialize.
logger – A logger object for logging debug info. [default: None]
- Returns:
a new list of Star instances
- interpolate(star, logger=None, inplace=False)[source]
Perform the interpolation to find the interpolated parameter vector at some position.
- Parameters:
star – A Star instance to which one wants to interpolate
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned star is the same object as the input star. [default: False]
- Returns:
a Star instance holding the interpolated parameters
- interpolateList(stars, logger=None, inplace=False)[source]
Perform the interpolation for a list of stars.
The base class just calls interpolate(star) for each star in the list, but in many cases, this may be more efficiently done with a matrix operation, so we make it available for derived classes to override.
- Parameters:
stars – A list of Star instances to interpolate.
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned stars are the same objects as the input stars. [default: False]
- Returns:
a list of Star instances with interpolated parameters
- classmethod parseKwargs(config_interp, logger=None)[source]
Parse the interp field of a configuration dict and return the kwargs to use for initializing an instance of the class.
The base class implementation just returns the kwargs as they are, but derived classes might want to override this if they need to do something more sophisticated with them.
- Parameters:
config_interp – The interpolator field of the configuration dict, config[‘interp’]
logger – A logger object for logging debug info. [default: None]
- Returns:
a kwargs dict to pass to the initializer
- classmethod process(config_interp, logger=None)[source]
Parse the interp field of the config dict.
- Parameters:
config_interp – The configuration dict for the interp field.
logger – A logger object for logging debug info. [default: None]
- Returns:
an Interp instance
- property property_names
List of properties used by this interpolant.
- classmethod read(reader, name)[source]
Read an Interp via a reader object.
- Parameters:
reader – A reader object that encapsulates the serialization format.
name – Name associated with this interpolator in the serialized output.
- Returns:
an interpolator built from serialized information.
- set_num(num)[source]
If there are multiple components involved in the fit, set the number to use for this model.
- solve(stars, logger=None)[source]
Solve for the interpolation coefficients given some data.
- Parameters:
stars – A list of Star instances to interpolate between
logger – A logger object for logging debug info. [default: None]
- write(writer, name)[source]
Write an Interp via a writer object.
Note: this only writes the initialization kwargs to the fits extension, not the parameters.
The base class implementation works if the class has a self.kwargs attribute and these are all simple values (str, float, or int).
However, the derived class will need to implement _finish_write to write the solution parameters to a binary table.
- Parameters:
writer – A writer object that encapsulates the serialization format.
name – A name to associate with this interpolator in the serialized output.
Mean interpolation
- class piff.Mean(logger=None)[source]
The simplest possible interpolation scheme. It just finds the mean of the parameter vectors and uses that at every position.
Use type name “Mean” in a config field to use this interpolant.
- interpolate(star, logger=None, inplace=False)[source]
Perform the interpolation to find the interpolated parameter vector at some position.
- Parameters:
star – A Star instance to which one wants to interpolate
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned star is the same object as the input star. [default: False]
- Returns:
a Star instance holding the interpolated parameters
Polynomial interpolation
- class piff.Polynomial(order=None, orders=None, poly_type='poly', logger=None)[source]
An interpolator that uses scipy curve_fit command to fit a polynomial surface to each parameter passed in independently.
Use type name “Polynomial” in a config field to use this interpolant.
- Parameters:
order – The maximum order in the polynomial. i.e. the maximum value of i+j where p(u,v) = sum c_{ij} x^i y^j. [required, unless orders is given]
orders – Optionally, a list of orders, one for each parameter to be interpolated. This list should be the same length as the number of parameters that will be given to interpolate.
poly_type –
- A string, one of the keys in the polynomial_types
dictionary. By default these are “poly” (ordinary polynomials), “chebyshev”, “legendre”, “laguerre”, “hermite”. To add more you can add a key to polynomial_types with the value of a function with the signature of numpy.polynomial.polynomial.polyval2d
- _setup_indices(nparam)[source]
An internal function that sets up the indices, given the number of parameters
- _set_function(poly_type)[source]
An internal function that sets the type of the polynomial interpolation used. The options are the keys in polynomial_types.
- Parameters:
poly_type – A string value, one of the keys from polynomial_types
- _generate_indices(order)[source]
Generate, for internal use, the exponents i,j used in the polynomial model p(u,v) = sum c_{ij} u^i v^j
This needs to be called whenever the order of the polynomial fit is changed. At the moment that is just when an object is initialized or updated from file.
- Parameters:
order – The maximum order of the polynomial; the max value of i+j where p(x,y) ~ x^i y^j
- _pack_coefficients(parameter_index, C)[source]
Pack the 2D matrix of coefficients used as the model fit parameters into a vector of coefficients, either so this can be passed as a starting point into the curve_fit routine or for serialization to file.
For subclasses, the 2D matrix format could be whatever you wanted as long as _initialGuess, _interpolationModel, and the pack and unpack functions are consistent. The intialGuess method can return and the _interpolationModel can accept parameters in whatever form you like (e.g. could be a dict if you want) as long as _pack_coefficients can convert this into a 1D array and _unpack_coefficients convert it the other way.
- Parameters:
parameter_index – The integer index of the parameter; this lets us find the order of the parameter from self.
C – A 2D matrix of polynomial coefficients in the form that the numpy polynomial form is expecting: p(x,y,c) = sum_{i,j} c_{ij} x^i y^j
- Returns coeffs:
A 1D numpy array of coefficients of length self.nvariable
- _unpack_coefficients(parameter_index, coeffs)[source]
Unpack a sequence of parameters into the 2D matrix for the given parameter_index (which determines the order of the matrix)
This function is the inverse of _pack_coefficients
- Parameters:
parameter_index – The integer index of the parameter being used
coeffs – A 1D numpy array of coefficients of length self.nvariable
- Returns:
A 2D matrix of polynomial coefficients in the form that the numpy polynomial form is expecting: p(x,y,c) = sum_{i,j} c_{ij} x^i y^j
- _interpolationModel(pos, C)[source]
Generate the polynomial variation of some quantity at x and y coordinates for a given coefficient matrix.
TODO: At the moment this function is expecting a numpy array of shape (2,nstar) for the positions. We might want to use a galsim position object instead since I think some code elsewhere in Piff is expecting this.
This is an internal method used during the fitting.
- Parameters:
pos – A numpy array of the u,v positions at which to build the model
C – A 2D matrix of polynomial coefficients in the form that the numpy polynomial form is expecting: p(x,y,c) = sum_{i,j} c_{ij} x^i y^j
- Returns:
A numpy array of the calculated p_x(x)*p_y(y) where the p functions are polynomials.
- _initialGuess(positions, parameter, parameter_index)[source]
Make an initial guess for a set of parameters to use in the fit for your model. This is passed to curve_fit as a starting point.
- Parameters:
positions – A list of positions ((u,v) in this case) of stars.
parameter – A numpy array of the measured values of a parameter for each star
parameter_index – The integer index of the parameter being used
- Returns:
A guess for the parameters. In this case a 2D matrix which is zero everywhere except for (0,0). This should correspond to a flat function of the parameters with value given by the mean.
- initialize(stars, logger=None)[source]
Initialization is just solving the interpolator with current stars. This then calls interpolateList, which will set the stars to have the right type of object in its star.fit.params attribute.
- Parameters:
stars – A list of Star instances to use to initialize.
logger – A logger object for logging debug info. [default: None]
- Returns:
a new list of Star instances
- interpolate(star, logger=None, inplace=False)[source]
Perform the interpolation to find the interpolated parameter vector at some position.
- Parameters:
star – A Star instance to which one wants to interpolate
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned star is the same object as the input star. [default: False]
- Returns:
a Star instance holding the interpolated parameters
- solve(stars, logger=None)[source]
Solve for the interpolation coefficients given some data, using the scipy.optimize.curve_fit routine, which uses Levenberg-Marquardt to find the least-squares solution.
This currently assumes that our positions pos are just u and v.
- Parameters:
stars – A list of Star instances to use for the interpolation.
logger – A logger object for logging debug info. [default: None]
Interpolation using basis functions
- class piff.BasisInterp[source]
An Interp class that works whenever the interpolating functions are linear sums of basis functions. Does things the “slow way” to be stable to degenerate fits to individual stars, instead of fitting to parameter sets produced by single stars.
First time coding this we will assume that each element of the PSF parameter vector p is a linear combination of the same set of basis functions across the focal plane,
\[p_i = \sum_{j} q_{ij} K_j(u,v,other stellar params).\]The property degenerate_points is set to True to indicate that this interpolator uses the alpha/beta quadratic form of chisq for each sample, rather than assuming that a best-fit parameter vector is available at every sample.
Internally we’ll store the interpolation coefficients in a 2d array of dimensions (nparams, nbases)
Note: This is an abstract base class. The concrete class you probably want to use is BasisPolynomial.
- basis(star)[source]
Return 1d array of polynomial basis values for this star
- Parameters:
star – A Star instance
- Returns:
1d numpy array with values of u^i v^j for 0<i+j<=order
- constant(value=1.0)[source]
Return 1d array of coefficients that represent a polynomial with constant value.
- Parameters:
value – The value to use as the constant term. [default: 1.]
- Returns:
1d numpy array with values of u^i v^j for 0<i+j<=order
- initialize(stars, logger=None)[source]
Initialize both the interpolator to some state prefatory to any solve iterations and initialize the stars for use with this interpolator.
This class will initialize everything to have constant PSF parameter vector taken from the first Star in the list.
- Parameters:
stars – A list of Star instances to use to initialize.
logger – A logger object for logging debug info. [default: None]
- Returns:
A new list of Stars which have their parameters initialized.
- interpolate(star, logger=None, inplace=False)[source]
Perform the interpolation to find the interpolated parameter vector at some position.
- Parameters:
star – A Star instance to which one wants to interpolate
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned star is the same object as the input star. [default: False]
- Returns:
a Star instance holding the interpolated parameters
- solve(stars, logger=None)[source]
Solve for the interpolation coefficients given some data. The StarFit element of each Star in the list is assumed to hold valid alpha and beta members specifying depending of chisq on differential changes to its parameter vector.
- Parameters:
stars – A list of Star instances to interpolate between
logger – A logger object for logging debug info. [default: None]
- class piff.BasisPolynomial(order, keys=('u', 'v'), max_order=None, solver='scipy', use_qr=False, logger=None)[source]
A version of the Polynomial interpolator that works with BasisModels and can use the quadratic form of the chisq information it calculates. It works better than the regular Polynomial interpolator when there is missing or degenerate information.
The order is the highest power of a key to be used. This can be the same for all keys or you may provide a list of separate order values to be used for each key. (e.g. you may want to use 2nd order in the positions, but only 1st order in the color).
All combinations of powers of keys that have total order <= max_order are used. The maximum order is normally the maximum order of any given key’s order, but you may specify a larger value. (e.g. to use 1, x, y, xy, you would specify order=1, max_order=2.)
There are several options for what code to use for doing the linear algebra, controlled by the
solverparameter, which can take one of the following values:“scipy” uses regular numpy array functionality to build the matrices, and then uses
scipy.linalg.solveto find the solution. It starts by assuming the matrix is positive definite, and it falls back to an SVD solution when that is not the case.“qr” will also use numpy to build the matrices, but then it uses QR decomposition for the solution rather than the more direct least squares solution. QR decomposition requires more memory than the default and is somewhat slower (nearly a factor of 2); however, it is significantly less susceptible to numerical errors from high condition matrices.
“jax” uses the JAX module for building and solving the linear algebra equations rather than numpy/scipy. It should be equivalent in its results to the “scipy” option, but it may be faster if a multi-core cpu or gpu is available.
“cpp” uses the Eigen linear algebra package in C++ to build and solve the linear algebra equations rather than numpy/scipy. On a single core cpu (and more), it will be faster than default “scipy” solver if the number of training stars is more than ~30. With a Piff config using
PixelGridwithsize=25andinterp="Lanczos(11)"and a second order polynomial for interpolation, and running on ~O(100) PSFs reserved stars on a 4GB single core CPU, “cpp” solver is 60% faster than “scipy” solver.
Use type name “BasisPolynomial” in a config field to use this interpolant.
- Parameters:
order – The order to use for each key. Can be a single value (applied to all keys) or an array matching number of keys.
keys – List of keys for properties that will be used as the polynomial arguments. [default: (‘u’,’v’)]
max_order – The maximum total order to use for cross terms between keys. [default: None, which uses the maximum value of any individual key’s order] If this is an integer, it applies to all pairs, but you may also specify a dict mapping pairs of keys to an integer. E.g. {(‘u’,’v’):3, (‘u’,’z’):0, (‘v’,’z’):0}. This sets the maximum order for cross terms between these pairs. Furthermore, any pairs for which you want to skip cross terms (max=0) may be omitted from the dict.
solver – Which solver to use. Solvers available are “scipy”, “qr”, “jax”, “cpp”. See above for details. [default: ‘scipy’]
logger –
A logger object for logging debug info. [default: None]
- basis(star)[source]
Return 1d array of polynomial basis values for this star
- Parameters:
star – A Star instance
- Returns:
1d numpy array with values of u^i v^j for 0<i+j<=order
- constant(value=1.0)[source]
Return 1d array of coefficients that represent a polynomial with constant value.
- Parameters:
value – The value to use as the constant term. [default: 1.]
- Returns:
1d numpy array with values of u^i v^j for 0<i+j<=order
- getProperties(star)[source]
Extract the appropriate properties to use as the independent variables for the interpolation.
The base class implementation returns the field position (u,v) as a 1d numpy array.
- Parameters:
star – A Star instance from which to extract the properties to use.
- Returns:
A numpy vector of these properties.
- property property_names
List of properties used by this interpolant.
K-Nearest Neighbors
- class piff.KNNInterp(keys=('u', 'v'), n_neighbors=15, weights='uniform', algorithm='auto', p=2, logger=None)[source]
An interpolator that uses sklearn KNeighborsRegressor to interpolate a single surface
Use type name “KNN” or “KNearestNeighbors” in a config field to use this interpolant.
- Parameters:
keys – A list of star attributes to interpolate from [default: (‘u’, ‘v’)]
n_neighbors – Number of neighbors used for interpolation. [default: 15]
weights – Weight function used in prediction. Possible values are ‘uniform’, ‘distance’, and a callable function which accepts an array of distances and returns an array of the same shape containing the weights. [default: ‘uniform’]
algorithm – Algorithm used to compute nearest neighbors. Possible values are ‘ball_tree’, ‘kd_tree’, ‘brute’, and ‘auto’, which tries to determine the best choice. [default: ‘auto’]
p – Power parameter of distance metrice. p=2 is default euclidean distance, p=1 is manhattan. [default: 2]
logger –
A logger object for logging debug info. [default: None]
- _fit(locations, targets, logger=None)[source]
Update the Neighbors Regressor with data
- Parameters:
locations – The locations for interpolating. (n_samples, n_features). (In sklearn parlance, this is ‘X’.)
targets – The target values. (n_samples, n_targets). (In sklearn parlance, this is ‘y’.)
- _predict(locations, logger=None)[source]
Predict from knn.
- Parameters:
locations – The locations for interpolating. (n_samples, n_features). In sklearn parlance, this is ‘X’
- Returns:
Regressed parameters y (n_samples, n_targets)
- getProperties(star, logger=None)[source]
Extract the appropriate properties to use as the independent variables for the interpolation.
Take self.keys from star.data
- Parameters:
star – A Star instances from which to extract the properties to use.
- Returns:
A np vector of these properties.
- initialize(stars, logger=None)[source]
Initialize both the interpolator to some state prefatory to any solve iterations and initialize the stars for use with this interpolator.
- Parameters:
stars – A list of Star instances to interpolate between
logger – A logger object for logging debug info. [default: None]
- interpolate(star, logger=None, inplace=False)[source]
Perform the interpolation to find the interpolated parameter vector at some position.
- Parameters:
star – A Star instance to which one wants to interpolate
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned star is the same object as the input star. [default: False]
- Returns:
a Star instance with its StarFit member holding the interpolated parameters
- interpolateList(stars, logger=None, inplace=False)[source]
Perform the interpolation for a list of stars.
- Parameters:
stars – A list of Star instances to which to interpolate.
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned stars are the same objects as the input stars. [default: False]
- Returns:
a list of Star instances with interpolated parameters
- property property_names
List of properties used by this interpolant.
Gaussian process interpolation
- class piff.GPInterp(keys=('u', 'v'), kernel='RBF(1)', optimizer='isotropic', normalize=True, l0=3000.0, white_noise=0.0, n_neighbors=4, average_fits=None, nbins=20, min_sep=None, max_sep=None, rows=None, logger=None)[source]
An interpolator that models the underlying field as a Gaussian process.
Gaussian process regression, also known as “Kriging,” assumes that the parameters are drawn from a multi-dimensional Gaussian random field across the (u, v) space. It requires an estimate of the spatial covariance function of p, commonly referred to as the kernel.
The interpolation estimate at an arbitrary location (u, v) is the minimum variance unbiased estimate from the Gaussian distribution at that location conditioned on the values of the parameters measured at all the PSF stars.
The implemention of this class use the treegp module. (https://github.com/PFLeget/treegp) It can use any of the kernels defined in scikit-learn (https://scikit-learn.org/) to define the covariance matrix as well as a custom VonKarman kernel defined in treegp. The default kernel is the so-called “squared exponential” or “radial basis function” kernel, known as RBF in scikit-learn.
The default behavior involves measuring the radial two-point correlation function with TreeCorr (https://github.com/rmjarvis/TreeCorr) and then fitting the hyper-parameters of the kernel that best fits this measurement. This can be done either isotropically or anisotropically. There are also options to use the traditional maximum likelihood optimization or no optimization if preferred. See the
optimizerparameter below.Use type name “GP” or “GaussianProcess” in a config field to use this interpolant.
- Parameters:
keys – A list of keys for properties that will be interpolated. Must be 2 properties, which can be used to calculate a 2-point correlation function. [default: (‘u’,’v’)]
kernel – A string that can be evaled to make a
sklearn.gaussian_process.kernels.Kernelobject. Could be also a list of Kernels objects (one per PSF param). The reprs of sklear Kernels will work, as well as the repr of a custom treegp VonKarman object. [default: ‘RBF(1)’]optimizer –
Indicates which techniques to use for optimizing the kernel. Four options are available:
”isotropic” = use an isotropic radial 2-point correlation function estimated by TreeCorr.
”anisotropic” = use an anisotropic two-dimensional 2-point correlation function estimated by TreeCorr.
”likelihood” = use the classical Gaussian process maximum likelihood method.
”none” = don’t do any kernal optimization.
[default: “isotropic”]
rows – A list of integer which indicates which rows of Star.fit.param need to be interpolated using GPs. [default: None, which means all rows]
normalize – Whether to normalize the interpolation parameters to have a mean of 0. Normally, the parameters being interpolated are not mean 0, so you would want this to be True, but if your parameters have an a priori mean of 0, then subtracting off the realized mean would be invalid. [default: True]
white_noise – A float value that indicate the ammount of white noise that you want to use during the gp interpolation. This is an additional uncorrelated noise added in quadrature to the measured uncertainties of the PSF parameters. This should be given as a “sigma” value, not a variance. [default: 0.]
nbins – Number of bins (in each direction using a 2D correlation function) used in TreeCorr to compute the 2-point correlation function. Used only if optimizer is “isotropic” or “anisotropic”. [default: 20]
min_sep – Minimum separation between pairs when computing 2-point correlation function in arcsec (or more generally the same units as the keys). Compute automaticaly if it is not given. Used only if optimizer is “isotropic” or “anisotropic”. [default: None]
max_sep – Maximum separation between pairs when computing 2-point correlation function in arcsec (or more generally the same units as the keys). Compute automaticaly if it is not given. Used only if optimizer is “isotropic” or “anisotropic”. [default: None]
l0 – Initial guess for correlation length when optimzer is “anisotropic” in arcsec (or more generally the same units as the keys). [default: 3000.]
average_fits – A fits file that have the spatial average functions of PSF parameters build in it. Build using meanify and piff output across different exposures. See meanify documentation for details. [default: None]
n_neighbors – Number of neighbors to use for interpolating the spatial average using a KNeighbors interpolation. Used only if average_fits is not None. [defaulf: 4]
logger –
A logger object for logging debug info. [default: None]
- _fit(X, y, y_err=None, logger=None)[source]
Update the GaussianProcess with data
- Parameters:
X – The independent covariates. (n_samples, n_features)
y – The dependent responses. (n_samples, n_targets)
y_err – Error of y. (n_samples, n_targets)
logger – A logger object for logging debug info. [default: None]
- _predict(Xstar)[source]
Predict responses given covariates. :param X: The independent covariates at which to interpolate. (n_samples, n_features). :returns: Regressed parameters (n_samples, n_targets)
- getProperties(star, logger=None)[source]
Extract the appropriate properties to use as the independent variables for the interpolation.
Take self.keys from star.data
- Parameters:
star – A Star instances from which to extract the properties to use.
- Returns:
A np vector of these properties.
- initialize(stars, logger=None)[source]
Initialize both the interpolator to some state prefatory to any solve iterations and initialize the stars for use with this interpolator.
- Parameters:
stars – A list of Star instances to interpolate between
logger – A logger object for logging debug info. [default: None]
- interpolate(star, logger=None, inplace=False)[source]
Perform the interpolation to find the interpolated parameter vector at some position.
- Parameters:
star – A Star instance to which one wants to interpolate
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned star is the same object as the input star. [default: False]
- Returns:
a Star instance with its StarFit member holding the interpolated parameters
- interpolateList(stars, logger=None, inplace=False)[source]
Perform the interpolation for a list of stars.
- Parameters:
star_list – A list of Star instances to which to interpolate.
logger – A logger object for logging debug info. [default: None]
inplace – Whether to update the parameters in place, in which case the returned stars are the same objects as the input stars. [default: False]
- Returns:
a list of Star instances with interpolated parameters
- property property_names
List of properties used by this interpolant.