Ticket #8018: trac_8018-numerical-eigenforms.patch

File trac_8018-numerical-eigenforms.patch, 11.2 KB (added by rbeezer, 13 years ago)

  • sage/modular/modform/numerical.py 

    # HG changeset patch
# User Rob Beezer <beezer@ups.edu>
# Date 1264309004 28800
# Node ID fd30b4d4509853ce1580103dd09f634e3ef79c3a
# Parent  6c4190c4ba571ea4615296cce5a44a3d87c2cd0c
[mq]: 8018-modularform-eigen

diff -r 6c4190c4ba57 -r fd30b4d45098 sage/modular/modform/numerical.py
    a b  
    1818from sage.misc.misc              import verbose
    1919from sage.rings.all              import CDF, Integer, QQ, next_prime, prime_range
    2020from sage.misc.prandom           import randint
     21from sage.matrix.constructor     import matrix
     22
     23# This variable controls importing the SciPy library sparingly
     24scipy=None
    2125
    2226class NumericalEigenforms(SageObject):
    2327    """
    … …  
    2731
    2832    - ``group`` - a congruence subgroup of a Dirichlet character of
    2933      order 1 or 2
    30      
     34
    3135    - ``weight`` - an integer >= 2
    32    
     36
    3337    - ``eps`` - a small float; abs( ) < eps is what "equal to zero" is
    3438      interpreted as for floating point numbers.
    35      
     39
    3640    - ``delta`` - a small-ish float; eigenvalues are considered distinct
    3741      if their difference has absolute value at least delta
    38      
     42
    3943    - ``tp`` - use the Hecke operators T_p for p in tp when searching
    4044      for a random Hecke operator with distinct Hecke eigenvalues.
    41              
     45
    4246    OUTPUT:
    4347
    4448    A numerical eigenforms object, with the following useful methods:
    … …  
    5155          [[eigenvalues of T_2],
    5256           [eigenvalues of T_3],
    5357           [eigenvalues of T_5], ...]
    54            
     58
    5559    - :meth:`systems_of_eigenvalues` - a list of the systems of
    5660      eigenvalues of eigenforms such that the chosen random linear
    5761      combination of Hecke operators has multiplicity 1 eigenvalues.
    5862
    5963    EXAMPLES::
    60    
     64
    6165        sage: n = numerical_eigenforms(23)
    6266        sage: n == loads(dumps(n))
    6367        True
    … …  
    8690        Create a new space of numerical eigenforms.
    8791
    8892        EXAMPLES::
    89        
     93
    9094            sage: numerical_eigenforms(61) # indirect doctest
    9195            Numerical Hecke eigenvalues for Congruence Subgroup Gamma0(61) of weight 2
    9296        """
    … …  
    107111        symbols, and -1 otherwise.
    108112
    109113        EXAMPLES::
    110        
     114
    111115            sage: n = numerical_eigenforms(23)
    112116            sage: n.__cmp__(loads(dumps(n)))
    113117            0
    … …  
    124128        Return the level of this set of modular eigenforms.
    125129
    126130        EXAMPLES::
    127        
     131
    128132            sage: n = numerical_eigenforms(61) ; n.level()
    129133            61
    130134        """
    … …  
    135139        Return the weight of this set of modular eigenforms.
    136140
    137141        EXAMPLES::
    138        
     142
    139143            sage: n = numerical_eigenforms(61) ; n.weight()
    140144            2
    141145        """
    … …  
    146150        Print string representation of self.
    147151
    148152        EXAMPLES::
    149        
     153
    150154            sage: n = numerical_eigenforms(61) ; n
    151155            Numerical Hecke eigenvalues for Congruence Subgroup Gamma0(61) of weight 2
    152156
    … …  
    162166        set of modular eigenforms.
    163167
    164168        EXAMPLES::
    165        
     169
    166170            sage: n = numerical_eigenforms(61) ; n.modular_symbols()
    167171            Modular Symbols space of dimension 5 for Gamma_0(61) of weight 2 with sign 1 over Rational Field
    168172        """
    … …  
    177181            return M
    178182
    179183    def _eigenvectors(self):
    180         """
     184        r"""
    181185        Find numerical approximations to simultaneous eigenvectors in
    182186        self.modular_symbols() for all T_p in self._tp.
    183187
    184188        EXAMPLES::
    185        
     189
    186190            sage: n = numerical_eigenforms(61)
    187191            sage: n._eigenvectors() # random order
    188192            [              1.0    0.289473640239    0.176788851952    0.336707726757  2.4182243084e-16]
    … …  
    190194            [                0    0.413171180356    0.141163094698   0.0923242547901    0.707106781187]
    191195            [                0    0.826342360711    0.282326189397     0.18464850958 6.79812569682e-16]
    192196            [                0      0.2402380858    0.792225196393    0.905370774276 4.70805946682e-16]
     197
     198        TESTS:
     199
     200        This tests if this routine selects only eigenvectors with
     201        multiplicity one.  Two of the eigenvalues are
     202        (roughly) -92.21 and -90.30 so if we set ``eps = 2.0``
     203        then they should compare as equal, causing both eigenvectors
     204        to be absent from the matrix returned.  The remaining eigenvalues
     205        (ostensibly unique) are visible in the test, which should be
     206        indepedent of which eigenvectors are returned, but it does presume
     207        an ordering of these eigenvectors for the test to succeed.
     208        This exercises a correction in Trac 8018. ::
     209
     210            sage: n = numerical_eigenforms(61, eps=2.0)
     211            sage: evectors = n._eigenvectors()
     212            sage: evalues = diagonal_matrix(CDF, [-283.0, 108.522012456, 142.0])
     213            sage: diff = n._hecke_matrix*evectors - evectors*evalues
     214            sage: sum([abs(diff[i,j]) for i in range(5) for j in range(3)]) < 1.0e-9
     215            True
    193216        """
    194217        try:
    195218            return self.__eigenvectors
    … …  
    206229            t += randint(-50,50)*M.T(p).matrix()
    207230
    208231        self._hecke_matrix = t
    209         # evals, B = t.change_ring(CDF).right_eigenvectors()
    210         from sage.matrix.constructor import matrix
    211         spectrum = t.change_ring(CDF).right_eigenvectors()
    212         evals = [spectrum[i][0] for i in range(len(spectrum))]
    213         B = matrix(CDF, [spectrum[i][1][0] for i in range(len(spectrum))]).transpose()
    214        
    215         # Find the eigenvalues that occur with multiplicity 1 up
    216         # to the given eps.
     232
     233        global scipy
     234        if scipy is None:
     235            import scipy
     236        import scipy.linalg
     237        evals,eig = scipy.linalg.eig(self._hecke_matrix.numpy(), right=True, left=False)
     238        B = matrix(eig)
     239        v = [CDF(evals[i]) for i in range(len(evals))]
     240
     241        # Determine the eigenvectors with eigenvalues of multiplicity
     242        # one, with equality controlled by the value of eps
     243        # Keep just these eigenvectors
    217244        eps = self._eps
    218         v = list(evals)
    219         v.sort()
    220245        w = []
    221246        for i in range(len(v)):
    222247            e = v[i]
    223248            uniq = True
    224249            for j in range(len(v)):
    225                 if i != j and abs(e-v[j]) < eps:
     250                if uniq and i != j and abs(e-v[j]) < eps:
    226251                    uniq = False
    227252            if uniq:
    228253                w.append(i)
    229254        self.__eigenvectors = B.matrix_from_columns(w)
    230         return B
     255        return self.__eigenvectors
    231256
    232257    def _easy_vector(self):
    233258        """
    … …  
    245270           appropriate multiple.  Repeat.
    246271
    247272        EXAMPLES::
    248        
     273
    249274            sage: n = numerical_eigenforms(37)
    250275            sage: n._easy_vector()                 # slightly random output
    251276            (1.0, 1.0, 0)
    … …  
    265290        x = (CDF**E.nrows()).zero_vector()
    266291        if E.nrows() == 0:
    267292            return x
    268        
    269        
     293
     294
    270295
    271296        def best_row(M):
    272297            """
    … …  
    274299            with the most entries supported outside [-delta, delta].
    275300
    276301            EXAMPLES::
    277            
     302
    278303                sage: numerical_eigenforms(61)._easy_vector() # indirect doctest
    279304                (1.0, 1.0, 0, 0, 0)
    280305            """
    … …  
    298323            i, f = best_row(C)
    299324            x[i] += 1   # simplistic
    300325            e = x*E
    301        
     326
    302327        self.__easy_vector = x
    303328        return x
    304329
    … …  
    307332        Return all eigendata for self._easy_vector().
    308333
    309334        EXAMPLES::
    310        
     335
    311336            sage: numerical_eigenforms(61)._eigendata() # random order
    312337            ((1.0, 0.668205013164, 0.219198805797, 0.49263343893, 0.707106781187), (1.0, 1.49654668896, 4.5620686498, 2.02990686579, 1.41421356237), [0, 1], (1.0, 1.0))
    313338        """
    … …  
    316341        except AttributeError:
    317342            pass
    318343        x = self._easy_vector()
    319        
     344
    320345        B = self._eigenvectors()
    321346        def phi(y):
    322347            """
    … …  
    324349            linear combination of basis vectors.
    325350
    326351            EXAMPLES::
    327            
     352
    328353                sage: n = numerical_eigenforms(61) # indirect doctest
    329354                sage: n._eigendata() # random order
    330                 ((1.0, 0.668205013164, 0.219198805797, 0.49263343893, 0.707106781187), (1.0, 1.49654668896, 4.5620686498, 2.02990686579, 1.41421356237), [0, 1], (1.0, 1.0))               
     355                ((1.0, 0.668205013164, 0.219198805797, 0.49263343893, 0.707106781187), (1.0, 1.49654668896, 4.5620686498, 2.02990686579, 1.41421356237), [0, 1], (1.0, 1.0))
    331356            """
    332357            return y.element() * B
    333        
     358
    334359        phi_x = phi(x)
    335360        V = phi_x.parent()
    336361        phi_x_inv = V([a**(-1) for a in phi_x])
    … …  
    355380        - ``list`` - a list of double precision complex numbers
    356381
    357382        EXAMPLES::
    358        
     383
    359384            sage: n = numerical_eigenforms(11,4)
    360385            sage: n.ap(2) # random order
    361386            [9.0, 9.0, 2.73205080757, -0.732050807569]
    … …  
    380405        a = Sequence(self.eigenvalues([p])[0], immutable=True)
    381406        self._ap[p] = a
    382407        return a
    383            
     408
    384409    def eigenvalues(self, primes):
    385410        """
    386411        Return the eigenvalues of the Hecke operators corresponding
    387412        to the primes in the input list of primes.   The eigenvalues
    388         match up between one prime and the next. 
     413        match up between one prime and the next.
    389414
    390415        INPUT:
    391416
    392417        - ``primes`` - a list of primes
    393418
    394419        OUTPUT:
    395        
     420
    396421        list of lists of eigenvalues.
    397422
    398423        EXAMPLES::
    399        
     424
    400425            sage: n = numerical_eigenforms(1,12)
    401426            sage: n.eigenvalues([3,5,13])
    402427            [[177148.0, 252.0], [48828126.0, 4830.0], [1.79216039404e+12, -577737.999...]]
    … …  
    413438            linear combination of basis vectors.
    414439
    415440            EXAMPLES::
    416            
     441
    417442                sage: n = numerical_eigenforms(1,12)  # indirect doctest
    418443                sage: n.eigenvalues([3,5,13])
    419444                [[177148.0, 252.0], [48828126.0, 4830.0], [1.79216039404e+12, -577737.999...]]
    … …  
    434459        up to bound.
    435460
    436461        EXAMPLES::
    437        
     462
    438463            sage: numerical_eigenforms(61).systems_of_eigenvalues(10)
    439464            [
    440465            [-1.48119430409..., 0.806063433525..., 3.15632517466..., 0.675130870567...],
    … …  
    454479        v.sort()
    455480        v.set_immutable()
    456481        return v
    457        
     482
    458483    def systems_of_abs(self, bound):
    459484        """
    460485        Return the absolute values of all systems of eigenvalues for
    461486        self for primes up to bound.
    462487
    463488        EXAMPLES::
    464        
     489
    465490            sage: numerical_eigenforms(61).systems_of_abs(10)
    466491            [
    467492            [0.311107817466, 2.90321192591, 2.52542756084, 3.21431974338],
    … …  
    488513    indices where `|v|` is larger than eps.
    489514
    490515    EXAMPLES::
    491    
     516
    492517        sage: sage.modular.modform.numerical.support( numerical_eigenforms(61)._easy_vector(), 1.0 )
    493518        []
    494        
     519
    495520        sage: sage.modular.modform.numerical.support( numerical_eigenforms(61)._easy_vector(), 0.5 )
    496521        [0, 1]
    497522
    … …  
    499524    return [i for i in range(v.degree()) if abs(v[i]) > eps]
    500525
    501526
    502    
     527