
:>"^J\                 @   s  d  Z  d d l Z y e j Z Wn> e k
 r_ d d l Z e j d e j  d d   Z Yn Xd d   Z e j j	   d Z
 e j e
  Z Gd	 d
   d
 e  Z d d   Z d d   Z d d   Z d d d d d d  Z d Z d d d d d d d d d  Z d d   Z d d   Z d d   Z d d d d d  Z d d   Z d d  d!  Z d" d#   Z d$ d%   Z d& d'   Z d( d)   Z d* d+   Z d, d-   Z  d. d/   Z! d0 d1   Z" d2 d3   Z# d S)4a  A state-emitting MarkovModel.

Note terminology similar to Manning and Schutze is used.


Functions:
train_bw        Train a markov model using the Baum-Welch algorithm.
train_visible   Train a visible markov model using MLE.
find_states     Find the a state sequence that explains some observations.

load            Load a MarkovModel.
save            Save a MarkovModel.

Classes:
MarkovModel     Holds the description of a markov model
    NzVFor optimal speed, please update to Numpy version 1.3 or later (current version is %s)c             C   sf   | |  d k r | S|  | d k r( |  St  |  |  } | t j t j |  |  t j | |   S)z>Implement logaddexp method if Numpy version is older than 1.3.d   )minnumpylogexp)ZlogxZlogyZminxy r   4/tmp/pip-build-ww9dw3qa/biopython/Bio/MarkovModel.py	logaddexp&   s    r	   c             C   sX   i  } t  |  d d d   } t |   d } x" | D] \ } } | | | | <q6 W| S)zAReturn a dictionary of values with their sequence offset as keys.N   )	enumeratelen)valuesdentriesnindexkeyr   r   r   	itemindex0   s    r   gYnc               @   s7   e  Z d  Z d Z d d d d d  Z d d   Z d S)MarkovModelz+Create a state-emitting MarkovModel object.Nc             C   s1   | |  _  | |  _ | |  _ | |  _ | |  _ d S)zInitialize the class.N)statesalphabet	p_initialp_transition
p_emission)selfr   r   r   r   r   r   r   r   __init__C   s
    				zMarkovModel.__init__c             C   s=   d d l  m } |   } t |  |  | j d  | j   S)z9Create a string representation of the MarkovModel object.r   )StringIO)Z	Bio._py3kr   saveseekread)r   r   handler   r   r   __str__L   s
    	zMarkovModel.__str__)__name__
__module____qualname____doc__r   r"   r   r   r   r   r   @   s   r   c             C   s5   |  j    } | j |  s1 t d | | f   | S)zNRead the first line and evaluate that begisn with the correct start (PRIVATE).zI expected %r but got %r)readline
startswith
ValueError)r!   startliner   r   r   _readline_and_check_startU   s    r,   c             C   s
  t  |  d  } | j   d d  } t  |  d  } | j   d d  } t | |  } t |  t |  } } t j |  | _ t  |  d  } xN t t |   D]: } t  |  d | |  } t | j   d  | j | <q Wt j | | f  | _	 t  |  d  } xg t t |   D]S } t  |  d | |  } d d	   | j   d d  D | j	 | d d  f <qWt j | | f  | _
 t  |  d
  } xg t t |   D]S } t  |  d | |  } d d	   | j   d d  D | j
 | d d  f <qW| S)z.Parse a file handle into a MarkovModel object.zSTATES:r
   Nz	ALPHABET:zINITIAL:z  %s:zTRANSITION:c             S   s   g  |  ] } t  |   q Sr   )float).0vr   r   r   
<listcomp>v   s   	 zload.<locals>.<listcomp>z	EMISSION:c             S   s   g  |  ] } t  |   q Sr   )r-   )r.   r/   r   r   r   r0   }   s   	 r   )r,   splitr   r   r   zerosr   ranger-   r   r   )r!   r+   r   r   mmNMir   r   r   load]   s,    !::r8   c             C   sQ  | j  } | d d j |  j   | d d j |  j   | d  x? t t |  j   D]( } | d |  j | |  j | f  q] W| d  xU t t |  j   D]> } | d |  j | d j d d	   |  j | D  f  q W| d
  xU t t |  j   D]> } | d |  j | d j d d	   |  j | D  f  qWd S)z$Save MarkovModel object into handle.zSTATES: %s
 zALPHABET: %s
