Fitting a warp models errors

[dan-adi]

Hello,
I tried to model a galaxy (MCG +07-16-015) captured using a Sloan g' filter with my telescope and a CCD camera(KAF 16200 chip, gain 0.6 e/ADU, RN 10e), but every warp model fails. Perhaps someone can point me to a mistake.
The fits image has been calibrated with darks,flats,bias. I attached the compressed fits file
Also, how do I compute the variance?

Here are my steps:

1. Make a func to substract the background

def subtract_background(image):
        
        sigma_clip = SigmaClip(sigma=3.0)
        bkg_estimator = MedianBackground()
        bkg = Background2D(image, (50, 50), filter_size=(3, 3),sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)
        background = bkg.background
        image = image - background
        return image

2. Open the image, crop, substract background, make some masks for the saturated stars

hdu =  fits.open('_g.fit')
x_center, y_center = 2129, 2247
width, height = 150, 150  # box, center [75,75]

# Calculate the slicing indices
x1, x2 = x_center - width // 2, x_center + width // 2
y1, y2 = y_center - height // 2, y_center + height // 2

# Crop the image and substract the bkg

crop_data = hdu[0].data[y1:y2, x1:x2]
substracted_bkg = subtract_background(crop_data)
final_image = np.array(substracted_bkg, dtype=np.float64)

#make masks for stars
mask = np.zeros(final_image.shape, dtype = bool)
mask[112:124,30:42] = True
mask[8:20,24:34] = True
mask[55:79,98:118] = True
mask[135:150,135:148] = True
mask[52:60,120:128] = True
mask[76:82,126:130] = True

# Create a target object
mytarget = ap.image.Target_Image(
    data=final_image,
    pixelscale=0.92, 
    zeropoint=26.7,
    mask=mask,
    variance='auto',)

fig1, ax1 = plt.subplots(figsize=(8, 8))
ap.plots.target_image(fig1, ax1, mytarget)
plt.show()

3. Make a simple model

model2 =ap.models.AstroPhot_Model(
    name="Model 2",
    model_type="sersic warp galaxy model",
    target=mytarget,  
    window=[[45, 97], [35, 120]],
    parameters = {"center": [75,75]} 
)

model2.initialize()

4. Fit error

result = ap.fit.LM(model2, verbose=1).fit()

print("Fit message:", result.message

Error

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[318], line 1
----> 1 result = ap.fit.LM(model2, verbose=1).fit()
      3 print("Fit message:", result.message)  # the fitter will return a message about its convergence

File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astrophot/fit/lm.py:165, in LM.__init__(self, model, initial_state, max_iter, relative_tolerance, ndf, **kwargs)
    155 def __init__(
    156     self,
    157     model,
   (...)
    162     **kwargs,
    163 ):
--> 165     super().__init__(
    166         model,
    167         initial_state,
    168         max_iter=max_iter,
    169         relative_tolerance=relative_tolerance,
    170         **kwargs,
    171     )
    172     # The forward model which computes the output image given input parameters
    173     self.forward = partial(model, as_representation=True)

File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astrophot/fit/base.py:68, in BaseOptimizer.__init__(self, model, initial_state, relative_tolerance, fit_window, **kwargs)
     66 if initial_state is None:
     67     self.model.initialize()
---> 68     initial_state = self.model.parameters.vector_representation()
     69 else:
     70     initial_state = torch.as_tensor(
     71         initial_state, dtype=AP_config.ap_dtype, device=AP_config.ap_device
     72     )

File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astrophot/param/parameter.py:214, in Parameter_Node.vector_representation(self)
    207 def vector_representation(self):
    208     """This returns a vector (see vector_values) with the representation
    209     for each leaf node. The representation is an alternative view
    210     of each value which is mapped into the (-inf, inf) range where
    211     optimization is more stable.
    212 
    213     """
--> 214     return self.vector_transform_val_to_rep(self.vector_values())

File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astrophot/param/parameter.py:413, in Parameter_Node.vector_transform_val_to_rep(self, val)
    410 reps = []
    411 for node in flat.values():
    412     reps.append(
--> 413         node.vector_transform_val_to_rep(
    414             val[mask[:loc].sum().int() : mask[: loc + node.size].sum().int()]
    415         )
    416     )
    417     loc += node.size
    418 if len(reps) > 0:

File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/astrophot/param/parameter.py:400, in Parameter_Node.vector_transform_val_to_rep(self, val)
    395         rep = val
    396     else:
    397         rep = boundaries(
    398             val,
    399             (
--> 400                 None if self.limits[0] is None else self.limits[0][self.mask],
    401                 None if self.limits[1] is None else self.limits[1][self.mask],
    402             ),
    403         )
    404     return rep
    406 mask = self.vector_mask()

IndexError: too many indices for tensor of dimension 0

_g.fit.zip

[ConnorStoneAstro]

@dan-adi looks like you've done everything right here. There is a bug in how the parameter limits are treated for parameters like the warp that arent scalars. Not sure how this got through my tests, but its an easy fix. I'll update AstroPhot with a patch today :)

[dan-adi]

Thank you Connor. What were the chances :). For the record I'm just an amateur astronomer. Got a little bored just taking pretty pictures with the telescope, so I thought I should learn some photometry since I know a little bit of python.
Thank you very much for helping! I'll update later this evening

[ConnorStoneAstro]

Also @dan-adi about computing the variance. This is often passed from the data reduction pipeline. In a single frame which is converted to units of counts is actually its own variance map (since for a Poisson distribution the variance and mean are the same). When a data reduction pipeline corrects for vignetting, darks, flats, hot pixels, and does drizzle/stacking, it usually tracks the variance through all those processes and returns a "Weight map". You can actually pass a weight map to an AstroPhot Target_Image objects just as weight=weight_map, in the same way you would pass variance=variance_map,. I see you're using variance="auto", right now, this is a good start for testing, but you'll want the true weight map to get the errors right on the final values. I hope that helps!

