def compareOutputGeometry(self, orientedImageData, spacing, origin, directions):
   if orientedImageData is None:
     logging.error('Invalid input oriented image data')
     return False
   if (not isinstance(spacing, list) and not isinstance(spacing, tuple)) \
       or (not isinstance(origin, list) and not isinstance(origin, tuple)) \
       or not isinstance(directions, list):
     logging.error('Invalid baseline object types - need lists')
     return False
   if len(spacing) != 3 or len(origin) != 3 or len(directions) != 3 \
       or len(directions[0]) != 3 or len(directions[1]) != 3 or len(directions[2]) != 3:
     logging.error('Baseline lists need to contain 3 elements each, the directions 3 lists of 3')
     return False
   import numpy
   tolerance = 0.0001
   actualSpacing = orientedImageData.GetSpacing()
   actualOrigin = orientedImageData.GetOrigin()
   actualDirections = [[0]*3,[0]*3,[0]*3]
   orientedImageData.GetDirections(actualDirections)
   for i in [0,1,2]:
     if not numpy.isclose(spacing[i], actualSpacing[i], tolerance):
       logging.warning('Spacing discrepancy: ' + str(spacing) + ' != ' + str(actualSpacing))
       return False
     if not numpy.isclose(origin[i], actualOrigin[i], tolerance):
       logging.warning('Origin discrepancy: ' + str(origin) + ' != ' + str(actualOrigin))
       return False
     for j in [0,1,2]:
       if not numpy.isclose(directions[i][j], actualDirections[i][j], tolerance):
         logging.warning('Directions discrepancy: ' + str(directions) + ' != ' + str(actualDirections))
         return False
   return True

def test_point_setitem():
    p = Point()

    p[0] = 6.0
    assert p[0] == 6.0

    p[1] = 16.0
    p[1] += 600.0
    assert np.isclose(p[1], 616.0)

    p[2] = 111.0
    p[2] *= 12.0
    p[2] /= 2
    assert np.isclose(p[2], 666.0)

    with pytest.raises(IndexError):
        p[3] = 6666.0

    p[:] = (0, 0, 0)
    assert np.all(p[:] == 0)

    p[:] = (1, 2, 3)
    assert np.all(p[:] == (1, 2, 3))

    p[:] += np.array((1, 2, 3))
    assert np.all(p[:] == (2, 4, 6))

    p[:] /= 2
    assert np.all(p[:] == (1, 2, 3))

    p[:] *= np.array((2., 2., 2.))
    assert np.all(p[:] == (2, 4, 6))

def test_init():
  alm = {
    'gps': {
        'a': 8,
        'af0': 9,
        'af1': 10,
        'argp': 7,
        'ecc': 8,
        'inc': 3,
        'ma': 8,
        'raaw': 6,
        'rora': 4,
        'toa': 2,
        'week': 11
    },
    'healthy': 1,
    'sid': {
        'band': 0,
        'constellation': 0,
        'sat': 1
    },
    'valid': 1,
  }

  satAlmanac = swiftnav.almanac.Almanac(**alm)
  assert np.isclose(alm['gps']['a'], satAlmanac.gps['a'])
  assert np.isclose(alm['gps']['ecc'], satAlmanac.gps['ecc'])

