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.

Sign up for GitHub

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

Jump to bottom

Simplified implementing new models #980 #1203

Merged
merged 7 commits into from
Jul 8, 2013
Merged

Simplified implementing new models #980 #1203

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8f61867
Simplified creating new models #980
adonath Jun 25, 2013
59cd4a6
WIP
adonath Jun 28, 2013
b00b6f9
WIP
adonath Jul 5, 2013
46cc8ac
Fixed the points discussed in #1203
adonath Jul 5, 2013
fb4b5c7
Removed **kwargs in LSQFitter constructor
adonath Jul 5, 2013
38aaa50
Fixed docstring error in Custom1DModel
adonath Jul 5, 2013
a6f74dc
Removed jacobian_func and added estimate_jacobian to docstring
adonath Jul 6, 2013
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 49 additions & 2 deletions astropy/modeling/core.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,10 @@
from . import constraints
from .utils import InputParameterError


__all__ = ['Model', 'ParametricModel', 'PCompositeModel', 'SCompositeModel',
'LabeledInput', '_convert_input', '_convert_output']
'LabeledInput', '_convert_input', '_convert_output',
'Parametric1DModel']


def _convert_input(x, pdim):
Expand Down Expand Up @@ -705,7 +707,7 @@ def __init__(self, transforms, inmap=None, outmap=None, n_inputs=None, n_outputs
def inverse(self):
try:
transforms = [tr.inverse() for tr in self._transforms[::-1]]
except NotIMplementedError:
except NotImplementedError:
raise
if self._inmap is not None:
inmap = self._inmap[::-1]
Expand Down Expand Up @@ -824,3 +826,48 @@ def __call__(self, *data):
for tr in self._transforms:
result += tr(*data)
return result


class Parametric1DModel(ParametricModel):
"""
Base class for one dimensional parametric models
This class provides an easier interface to defining new models.
Examples can be found in functional_models.py
Parameters
----------
parameter_dict : dictionary
Dictionary of model parameters with initialisation values
{'parameter_name': 'parameter_value'}
"""
deriv = None
linear = False

def __init__(self, param_dict, **cons):
# Get parameter dimension
param_dim = np.size(param_dict[self.param_names[0]])

# Initialize model parameters. This is preliminary as long there is
# no new parameter class. It may be more reasonable and clear to init
# the parameters in the model constructor itself, with constraints etc.
for param_name in self.param_names:
setattr(self, "_" + param_name, parameters.Parameter(name=param_name,
val=param_dict[param_name], mclass=self, param_dim=param_dim))

super(Parametric1DModel, self).__init__(self.param_names, n_inputs=1,
n_outputs=1, param_dim=param_dim, **cons)

def __call__(self, x):
"""
Transforms data using this model.
Parameters
----------
x : array like or a number
input
"""
x, fmt = _convert_input(x, self.param_dim)
result = self.eval(x, *self.param_sets)
Copy link
Contributor

Choose a reason for hiding this comment

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

Initially the parameters list attribute was used for evaluating the model but because of the multiple parameter sets the evaluation was changed to use param_sets, which was a mistake since constructing param_sets uses properties to look up the values. I would like to go back to using a list or as close as possible, something like this:

p = np.asarray(self.parameters)
p.shape = (len(self.param_names), self.param_dim)
result = self.eval(x, *p)

Timing shows this is 1.5 x faster than using param_sets. If you agree with this, would you make the change?
There are other improvements that can be done but they are related also to constraints handling and deserve a separate PR.

Copy link
Contributor

Choose a reason for hiding this comment

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

Oh, I take this back. The reason it wasn't done this way to start is that this only works with scalar parameters.

return _convert_output(result, fmt)
33 changes: 19 additions & 14 deletions astropy/modeling/fitting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,10 +81,6 @@ def __init__(self, model):
self._fitpars = self.model.constraints.fitpars[:]
else:
self._fitpars = self.model._parameters[:]
if self._model.deriv is None:
self.dfunc = None
else:
self.dfunc = self._wrap_deriv
self._weights = None

@property
Expand Down Expand Up @@ -172,22 +168,22 @@ def _wrap_deriv(self, p, x, y, z=None):
if fixed_and_tied:
pars = self.model.constraints.modelpars
if z is None:
fullderiv = self.model.deriv(pars, x, y)
fullderiv = np.array(self.model.deriv(x, *pars))
else:
fullderiv = self.model.deriv(pars, x, y, z)
fullderiv = np.array(self.model.deriv(x, y, *pars))
ind = range(len(self.model.param_names))
for name in fixed_and_tied:
index = self.model.param_names.index(name)
ind.remove(index)
res = np.empty((fullderiv.shape[0], fullderiv.shape[1] - len(ind)))
res = fullderiv[:, ind]
res = np.empty((fullderiv.shape[1] - len(ind), fullderiv.shape[0]))
res = fullderiv[ind, :]
return res
else:
pars = p[:]
if z is None:
return self.model.deriv(pars, x, y)
return self.model.deriv(x, *pars)
else:
return self.model.deriv(pars, x, y, z)
return self.model.deriv(x, y, *pars)

def _validate_constraints(self):
fname = self.__class__.__name__
Expand Down Expand Up @@ -483,7 +479,8 @@ def covar(self):
log.info("Could not construct a covariance matrix")
return None

def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER, epsilon=EPS):
def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER,
epsilon=EPS, estimate_jacobian=False):
"""
Fit data to this model.
Expand All @@ -505,7 +502,10 @@ def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER, epsilon=EPS):
epsfcn is less than the machine precision, it is
assumed that the relative errors in the functions are
of the order of the machine precision.
estimate_jacobian : bool
If False (default) and if the model has a deriv method,
it will be used. Otherwise the Jacobian will be estimated.
If True, the Jacobian will be estimated in any case.
"""
from scipy import optimize
x = np.asarray(x, dtype=np.float)
Expand All @@ -527,9 +527,14 @@ def __call__(self, x, y, z=None, weights=None, maxiter=MAXITER, epsilon=EPS):
meas = np.asarray(z, dtype=np.float)
farg = (x, y, meas)

if self._model.deriv is None or estimate_jacobian:
self.dfunc = None
else:
self.dfunc = self._wrap_deriv

self.fitpars, status, dinfo, mess, ierr = optimize.leastsq(
self.errorfunc, self.fitpars, args=farg, Dfun=self.dfunc,
maxfev=maxiter, epsfcn=epsilon, full_output=True)
col_deriv=1, maxfev=maxiter, epsfcn=epsilon, full_output=True)
self.fit_info.update(dinfo)
self.fit_info['status'] = status
self.fit_info['message'] = mess
Expand Down Expand Up @@ -737,7 +742,7 @@ def errorfunc(self, fps, *args):
plen = sl.stop - sl.start
mpars.extend(mfpars[:plen])
del mfpars[:plen]
modelfit = model.eval(margs[:-1], mpars)
modelfit = model.eval(margs[:-1], *mpars)
fitted.extend(modelfit - margs[-1])
return np.ravel(fitted)

Expand Down
Loading