← Back to branch summary

~ubuntu-branches/ubuntu/natty/python-cogent/natty

« back to all changes in this revision

Viewing changes to cogent/evolve/parameter_controller.py

  • Committer: Bazaar Package Importer
  • Author(s): Steffen Moeller
  • Date: 2010-12-04 22:30:35 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20101204223035-j11kinhcrrdgg2p2
Tags: 1.5-1
* Bumped standard to 3.9.1, no changes required.
* New upstream version.
  - major additions to Cookbook
  - added AlleleFreqs attribute to ensembl Variation objects.
  - added getGeneByStableId method to genome objects.
  - added Introns attribute to Transcript objects and an Intron class.
  - added Mann-Whitney test and a Monte-Carlo version
  - exploratory and confirmatory period estimation techniques (suitable for
    symbolic and continuous data)
  - Information theoretic measures (AIC and BIC) added
  - drawing of trees with collapsed nodes
  - progress display indicator support for terminal and GUI apps
  - added parser for illumina HiSeq2000 and GAiix sequence files as 
    cogent.parse.illumina_sequence.MinimalIlluminaSequenceParser.
  - added parser to FASTQ files, one of the output options for illumina's
    workflow, also added cookbook demo.
  - added functionality for parsing of SFF files without the Roche tools in
    cogent.parse.binary_sff
  - thousand fold performance improvement to nmds
  - >10-fold performance improvements to some Table operations

Show diffs side-by-side

added added

removed removed

Lines of Context:
cogent/evolve/parameter_controller.py
5
5
"""
6
6
 
7
7
from __future__ import with_statement
 
8
 
 
9
import numpy
 
10
import pickle, warnings
 
11
 
8
12
from cogent.core.tree import TreeError
9
13
from cogent.evolve import likelihood_calculation
10
14
from cogent.align import dp_calculation
11
15
from cogent.evolve.likelihood_function import LikelihoodFunction as _LF
12
16
from cogent.recalculation.scope import _indexed
 
17
from cogent.maths.stats.information_criteria import aic, bic
13
18
 
14
 
import numpy
15
 
import pickle
16
19
 
17
20
from cogent.align.pairwise import AlignableSeq
 
21
from cogent.util.warning import discontinued
18
22
 
19
 
__author__ = "Andrew Butterfield"
 
23
__author__ = "Peter Maxwell"
20
24
__copyright__ = "Copyright 2007-2009, The Cogent Project"
21
25
__credits__ = ["Andrew Butterfield", "Peter Maxwell", "Gavin Huttley",
22
26
                    "Helen Lindsay"]
23
27
__license__ = "GPL"
24
 
__version__ = "1.4.1"
 
28
__version__ = "1.5.0"
25
29
__maintainer__ = "Gavin Huttley"
26
30
__email__ = "gavin.huttley@anu.ed.au"
27
31
__status__ = "Production"
64
68
        self.motifs = self._motifs = model.getMotifs()
65
69
        self._mprob_motifs = list(model.getMprobAlphabet())
66
70
        defn = self.makeLikelihoodDefn(**kw)
67
 
        self.real_par_controller = defn.makeParamController()
 
71
        super(_LF, self).__init__(defn)
68
72
        self.setDefaultParamRules()
69
73
        self.setDefaultTreeParameterRules()
70
74
        self.mprobs_from_alignment = motif_probs_from_align
76
80
        f = open(filename, 'w')
77
81
        temp = {}
78
82
        try:
79
 
            for d in self.real_par_controller.defns:
 
83
            for d in self.defns:
80
84
                temp[id(d)] = d.values
81
85
                del d.values
82
86
            pickle.dump((1.0, None, self), f)
83
87
        finally:
84
 
            for d in self.real_par_controller.defns:
 
88
            for d in self.defns:
85
89
                if id(d) in temp:
86
90
                    d.values = temp[id(d)]
87
91
    
88
 
    def updateIntermediateValues(self):
89
 
        self.real_par_controller.update()
90
 
    
91
 
    def optimise(self, *args, **kw):
92
 
        return_calculator = kw.pop('return_calculator', False)
93
 
        lc = self.real_par_controller.makeCalculator()
94
 
        lc.optimise(*args, **kw)
95
 
        self.real_par_controller.updateFromCalculator(lc)
96
 
        if return_calculator:
97
 
            return lc
98
 
    
99
 
    def graphviz(self, **kw):
100
 
        lc = self.real_par_controller.makeCalculator()
101
 
        return lc.graphviz(**kw)
102
 
    
103
 
    def __repr__(self):
104
 
        return repr(self.real_par_controller)
105
 
    
106
92
    def setDefaultTreeParameterRules(self):
107
93
        """Lengths are set to the values found in the tree (if any), and
108
94
        free to be optimised independently.
109
95
        Other parameters are scoped based on the unique values found in the
110
96
        tree (if any) or default to having one value shared across the whole