def _isclosemod(a, b, atol=1E-5, mod=2*pi):
    """
    Return whether two numbers (or arrays) are within atol of each other
    in the modulo space determined by mod.
    """
    return (isclose(a%mod, b%mod, atol=atol) 
            or isclose((a+atol)%mod, (b+atol)%mod, atol=atol))

 def test_fill_empty(self):
     empty = HistogramData(bins = np.linspace(0, 10, 11))
     empty.fill(5)
     assert empty.total == 1
     empty.fill(5, 0.4)
     assert np.isclose(empty.total_weight, 1.4)
     assert np.isclose(empty.total, 2)

    def test_monopole_fluxpoints(self):
        """Tests monopole flux points."""

        field = ElectricField([PointCharge(2, [0, 0])])
        circle = GaussianCircle([0, 0], 10)

        fluxpoints = circle.fluxpoints(field, 4)
        self.assertEqual(len(fluxpoints), 4)
        self.assertTrue(isclose(fluxpoints,
                                [[10, 0], [0, 10], [-10, 0], [0, -10]]).all())

        fluxpoints = circle.fluxpoints(field, 14)
        self.assertEqual(len(fluxpoints), 14)
        self.assertTrue(isclose(fluxpoints[0], [10, 0]).all())
        self.assertTrue(isclose(fluxpoints[7], [-10, 0]).all())

        x1 = fluxpoints[1:7]
        x2 = fluxpoints[-1:7:-1]
        x2[:, 1] = fabs(x2[:, 1])
        self.assertTrue(isclose(x1, x2).all())

        x1 = append(fluxpoints[-3:], fluxpoints[:4], axis=0)
        x2 = fluxpoints[-4:3:-1]
        x2[:, 0] = fabs(x2[:, 0])
        self.assertEqual(len(x1), len(x2))
        self.assertTrue(isclose(x1, x2).all())

    def test_moma_sanity(self, model):
        """Test optimization criterion and optimality."""
        try:
            solver = sutil.get_solver_name(qp=True)
            model.solver = solver
        except sutil.SolverNotFound:
            pytest.skip("no qp support")

        sol = model.optimize()
        with model:
            model.reactions.PFK.knock_out()
            knock_sol = model.optimize()
            ssq = (knock_sol.fluxes - sol.fluxes).pow(2).sum()

        with model:
            add_moma(model)
            model.reactions.PFK.knock_out()
            moma_sol = model.optimize()
            moma_ssq = (moma_sol.fluxes - sol.fluxes).pow(2).sum()

        # Use normal FBA as reference solution.
        with model:
            add_moma(model, solution=sol)
            model.reactions.PFK.knock_out()
            moma_ref_sol = model.optimize()
            moma_ref_ssq = (moma_ref_sol.fluxes - sol.fluxes).pow(2).sum()

        assert numpy.isclose(moma_sol.objective_value, moma_ssq)
        assert moma_ssq < ssq
        assert numpy.isclose(moma_sol.objective_value,
                             moma_ref_sol.objective_value)
        assert numpy.isclose(moma_ssq, moma_ref_ssq)

    def get_vanadium(self, detector_mask, m1, colltrans, exp, indir):
        """
        This function returns either (vanadium_count, vanadium_monitor, None) or
        (None, None, vcorr) depending what type of file is provided by getProperty("Vanadium")
        """
        if not self.getProperty("Normalise").value:
            return None, None, np.ones(44)[detector_mask]

        vanadium_filename = self.getProperty("Vanadium").value
        if vanadium_filename:
            if vanadium_filename.split('.')[-1] == 'dat':
                vanadium = np.genfromtxt(vanadium_filename)
                vanadium_count = vanadium[:, 5:49].sum(axis=0)[detector_mask]
                vanadium_monitor = vanadium[:, 3].sum()
                logger.notice("Using vanadium data file: {}".format(vanadium_filename))
                return vanadium_count, vanadium_monitor, None
            else:
                vcorr_filename = vanadium_filename
        else: # Find adjacent vcorr file
            # m1 is the monochromator angle
            # m1 = 0 -> Ge 115, 1.54A
            # m1 = 9.45 -> Ge 113, 2.41A
            # colltrans is the collimator position, whether in or out of the beam
            # colltrans = 0 -> IN
            # colltrans = +/-80 -> OUT
            vcorr_filename = 'HB2A_{}__Ge_{}_{}_vcorr.txt'.format(exp,
                                                                  115 if np.isclose(m1, 0, atol=0.1) else 113,
                                                                  "IN" if np.isclose(colltrans, 0, atol=0.1) else "OUT")
        vcorr_filename = os.path.join(indir, vcorr_filename)
        logger.notice("Using vcorr file: {}".format(vcorr_filename))
        if not os.path.isfile(vcorr_filename):
            raise RuntimeError("Vanadium file {} does not exist".format(vcorr_filename))

        return None, None, np.genfromtxt(vcorr_filename)[detector_mask]

    def test(self):
        global filename
        if filename is None: filename = corsika.example_data_dir + '/DAT000002-32'
        assert os.path.exists(filename)

        raw = corsika.RawStream(filename)
        block = corsika.Block()

        # get the run header, event header and first particle block
        raw.get_next_block(block)
        assert block.ID == 'RUNH'
        raw.get_next_block(block)
        assert block.ID == 'EVTH'
        raw.get_next_block(block)
        assert numpy.all(numpy.isclose(reference, block.data.reshape((-1,7))))
        n_blocks = 3
        while raw.get_next_block(block):
            n_blocks += 1

        # check total number of blocks
        assert n_blocks == 4725
        raw.close()

        # check particle iteration
        raw = corsika.RawStream(filename)
        raw.get_next_block(block)
        raw.get_next_block(block)
        particles = raw.particles() # this works because it is positioned right before the particle block
        for i,p in enumerate(raw.particles()):
            if i >= reference.shape[0]: break
            assert numpy.all(numpy.isclose([p.px, p.py, p.pz, p.x, p.y, p.t_or_z], reference[i][1:], 3))

def likelihood_check(obs_distns,trans_matrix,init_distn,data,target_val):
    for cls in [m.HMMPython, m.HMM]:
        hmm = cls(alpha=6.,init_state_concentration=1, # placeholders
                obs_distns=obs_distns)
        hmm.trans_distn.trans_matrix = trans_matrix
        hmm.init_state_distn.weights = init_distn
        hmm.add_data(data)

        # test default log_likelihood method

        assert np.isclose(target_val, hmm.log_likelihood())

        # manual tests of the several message passing methods

        states = hmm.states_list[-1]

        states.clear_caches()
        states.messages_forwards_normalized()
        assert np.isclose(target_val,states._normalizer)

        states.clear_caches()
        states.messages_forwards_log()
        assert np.isinf(target_val) or np.isclose(target_val,states._normalizer)

        states.clear_caches()
        states.messages_backwards_log()
        assert np.isinf(target_val) or np.isclose(target_val,states._normalizer)

        # test held-out vs in-model

        assert np.isclose(target_val, hmm.log_likelihood(data))

def test_colorbar_renorm():
    x, y = np.ogrid[-4:4:31j, -4:4:31j]
    z = 120000*np.exp(-x**2 - y**2)

    fig, ax = plt.subplots()
    im = ax.imshow(z)
    cbar = fig.colorbar(im)
    assert np.allclose(cbar.ax.yaxis.get_majorticklocs(),
                       np.arange(0, 120000.1, 15000))

    cbar.set_ticks([1, 2, 3])
    assert isinstance(cbar.locator, FixedLocator)

    norm = LogNorm(z.min(), z.max())
    im.set_norm(norm)
    assert isinstance(cbar.locator, _ColorbarLogLocator)
    assert np.allclose(cbar.ax.yaxis.get_majorticklocs(),
                       np.logspace(-8, 5, 14))
    # note that set_norm removes the FixedLocator...
    assert np.isclose(cbar.vmin, z.min())
    cbar.set_ticks([1, 2, 3])
    assert isinstance(cbar.locator, FixedLocator)
    assert np.allclose(cbar.ax.yaxis.get_majorticklocs(),
                       [1.0, 2.0, 3.0])

    norm = LogNorm(z.min() * 1000, z.max() * 1000)
    im.set_norm(norm)
    assert np.isclose(cbar.vmin, z.min() * 1000)
    assert np.isclose(cbar.vmax, z.max() * 1000)

