def calc_B_1s_diss(self,op,n):
        """Applies a single-site operator to a single site and returns
        the parameter tensor for that site after the change with the
        change in the norm of the state projected out.
        
        Parameters
        ----------
        op : ndarray or callable
            The single-site operator. See self.expect_1s().
        n: int
            The site to apply the operator to.
        """
        if callable(op):
            op = sp.vectorize(op, otypes=[sp.complex128])
            op = sp.fromfunction(op, (self.q[n], self.q[n]))
            
        newAn = sp.zeros_like(self.A[n])
        
        for s in range(self.q[n]):
            for t in range(self.q[n]):
                newAn[s] += self.A[n][t] * op[s, t]
                
        r_nm1 = TDVP.tm.eps_r_noop(self.r[n], newAn, newAn)
        ev = mm.adot(self.l[n - 1], r_nm1)

        newAn -= ev * self.A[n] #norm-fixing
        
        return newAn

def interp2d(qx, qy, qz):
    Vandermonde = sp.zeros((4,4))
    Vandermonde[:,0] = 1
    Vandermonde[:,1] = qx
    Vandermonde[:,2] = qy
    Vandermonde[:,3] = qx*qy
    Vinv = sp.linalg.inv(Vandermonde)
    
    print 'Vandermonde\n', Vandermonde
    print
    print 'Vandermonde inverse official \n', Vinv    
    Vinv = inverse(Vandermonde, 4)
    print 'Vandermonde inverse Gauss \n', Vinv
    V22 = sp.copy(Vinv.T)
    print 'Identity check'
    print sp.dot(Vinv,Vandermonde)
    print 'Transpose official'
    print V22
    for i in range(3):
        for j in range(i+1,4):
            d = Vinv[i,j]
            Vinv[i,j]= Vinv[j,i]
            Vinv[j,i]= d
    print 'Index ranspose\n', Vinv
    print 'Check transpose\n', Vinv-V22
    
    
    def SU2(x,y):
        RHS = sp.array([1,x,y,x*y])
        b = sp.dot(Vinv,RHS)
        return sp.dot(b,qz.T)
    SU2 = sp.vectorize(SU2)
    return SU2

def relabel_map(label_image, mapping, key=lambda x, y: x[y]):
    """
    Relabel an image using the supplied mapping.
    
    The mapping can be any kind of subscriptable object. The respective region id is used
    to access the new value from the mapping. The key keyword parameter can be used to
    supply another access function. The key function must have the signature
    key(mapping, region-id) and return the new region-id to assign.
    
    @param label_image A nD label_image
    @type label_image sequence
    @param mapping A mapping object
    @type mapping subscriptable object
    
    @return A binary image
    @rtype numpy.ndarray
    
    @raise ArgumentError If a region id is missing in the supplied mapping
    """    
    label_image = scipy.array(label_image)
    
    def _map(x):
        try:
            return key(mapping, x)
        except Exception as e:
            raise ArgumentError('No conversion for region id {} found in the supplied mapping. Error: {}'.format(x, e))
    
    vmap = scipy.vectorize(_map, otypes=[label_image.dtype])
         
    return vmap(label_image)