z	INITIAL:
z	  %s: %g
zTRANSITION:
z	  %s: %s
c             s   s   |  ] } t  |  Vq d  S)N)str)r.   xr   r   r   	<genexpr>   s    zsave.<locals>.<genexpr>z
EMISSION:
c             s   s   |  ] } t  |  Vq d  S)N)r:   )r.   r;   r   r   r   r<      s    N)	writejoinr   r   r3   r   r   r   r   )r4   r!   wr7   r   r   r   r      s    	
&
<
r   c                s  t  |   t  |  } } | s+ t d   | d k	 rd t j |  } | j | f k rd t d   | d k	 r t j |  } | j | | f k r t d   | d k	 r t j |  } | j | | f k r t d   g  }	 t |    x. | D]& }
 |	 j   f d d   |
 D  q Wd d   |	 D } t |  d	 k rPt d
   t | | |	 d | d | d | d | } | \ } } } t	 |  | | | |  S)a  Train a MarkovModel using the Baum-Welch algorithm.

    Train a MarkovModel using the Baum-Welch algorithm.  states is a list
    of strings that describe the names of each state.  alphabet is a
    list of objects that indicate the allowed outputs.  training_data
    is a list of observations.  Each observation is a list of objects
    from the alphabet.

    pseudo_initial, pseudo_transition, and pseudo_emission are
    optional parameters that you can use to assign pseudo-counts to
    different matrices.  They should be matrices of the appropriate
    size that contain numbers to add to each parameter matrix, before
    normalization.

    update_fn is an optional callback that takes parameters
    (iteration, log_likelihood).  It is called once per iteration.
    zNo training data given.Nz$pseudo_initial not shape len(states)z5pseudo_transition not shape len(states) X len(states)z5pseudo_emission not shape len(states) X len(alphabet)c                s   g  |  ] }   |  q Sr   r   )r.   r;   )indexesr   r   r0      s   	 ztrain_bw.<locals>.<listcomp>c             S   s   g  |  ] } t  |   q Sr   )r   )r.   r;   r   r   r   r0      s   	 r   z,I got training data with outputs of length 0pseudo_initialpseudo_transitionpseudo_emission	update_fn)
r   r)   r   asarrayshaper   appendr   _baum_welchr   )r   r   training_datarA   rB   rC   rD   r5   r6   training_outputsoutputsZlengthsr;   r   r   r   r   )r@   r   train_bw   s:    $	rL   i  c
             C   s  | d k r t  |   } n t | |  f  } | d k rN t  |  |  f  } n t | |  |  f  } | d k r t  |  | f  } n t | |  | f  } t j |  }
 t j |  } t j |  } | d k	 r t j |  } n d } | d k	 rt j |  } n d } | d k	 r,t j |  } n d } d } x t t  D] } t } x6 | D]. } | t |  | | |
 | | | | | 	 7} qXW|	 d k	 r|	 | |  | d k	 rt j | |  d k  rP| } qEWt	 d t   d d   |
 | | f D S)zfImplement the Baum-Welch algorithm to evaluate unknown parameters in the MarkovModel object (PRIVATE).Ng?z%HMM did not converge in %d iterationsc             S   s   g  |  ] } t  j |   q Sr   )r   r   )r.   _r   r   r   r0     s   	 z_baum_welch.<locals>.<listcomp>)
_random_norm_copy_and_checkr   r   r3   MAX_ITERATIONSLOG0_baum_welch_onefabsRuntimeError)r5   r6   rJ   r   r   r   rA   rB   rC   rD   
lp_initiallp_transitionlp_emissionlpseudo_initiallpseudo_transitionlpseudo_emissionZ	prev_llikr7   ZllikrK   r   r   r   rH      sJ    		%

rH   c	          	   C   s  t  |  }	 t |  |	 | | | |  }
 t |  |	 | | |  } t j |  |  |	 f  } x t |	  D] } | | } t j |  |  f  } xo t |   D]a } xX t |   D]J } |
 | | | | | | | | | | | d } | | | | <q Wq W| t |  | d d  d d  | f <qd Wt j |  |	 f  } xR t |	  D]D } x; t |   D]- } t | | d d  | f  | | | <q]WqJWt j |   } x4 t |   D]& } t | | d d  f  | | <qW| d d  d f } | d k	 rt | |  } | t |  } x t |   D] } xC t |   D]5 } t | | | d d  f  | | | | | <q9W| d k	 r&t | | |  | | <| | t | |  | | <q&Wx t |   D] } t j |  t } xY t |	  D]K } | | } x8 t |   D]* } t	 | | | | | | f  | | <qWqW| t |  } | d k	 rvt | | |  } | t |  } | | | d d  f <qWt |
 d d  |	 f  S)zExecute one step for Baum-Welch algorithm (PRIVATE).

    Do one iteration of Baum-Welch based on a sequence of output.
    Changes the value for lp_initial, lp_transition and lp_emission in place.
    r
   Nr   )