def test_gen_design():
    from brainiak.utils.utils import gen_design
    import numpy as np
    import os.path
    files = {'FSL1': 'example_stimtime_1_FSL.txt',
             'FSL2': 'example_stimtime_2_FSL.txt',
             'AFNI1': 'example_stimtime_1_AFNI.txt'}
    for key in files.keys():
        files[key] = os.path.join(os.path.dirname(__file__), files[key])
    design1 = gen_design(stimtime_files=files['FSL1'], scan_duration=[48, 20],
                         TR=2, style='FSL')
    assert design1.shape == (34, 1), 'Returned design matrix has wrong shape'
    assert design1[24] == 0, (
        "gen_design should generated design matrix for each run separately "
        "and concatenate them.")
    design2 = gen_design(stimtime_files=[files['FSL1'], files['FSL2']],
                         scan_duration=[48, 20], TR=2, style='FSL')
    assert design2.shape == (34, 2), 'Returned design matrix has wrong shape'
    design3 = gen_design(stimtime_files=files['FSL1'], scan_duration=68, TR=2,
                         style='FSL')
    assert design3[24] != 0, (
        'design matrix should be non-zero 8 seconds after an event onset.')
    design4 = gen_design(stimtime_files=[files['FSL2']],
                         scan_duration=[48, 20], TR=2, style='FSL')
    assert np.all(np.isclose(design1 * 0.5, design4)), (
        'gen_design does not treat missing values correctly')
    design5 = gen_design(stimtime_files=[files['FSL2']],
                         scan_duration=[48, 20], TR=1)
    assert np.all(np.isclose(design4, design5[::2])), (
        'design matrices sampled at different frequency do not match'
        ' at corresponding time points')
    design6 = gen_design(stimtime_files=[files['AFNI1']],
                         scan_duration=[48, 20], TR=2, style='AFNI')
    assert np.all(np.isclose(design1, design6)), (
        'design matrices generated from AFNI style and FSL style do not match')

    def test_case1_vs_npss(self):


        component  =  pod_mach.PodMach()

        prob = create_problem(component)

        prob.setup()

        prob['comp.gam'] = 1.4
        prob['comp.R'] = 287.0
        prob['comp.BF'] = .9
        prob['comp.A_pod'] = 1.4
        prob['comp.L'] = 22.0
        prob['comp.prc'] = 12.5
        prob['comp.p_tube'] = 850.0
        prob['comp.T_ambient'] = 298.0
        prob['comp.mu'] = 1.846e-5
        prob['comp.M_duct'] = .95
        prob['comp.M_diff'] = .6
        prob['comp.cp'] = 1009.0
        prob['comp.delta_star'] = .07
        prob['comp.M_pod'] = .8

        prob.run()

        assert np.isclose(prob['comp.Re'], 3278799.304354, rtol=0.1)
        assert np.isclose(prob['comp.A_tube'], 18.600833, rtol=0.1)

def test_ISFC():
    curr_dir = os.path.dirname(__file__)

    mask_fname = os.path.join(curr_dir, 'mask.nii.gz')
    mask = io.load_boolean_mask(mask_fname)
    fnames = [os.path.join(curr_dir, 'subj1.nii.gz'),
              os.path.join(curr_dir, 'subj2.nii.gz')]
    masked_images = image.mask_images(io.load_images(fnames), mask)

    D = image.MaskedMultiSubjectData.from_masked_images(masked_images,
                                                        len(fnames))

    assert D.shape == (4, 5, 2), "Loaded data has incorrect shape"

    (ISFC, p) = brainiak.isfc.isfc(D, return_p=True, num_perm=100,
                                   two_sided=True, random_state=0)

    ground_truth = \
        [[1, 1, 0, -1],
         [1, 1, 0, -1],
         [0, 0, 1,  0],
         [-1, -1, 0, 1]]

    ground_truth_p = 1 - np.abs(ground_truth)

    assert np.isclose(ISFC, ground_truth).all(), \
        "Calculated ISFC does not match ground truth"

    assert np.isclose(p, ground_truth_p).all(), \
        "Calculated p values do not match ground truth"

def _get_intersection_bound_vector_plane(bound_vector, plane):
    distance_to_plane = dot(
        plane.point_in_plane - bound_vector.initial_point,
        plane.normal_vector)
    projected_vector_length = dot(
        bound_vector.free_vector,
        plane.normal_vector)

    distance_to_plane_close_to_zero = isclose(
        distance_to_plane,
        0,
        **config['numbers_close_kwargs'])
    projected_vector_length_close_to_zero = isclose(
        projected_vector_length,
        0,
        **config['numbers_close_kwargs'])
    if (
            distance_to_plane_close_to_zero and
            projected_vector_length_close_to_zero):
        return bound_vector

    with errstate(divide='ignore'):
        param = nan_to_num(distance_to_plane / projected_vector_length)

    # TODO: add distinction for included and excluded initial and terminal
    # points
    if 0 <= param <= 1:
        intersection = (
            bound_vector.initial_point +
            param*(bound_vector.terminal_point - bound_vector.initial_point))
    else:
        intersection = None
    return intersection