[dan-adi]

I use PixInsight for stacking and correcting with flats, darks, etc. I'll check and see if the software generates the Weight map.

[ConnorStoneAstro]

Ah, pixinsight might not compute a weight map, thats usually the domain of astronomer focused pipelines. Maybe check out SWARP (https://github.com/astromatic/swarp and https://www.astromatic.net/software/swarp/) to learn more!

[dan-adi]

Apparently Pixinsight makes a weight map, but it's something different. If I used their generated map, I got enormous Chi^2/Dof values, in the millions, without converging to a solution.
So after a bit of searching I made a map from the background RMS image like this:

# Load a FITS file
hdu = fits.open('_g.fit')
data = hdu[0].data

# Estimate the background
bkg_estimator = MedianBackground()
bkg = Background2D(data, (50, 50), filter_size=(3, 3), bkg_estimator=bkg_estimator)

# Background image
background = bkg.background

# Estimate the background RMS
rms_estimator = MADStdBackgroundRMS()
bkg_rms = Background2D(data, (50, 50), filter_size=(3, 3), bkg_estimator=bkg_estimator, bkgrms_estimator=rms_estimator)

# Background RMS image
background_rms = bkg_rms.background_rms

# Calculate the variance from the RMS (variance = RMS^2)
variance = background_rms**2

Using this map and a nuke model it worked much better:

#open images
image =  fits.open('_g_drizzle_r_ABE.fit')
varianceMap = fits.open('variance_from_background_rms.fits')

x_center, y_center = 2129, 2247
width, height = 201, 201  # box, center [75,75]

# Calculate the slicing indices
x1, x2 = x_center - width // 2, x_center + width // 2
y1, y2 = y_center - height // 2, y_center + height // 2

# Crop the image and the weight map image
crop_data = image[0].data[y1:y2, x1:x2]
varianceMapCrop = varianceMap[0].data[y1:y2, x1:x2]


final_image = np.array(crop_data, dtype=np.float64)

mask = np.zeros(final_image.shape, dtype = bool)
mask[75:85,5:18] = True
mask[26:45,45:65] = True
mask[80:112,124:145] = True

# Create a target object
mytarget = ap.image.Target_Image(
    data=final_image,
    pixelscale=0.92, 
    zeropoint=26.856,
    mask=mask,
    variance = varianceMapCrop
)

model =ap.models.AstroPhot_Model(
    name="Model 1",
    model_type="nuker galaxy model",
    target=mytarget,
    parameters = {"center":[91,94],
                  
                 },
    window=[[60, 140], [50, 140]],
)

model.initialize()

Result

result = ap.fit.LM(model, verbose=1).fit()
print("Fit message:", result.message)

Chi^2/DoF: 120.65020806472212, L: 1.0
Chi^2/DoF: 26.651451561484322, L: 0.1111111111111111
Chi^2/DoF: 13.518538907422496, L: 0.012345679012345678
Chi^2/DoF: 10.59589260468895, L: 0.0013717421124828531
Chi^2/DoF: 10.477631380487722, L: 0.0016765736930345982
Chi^2/DoF: 10.40660630524753, L: 0.0020491456248200642
Chi^2/DoF: 10.34659758289024, L: 0.0025045113192245227
Chi^2/DoF: 10.32591635505707, L: 0.0030610693901633058
Chi^2/DoF: 10.158265707382416, L: 0.0037413070324218184
Chi^2/DoF: 10.11480245461045, L: 0.004572708595182222
Chi^2/DoF: 10.024429537684666, L: 0.0005080787327980246
Chi^2/DoF: 9.999117661169763, L: 0.0006209851178642523
Chi^2/DoF: 9.995032680270635, L: 0.000758981810722975
Chi^2/DoF: 9.992223587707235, L: 0.11224497667469774
Chi^2/DoF: 9.986625361833337, L: 0.0013857404527740462
Chi^2/DoF: 9.985328253982487, L: 0.001693682775612723
Chi^2/DoF: 9.981969695349255, L: 3.1869629471808213e-09
Chi^2/DoF: 9.979900141518465, L: 0.07590597901623809
Chi^2/DoF: 9.97889607161929, L: 1.4283049109165097e-07
Chi^2/DoF: 9.978714213094124, L: 1.7633393961932218e-09
Chi^2/DoF: 9.97865368244489, L: 1e-09
Chi^2/DoF: 9.978607691672263, L: 0.0021652412222222224
Chi^2/DoF: 9.978516607530377, L: 4.0742833583073624e-09
Final Chi^2/DoF: 9.978512451802006, L: 4.0742833583073624e-09. Converged: success
Fit message: success

Plot

I have to see if I can make it better

[ConnorStoneAstro]

Hi @dan-adi I have just pushed the patch to fix this bug, check out v0.16.4. Thanks for bringing it to my attention! Stuff like this makes the code way better as we root out all the bugs :)

[dan-adi]

It worked thank you.

Is there an easy way to retrieve the surface brightness at a corresponding radii?
I'm thinking of determining the surface brightness of the galaxy, besides the profile, so I can compare with SDSS catalogue

[ConnorStoneAstro]

Hi @dan-adi yay! glad it worked :) There are a couple ways to do this. For one, you can take the parameter values and then look up the nuker radial intensity equation, then you can write it out yourself. An easier way is to call the model.radial_model(R) where R is a PyTorch tensor of radius values (in arcsec); that will return intensities I, to get SB you will need to do zeropoint - 2.5*np.log10(I). You might need to multiply I by pixelscale^2 as well to get the right units