<jupyter_start><jupyter_text>Overview of core astropy.modeling features- simple models- fitters- combined models- creating new models<jupyter_code>import numpy as np
%matplotlib inline
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>Simple models<jupyter_code>from astropy.modeling import models<jupyter_output><empty_output><jupyter_text>Models are defined with their parameters.<jupyter_code>models.Gaussian1D.param_names
gauss = models.Gaussian1D(amplitude=5.3, mean=3.1, stddev=1.1)<jupyter_output><empty_output><jupyter_text>Parameters are objects with attributes:<jupyter_code>gauss.mean
print(gauss.mean.name)
print(gauss.mean.value)

print(gauss.mean.default)
models.Gaussian1D(amplitude=1, stddev=1)<jupyter_output><empty_output><jupyter_text>astropy.modeling supports parameter constraints:<jupyter_code>print(gauss.mean.constraints) # more on this later
print(gauss.parameters)<jupyter_output><empty_output><jupyter_text>Parameter values can be changed by assignment:<jupyter_code>gauss.mean = 1.4
print(gauss.mean)
print(gauss.parameters)<jupyter_output><empty_output><jupyter_text>Models are evaluated like functions, by passing the inputs. Each dimension is a separate input.<jupyter_code>x = np.linspace(-2, 5, 100)
data = gauss(x)
plt.plot(x, data)<jupyter_output><empty_output><jupyter_text>Add noise to the data and fit a gaussian.<jupyter_code>data = gauss(x) + np.random.normal(0, 0.4, x.shape)
plt.plot(x, data)<jupyter_output><empty_output><jupyter_text>Fitters<jupyter_code>from astropy.modeling import fitting<jupyter_output><empty_output><jupyter_text>Create a fitter which uses the Levenberg-Marquardt optimization algorithm and least squares statistics.Then pass the model and the data to the fitter. The output is a new model with fitted parameters.<jupyter_code>fitter = fitting.LevMarLSQFitter()
model = fitter(gauss, x, data)
print(model.parameters)
#print(gauss.parameters)
plt.plot(x, data, label='data')
plt.plot(x, model(x), 'r', label='gaussian')
plt.legend(loc=2)<jupyter_output><empty_output><jupyter_text>Fitters support parameter constraints of type *fixed*, *tied* and *bounds*.Let's constrain the amplitude to its initial value of 5.3.<jupyter_code>gauss.amplitude.fixed = True
model = fitter(gauss, x, data)
print(model.parameters)<jupyter_output><empty_output><jupyter_text>Bounds can be set on a parameter either by using *min* and *max* or the *bounds* attribute.<jupyter_code>gauss.stddev.min = .1
gauss.stddev.max = .6
print(gauss.stddev.bounds)<jupyter_output><empty_output><jupyter_text>To show the mechanism of tieing (or linking) two parameters in the next example the stddev parameter is tied to the amplitude.<jupyter_code>def tie_stddev_ampl(model):
    return model.amplitude / 3.78

gauss.stddev.tied = tie_stddev_ampl
model_tied = fitter(gauss, x, data)
print(model_tied)
plt.plot(x, data)
plt.plot(x, model_tied(x))<jupyter_output><empty_output><jupyter_text>Support for parameter constraints varies between fitters. To see what is supported:<jupyter_code>print(fitting.LevMarLSQFitter.supported_constraints)<jupyter_output><empty_output><jupyter_text>astropy.modeling has several other fitters: *SimplexLSQFitter*, *SLSQPLSQFitter* and *LinearLSQFitter*.The last one can be used only with linear models and provides an exact solution.Create a Chebyshev model, evaluate it and add noise to the data.<jupyter_code>cheb1 = models.Chebyshev1D(degree=3, c0=1, c2=1, c3=1)

data = cheb1(x) + np.random.normal(0, 10, x.shape)
plt.plot(x, data, 'r')<jupyter_output><empty_output><jupyter_text>Fit a *Chebyshev1D* polynomial using the *LinearLSQFitter*.<jupyter_code>linfitter = fitting.LinearLSQFitter()
model = linfitter(cheb1, x, data)

plt.plot(x, data, label='data')
plt.plot(x, model(x), 'r', label='fitted Chebyshev1D')
plt.legend(loc=2)<jupyter_output><empty_output><jupyter_text>Exercise 1:Generate fake data```np.random.seed(0)x = np.linspace(-5., 5., 200)y = 3 * np.exp(-0.5 * (x - 1.3)**2 / 0.8**2)y += np.random.normal(0., 0.2, x.shape)```- Fit the data with a Trapezoid1D model.- Fit a Gaussian1D model to it.- Display the results.  Model Sets There are cases when it's useful to describe many models of the same type but with different parameter values.This could be done by passing *n_models* to the model with an integer value indicating the number of models.This is especially useful in the context of simultaneously fitting many linear models using the *LinearLSQFittter*.Evaluation of sets of models works for all models while fitting of sets of models is currently supported only for linear models.<jupyter_code>poly = models.Polynomial1D(degree=2, c0=0.5, c1=0, c2=1) # one model

poly10 = models.Polynomial1D(degree=2, c0=0.5*np.ones(10), c2=np.ones(10), n_models=10) # 10 models
print(poly10)
poly10real = models.Polynomial1D(degree=2, c0=0.5*np.ones(10) * np.random.normal(0, .1, 10),
                             c2=np.ones(10) * np.random.normal(1, .1, 10), n_models=10)
x = np.linspace(-1, 1, 21)
data = poly10real(x, model_set_axis=False)
print(poly10real)