def test_Cacciato09Sats2():
    """
    Check that the model behavior is altered in the expected way by changing
    param_dict values.
    """
    model = Cacciato09Sats(threshold=9.5)
    nsat_exp = model.mean_occupation(prim_haloprop=5e13)
    # Increasing b_0 by x should increase the occupation by exactly 10**x.
    model.param_dict['b_0'] += 0.1
    nsat_exp_new = model.mean_occupation(prim_haloprop=5e13)
    assert np.isclose(nsat_exp_new, nsat_exp * 10**0.1, rtol=1e-2, atol=1.e-2)

    # Increasing b_1 by x should increase the occupation by exactly
    # 10**(x * (log prim_haloprop - 12.0)).
    model.param_dict['b_0'] -= 0.1
    model.param_dict['b_1'] += 0.1
    nsat_exp_new = model.mean_occupation(prim_haloprop=5e13)
    assert np.isclose(nsat_exp_new, nsat_exp * 10**(
        0.1 * (np.log10(5e13) - 12.0)), rtol=1e-2, atol=1.e-2)

    # Increasing b_2 by x should increase the occupation by exactly
    # 10**(x * (log prim_haloprop - 12.0)**2).
    model.param_dict['b_1'] -= 0.1
    model.param_dict['b_2'] += 0.1
    nsat_exp_new = model.mean_occupation(prim_haloprop=5e13)
    assert np.isclose(nsat_exp_new, nsat_exp * 10 ** (
        0.1 * (np.log10(5e13) - 12.0)**2), rtol=1e-2, atol=1.e-2)

def peak_effects(energies):
    """
    This module checks to see if a given list of peak energies contains any
    possible single escape peaks, double escape peaks, and sum peaks.
    """

    single_escape_peak = []
    single_escape_peak_index = []
    double_escape_peak = []
    double_escape_peak_index = []

    for i in range(len(energies)):
        #checks to see if condition for escape peaks is fulfilled.
        if energies[i] >= float(1022):
            for j in range(len(energies)):
                #checks to see if peak is a single escape peak.
                if np.isclose(energies[i],energies[j]+511,atol=1) == True:
                    single_escape_peak.append(energies[j])
                    single_escape_peak_index.append(j)
                #checks to see if peak is a double escape peak.
                if np.isclose(energies[i],energies[j]+1022,atol=1) == True:
                    double_escape_peak.append(energies[j])
                    double_escape_peak_index.append(j)
                    """
                    for m in range(len(energies)):
                        if np.isclose(energies[j],energies[m]+511,atol=1) == True:
                            double_escape_peak.append(energies[m])
                            double_escape_peak_index.append(m)
                    """
    #gets rid of the single escape peak that has the same value as a double
    #escape peak(s).
    double_count = []
    double_count_index = []

    for i in range(len(single_escape_peak)):
        for j in range(len(double_escape_peak)):
            if single_escape_peak[i] == double_escape_peak[j]:
                double_count.append(single_escape_peak[i])
                double_count_index.append(energies.index(single_escape_peak[i]))

    single_escape_peak = [x for x in single_escape_peak if x not in double_count]
    single_escape_peak_index = [x for x in single_escape_peak_index if x not in double_count_index]

    #finds the index of the original energy that is responsible for the single
    #escape peaks and the double escape peaks.
    origin_index_se = []
    origin_index_de = []
    for i in range(len(energies)):
        for j in range(len(single_escape_peak)):
            if np.isclose(energies[i],single_escape_peak[j]+511,atol=1) == True:
                origin_index_se.append(i)
        for k in range(len(double_escape_peak)):
            if np.isclose(energies[i],double_escape_peak[k]+1022,atol=1) == True:
                origin_index_de.append(i)

    peak_effects = {'single_escape_peak':single_escape_peak,'double_escape_peak':double_escape_peak,
                    'single_escape_peak_index':single_escape_peak_index,'double_escape_peak_index':double_escape_peak_index,
                    'origin_index_se':origin_index_se,'origin_index_de':origin_index_de}

    return(peak_effects)

