FitModel including disk_teff and disk_radius: TypeCheckError

[gabrielastro]

Following the great (=useful and clear) tutorial at https://species.readthedocs.io/en/latest/tutorials/fitting_model_spectra.html but modifying it a bit, I am doing:

fit = FitModel(object_name='beta Pic b', model='atmo-neq-weak',
  inc_phot=False, inc_spec=True, fit_corr=None, apply_weights=False,
  bounds={'teff': (1500., 1800.),
          'radius': (0.5, 2.),
          'GPI_Y': ((0.5, 1.5), None),
          'disk_teff': (100,1e4),
          'disk_radius': (0.1,1e3)} )

i.e., setting up an atmospheric model and adding a (very free!) blackbody. Then, I do (purposefully leaving the mass very free: 10±9 MJ, i.e., 1–19 MJ to 1 sigma):

fit.run_ultranest(tag='betapic', min_num_live_points=500,
    output='betapic/', prior={'mass': (10., 9.)})

or

fit.run_multinest(tag='betapic', n_live_points=500,
    output='betapic_2', prior={'mass': (10., 9.)})

but both lead to

TypeCheckError: the return value (float) is not an instance of numpy.float64

If using ultranest, the program stops at the error, whereas multinest (fortunately!) ignores it and prints it sixteen bazillion times (I counted them; doing Ctrl+C does not interrupt the run either…) but otherwise works.

The problem does not occur (the fitting runs as advertised) if not including disk_teff and disk_radius.

[tomasstolker]

Thanks a lot for opening this issue Gabriel!

Would you mind posting the backtrace of the error (i.e. including the file names and line numbers)?

Probably it will be an easy fix, since the type checking is too strict somewhere. It seems to have returned a regular float somewhere instead of a float produced by NumPy, so both should be supported.

[tomasstolker]

For aborting from multinest, it may work to use Ctrl+Z instead of Ctrl-C.

[gabrielastro]

With pleasure! By the way, with Crtl+Z one only suspends ipython, unfortunately, but still this is fortunate because after a number of these warnings, the chain sort of works (the posteriors are not smooth but some non-pathological distribution is obtained).

Here is the backtrace (left me know if you need more):

~$ ipython3 
Python 3.9.2 (default, Feb 28 2021, 17:03:44) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.20.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: from species import SpeciesInit
   ...: from species.data.database import Database
   ...: from species.fit.fit_model import FitModel
   ...: from species.read.read_model import ReadModel

In [2]: SpeciesInit()
   ...: database = Database()
==============
species v0.7.4
==============

Working folder: […]/Test
Configuration settings:
   - Database: […]/Test/species_database.hdf5
   - Data folder: […]/Test/data
   - Interpolation method: linear
   - Magnitude of Vega: 0.03
Creating species_database.hdf5... [DONE]

In [3]: database.add_object('beta Pic b', parallax=(50.93, 0.15), app_mag=None, flux_density=None,deredden=None, spectrum={'GPI_Y': ('betapicb_gpi_y.dat', None, 40.)})  # downloaded previously from Tomas' website
Adding object: beta Pic b
   - Parallax (mas) = 50.93 +/- 0.15
   - Spectrum:
      - Database tag: GPI_Y
      - Filename: betapicb_gpi_y.dat
      - Data shape: (29, 3)
      - Wavelength range (um): 0.98 - 1.13
      - Mean flux (W m-2 um-1): 5.40e-15
      - Mean error (W m-2 um-1): 1.87e-15
   - Spectral resolution:
      - GPI_Y: 40.0

In [4]: database.add_model(model='atmo-neq-weak', teff_range=(1500., 1800.))
Unpacking ATMO NEQ weak model spectra (276 MB)... [DONE]
Please cite Phillips et al. (2020) when using ATMO NEQ weak in a publication
Reference URL: https://ui.adsabs.harvard.edu/abs/2020A%26A...637A..38P/abstract
Wavelength range (um) = 0.2 - 6000
Spectral resolution = 10000
Teff range (K) = 1500.0 - 1800.0
Adding ATMO NEQ weak model spectra... [DONE]                                   
Grid points stored in the database:
   - Teff = [1500. 1600. 1700. 1800.]
   - log(g) = [2.5 3.  3.5 4.  4.5 5.  5.5]
Number of grid points per parameter:
   - teff: 4
   - logg: 7
Fix missing grid points with a linear interpolation:
Number of stored grid points: 28
Number of interpolated grid points: 0
Number of missing grid points: 0