def simple_wf():
    wf = WorkFlow()

    
    # declare docs:
    wf.docs = [
        Doc("Testing advanced 1",  "a b c c c d d d d e"),
        Doc("Testing advanced 2", "a a a a a b c c c d d d e e"),
        Doc("Testing advanced 3", "b b b b f f f f"),
    ]
    
    # declare preprocessing
    wf.title_index = [ doc.title for doc in wf.docs] 
    wf.word_index =  ['a', 'b', 'c', 'd', 'e', 'f']
    
    # prepare tf matrix     
    wf.count_matrix = array([ 
        [1, 5, 0], #a 
        [1, 1, 4], #b 
        [3, 3, 0], #c 
        [4, 3, 0], #d 
        [1, 2, 0], #e 
        [0, 0, 4], #f 
    ])      
    count_to_tf = vectorize(math_utils.count_to_tf)
    wf.tf_mat = count_to_tf(wf.count_matrix)

    # prepare lidf vector: log of inverted df
    wf.n_docs       = float(len(wf.docs))                
    wf.df_vec       = array([2.0, 3.0, 2.0, 2.0, 2.0, 1.0])
    wf.idf_vec      = wf.n_docs / wf.df_vec
    wf.log_idf_vec  = log(wf.idf_vec)
    wf.wieghts_mat  = matrix(wf.tf_mat * wf.log_idf_vec[:,None])
    
    return wf

    def expect_2s(self, op, n):
        """Computes the expectation value of a nearest-neighbour two-site operator.

        The operator should be a q[n] x q[n + 1] x q[n] x q[n + 1] array
        such that op[s, t, u, v] = <st|op|uv> or a function of the form
        op(s, t, u, v) = <st|op|uv>.

        Parameters
        ----------
        o : ndarray or callable
            The operator array or function.
        n : int
            The leftmost site number (operator acts on n, n + 1).
        """
        A = self.get_A(n)
        Ap1 = self.get_A(n + 1)
        AA = tm.calc_AA(A, Ap1)

        if callable(op):
            op = sp.vectorize(op, otypes=[sp.complex128])
            op = sp.fromfunction(op, (A.shape[0], Ap1.shape[0], A.shape[0], Ap1.shape[0]))

        C = tm.calc_C_mat_op_AA(op, AA)
        res = tm.eps_r_op_2s_C12_AA34(self.get_r(n + 1), C, AA)
        return mm.adot(self.get_l(n - 1), res)

    def expect_3s(self, op, n):
        """Computes the expectation value of a nearest-neighbour three-site operator.

        The operator should be a q[n] x q[n + 1] x q[n + 2] x q[n] x
        q[n + 1] x q[n + 2] array such that op[s, t, u, v, w, x] =
        <stu|op|vwx> or a function of the form op(s, t, u, v, w, x) =
        <stu|op|vwx>.

        The state must be up-to-date -- see self.update()!

        Parameters
        ----------
        op : ndarray or callable
            The operator array or function.
        n : int
            The leftmost site number (operator acts on n, n + 1, n + 2).

        Returns
        -------
        expval : floating point number
            The expectation value (data type may be complex)
        """
        A = self.A[n]
        Ap1 = self.A[n + 1]
        Ap2 = self.A[n + 2]
        AAA = tm.calc_AAA(A, Ap1, Ap2)

        if callable(op):
            op = sp.vectorize(op, otypes=[sp.complex128])
            op = sp.fromfunction(op, (A.shape[0], Ap1.shape[0], Ap2.shape[0],
                                      A.shape[0], Ap1.shape[0], Ap2.shape[0]))

        C = tm.calc_C_3s_mat_op_AAA(op, AAA)
        res = tm.eps_r_op_3s_C123_AAA456(self.r[n + 2], C, AAA)
        return m.adot(self.l[n - 1], res)

 def expect_1s_diss(self,op,n):
     """Applies a single-site operator to a single site and returns
     the value after the change. In contrast to
     mps_gen.apply_op_1s, this routine does not change the state itself.
     
     Also, this does not perform self.update().
     
     Parameters
     ----------
     op : ndarray or callable
         The single-site operator. See self.expect_1s().
     n: int
         The site to apply the operator to.
     """
     if callable(op):
         op = sp.vectorize(op, otypes=[sp.complex128])
         op = sp.fromfunction(op, (self.q[n], self.q[n]))
         
     newAn = sp.zeros_like(self.A[n])
     
     for s in xrange(self.q[n]):
         for t in xrange(self.q[n]):
             newAn[s] += self.A[n][t] * op[s, t]
             
     return newAn

    def compute_estimator(self, svm, validation_set):
        """Compute least-square estimator for validation set"""
        val_set_dp = validation_set['dtp']
        val_set_cl = validation_set['cl']
        estimator = 0.0
        validation_size = len(val_set_dp)
        print "Computing estimator with validation part..."

        val_output = svm.get_output_2d(val_set_dp)

        # Compute estimator
        diff = val_set_cl - val_output
        estimator = diff.dot(diff).sum()
        # See p. 183
        estimator *= 1.0 * self.M / self.n

        classify_vect = s.vectorize(svm.classify_output)
        output_classes = classify_vect(val_output)

        diff_classes = output_classes - val_set_cl
        errors = s.count_nonzero(diff_classes)
        classified_correctly = validation_size - errors

        print "Classified correctly: %u/%u (%.2f%%)" % \
              (classified_correctly, validation_size,
               100.0 * classified_correctly / validation_size)
        return estimator

 def classify_2d(self, xs):
     """Classify array of data points, return an array of 1 and -1"""
     size = len(xs)
     output = self.get_output_2d(xs)
     classify_vect = s.vectorize(self.classify_output)
     output_classes = classify_vect(output)
     return output_classes