r   _forward	_backwardr   r2   r3   _logsum
_logvecaddrQ   r	   )r5   r6   rK   rU   rV   rW   rX   rY   rZ   TZfmatZbmatZlp_arctkZlp_traverser7   jlpZlp_arcout_tZ	lp_arcoutZksumr   r   r   rR     sT    
#-/$3 
,	rR   c             C   s   t  j |  | d f  } | | d d  d f <x t d | d  D] } | | d } xx t |   D]j }	 t }
 xM t |   D]? } | | | d | | |	 | | | } t |
 |  }
 q} W|
 | |	 | <qd WqC W| S)zImplement forward algorithm (PRIVATE).

    Calculate a Nx(T+1) matrix, where the last column is the total
    probability of the output.
    r
   Nr   )r   r2   r3   rQ   r	   )r5   r_   rU   rV   rW   rK   matrixr`   ra   rb   lprobr7   rc   r   r   r   r[   c  s    *r[   c             C   s   t  j |  | d f  } x t | d d d  D] } | | } xx t |   D]j } t }	 xM t |   D]? }
 | |
 | d | | |
 | | | } t |	 |  }	 qf W|	 | | | <qM Wq0 W| S)z'Implement backward algorithm (PRIVATE).r
   r   r   )r   r2   r3   rQ   r	   )r5   r_   rV   rW   rK   rd   r`   ra   r7   re   rb   rc   r   r   r   r\   z  s    
r\   c                s  t  |   t  |  } } | d k	 rR t j |  } | j | f k rR t d   | d k	 r t j |  } | j | | f k r t d   | d k	 r t j |  } | j | | f k r t d   g  g  } }	 t |    t |    xx | D]p \ }
 } t  |  t  |
  k r&t d   | j  f d d   | D  |	 j   f d d   |
 D  q Wt | | |	 | | | |  } | \ } } } t |  | | | |  S)	a  Train a visible MarkovModel using maximum likelihoood estimates for each of the parameters.

    Train a visible MarkovModel using maximum likelihoood estimates
    for each of the parameters.  states is a list of strings that
    describe the names of each state.  alphabet is a list of objects
    that indicate the allowed outputs.  training_data is a list of
    (outputs, observed states) where outputs is a list of the emission
    from the alphabet, and observed states is a list of states from
    states.

    pseudo_initial, pseudo_transition, and pseudo_emission are
    optional parameters that you can use to assign pseudo-counts to
    different matrices.  They should be matrices of the appropriate
    size that contain numbers to add to each parameter matrix.
    Nz$pseudo_initial not shape len(states)z5pseudo_transition not shape len(states) X len(states)z5pseudo_emission not shape len(states) X len(alphabet)zstates and outputs not alignedc                s   g  |  ] }   |  q Sr   r   )r.   r;   )states_indexesr   r   r0     s   	 z!train_visible.<locals>.<listcomp>c                s   g  |  ] }   |  q Sr   r   )r.   r;   )outputs_indexesr   r   r0     s   	 )	r   r   rE   rF   r)   r   rG   _mler   )r   r   rI   rA   rB   rC   r5   r6   training_statesrJ   ZtoutputsZtstatesr;   r   r   r   r   )rg   rf   r   train_visible  s2     $rj   c             C   s'  t  j |   } | r | | } x  | D] } | | d d 7<q& Wt |  } t  j |  |  f  }	 | rs |	 | }	 x\ | D]T } xK t t |  d  D]3 }
 | |
 | |
 d } } |	 | | f d 7<q Wqz WxZ t t |	   D]F } |	 | d d  f t |	 | d d  f  |	 | d d  f <q Wt  j |  | f  } | rT| | } t  j |  | f  } xQ t | |  D]@ \ } } x1 t | |  D]  \ } } | | | f d 7<qWqyWxZ t t |   D]F } | | d d  f t | | d d  f  | | d d  f <qW| |	 | f S)z<Implement Maximum likelihood estimation algorithm (PRIVATE).r   r
   N)r   r2   
_normalizer3   r   sumoneszip)r5   r6   rJ   ri   rA   rB   rC   r   r   r   r   r7   rb   r   rK   osr   r   r   rh     s2    