In [5]: fit = FitModel(object_name='beta Pic b', model='atmo-neq-weak', inc_phot=False, inc_spec=True, fit_corr=None, apply_weights=False, bounds={'teff': (1500., 1800.), 'radius': (0.5, 2.), 'GPI_Y': ((0.5, 1.5), None), 'disk
   ...: _teff': (100,1e4), 'disk_radius': (0.1,1e3)} )
Interpolating GPI_Y... [DONE]
Interpolating GPI_Y...[…]/species/species/read/read_model.py:143: UserWarning: The 'blackbody' model spectra are not present in the database. Will try to add the model grid. If this does not work (e.g. currently without an internet connection) then please use the 'add_model' method of 'Database' to add the grid of spectra at a later moment.
  warnings.warn(
Unpacking blackbody model spectra (56 MB)... [DONE]
Wavelength range (um) = 0.1 - 5000
Spectral resolution = 1000
Teff range (K) = 10 - 20000
Adding blackbody model spectra... [DONE]                       
Grid points stored in the database:
   - Teff = [1.00e+01 2.00e+01 3.00e+01 4.00e+01 5.00e+01 6.00e+01 7.00e+01 8.00e+01
 9.00e+01 1.00e+02 1.50e+02 2.00e+02 2.50e+02 3.00e+02 3.50e+02 4.00e+02
 4.50e+02 5.00e+02 5.50e+02 6.00e+02 6.50e+02 7.00e+02 7.50e+02 8.00e+02
 8.50e+02 9.00e+02 9.50e+02 1.00e+03 1.05e+03 1.10e+03 1.15e+03 1.20e+03
 1.25e+03 1.30e+03 1.35e+03 1.40e+03 1.45e+03 1.50e+03 1.55e+03 1.60e+03
 1.65e+03 1.70e+03 1.75e+03 1.80e+03 1.85e+03 1.90e+03 1.95e+03 2.00e+03
 2.05e+03 2.10e+03 2.15e+03 2.20e+03 2.25e+03 2.30e+03 2.35e+03 2.40e+03
 2.45e+03 2.50e+03 2.55e+03 2.60e+03 2.65e+03 2.70e+03 2.75e+03 2.80e+03
 2.85e+03 2.90e+03 2.95e+03 3.00e+03 3.05e+03 3.10e+03 3.15e+03 3.20e+03
 3.25e+03 3.30e+03 3.35e+03 3.40e+03 3.45e+03 3.50e+03 3.55e+03 3.60e+03
 3.65e+03 3.70e+03 3.75e+03 3.80e+03 3.85e+03 3.90e+03 3.95e+03 4.00e+03
 4.05e+03 4.10e+03 4.15e+03 4.20e+03 4.25e+03 4.30e+03 4.35e+03 4.40e+03
 4.45e+03 4.50e+03 4.55e+03 4.60e+03 4.65e+03 4.70e+03 4.75e+03 4.80e+03
 4.85e+03 4.90e+03 4.95e+03 5.00e+03 5.10e+03 5.20e+03 5.30e+03 5.40e+03
 5.50e+03 5.60e+03 5.70e+03 5.80e+03 5.90e+03 6.00e+03 6.10e+03 6.20e+03
 6.30e+03 6.40e+03 6.50e+03 6.60e+03 6.70e+03 6.80e+03 6.90e+03 7.00e+03
 7.10e+03 7.20e+03 7.30e+03 7.40e+03 7.50e+03 7.60e+03 7.70e+03 7.80e+03
 7.90e+03 8.00e+03 8.10e+03 8.20e+03 8.30e+03 8.40e+03 8.50e+03 8.60e+03
 8.70e+03 8.80e+03 8.90e+03 9.00e+03 9.10e+03 9.20e+03 9.30e+03 9.40e+03
 9.50e+03 9.60e+03 9.70e+03 9.80e+03 9.90e+03 1.00e+04 1.01e+04 1.02e+04
 1.03e+04 1.04e+04 1.05e+04 1.06e+04 1.07e+04 1.08e+04 1.09e+04 1.10e+04
 1.11e+04 1.12e+04 1.13e+04 1.14e+04 1.15e+04 1.16e+04 1.17e+04 1.18e+04
 1.19e+04 1.20e+04 1.21e+04 1.22e+04 1.23e+04 1.24e+04 1.25e+04 1.26e+04
 1.27e+04 1.28e+04 1.29e+04 1.30e+04 1.31e+04 1.32e+04 1.33e+04 1.34e+04
 1.35e+04 1.36e+04 1.37e+04 1.38e+04 1.39e+04 1.40e+04 1.41e+04 1.42e+04
 1.43e+04 1.44e+04 1.45e+04 1.46e+04 1.47e+04 1.48e+04 1.49e+04 1.50e+04
 1.51e+04 1.52e+04 1.53e+04 1.54e+04 1.55e+04 1.56e+04 1.57e+04 1.58e+04
 1.59e+04 1.60e+04 1.61e+04 1.62e+04 1.63e+04 1.64e+04 1.65e+04 1.66e+04
 1.67e+04 1.68e+04 1.69e+04 1.70e+04 1.71e+04 1.72e+04 1.73e+04 1.74e+04
 1.75e+04 1.76e+04 1.77e+04 1.78e+04 1.79e+04 1.80e+04 1.81e+04 1.82e+04
 1.83e+04 1.84e+04 1.85e+04 1.86e+04 1.87e+04 1.88e+04 1.89e+04 1.90e+04
 1.91e+04 1.92e+04 1.93e+04 1.94e+04 1.95e+04 1.96e+04 1.97e+04 1.98e+04
 1.99e+04 2.00e+04]
Number of grid points per parameter:
   - teff: 258
[…]/species/species/util/data_util.py:379: RuntimeWarning: divide by zero encountered in log10
  flux = np.log10(flux)
Number of stored grid points: 0
Number of interpolated grid points: 0
Number of missing grid points: 0
 [DONE]
Fitting 7 parameters:
   - teff
   - logg
   - radius
   - parallax
   - disk_teff
   - disk_radius
   - scaling_GPI_Y
Prior boundaries:
   - teff = (1500.0, 1800.0)
   - radius = (0.5, 2.0)
   - disk_teff = (100, 10000.0)
   - disk_radius = (0.1, 1000.0)
   - logg = (2.5, 5.5)
   - scaling_GPI_Y = (0.5, 1.5)
Weights for the log-likelihood function:
   - GPI_Y = 1.00

In [6]: fit.run_ultranest(tag='betapic', min_num_live_points=500, output='betapic/', prior={'mass': (10., 9.)})
Running nested sampling with UltraNest...
Creating directory for new run betapic/run2
---------------------------------------------------------------------------
TypeCheckError                            Traceback (most recent call last)
<ipython-input-7-ca466f5e96fd> in <module>
----> 1 fit.run_ultranest(tag='betapic', min_num_live_points=500, output='betapic/', prior={'mass': (10., 9.)})

[…]/species/species/fit/fit_model.py in run_ultranest(self, tag, min_num_live_points, output, prior)
   2088             return self.lnlike_func(params, prior=prior)
   2089 
-> 2090         sampler = ultranest.ReactiveNestedSampler(
   2091             self.modelpar,
   2092             lnlike_ultranest,

~/.local/lib/python3.9/site-packages/ultranest/integrator.py in __init__(self, param_names, loglike, transform, derived_param_names, wrapped_params, resume, run_num, log_dir, num_test_samples, draw_multiple, num_bootstraps, vectorized, ndraw_min, ndraw_max, storage_backend, warmstart_max_tau)                                                                                                                                                               
   1194         self.ndraw_max = ndraw_max
   1195         self.build_tregion = transform is not None
-> 1196         if not self._check_likelihood_function(transform, loglike, num_test_samples):
   1197             assert self.log_to_disk
   1198             if resume_similar and self.log_to_disk:

~/.local/lib/python3.9/site-packages/ultranest/integrator.py in _check_likelihood_function(self, transform, loglike, num_test_samples)
   1256                 "Error in transform function: returned shape is %s, expected %s" % (
   1257                     np.shape(p), (num_test_samples, self.num_params)))
-> 1258             logl = loglike(p)
   1259             assert np.logical_and(u > 0, u < 1).all(), (
   1260                 "Error in transform function: u was modified!")

~/.local/lib/python3.9/site-packages/ultranest/utils.py in vectorized(args)
    132     def vectorized(args):
    133         """Vectorized version of function."""
--> 134         return np.asarray([function(arg) for arg in args])
    135 
    136     # give a user-friendly name to the vectorized version of the function

~/.local/lib/python3.9/site-packages/ultranest/utils.py in <listcomp>(.0)
    132     def vectorized(args):
    133         """Vectorized version of function."""
--> 134         return np.asarray([function(arg) for arg in args])
    135 
    136     # give a user-friendly name to the vectorized version of the function

[…]/species/species/fit/fit_model.py in lnlike_ultranest(params)
   2086             """
   2087 
-> 2088             return self.lnlike_func(params, prior=prior)
   2089 
   2090         sampler = ultranest.ReactiveNestedSampler(

~/.local/lib/python3.9/site-packages/typeguard/_functions.py in check_return_type(func_name, retval, annotation, memo)
    165 
    166     try:
--> 167         check_type_internal(retval, annotation, memo)
    168     except TypeCheckError as exc:
    169         # Allow NotImplemented if this is a binary magic method (__eq__() et al)

~/.local/lib/python3.9/site-packages/typeguard/_checkers.py in check_type_internal(value, annotation, memo)
    762     if isclass(origin_type):
    763         if not isinstance(value, origin_type):
--> 764             raise TypeCheckError(f"is not an instance of {qualified_name(origin_type)}")
    765     elif type(origin_type) is str:  # noqa: E721
    766         warnings.warn(

TypeCheckError: the return value (float) is not an instance of numpy.float64

I tried slightly different numbers and got the same problem.

On a different machine (with a slightly different python distribution, I guess), with the exact same commands, I get:

In [6]: fit.run_ultranest(tag='betapic', min_num_live_points=500, output='betapic/', prior={'mass': (10., 9.)})                                                                                       
Running nested sampling with UltraNest...
Creating directory for new run betapic/run4
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-5-ca466f5e96fd> in <module>
----> 1 fit.run_ultranest(tag='betapic', min_num_live_points=500, output='betapic/', prior={'mass': (10., 9.)})

~/.local/lib/python3.8/site-packages/typeguard/__init__.py in wrapper(*args, **kwargs)
   1031         memo = _CallMemo(python_func, _localns, args=args, kwargs=kwargs)
   1032         check_argument_types(memo)
-> 1033         retval = func(*args, **kwargs)
   1034         try:
   1035             check_return_type(retval, memo)

[…]/species/species/fit/fit_model.py in run_ultranest(self, tag, min_num_live_points, output, prior)
   2088             return self.lnlike_func(params, prior=prior)
   2089 
-> 2090         sampler = ultranest.ReactiveNestedSampler(
   2091             self.modelpar,
   2092             lnlike_ultranest,

~/.local/lib/python3.8/site-packages/ultranest/integrator.py in __init__(self, param_names, loglike, transform, derived_param_names, wrapped_params, resume, run_num, log_dir, num_test_samples, draw_multiple, num_bootstraps, vectorized, ndraw_min, ndraw_max, storage_backend, warmstart_max_tau)
   1187         self.ndraw_max = ndraw_max
   1188         self.build_tregion = transform is not None
-> 1189         if not self._check_likelihood_function(transform, loglike, num_test_samples):
   1190             assert self.log_to_disk
   1191             if resume_similar and self.log_to_disk:

~/.local/lib/python3.8/site-packages/ultranest/integrator.py in _check_likelihood_function(self, transform, loglike, num_test_samples)
   1254             assert np.shape(logl) == (num_test_samples,), (
   1255                 "Error in loglikelihood function: returned shape is %s, expected %s" % (np.shape(logl), (num_test_samples,)))
-> 1256             assert np.isfinite(logl).all(), (
   1257                 "Error in loglikelihood function: returned non-finite number: %s for input u=%s p=%s" % (logl, u, p))
   1258 

AssertionError: Error in loglikelihood function: returned non-finite number: [           -inf -45780.27996034] for input u=[[0.08830236 0.06385349 0.33097212 0.02408574 0.81923631 0.60554687
  0.20954194]
 [0.31167665 0.90636524 0.27903744 0.30300557 0.07012407 0.68003549
  0.78115951]] p=[[1.52649071e+03 2.69156047e+00 9.96458175e-01 5.06336221e+01
  8.21043949e+03 6.05586319e+02 7.09541940e-01]
 [1.59350300e+03 5.21909573e+00 9.18556156e-01 5.08526337e+01
  7.94228337e+02 6.80067485e+02 1.28115951e+00]]

This might be a more useful output. Something about an inf in the blackbody part?

[tomasstolker]

It seems that UltraNest can not handle inf. Other samplers typically do. I think this has been fixed in commit bc6a4d3a4129fb792a57c7ce0b804232e1d8c9ec but please give it a try!

[tomasstolker]

The first issue with the returned float type has also been fixed.

[gabrielastro]

Perfect, bedankt! Indeed it works.