111
97
        tree""" 
112
 
        with self.real_par_controller.updatesPostponed():
 
98
        with self.updatesPostponed():
113
99
            edges = self.tree.getEdgeVector()
114
100
            for par_name in self.model.getParamList():
115
101
                try:
157
143
        assert expm in ['pade', 'either', 'eigen', 'checked'], expm
158
144
        self.setParamRule('expm', is_const=True, value=expm)
159
145
 
160
 
    def makeCalculator(self, aligns):
161
 
        # deprecate
162
 
        self.setAlignment(aligns)
163
 
        if getattr(self, 'used_as_calculator', False):
164
 
            warnings.warn('PC used as two different calculators', stacklevel=2)
165
 
        self.used_as_calculator = True
166
 
        return self
 
146
    def makeCalculator(self, *args, **kw):
 
147
        if args:
 
148
            discontinued('method', "makeCalculator(aligns)", '1.6')
 
149
            # and shadowing a quite different superclass method.
 
150
            self.setAlignment(*args)
 
151
            if getattr(self, 'used_as_calculator', False):
 
152
                warnings.warn('PC used as two different calculators', stacklevel=2)
 
153
            self.used_as_calculator = True
 
154
            return self
 
155
        else:
 
156
            return super(_LF, self).makeCalculator(**kw)
167
157
    
168
158
    def _process_scope_info(self, edge=None, tip_names=None, edges=None,
169
159
            is_clade=None, is_stem=None, outgroup_name=None):
250
240
        elif init is not None:
251
241
            assert not value
252
242
            value = init
253
 
        self.real_par_controller.assignAll(par_name, scopes, value, lower,
254
 
                upper, is_const, is_independent)
 
243
        self.assignAll(par_name, scopes, value, lower, upper, is_const, 
 
244
                is_independent)
255
245
    
256
246
    def setLocalClock(self, tip1name, tip2name):
257
247
        """Constrain branch lengths for tip1name and tip2name to be equal.
275
265
        if tree is None:
276
266
            tree = self.tree
277
267
        
278
 
        with self.real_par_controller.updatesPostponed():
 
268
        with self.updatesPostponed():
279
269
            for edge in tree.getEdgeVector():
280
270
                if edge.Length is None or edge.Name in exclude_list:
281
271
                    continue
282
272
                self.setParamRule("length", edge=edge.Name, is_const=1,
283
273
                                        value=edge.Length)
284
274
    
 
275
    def getAic(self, second_order=False):
 
276
        """returns Aikake Information Criteria
 
277
        
 
278
        Arguments:
 
279
            - second_order: if true, the second-order AIC is returned,
 
280
              adjusted by the alignment length"""
 
281
        if second_order:
 
282
            sequence_length = sum(len(self.getParamValue('lht', locus=l).index)
 
283
                                    for l in self.locus_names)
 
284
        else:
 
285
            sequence_length = None
 
286
        
 
287
        lnL = self.getLogLikelihood()
 
288
        nfp = self.getNumFreeParams()
 
289
        return aic(lnL, nfp, sequence_length)
 
290
    
 
291
    def getBic(self):
 
292
        """returns the Bayesian Information Criteria"""
 
293
        sequence_length = sum(len(self.getParamValue('lht', locus=l).index)
 
294
                                for l in self.locus_names)
 
295
        lnL = self.getLogLikelihood()
 
296
        nfp = self.getNumFreeParams()
 
297
        return bic(lnL, nfp, sequence_length)
285
298
 
286
299
class AlignmentLikelihoodFunction(_LikelihoodParameterController):
287
300
    
288
301
    def setDefaultParamRules(self):
289
302
        try:
290
 
            self.real_par_controller.assignAll(
 
303
            self.assignAll(
291
304
                'fixed_motif', None, value=-1, const=True, independent=True)
292
305
        except KeyError:
293
306
            pass
294
307
    
295
 
    def makeLikelihoodDefn(self, sites_independent=True):
296
 
        defns = self.model.makeParamControllerDefns(
297
 
                bin_names=self.bin_names)
 
308
    def makeLikelihoodDefn(self, sites_independent=True, discrete_edges=None):
 
309
        defns = self.model.makeParamControllerDefns(bin_names=self.bin_names)
 
310
        if discrete_edges is not None:
 
311
            from discrete_markov import PartialyDiscretePsubsDefn
 
312
            defns['psubs'] = PartialyDiscretePsubsDefn(
 
313
                    self.motifs, defns['psubs'], discrete_edges)
298
314
        return likelihood_calculation.makeTotalLogLikelihoodDefn(
299
315
            self.tree, defns['align'], defns['psubs'], defns['word_probs'],
300
316
            defns['bprobs'], self.bin_names, self.locus_names,
316
332
                                (self.tree.getTipNames(), aln.getSeqNames(),
317
333
                                locus_name)
318
334
            assert not "root" in aln.getSeqNames(), "'root' is a reserved name."
319
 
        with self.real_par_controller.updatesPostponed():
 
335
        with self.updatesPostponed():
320
336
            for (locus_name, align) in zip(self.locus_names, aligns):
321
 
                self.real_par_controller.assignAll(
 
337
                self.assignAll(
322
338
                        'alignment', {'locus':[locus_name]},
323
339
                        value=align, const=True)
324
340
                if self.mprobs_from_alignment:
352
368
        self.setPogs(leaves, locus=locus)
353
369
    
354
370
    def setPogs(self, leaves, locus=None):
355
 
        with self.real_par_controller.updatesPostponed():
 
371
        with self.updatesPostponed():
356
372
            for (name, pog) in leaves.items():
357
373
                self.setParamRule('leaf', edge=name, value=pog, is_const=True)
358
374
            if self.mprobs_from_alignment: