Issue with LogL values changing for different parameter values using JointLikelihood

[ghost]

I tried to develop a morphology lookup table very similar to template_model.py. In my case, I adapted the template model to read in FITS files to build a data frame of flux values that is multi-index by morphology parameters and information of energy, RA, and Dec. I first interpolate over the morphology parameters and then over energy, RA, and dec. The full program is here:

My problem comes to running the fit for my model. When I try to make a fit using JointLikelihood, the parameters seem to be sent to the fitter, but the logL values do not change

The first interpolation is done in this function

def _prepare_interpolators(self):

        #Figure out the shape of the data matrices
        #para_shape = [x.shape[0] for x in list(self._parameters_grids.values())]
        para_shape = map(lambda x: len(x), list(self._parameters_grids.values()))

        #interpolate over the parameters
        self._interpolators = []

        for e in self._E:
            for b in self._B:
                for l in self._L:

                    this_data = np.array(self._data_frame[tuple([e, l, b])].values.reshape(*para_shape), dtype=float)

                    #ti = GridInterpolator(self._parameters_grids.values(), this_data, bounds_error=False, fill_value=0.0 )
                    ti = GridInterpolator(tuple(map(lambda x: np.array(x), list(self._parameters_grids.values()))),
                                        this_data, bounds_error=False, fill_value=0.0)
                    self._interpolators.append(ti)

the parameter values are to evaluated to give an interpolated map and be used for interpolation over energy, RA, and dec

    def _interpolate(self, parameter_values):

        #create an interpolated map
        self._interpolated_map = np.array(map(lambda j: self._interpolators[j](np.atleast_1d(parameter_values)),
                                        range(len(self._interpolators))))

        #figure out the shape of the map
        map_shape = map(lambda x: len(x), self._map_parameters_grids.values())

        #interpolate over map energy, ra, and dec
        self._interpolator = GridInterpolator((self._E, self._L, self._B),
                                            self._interpolated_map.reshape(*map_shape), bounds_error=False, fill_value=0.0)

the function _interpolate gets called inside the evaluate function, and the self._interpolator should evalute the different values of energy, ra, and dec

    def evaluate(self, x, y, z, K, hash, lon0, lat0, *args):

        lon = x
        lat = y

        angsep = angular_distance_fast(lon0, lat0, lon, lat)

        #transform energy from keV to MeV
        #galprop likes MeV, 3ML likes keV
        convert_val = np.log10((u.MeV.to('keV')/u.keV).value)

        #if only one energy is passed, make sure we can iterate just once
        if not isinstance(z, np.ndarray):        
            z = np.array(z)

        energy = np.log10(z) - convert_val

        #Interpolated values can be cached-
        if self._interpmap is None:

            if lon.size != lat.size:
                raise AttributeError("Lons and lats should have the same size!")

            print("Building interpolation map....")

            #interpolate over  the spatial parameters
            print(args)
            self._interpolate(args)

            #create the array for the interpolated values
            f_new = np.zeros([energy.size, lat.size]) #Nicola : The image is a one dimensional array
            num = lat.size

            print("Interpolating over energy, RA, and Dec....")

            #evaluate interpolation function
            for i,e in enumerate(energy):

                engs = np.repeat(e, num)
                slice_points = zip(engs, lon, lat)
                interpolated_slice = self._interpolator(slice_points)

                f_new[i] = interpolated_slice

            print("Finished building interpolation map!")

            assert np.all(np.isfinite(f_new)), "some interpolated values are wrong"

            self._interpmap = f_new

        A = np.multiply(K, self._interpmap/(10**convert_val))
        #A = np.multiply(K, self._interpmap)

        return A.T

here is an example of the problem using the diffusion coefficient as a fitting parameter

try:

    jl.set_minimizer("minuit")

    param_df, like_df = jl.fit(quiet=False, compute_covariance=True)

except:

    jl.set_minimizer("ROOT")

    param_df, like_df = jl.fit(quiet=False, compute_covariance=True)

trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.5187e+26 -> logL = -11079569737.757 trial values: 3.3807e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.5187e+26 -> logL = -11079569737.757 trial values: 3.3807e+26 -> logL = -11079569737.757 trial values: 4.0881e+26 -> logL = -11079569737.757 trial values: 2.7351e+26 -> logL = -11079569737.757 trial values: 4.3404e+26 -> logL = -11079569737.757 trial values: 2.395e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.4501e+26 -> logL = -11079569737.757 trial values: 3.45e+26 -> logL = -11079569737.757 trial values: 3.4503e+26 -> logL = -11079569737.757 trial values: 3.4497e+26 -> logL = -11079569737.757 trial values: 3.453e+26 -> logL = -11079569737.757 trial values: 3.4471e+26 -> logL = -11079569737.757 trial values: 3.4798e+26 -> logL = -11079569737.757

[henrikef]

Hi Ramiro, I think it would be helpful if you included the full function definition here.

[henrikef]

Btw, were you able to verify in the mean time that your function does what you want it to do? I would start with that before trying to fit anything

First step could be to just verify the function values (define a source with this morphology, set the parameters, print the values at a few different energies/positions, change the parameter(s), print again etc).

Next step could be to make model maps for a few different values of your parameter(s).