def run_checks(self, **kwargs):
        """
        Check to see whether the system is expected to be computable.

        This is called by default for each set_value but will only raise a
        logger warning if fails.  This is also called immediately when calling
        :meth:`run_compute`.

        kwargs are passed to override currently set values as if they were
        sent to :meth:`run_compute`.

        :return: True if passed, False if failed and a message
        """

        # make sure all constraints have been run
        changed_params = self.run_delayed_constraints()

        hier = self.hierarchy
        if hier is None:
            return True, ''
        for component in hier.get_stars():
            kind = hier.get_kind_of(component)
            comp_ps = self.get_component(component)

            if not len(comp_ps):
                return False, "component '{}' in the hierarchy is not in the bundle".format(component)

            parent = hier.get_parent_of(component)
            parent_ps = self.get_component(parent)
            if kind in ['star']:
                # ignore the single star case
                if parent:
                    # contact systems MUST by synchronous
                    if hier.is_contact_binary(component):
                        if self.get_value(qualifier='syncpar', component=component, context='component', **kwargs) != 1.0:
                            return False,\
                                'contact binaries must by synchronous, but syncpar@{}!=1'.format(component)

                        if self.get_value(qualifier='ecc', component=parent, context='component', **kwargs) != 0.0:
                            return False,\
                                'contact binaries must by circular, but ecc@{}!=0'.format(component)

                        if self.get_value(qualifier='pitch', component=component, context='component', **kwargs) != 0.0:
                            return False,\
                                'contact binaries must be aligned, but pitch@{}!=0'.format(component)

                        if self.get_value(qualifier='yaw', component=component, context='component', **kwargs) != 0.0:
                            return False,\
                                'contact binaries must be aligned, but yaw@{}!=0'.format(component)

                    # MUST NOT be overflowing at PERIASTRON (d=1-ecc, etheta=0)

                    requiv = comp_ps.get_value('requiv', unit=u.solRad, **kwargs)
                    requiv_max = comp_ps.get_value('requiv_max', unit=u.solRad, **kwargs)



                    if hier.is_contact_binary(component):
                        if np.isnan(requiv) or requiv > requiv_max:
                            return False,\
                                '{} is overflowing at L2/L3 (requiv={}, requiv_max={})'.format(component, requiv, requiv_max)


                        requiv_min = comp_ps.get_value('requiv_min')

                        if np.isnan(requiv) or requiv <= requiv_min:
                            return False,\
                                 '{} is underflowing at L1 and not a contact system (requiv={}, requiv_min={})'.format(component, requiv, requiv_min)
                        elif requiv <= requiv_min * 1.001:
                            return False,\
                                'requiv@{} is too close to requiv_min (within 0.1% of critical).  Use detached/semidetached model instead.'.format(component)

                    else:
                        if requiv > requiv_max:
                            return False,\
                                '{} is overflowing at periastron (requiv={}, requiv_max={})'.format(component, requiv, requiv_max)

            else:
                raise NotImplementedError("checks not implemented for type '{}'".format(kind))

        # we also need to make sure that stars don't overlap each other
        # so we'll check for each pair of stars (see issue #70 on github)
        # TODO: rewrite overlap checks
        for orbitref in []: #hier.get_orbits():
            if len(hier.get_children_of(orbitref)) == 2:
                q = self.get_value(qualifier='q', component=orbitref, context='component', **kwargs)
                ecc = self.get_value(qualifier='ecc', component=orbitref, context='component', **kwargs)

                starrefs = hier.get_children_of(orbitref)
                if hier.get_kind_of(starrefs[0]) != 'star' or hier.get_kind_of(starrefs[1]) != 'star':
                    # print "***", hier.get_kind_of(starrefs[0]), hier.get_kind_of(starrefs[1])
                    continue
                if self.get_value(qualifier='pitch', component=starrefs[0])!=0.0 or \
                        self.get_value(qualifier='pitch', component=starrefs[1])!=0.0 or \
                        self.get_value(qualifier='yaw', component=starrefs[0])!=0.0 or \
                        self.get_value(qualifier='yaw', component=starrefs[1])!=0.0:

                    # we cannot run this test for misaligned cases
                   continue

                comp0 = hier.get_primary_or_secondary(starrefs[0], return_ind=True)
                comp1 = hier.get_primary_or_secondary(starrefs[1], return_ind=True)
                q0 = roche.q_for_component(q, comp0)
                q1 = roche.q_for_component(q, comp1)

                F0 = self.get_value(qualifier='syncpar', component=starrefs[0], context='component', **kwargs)
                F1 = self.get_value(qualifier='syncpar', component=starrefs[1], context='component', **kwargs)

                pot0 = self.get_value(qualifier='pot', component=starrefs[0], context='component', **kwargs)
                pot0 = roche.pot_for_component(pot0, q0, comp0)

                pot1 = self.get_value(qualifier='pot', component=starrefs[1], context='component', **kwargs)
                pot1 = roche.pot_for_component(pot1, q1, comp1)

                xrange0 = libphoebe.roche_xrange(q0, F0, 1.0-ecc, pot0+1e-6, choice=0)
                xrange1 = libphoebe.roche_xrange(q1, F1, 1.0-ecc, pot1+1e-6, choice=0)

                if xrange0[1]+xrange1[1] > 1.0-ecc:
                    return False,\
                        'components in {} are overlapping at periastron (change ecc@{}, syncpar@{}, or syncpar@{}).'.format(orbitref, orbitref, starrefs[0], starrefs[1])


        # check to make sure passband supports the selected atm
        for pbparam in self.filter(qualifier='passband').to_list():
            pb = pbparam.get_value()
            pbatms = _pbtable[pb]['atms']
            # NOTE: atms are not attached to datasets, but per-compute and per-component
            for atmparam in self.filter(qualifier='atm', kind='phoebe').to_list():
                atm = atmparam.get_value()
                if atm not in pbatms:
                    return False, "'{}' passband ({}) does not support atm='{}' ({}).".format(pb, pbparam.twig, atm, atmparam.twig)

        # check length of ld_coeffs vs ld_func and ld_func vs atm
        def ld_coeffs_len(ld_func, ld_coeffs):
            # current choices for ld_func are:
            # ['interp', 'uniform', 'linear', 'logarithmic', 'quadratic', 'square_root', 'power', 'claret', 'hillen', 'prsa']
            if ld_func == 'interp':
                return True,
            elif ld_func in ['linear'] and len(ld_coeffs)==1:
                return True,
            elif ld_func in ['logarithmic', 'square_root', 'quadratic'] and len(ld_coeffs)==2:
                return True,
            elif ld_func in ['power'] and len(ld_coeffs)==4:
                return True,
            else:
                return False, "ld_coeffs={} wrong length for ld_func='{}'.".format(ld_coeffs, ld_func)

        for component in self.hierarchy.get_stars():
            # first check ld_coeffs_bol vs ld_func_bol
            ld_func = str(self.get_value(qualifier='ld_func_bol', component=component, context='component', check_visible=False, **kwargs))
            ld_coeffs = np.asarray(self.get_value(qualifier='ld_coeffs_bol', component=component, context='component', check_visible=False, **kwargs))
            check = ld_coeffs_len(ld_func, ld_coeffs)
            if not check[0]:
                return check

            if ld_func != 'interp':
                check = libphoebe.ld_check(ld_func, ld_coeffs)
                if not check:
                    return False, 'ld_coeffs_bol={} not compatible for ld_func_bol=\'{}\'.'.format(ld_coeffs, ld_func)

            for dataset in self.datasets:
                if dataset=='_default' or self.get_dataset(dataset=dataset, kind='*dep').kind not in ['lc_dep', 'rv_dep']:
                    continue
                ld_func = str(self.get_value(qualifier='ld_func', dataset=dataset, component=component, context='dataset', **kwargs))
                ld_coeffs = np.asarray(self.get_value(qualifier='ld_coeffs', dataset=dataset, component=component, context='dataset', check_visible=False, **kwargs))
                if ld_coeffs is not None:
                    check = ld_coeffs_len(ld_func, ld_coeffs)
                    if not check[0]:
                        return check

                if ld_func != 'interp':
                    check = libphoebe.ld_check(ld_func, ld_coeffs)
                    if not check:
                        return False, 'ld_coeffs={} not compatible for ld_func=\'{}\'.'.format(ld_coeffs, ld_func)

                if ld_func=='interp':
                    for compute in kwargs.get('computes', self.computes):
                        atm = self.get_value(qualifier='atm', component=component, compute=compute, context='compute', **kwargs)
                        if atm != 'ck2004':
                            return False, "ld_func='interp' only supported by atm='ck2004'.  Either change atm@{} or ld_func@{}@{}".format(component, component, dataset)

        # mesh-consistency checks
        for compute in self.computes:
            mesh_methods = [p.get_value() for p in self.filter(qualifier='mesh_method', compute=compute, force_ps=True).to_list()]
            if 'wd' in mesh_methods:
                if len(set(mesh_methods)) > 1:
                    return False, "all (or none) components must use mesh_method='wd'."

        #### WARNINGS ONLY ####
        # let's check teff vs gravb_bol and irrad_frac_refl_bol
        for component in self.hierarchy.get_stars():
            teff = self.get_value(qualifier='teff', component=component, context='component', unit=u.K, **kwargs)
            gravb_bol = self.get_value(qualifier='gravb_bol', component=component, context='component', **kwargs)

            if teff >= 8000. and gravb_bol < 0.9:
                return None, "'{}' probably has a radiative atm (teff={:.0f}K>8000K), for which gravb_bol=1.00 might be a better approx than gravb_bol={:.2f}.".format(component, teff, gravb_bol)
            elif teff <= 6600. and gravb_bol >= 0.9:
                return None, "'{}' probably has a convective atm (teff={:.0f}K<6600K), for which gravb_bol=0.32 might be a better approx than gravb_bol={:.2f}.".format(component, teff, gravb_bol)
            elif gravb_bol < 0.32 or gravb_bol > 1.00:
                return None, "'{}' has intermittent temperature (6600K<teff={:.0f}K<8000K), gravb_bol might be better between 0.32-1.00 than gravb_bol={:.2f}.".format(component, teff, gravb_bol)

        for component in self.hierarchy.get_stars():
            teff = self.get_value(qualifier='teff', component=component, context='component', unit=u.K, **kwargs)
            irrad_frac_refl_bol = self.get_value(qualifier='irrad_frac_refl_bol', component=component, context='component', **kwargs)

            if teff >= 8000. and irrad_frac_refl_bol < 0.8:
                return None, "'{}' probably has a radiative atm (teff={:.0f}K>8000K), for which irrad_frac_refl_bol=1.00 might be a better approx than irrad_frac_refl_bol={:.2f}.".format(component, teff, irrad_frac_refl_bol)
            elif teff <= 6600. and irrad_frac_refl_bol >= 0.75:
                return None, "'{}' probably has a convective atm (teff={:.0f}K<6600K), for which irrad_frac_refl_bol=0.6 might be a better approx than irrad_frac_refl_bol={:.2f}.".format(component, teff, irrad_frac_refl_bol)
            elif irrad_frac_refl_bol < 0.6:
                return None, "'{}' has intermittent temperature (6600K<teff={:.0f}K<8000K), irrad_frac_refl_bol might be better between 0.6-1.00 than irrad_frac_refl_bol={:.2f}.".format(component, teff, irrad_frac_refl_bol)

        # TODO: add other checks
        # - make sure all ETV components are legal
        # - check for conflict between dynamics_method and mesh_method (?)

        # we've survived all tests
        return True, ''