Q3D / q3dfit

Python software for fitting integral field spectroscopic data
GNU General Public License v3.0
11 stars 3 forks source link

lines not being fit #15

Closed jerome-seebeck closed 1 year ago

jerome-seebeck commented 1 year ago

I was testing some of the example notebooks (PG1411 and J1652) and had the same issue in both. After the data was fit and saved as I attempted to reload it and plot the data with "q3do.plot_line(q3di,plotargs=argsplotline)" in cell [31] I get the result "plot_line: no lines to plot!". The notebook was not altered at all and I am using v1.1.1 installed via pip.

I can see at some point in the fitting procedure from the output that FIT STATUS changes from 1 to None but I am not sure what causes this. This is supported by the change in "dolinefit" flag from the q3di input to the q3do output from True to False.

Thanks so much for any help!

Update: I think what's happening is that when checkcomp is run on Halpha during one of the iterations it doesn't find any good components and thus leaves all of them at zero. Thus the fit gets rerun with ncomp = 0 (or an array of zeros, not exactly sure how the functions work). This is what causes the output line in the above file of FITLOOP: Repeating the fit of Halpha with 0 components.. Then when the fit is rerun for Halpha with no nonzero components it is killed by line 690 in fitspec.py: q3do.dolinefit = False (when len(nonzerocomp) <= 0).

I'm not sure if this is an intentional feature and there's a separate bug causing it for my versions of packages and the example notebook or if this is even the root cause but I thought it was important to note.

jerome-seebeck commented 1 year ago

Ultimately was an error with not having the proper Python environment set up I was using 3.10.9 and moved to 3.9.11

drupke commented 1 year ago

Jerome, thanks for following up. However, if the error appears in Python 3.10 that is still an issue. Users should be able to run any version of Python >~3.8. E.g., I am running q3dfit on 3.10.

Is the only difference in your new environment the lower Python version?

jerome-seebeck commented 1 year ago

It was the only difference. Although I am using Anaconda so that likely changed some other package versions. While single-core fitting was fixed with downgrading I had trouble with the multicore fitting which just runs indefinitely. Which specific python 3.10 version are you running? I'm having errors with 3.10.9.

It seems I'm having errors like the original one described in this thread when I use python 3.10.9. But when I downgrade to 3.9.11, I am able to fit the notebooks correctly with single core but multicore breaks down.

jerome-seebeck commented 1 year ago

An update:

I have q3dfit working on python 3.9.11 with the multicore. However, while this works for fitting of the example cube PG1411 it's not fitting nearly the number of spaxels as were fit in the example notebook. For most spaxels it has the same issue as I was having while using python 3.10.6 where it exits the line fitting after the message FITLOOP: Repeating the fit of Halpha with 0 components.. This results in a fairly empty cube which you can see in the attached image (top) of the Halpha maps. The points that are there seem fairly accurate to me though. I have also included the maps from the example notebook (bottom).

image

image

I think possibly this might be an issue of fitting tolerances. In the current jnb all of the linefit tolerances are set to 1e-15. Thus I set them back to what was listed as the default of 1e-8 and ran the fitting again. This fits a few more spaxels but still far from the number fit in the example jnb. If this is the issue and I just need to drop the tolerance even lower let me know and then q3dfit is working as expected but if not then the core issue of not having lines fit still persists. I would be curious what tolerances the example PG1411 jnb cube was fit with.

Regardless, it's odd to me that the same data for the spaxel (14, 11) in PG1411 is able to be fit in python 3.9.11 but not in 3.10.6 if this is a tolerance issue.

gabrieloliveira2428 commented 1 year ago

Hi, I'm using the 1.1.3rc2 version (from pip) and the line fitting issue still persists, with the error "plot_line: no lines to plot!". I managed to find a work-around for the problem, removing the lines with 0 components from the newncomp dict, and removing the same thing from the ncomp dict, right after "if somedata" and before all the code runs in the fitloop.py file. Don't know if I solved the issue, because now the fitting runs, but the second loop just copies the guess values from the first one. The continuum fits smoothly, so I think the problem is only associated with the line-emission fitting part.

lines continuum qso total

drupke commented 1 year ago

@gabrieloliveira2428, you will also need to make some by-hand changes to the notebook. I sent these over Slack to @jerome-seebeck, but you may not have seen them. Once things seem stable again, I'll package all of these changes into a patch release and we can hopefully move on! Will you make these changes and then try again? I assume the plots you are showing are for the default (14,11) spaxel?

Short answer is the specifications for the optimizer in lmfit. TLDR; make the following changes to your PG1411 notebook. I will include these in the next patch release. In the options to lmfit under line fitting:

#q3di.argslinefit['method'] = 'leastsq'
#q3di.argslinefit['iter_cb'] = 'per_iteration'
# As an example, to change the criteria for fit convergence from the defaults of 1.e-8 to 1.e-15:
q3di.argslinefit['ftol'] = 1.e-15
q3di.argslinefit['gtol'] = 1.e-15
q3di.argslinefit['xtol'] = 1.e-15
q3di.argslinefit['x_scale'] = 'jac'
q3di.argslinefit['tr_solver'] = 'lsmr'
# .. and the "suitable step length for the forward- difference approximation of the Jacobian. Normally the actual step length will be sqrt(epsfcn)*x"
# this is only for scipy.optimize.leastsq
#q3di.argslinefit['epsfcn'] = 1.e-15

In the options to lmfit under continuum fitting:

#q3di.argscontfit['method'] = 'leastsq'

argslmfit = dict()
argslmfit['ftol'] = 1.e-15
argslmfit['gtol'] = 1.e-15
argslmfit['xtol'] = 1.e-15
#argslmfit['epsfcn'] = 1.e-15

q3di.argscontfit['argslmfit'] = argslmfit

I also recommend a slight change to the normalization:

q3di.argsreadcube = {'wmapext': None,
                     'waveunit_in': 'Angstrom',
                     'fluxunit_in': 'erg/s/cm2/Angstrom',
                     'fluxnorm': 1e3}
cube = q3di.load_cube()

_Slightly longer answer: The problem for which we are optimizing is nonlinear and has some degeneracies. After playing around with the leastsq method in scipy.optimize while fitting the [SIV] data for F2M1106, I realized that one could use this with lmfit despite the fact that leastsq doesn't support boundary conditions, because lmfit has its own boundary condition constraints implemented. leastsq has an advantage over leastsquares (what we had been using) in that it's a more direct descendent of the MINPACK algorithm that's also implemented in the IDL version, MPFIT. It's also faster. So, long story short, I changed the default method in the PG1411 notebook but did not extensively test it beyond the (14,11) spaxel. I.e., I did not run the full cube. It sounds like leastsq is not in fact working as well for this case. In general, the exact optimization method and settings we choose may depend on the situation, which could be the ultimate lesson here. I also continue to be a little dissatisfied with lmfit for these complex problems, but maybe that's because it's trying to be a one-size-fits-all. I'm unfortunately not enough of a whiz to write our own optimization approach a la Cappellari. Maybe someday. In fact, Michelle Cappellari's latest ppxf paper is a bit dense but has a good summary of the challenges of doing these sorts of fits: https://arxiv.org/pdf/2208.14974.pdf

drupke commented 1 year ago

OK, I've patch incremented to v1.1.4 so hopefully this issue is resolved. @gabrieloliveira2428, @jerome-seebeck, please reopen if you continue to have this particular issue.