fitpoly = linfitter(poly10, x, data)
print(fitpoly)
fitdata = fitpoly(x, model_set_axis=False)

for i in fitdata:
    plt.plot(x, i)<jupyter_output><empty_output><jupyter_text>Compound models astropy.modeling supports model combination using arithmetic operators and the specially defined **join (&)** and **composition ( | )** operators.<jupyter_code>g1 = models.Gaussian1D(2, 0, 0.2, name='g1') # models have names
g2 = models.Gaussian1D(3, 4, 1.5, name='g2')
p1 = models.Polynomial1D(1, c0=1, c1=0.2, name='poly')
model = g1 + g2 + p1<jupyter_output><empty_output><jupyter_text>Evaluating *model* is equivalent to evluating *g1(x) + g2(x) + p1(x)*.<jupyter_code>model(2.5) == g1(2.5) + g2(2.5) + p1(2.5)<jupyter_output><empty_output><jupyter_text>The other arithmetic operators work similarly. The only requirement is that all models accept the same input.Compound models support the *param_names* and *parameters* attributes the same way simple models do.<jupyter_code>print(model.param_names)
print(model.parameters)<jupyter_output><empty_output><jupyter_text>Parameter assignment works too.<jupyter_code>model.amplitude_0 = 10
print(model.parameters)<jupyter_output><empty_output><jupyter_text>If child models have names they can be printed, otherwise *None* is printed.<jupyter_code>print(model.submodel_names)<jupyter_output><empty_output><jupyter_text>Compound models can be sliced or indexed using their names or numerical indices like python lists.<jupyter_code>print(model['g2'])
print(model['g2' : ])<jupyter_output><empty_output><jupyter_text>The composition operator, *|*, combines models serially by chaining them one after the other. The number of outputs of a model must match the number of inputs of the next one.<jupyter_code>p2_1 = models.Polynomial2D(degree=1, c0_0=2, c0_1=.1, c1_0=2, name='Poly_X')
rot = models.Rotation2D(angle=23.1, name='Rotation')<jupyter_output><empty_output><jupyter_text>Chaining these two models as below is an error because the polynomial has one output, while the rotation has two inputs.<jupyter_code>broken_model = p2_1 | rot<jupyter_output><empty_output><jupyter_text>The join operator, *&*, evaluates the child models on independent inputs and the results are concatenated.The number of inputs passed to the combined model must equal the total number of inputs of all models.<jupyter_code>p2_2 = models.Polynomial2D(degree=1, c0_0=1, c1_0=2, name='Poly_Y')
model = p2_1 & p2_2 | rot
#x = np.linspace(1, 10, 23)
#y = np.linspace(-10, 10, 23)
#plt.imshow(model(x, y, x, y))
x = 1.2
y = 1.4
print(model(x, y, x, y))<jupyter_output><empty_output><jupyter_text>The **Mapping** model takes a tuple of indices into the inputs and returns the corresponding inputs. It is useful for changing the order of inputs, dropping or adding inputs.<jupyter_code>model = models.Mapping((0, 1, 0, 1)) | p2_1 & p2_2 | rot
print(model(x, y))<jupyter_output><empty_output><jupyter_text>Compound models also support parameter constraints. Constraints are defined on parameters of the compound model, not parameters of the submodels. Example 2:- read a spectrum from a text file (data/sample.txt).- Using the rest wavelengths as initial values, fit a gaussian to the H beta and OIII lines.Use the rest wavelengths as initial values for the locaiton of the lines.```Hbeta = 4862.721Halpha = 6564.614OIII_1 = 4958.911OIII_2 = 5008.239Na = 6549.86Nb = 6585.27Sa = 6718.29Sb = 6732.68```  Inverse of a model Models have their analytical inverse already defined. It is also possible to assign a "custom_inverse" by assigning a model to the *inverse* attribute.<jupyter_code>print(rot.inverse)
print(p2_2.inverse)
p2_2.inverse = p2_1
print(p2_2.inverse)<jupyter_output><empty_output><jupyter_text>Creating a new model Quite a few models are already defined in modeling.<jupyter_code>models.<jupyter_output><empty_output><jupyter_text>In most cases a new model can be easily defined following an existing model as an example.However, there's also a decorator, which works with user defined functions and turns them onto models.<jupyter_code>from astropy.modeling.models import custom_model

@custom_model
def sine_model(x, amplitude=1, frequency=1):
    return amplitude * np.sin(2 * np.pi * frequency * x)

model = sine_model(amplitude=3, frequency=2.1) # initialize the model
print(model(0.25))<jupyter_output><empty_output><jupyter_text>To supply also a derivative, *custom_model* can be used as a function.<jupyter_code>def sine_model(x, amplitude=1, frequency=1):
    return amplitude * np.sin(2 * np.pi * frequency * x)

def sine_deriv(x, amplitude=1, frequency=1):
    return 2 * np.pi * amplitude * np.cos(2 * np.pi * frequency *x)

SineModel = custom_model(sine_model, fit_deriv=sine_deriv) # create the class
model = SineModel(3, 2.1)# and initialize the model
print(model(0.25))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11117836965294592, 'fraction_numerical': 0.031376109765940276, 'mean_word_length': 6.538403041825095, 'pattern_counts': {'":': 0, '<': 125, '<?xml version=': 0, '>': 125, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '239564', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3602, 'n_tokens_neox': 3393, 'n_words': 929, 'path': 'AAS2016/notebooks/aas2016_Modeling.ipynb', 'repo_name': 'astropy/astropy-workshops', 'size': '9912'}