D
Drh   c             C   s   t  j |   g S)z?Return indeces of the maximum values aong the vector (PRIVATE).)r   Zargmax)Zvector	allowancer   r   r   	_argmaxes  s    rr   c       
         s   |   t   j  } t j  j t  } t j  j t  } t j  j t  } t  j	      f d d   | D } t
 | | | | |  } xV t t  |   D]B } | | \ } }	  f d d   | D t j |	  f | | <q W| S)zaFind states in the given Markov model output.

    Returns a list of (states, score) tuples.
    c                s   g  |  ] }   |  q Sr   r   )r.   r;   )r@   r   r   r0     s   	 zfind_states.<locals>.<listcomp>c                s   g  |  ] }   j  |  q Sr   )r   )r.   r;   )r4   r   r   r0     s   	 )r   r   r   r   r   VERY_SMALL_NUMBERr   r   r   r   _viterbir3   r   )
Zmarkov_modeloutputr5   rU   rV   rW   resultsr7   r   scorer   )r@   r4   r   find_states  s    0rx   c             C   s  t  |  } g  } x( t |   D] } | j d g |  q Wt j |  | f  } | | d d  | d f | d d  d f <x t d |  D] }	 | |	 }
 x t |   D]t } | d d  |	 d f | d d  | f | | |
 f } t |  } | | d | | |	 f <| | | |	 <q Wq Wg  } g  } t | d d  | d f  } x7 | D]/ } | j | d | g | | | d f  q\Wx | r| j   \ }	 } } |	 d k r| j | | f  q| | d |	 } x/ | D]' } | j |	 d | g | | f  qWqW| S)zSImplement Viterbi algorithm to find most likely states for a given input (PRIVATE).Nr   r
   )r   r3   rG   r   r2   rr   pop)r5   rU   rV   rW   ru   r_   Z	backtracer7   Zscoresr`   ra   rb   Zi_scoresr@   Z
in_processrv   r   rw   r   r   r   rt     s6    .
+ -	)rt   c             C   s   t  |  j  d k r. |  t t |    }  n t  |  j  d k r xi t t  |    D]F } |  | d d  f t |  | d d  f  |  | d d  f <qV Wn t d   |  S)z"Normalize matrix object (PRIVATE).r
      Nz&I cannot handle matrixes of that shape)r   rF   r-   rl   r3   r)   )rd   r7   r   r   r   rk   0  s    Grk   c             C   s   t  j |   } t |  S)z%Normalize a uniform matrix (PRIVATE).)r   rm   rk   )rF   rd   r   r   r   _uniform_norm=  s    r{   c             C   s   t  j j |   } t |  S)z$Normalize a random matrix (PRIVATE).)r   randomrk   )rF   rd   r   r   r   rN   C  s    rN   c             C   s   t  j |  d d }  |  j | k r0 t d   t |  j  d k rs t  j t |   d  d k r t d   nt t |  j  d k r x\ t t |    D]9 } t  j t |  |  d  d k r t d |   q Wn t d	   |  S)
zFCopy a matrix and check its dimension. Normalize at the end (PRIVATE).copyr
   zIncorrect dimensiong      ?g{Gz?zmatrix not normalized to 1.0rz   zmatrix %d not normalized to 1.0z&I don't handle matrices > 2 dimensions)r   arrayrF   r)   r   rS   rl   r3   )rd   Zdesired_shaper7   r   r   r   rO   I  s    #rO   c             C   si   t  |  j  d k r9 t j |  t j |  j  f  } n |  } t } x | D] } t | |  } qL W| S)z/Implement logsum for a matrix object (PRIVATE).r
   )r   rF   r   ZreshapeproductrQ   r	   )rd   Zvecrl   numr   r   r   r]   ]  s    $r]   c             C   su   t  |   t  |  k s$ t d   t j t  |    } x5 t t  |    D]! } t |  | | |  | | <qL W| S)z5Implement a log sum for two vector objects (PRIVATE).zvectors aren't the same length)r   AssertionErrorr   r2   r3   r	   )Zlogvec1Zlogvec2Zsumvecr7   r   r   r   r^   i  s
    $r^   c             C   s   t  |   } t j |  S)z-Return the exponential of a logsum (PRIVATE).)r]   r   r   )numbersrl   r   r   r   _exp_logsumr  s    r   )$r&   r   r	   AttributeErrorwarningswarn__version__r   r|   seedrs   r   rQ   objectr   r,   r8   r   rL   rP   rH   rR   r[   r\   rj   rh   rr   rx   rt   rk   r{   rN   rO   r]   r^   r   r   r   r   r   <module>   sJ   

%	:	9T1'+	