def normalizeHistogram(data):
    histogram = scipy.ndimage.histogram(data.astype("f"), 0, 255, 256)
    cumulatedHistogram = scipy.cumsum(histogram)
    nch = cumulatedHistogram.astype("f")/len(data.flat)
    inch = (nch*255).astype("i")
    normalize = scipy.vectorize(lambda i: inch[i])
    return normalize(data)

 def funapprox(self, f):
     """Polynomial approximation of a function"""
     c = self.funfitf(f)
     basis = self.basis
     def fhat(x):
         return sum([ cj * phi(self._trans(x))
                      for phi, cj in zip(basis, c)])
     return sp.vectorize(fhat)

    def expect_1s_1s(self, op1, op2, d):
        """Computes the expectation value of two single site operators acting 
        on two different sites.
        
        The result is < op1_n op2_n+d > with the operators acting on sites
        n and n + d.
        
        See expect_1s().
        
        Requires d > 0.
        
        The state must be up-to-date -- see self.update()!
        
        Parameters
        ----------
        op1 : ndarray or callable
            The first operator, acting on the first site.
        op2 : ndarray or callable
            The second operator, acting on the second site.
        d : int
            The distance (number of sites) between the two sites acted on non-trivially.
            
        Returns
        -------
        expval : floating point number
            The expectation value (data type may be complex)
        """        
        
        assert d > 0, 'd must be greater than 1'
        
        if callable(op1):
            op1 = sp.vectorize(op1, otypes=[sp.complex128])
            op1 = sp.fromfunction(op1, (self.q, self.q))
        
        if callable(op2):
            op2 = sp.vectorize(op2, otypes=[sp.complex128])
            op2 = sp.fromfunction(op2, (self.q, self.q)) 
        
        r_n = tm.eps_r_op_1s(self.r, self.A, self.A, op2)

        for n in xrange(d - 1):
            r_n = tm.eps_r_noop(r_n, self.A, self.A)

        r_n = tm.eps_r_op_1s(r_n, self.A, self.A, op1)
         
        return m.adot(self.l, r_n)

    def expect_1s_1s(self, op1, op2, n1, n2):
        """Computes the expectation value of two single site operators acting 
        on two different sites.
        
        The result is < op1 op2 >.
        
        See expect_1s().
        
        Requires n1 < n2.
        
        The state must be up-to-date -- see self.update()!
        
        Parameters
        ----------
        op1 : ndarray or callable
            The first operator, acting on the first site.
        op2 : ndarray or callable
            The second operator, acting on the second site.
        n1 : int
            The site number of the first site.
        n2 : int
            The site number of the second site (must be > n1).
            
        Returns
        -------
        expval : floating point number
            The expectation value (data type may be complex)
        """        
        if callable(op1):
            op1 = sp.vectorize(op1, otypes=[sp.complex128])
            op1 = sp.fromfunction(op1, (self.q[n1], self.q[n1]))
        
        if callable(op2):
            op2 = sp.vectorize(op2, otypes=[sp.complex128])
            op2 = sp.fromfunction(op2, (self.q[n2], self.q[n2])) 
        
        r_n = tm.eps_r_op_1s(self.r[n2], self.A[n2], self.A[n2], op2)

        for n in reversed(xrange(n1 + 1, n2)):
            r_n = tm.eps_r_noop(r_n, self.A[n], self.A[n])

        r_n = tm.eps_r_op_1s(r_n, self.A[n1], self.A[n1], op1)
         
        return m.adot(self.l[n1 - 1], r_n)

    def _convert_snps_to_binary(self,reference_ecotype='6909'):
        missingVal = 'NA'
        decoder = {missingVal:-1} #Might cause errors somewhere???!!!
        coding_fun = scipy.vectorize(lambda x: decoder[x], otypes=['int8'])

        if reference_ecotype in self.accessions:
            ref_i = self.accessions.index(reference_ecotype)
        else:
            ref_i = 0
            log.warning("Given reference ecotype %s wasn't found, using %s as 0-reference." % \
                    (reference_ecotype, self.accessions[ref_i]))
        snps = []
        num = len(self.positions)
        positions = []
        chr_region_removed_snps = [0] * len(self.chr_regions)
        for snp_i, (snp, pos) in enumerate(itertools.izip(self._snps, self.positions)):
            chr_region_ix,chr_region = self.get_chr_region_from_index(snp_i)
            if snp_i % (num / 10) == 0:
                log.info('Converted %s/%s of SNPs.' % (snp_i,num))
            unique_nts = scipy.unique(snp).tolist()
            if missingVal in unique_nts:
                if len(unique_nts) != 3:
                    chr_region_removed_snps[chr_region_ix] +=1
                    continue #Skipping non-binary SNP
                else:
                    unique_nts.remove(self.missingVal)
            else:
                if len(unique_nts) != 2:
                    chr_region_removed_snps[chr_region_ix] +=1
                    continue #Skipping non-binary SNP
            if snp[ref_i] != missingVal:
                ref_nt = snp[ref_i]
                decoder[ref_nt] = 0
                unique_nts.remove(ref_nt)
                decoder[unique_nts[0]] = 1
            else:
                decoder[unique_nts[0]] = 0
                decoder[unique_nts[1]] = 1
            snps.append(coding_fun(snp))
            positions.append(pos)
        log.info('Removed %d non-binary SNPs out of %d, when converting to binary SNPs.'\
            % (len(self.positions) - len(positions), len(self.positions)))
        assert len(snps) == len(positions), 'Somthing odd with the lengths.'
        chr_regions = self.chr_regions[:]
        sum_removed_snps = 0
        for i,num_snps in enumerate(chr_region_removed_snps):
            sum_removed_snps +=num_snps
            if i == 0:
                chr_regions[0] = (0,chr_regions[0][1] - num_snps)
            else:
                chr_regions[i] = (chr_regions[i-1][1],chr_regions[i][1] - sum_removed_snps)
        self._chr_regions = chr_regions

        self._snps = snps
        self._positions = positions
        self._data_format = 'binary'

    def __init__(self, type, slen=4, alen=1, lexsize=256):
        """
        Initialize an agent with an exhaustive lexicon made up of
        4-vowel stems (TODO: incorporate consonants,
            affixes & variable word-lengths)
        """
        # vowels {i,u,e,o} in articulatory features (hi, bk, rd) \in {-1,0,1}
        self.vowels = N.array(((1.0, 0.0, 0.0),
                               (1.0, 1.0, 0.0),
                               (0.0, 0.0, 0.0),
                               (0.0, 1.0, 0.0)))
        self.vf = {(1.0, 0.0, 0.0): "i",
                   (1.0, 1.0, 0.0): "u",
                   (0.0, 0.0, 0.0): "e",
                   (0.0, 1.0, 0.0): "o"}
        self.consonants = list("bcdfghjklmnpqrstvwxyz")
        # acoustic:articulatory mapping fxn for vowel prototypes
        # acoustic reps are F1,F2' pairs, articulatory reps are feature-based
        self.vowel_map = {}
        self.vowel_spread = 0
        self.memory = N.empty((lexsize, slen, 2))
        # each agent has its own articulatory variability
        #TODO: maybe this should be inferred by the learners
        #      on the basis of their data?
        self.alpha = N.random.normal(15, 2)
        self.beta = N.random.normal(2, 0.25)
        if self.beta < 1.0:
            self.beta = 1.1

        if type == "learner":
            self.stems = N.empty((lexsize, 4, 3), dtype=float)
            #self.affixes = N.empty((1,4))
        elif type == "speaker":
            tmp = [[x, y, 0.0] for x in [0.0, 1.0] for y in [0.0, 1.0]]
            self.stems = N.array([[a, b, c, d] for a in tmp for b in tmp
                                  for c in tmp for d in tmp])
        else:
            sys.exit("Undefined agent type. Aborting.")
        # vectorized versions of some fxns
        self.vec_perceive = vectorize(self.perceive)
        self.vec_articulate = vectorize(self.articulate)
        self.vec_acoustify = vectorize(self.acoustify)

 def pm(m, digits=5):
     def get_digits(x):
         if x == 0: return 0
         return int(math.log(abs(x) * 10.0, 10.0))
     vd = scipy.vectorize(get_digits)
     dm = vd(m)
     widths = dm.max(0)
     (l, c) = m.shape
     for i in xrange(l):
         EigenBase.pv(m[i,:],
                      digits=digits,
                      left_digits=widths)