def test_Cacciato09Cens6():
    """Check that the median primary galaxy property behave accordingly.
    """
    model = Cacciato09Cens(threshold=9.5)
    prim_galprop_1 = model.median_prim_galprop(prim_haloprop=1e14)
    model.param_dict['log_M_1'] += 0.1
    prim_galprop_2 = model.median_prim_galprop(prim_haloprop=1e14*10**0.1)
    assert np.isclose(prim_galprop_1, prim_galprop_2, rtol=1e-6, atol=1.e-2)

    model.param_dict['log_L_0'] += 0.1
    prim_galprop_3 = model.median_prim_galprop(prim_haloprop=1e14*10**0.1)
    assert np.isclose(prim_galprop_2 * 10**0.1, prim_galprop_3, rtol=1e-6,
                      atol=1.e-2)

    model.param_dict['gamma_1'] += 0.1
    prim_galprop_4 = model.median_prim_galprop(prim_haloprop=1e14*10**0.1)
    assert prim_galprop_3 != prim_galprop_4

    model.param_dict['gamma_2'] += 0.1
    prim_galprop_5 = model.median_prim_galprop(prim_haloprop=1e14 * 10 ** 0.1)
    assert prim_galprop_4 != prim_galprop_5

    model.param_dict['sigma'] += 0.1
    prim_galprop_6 = model.median_prim_galprop(prim_haloprop=1e14 * 10 ** 0.1)
    assert np.isclose(prim_galprop_5, prim_galprop_6, rtol=1e-6, atol=1.e-2)

    def _create_D_matrix(self):
        N = self.level_counts_cs[-1]
        D = np.zeros((N, N, len(self.laser_intensity)), dtype='object')
        bxrho = BxRho_Voigt if self.shape.lower() == 'voigt' else BxRho_Lorentzian

        self.indices = []
        for laser_index, laser in enumerate(self.laser_intensity):
            for i, j in itertools.combinations(range(len(self.level_counts)), 2):
                for k, (fe, mze) in enumerate(zip(self.Flist[i], self.MFlist[i])):
                    for l, (fg, mzg) in enumerate(zip(self.Flist[j], self.MFlist[j])):
                        x = self.level_counts_cs[i] - self.level_counts[i] + k
                        y = self.level_counts_cs[j] - self.level_counts[j] + l
                        if np.isclose(self.A_array[i, j], 0) or np.isclose(self.partial_A[x, y], 0):
                            continue
                        purity = self._params['Purity'].value
                        frac = purity if self.mode[laser_index] == (mze - mzg) else (1.0 - purity) if self.mode[laser_index] == -(mze - mzg) else 0
                        if frac == 0:
                            pass
                        else:
                            intensity = frac * self._params['Laser_intensity_' + str(laser_index)].value
                            A = self._params['Transition_strength_' + str(i) + '_to_' + str(j)].value
                            mu = (self.energies[k] + self.energy_change[k]) - (self.energies[l] + self.energy_change[l])
                            kwargs = {'A': A, 'mu': mu, 'laser': intensity}
                            if self.shape.lower() == 'voigt':
                                kwargs['fwhmG'] = self._params['FWHMG_' + str(i) + '_to_' + str(j)].value * 1e6
                                kwargs['fwhmL'] = self._params['FWHML_' + str(i) + '_to_' + str(j)].value * 1e6
                            else:
                                kwargs['fwhm'] = self._params['FWHML_' + str(i) + '_to_' + str(j)].value * 1e6
                            D[x, y, laser_index] = bxrho(**kwargs)
                            self.indices.append((x, y, laser_index, i, j, mze, mzg))

        self.D = D

    def inverse_covariance_matrix(self, alpha=0, beta=0):
        if 'mag' in self.covs:
            invcovmat = self.covs['mag'].copy()
        else:
            invcovmat = 0
        if self.alphabeta_covmat:
            if np.isclose(alpha, self._last_alpha) and np.isclose(beta, self._last_beta):
                return self.invcov
            self._last_alpha = alpha
            self._last_beta = beta

            alphasq = alpha * alpha
            betasq = beta * beta
            alphabeta = alpha * beta
            if 'stretch' in self.covs:
                invcovmat += alphasq * self.covs['stretch']
            if 'colour' in self.covs:
                invcovmat += betasq * self.covs['colour']
            if 'mag_stretch' in self.covs:
                invcovmat += 2 * alpha * self.covs['mag_stretch']
            if 'mag_colour' in self.covs:
                invcovmat -= 2 * beta * self.covs['mag_colour']
            if 'stretch_colour' in self.covs:
                invcovmat -= 2 * alphabeta * self.covs['stretch_colour']

            delta = self.pre_vars + alphasq * self.stretch_var + \
                    + betasq * self.colour_var + 2.0 * alpha * self.cov_mag_stretch \
                    - 2.0 * beta * self.cov_mag_colour \
                    - 2.0 * alphabeta * self.cov_stretch_colour
        else:
            delta = self.pre_vars
        np.fill_diagonal(invcovmat, invcovmat.diagonal() + delta)
        self.invcov = np.linalg.inv(invcovmat)
        return self.invcov