tingyuansen/The_Payne

ValueError when running fit_normalized_spectrum_single_star_model

Closed this issue · 4 comments

Ompha commented

Hi Yuan-sen,
I am a grad student at JHU. And I am trying to extract stellar labels from my IRTF SpeX spectra for the wavelength range of 0.5 to 2.5 microns, with R~2000.

When I run fit_normalized_spectrum_single_star_model, as shown in the tutorial example, I encountered the ValueError: fp and xp are not of the same length..

A couple of things that I have tried to debug:

I've made sure that norm_spec, spec_err, wavelength and mask all have the same array length.

I am only fitting my spectrum in the wavelength range of the APOGEE spectra (15168 - 16936 Angstrom). The wavelength range is the same as the range in the tutorial, and I roughly have 600 data points in that range.

Additionally, I am using conda 4.9.2, in which I have scipy 1.6.1 and astropy of 4.2.
Thank you for your help. Below I attach the detailed error message:


ValueError                                Traceback (most recent call last)
<ipython-input-31-b93152322912> in <module>
     17                                 wavelength = wavelength[idx_wav],
     18                                 mask=mask[idx_wav],
---> 19                                 p0=None)
     20 

/anaconda3/lib/python3.7/site-packages/The_Payne-1.1-py3.7.egg/The_Payne/fitting.py in fit_normalized_spectrum_single_star_model(norm_spec, spec_err, NN_coeffs, wavelength, mask, p0)
     58     # run the optimizer
     59     popt, pcov = curve_fit(fit_func, xdata=[], ydata = norm_spec, sigma = spec_err, p0 = p0,
---> 60                 bounds = bounds, ftol = tol, xtol = tol, absolute_sigma = True, method = 'trf')
     61     pstd = np.sqrt(np.diag(pcov))
     62     model_spec = fit_func([], *popt)

/anaconda3/lib/python3.7/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    794 
    795         res = least_squares(func, p0, jac=jac, bounds=bounds, method=method,
--> 796                             **kwargs)
    797 
    798         if not res.success:

/anaconda3/lib/python3.7/site-packages/scipy/optimize/_lsq/least_squares.py in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    815         x0 = make_strictly_feasible(x0, lb, ub)
    816 
--> 817     f0 = fun_wrapped(x0)
    818 
    819     if f0.ndim != 1:

/anaconda3/lib/python3.7/site-packages/scipy/optimize/_lsq/least_squares.py in fun_wrapped(x)
    810 
    811     def fun_wrapped(x):
--> 812         return np.atleast_1d(fun(x, *args, **kwargs))
    813 
    814     if method == 'trf':

/anaconda3/lib/python3.7/site-packages/scipy/optimize/minpack.py in func_wrapped(params)
    485     elif transform.ndim == 1:
    486         def func_wrapped(params):
--> 487             return transform * (func(xdata, *params) - ydata)
    488     else:
    489         # Chisq = (y - yd)^T C^{-1} (y-yd)

/anaconda3/lib/python3.7/site-packages/The_Payne-1.1-py3.7.egg/The_Payne/fitting.py in fit_func(dummy_variable, *labels)
     42         norm_spec = spectral_model.get_spectrum_from_neural_net(scaled_labels = labels[:-1],
     43             NN_coeffs = NN_coeffs)
---> 44         norm_spec = utils.doppler_shift(wavelength, norm_spec, labels[-1])
     45         return norm_spec
     46 

/anaconda3/lib/python3.7/site-packages/The_Payne-1.1-py3.7.egg/The_Payne/utils.py in doppler_shift(wavelength, flux, dv)
     92     doppler_factor = np.sqrt((1 - dv/c)/(1 + dv/c))
     93     new_wavelength = wavelength * doppler_factor
---> 94     new_flux = np.interp(new_wavelength, wavelength, flux)
     95     return new_flux
     96 

<__array_function__ internals> in interp(*args, **kwargs)

/anaconda3/lib/python3.7/site-packages/numpy/lib/function_base.py in interp(x, xp, fp, left, right, period)
   1421         fp = np.concatenate((fp[-1:], fp, fp[0:1]))
   1422 
-> 1423     return interp_func(x, xp, fp, left, right)
   1424 
   1425 

ValueError: fp and xp are not of the same length.

Hey @Ompha. This is very interesting. Have you tried to simply do np.interp on your spectra and wavelength? You might also want to print the wavelength and flux shape before "new_flux = np.interp(new_wavelength, wavelength, flux)", and see what are their sizes. If this gets too tricky, I will be more than happy to iterate this problem with you further through email -- e.g., if you send me the spectrum, wavelegth, etc, that you are trying to fit.

Ompha commented

Hi @tingyuansen , thank you very much for the debug solution! I tried it out and figured out a workaround. I am able to get the code to work by interpolating the input flux data to have the identical wavelength data points as the tutorial.

A quick question: for the input labels, what do Vmacro [km/s] and Vturb [km/s] stand for? Thank you!

Hi Ompha,

That would be macroturbulence and microturbulence.

Hope this is helpful.

Ompha commented

Thank you!