def chebyshev_vec(values, degree=1):
    """Calculate the Chebyshev Polynobials
    
    This implementation uses sp.vectorize to vectorize python's math functions)"""
    
    values = sp.asarray(values)
    A = sp.zeros((degree, len(values)))
    
    A[0,:]=1
    try:
        A[1,:]=values
    except IndexError:
        return A
    
    cos = sp.vectorize(sp.math.cos)
    acos = sp.vectorize(sp.math.acos)
    
    for i in range(2,degree): 
        A[i,:] = cos(i*acos(values))

    return A

    def expect_1s_cor(self, op1, op2, n1, n2):
        """Computes the correlation of two single site operators acting on two different sites.

        See expect_1S().

        n1 must be smaller than n2.

        Assumes that the state is normalized.

        Parameters
        ----------
        op1 : function
            The first operator, acting on the first site.
        op2 : function
            The second operator, acting on the second site.
        n1 : int
            The site number of the first site.
        n2 : int
            The site number of the second site (must be > n1).
        """
        A1 = self.get_A(n1)
        A2 = self.get_A(n2)

        if callable(op1):
            op1 = sp.vectorize(op1, otypes=[sp.complex128])
            op1 = sp.fromfunction(op1, (A1.shape[0], A1.shape[0]))

        if callable(op2):
            op2 = sp.vectorize(op2, otypes=[sp.complex128])
            op2 = sp.fromfunction(op2, (A2.shape[0], A2.shape[0]))

        r_n = tm.eps_r_op_1s(self.get_r(n2), A2, A2, op2)

        for n in reversed(xrange(n1 + 1, n2)):
            r_n = tm.eps_r_noop(r_n, self.get_A(n), self.get_A(n))

        r_n = tm.eps_r_op_1s(r_n, A1, A1, op1)

        return mm.adot(self.get_l(n1 - 1), r_n)

    def __init__(self,c,grid,thresh=0.0000001,maxiter=30,req=None):
        self.grid = grid
        self.thresh = thresh
        self.maxiter = maxiter
        self.lapmat = disclap.disclap(grid)
        if req == None:
            self.req = 0.9*self.grid.r
        else:
            self.req = req

        self.c = c
        self.cvec = vectorize(c)(self.grid.zv)
        self.csqvec = abs(self.cvec*self.cvec)
        self.u = self._iterate()

 def _iterate(self):
     self.uzero = vectorize(uinitfunc)(self.req,self.grid.zv,self.csqvec)
     u = np.copy(self.uzero)
     n = 0
     delta = max(abs(self.op_bvec(u)))
     sys.stderr.write('PDE: GOAL=%f\n' % self.thresh)
     while (n < self.maxiter) and (delta > self.thresh):
         udot = self.step(u)
         u = u + 0.9*udot
         delta = max(abs(self.op_bvec(u)))
         udotnorm = max(abs(udot))
         n = n + 1
         sys.stderr.write('PDE: error=%f  delta=%f\n' % (delta,udotnorm))
     return u