def formula(menage, period):
        depcom = menage('depcom', period.first_month)

        departements_idf = [b'75', b'77', b'78', b'91', b'92', b'93', b'94', b'95']
        in_idf = sum([startswith(depcom, departement_idf) for departement_idf in departements_idf])

        rfr_declarants_principaux_du_menage = menage.members.has_role(FoyerFiscal.DECLARANT_PRINCIPAL) * menage.members.foyer_fiscal('rfr', period.n_2)
        rfr = menage.sum(rfr_declarants_principaux_du_menage)

        nb_members = menage.nb_persons()

        bareme_idf = select(
            [nb_members == 1, nb_members == 2, nb_members == 3, nb_members == 4, nb_members >= 5],
            [select([rfr <= 19875, rfr <= 24194], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 29171, rfr <= 35510], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 35032, rfr <= 42648], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 40905, rfr <= 49799], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 46798 + ((nb_members - 5) * 5882), rfr <= 56970 + ((nb_members - 5) * 7162)], [2, 1], 0)])

        bareme_out = select(
            [nb_members == 1, nb_members == 2, nb_members == 3, nb_members == 4, nb_members >= 5],
            [select([rfr <= 14360, rfr <= 18409], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 21001, rfr <= 26923], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 25257, rfr <= 32377], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 29506, rfr <= 37826], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier),
             select([rfr <= 33774 + ((nb_members - 5) * 4257), rfr <= 43297 + ((nb_members - 5) * 5454)], [TypesEligibiliteANAH.tres_modeste, TypesEligibiliteANAH.modestes], TypesEligibiliteANAH.a_verifier)])

        return where(in_idf, bareme_idf, bareme_out)

    def function(self, simulation, period):
        period = period
        depcom_holder = simulation.compute('depcom', period)

        depcom = self.cast_from_entity_to_roles(depcom_holder)
        depcom = self.filter_role(depcom, role = CHEF)
        return period, startswith(depcom, '973')

    def function(self, simulation, period):
        period = period.start.period(u'year').offset('first-of')
        depcom = simulation.calculate('depcom', period)

        return period, np.logical_or(startswith(depcom, '97'), startswith(depcom, '98'))

 def function(self, depcom):
     return np.logical_or(startswith(depcom, '97'), startswith(depcom, '98'))

 def function(self, simulation, period):
     depcom = simulation.calculate('depcom', period)
     return period, startswith(depcom, '973')

 def formula(menage, period, parameters):
     depcom = menage('depcom', period)
     return startswith(depcom, b'971')

    def _setup_spec_cube(self):

        self.spec_cube = True
        self.geometry = "spectral_cube"

        end = min(self.dimensionality+1,4)
        if self.events_data:
            ctypes = self.axis_names
        else:
            ctypes = np.array([self.primary_header["CTYPE%d" % (i)]
                               for i in range(1,end)])

        log_str = "Detected these axes: "+"%s "*len(ctypes)
        mylog.info(log_str % tuple([ctype for ctype in ctypes]))

        self.lat_axis = np.zeros((end-1), dtype="bool")
        for p in lat_prefixes:
            self.lat_axis += np_char.startswith(ctypes, p)
        self.lat_axis = np.where(self.lat_axis)[0][0]
        self.lat_name = ctypes[self.lat_axis].split("-")[0].lower()

        self.lon_axis = np.zeros((end-1), dtype="bool")
        for p in lon_prefixes:
            self.lon_axis += np_char.startswith(ctypes, p)
        self.lon_axis = np.where(self.lon_axis)[0][0]
        self.lon_name = ctypes[self.lon_axis].split("-")[0].lower()

        if self.lat_axis == self.lon_axis and self.lat_name == self.lon_name:
            self.lat_axis = 1
            self.lon_axis = 0
            self.lat_name = "Y"
            self.lon_name = "X"

        if self.wcs.naxis > 2:

            self.spec_axis = np.zeros((end-1), dtype="bool")
            for p in spec_names.keys():
                self.spec_axis += np_char.startswith(ctypes, p)
            self.spec_axis = np.where(self.spec_axis)[0][0]
            self.spec_name = spec_names[ctypes[self.spec_axis].split("-")[0][0]]

            self.wcs_2d = _astropy.pywcs.WCS(naxis=2)
            self.wcs_2d.wcs.crpix = self.wcs.wcs.crpix[[self.lon_axis, self.lat_axis]]
            self.wcs_2d.wcs.cdelt = self.wcs.wcs.cdelt[[self.lon_axis, self.lat_axis]]
            self.wcs_2d.wcs.crval = self.wcs.wcs.crval[[self.lon_axis, self.lat_axis]]
            self.wcs_2d.wcs.cunit = [str(self.wcs.wcs.cunit[self.lon_axis]),
                                     str(self.wcs.wcs.cunit[self.lat_axis])]
            self.wcs_2d.wcs.ctype = [self.wcs.wcs.ctype[self.lon_axis],
                                     self.wcs.wcs.ctype[self.lat_axis]]

            self._p0 = self.wcs.wcs.crpix[self.spec_axis]
            self._dz = self.wcs.wcs.cdelt[self.spec_axis]
            self._z0 = self.wcs.wcs.crval[self.spec_axis]
            self.spec_unit = str(self.wcs.wcs.cunit[self.spec_axis])

            if self.spectral_factor == "auto":
                self.spectral_factor = float(max(self.domain_dimensions[[self.lon_axis,
                                                                         self.lat_axis]]))
                self.spectral_factor /= self.domain_dimensions[self.spec_axis]
                mylog.info("Setting the spectral factor to %f" % (self.spectral_factor))
            Dz = self.domain_right_edge[self.spec_axis]-self.domain_left_edge[self.spec_axis]
            self.domain_right_edge[self.spec_axis] = self.domain_left_edge[self.spec_axis] + \
                                                     self.spectral_factor*Dz
            self._dz /= self.spectral_factor
            self._p0 = (self._p0-0.5)*self.spectral_factor + 0.5

        else:

            self.wcs_2d = self.wcs
            self.spec_axis = 2
            self.spec_name = "z"
            self.spec_unit = "code_length"

    def _parse_parameter_file(self):

        if self.parameter_filename.startswith("InMemory"):
            self.unique_identifier = time.time()
        else:
            self.unique_identifier = \
                int(os.stat(self.parameter_filename)[stat.ST_CTIME])

        # Determine dimensionality

        self.dimensionality = self.naxis
        self.geometry = "cartesian"

        # Sometimes a FITS file has a 4D datacube, in which case
        # we take the 4th axis and assume it consists of different fields.
        if self.dimensionality == 4: self.dimensionality = 3

        self.domain_dimensions = np.array(self.dims)[:self.dimensionality]
        if self.dimensionality == 2:
            self.domain_dimensions = np.append(self.domain_dimensions,
                                               [int(1)])

        self.domain_left_edge = np.array([0.5]*3)
        self.domain_right_edge = np.array([float(dim)+0.5 for dim in self.domain_dimensions])

        if self.dimensionality == 2:
            self.domain_left_edge[-1] = 0.5
            self.domain_right_edge[-1] = 1.5

        # Get the simulation time
        try:
            self.current_time = self.parameters["time"]
        except:
            mylog.warning("Cannot find time")
            self.current_time = 0.0
            pass

        # For now we'll ignore these
        self.periodicity = (False,)*3
        self.current_redshift = self.omega_lambda = self.omega_matter = \
            self.hubble_constant = self.cosmological_simulation = 0.0

        if self.dimensionality == 2 and self.z_axis_decomp:
            mylog.warning("You asked to decompose along the z-axis, but this is a 2D dataset. " +
                          "Ignoring.")
            self.z_axis_decomp = False

        if self.events_data: self.specified_parameters["nprocs"] = 1

        # If nprocs is None, do some automatic decomposition of the domain
        if self.specified_parameters["nprocs"] is None:
            if self.z_axis_decomp:
                nprocs = np.around(self.domain_dimensions[2]/8).astype("int")
            else:
                nprocs = np.around(np.prod(self.domain_dimensions)/32**self.dimensionality).astype("int")
            self.parameters["nprocs"] = max(min(nprocs, 512), 1)
        else:
            self.parameters["nprocs"] = self.specified_parameters["nprocs"]

        # Check to see if this data is in some kind of (Lat,Lon,Vel) format
        self.spec_cube = False
        self.wcs_2d = None
        x = 0
        for p in lon_prefixes+lat_prefixes+list(spec_names.keys()):
            y = np_char.startswith(self.axis_names[:self.dimensionality], p)
            x += np.any(y)
        if x == self.dimensionality:
            if self.axis_names == ['LINEAR','LINEAR']:
                self.wcs_2d = self.wcs
                self.lat_axis = 1
                self.lon_axis = 0
                self.lat_name = "Y"
                self.lon_name = "X"
            else:
                self._setup_spec_cube()