def update_image(original_im, ci_red, ci_green, ci_blue):

    # diagnostics = dict()

    original_im = scipy.transpose(original_im)
    # diagnostics['original_im'] = original_im
    # diagnostics['ci_red'] = ci_red
    # diagnostics['ci_green'] = ci_green
    # diagnostics['ci_blue'] = ci_blue

    new_r = scipy.multiply(original_im[0], original_im[0] > ci_red)

    new_g = scipy.multiply(original_im[1], original_im[1] > ci_green)

    new_b = scipy.multiply(original_im[2], original_im[2] > ci_blue)

    new_im = (new_r, new_g, new_b)

    new_im = scipy.transpose(new_im)
    # diagnostics['new_im'] = new_im

    # with open('/Users/lages/Documents/sauceda/pictures_processed/diagnostics'
    #           '.p', 'wb') as f:
    #     pickle.dump(diagnostics, f)

    return new_im

def update_image(original_im, ci_red, ci_green, ci_blue):

    ci_vec = sp.array((ci_red, ci_green, ci_blue))
    ci_matrix = sp.multiply(sp.ones(original_im.shape), ci_vec)
    new_im = sp.multiply(original_im, original_im > ci_matrix)

    return new_im

  def generateGaborMotherWavelet(self):
    pitch = 440.0
    sigma = 6.
    NL = 48
    NU = 39
    print 'sampling rate:', self.fs, 'Hz'
    fs = float(self.fs)
    self.sample_duration = 10.
    #asigma = 0.3
    limit_t = 0.1
    #zurashi = 1.

    #NS = NL + NU + 1
    f = sp.array([2**(i/12.) for i in range(NL+NU+1)]) * pitch*2**(-NL/12.)
    f = f[:, sp.newaxis]
    sigmao = sigma*10**(-3)*sp.sqrt(fs/f)
    t = sp.arange(-limit_t, limit_t+1/fs, 1/fs)

    inv_sigmao = sp.power(sigmao, -1)
    inv_sigmao_t = inv_sigmao * t
    t_inv_sigmao2 = sp.multiply(inv_sigmao_t, inv_sigmao_t)
    omega_t = 2*sp.pi*f*t
    gabor = (1/sp.sqrt(2*sp.pi))
    gabor = sp.multiply(gabor, sp.diag(inv_sigmao))
    exps = -0.5*t_inv_sigmao2+sp.sqrt(-1)*omega_t
    self.gabor = gabor*sp.exp(exps)

def Au(U,GF,EpsArr,NX,NY,NZ):
    """Returns the result of matrix-vector multiplication
       by the system matrix A=I-GX
    """
    # reshaping input vector into 4-D array
    Uarr=sci.reshape(U,(NX,NY,NZ,3))
    # extended zero-padded arrays
    Uext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Vext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Jext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    JFext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Uext[0:NX,0:NY,0:NZ,:]=Uarr
    # contrast current array
    s=0
    while s<=2:
        Jext[0:NX,0:NY,0:NZ,s]=Uext[0:NX,0:NY,0:NZ,s]*(EpsArr[0:NX,0:NY,0:NZ]-1.0)
        JFext[:,:,:,s]=fft.fftn(sci.squeeze(Jext[:,:,:,s]))
        s=s+1
    Vext[:,:,:,0]=Uext[:,:,:,0]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,0,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,0,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,0,2],JFext[:,:,:,2])))
    Vext[:,:,:,1]=Uext[:,:,:,1]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,1,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,1,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,1,2],JFext[:,:,:,2])))
    Vext[:,:,:,2]=Uext[:,:,:,2]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,2,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,2,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,2,2],JFext[:,:,:,2])))
    # reshaping output into column vector
    V=sci.reshape(Vext[0:NX,0:NY,0:NZ,:],(NX*NY*NZ*3,1))

    return V

 def _calcVanillaOnlineGradient(self, sample, shapedfitnesses):
     invSigma = inv(self.sigma)
     phi = zeros(self.numDistrParams)
     phi[: self.numParameters] = self._logDerivX(sample, self.x, invSigma)
     logDerivSigma = self._logDerivFactorSigma(sample, self.x, invSigma, self.factorSigma)
     phi[self.numParameters :] = logDerivSigma.flatten()
     index = len(self.allSamples) % self.batchSize
     self.phiSquareWindow[index] = multiply(phi, phi)
     baseline = self._calcBaseline(shapedfitnesses)
     gradient = multiply((ones(self.numDistrParams) * shapedfitnesses[-1] - baseline), phi)
     return gradient

 def fit(self, train_pairs, verbose=False):
     n = len(train_pairs)
     if n == 0:
         raise NameError('Error: Train set is empty')
     else:
         if verbose:
             print('fit: Fitting a multiplicative model on %d pairs' % n)
         bases = [w for w, _ in train_pairs]
         derivs = [w for _, w in train_pairs]
         B = self.space.get_rows(bases).mat
         D = self.space.get_rows(derivs).mat
         DB = sp.multiply(B, D)
         BB = sp.multiply(B, B)
         self.mul_vector = DB.sum(axis=0) / (n * BB.sum(axis=0))

    def _calcVanillaBatchGradient(self, samples, shapedfitnesses):
        invSigma = inv(self.sigma)

        phi = zeros((len(samples), self.numDistrParams))
        phi[:, : self.numParameters] = self._logDerivsX(samples, self.x, invSigma)
        logDerivFactorSigma = self._logDerivsFactorSigma(samples, self.x, invSigma, self.factorSigma)
        phi[:, self.numParameters :] = array(logDerivFactorSigma)
        Rmat = outer(shapedfitnesses, ones(self.numDistrParams))

        # optimal baseline
        self.phiSquareWindow = multiply(phi, phi)
        baselineMatrix = self._calcBaseline(shapedfitnesses)

        gradient = sum(multiply(phi, (Rmat - baselineMatrix)), 0)
        return gradient

def costFunctionReg(flattendTheta, X, y, lmbda):
    """
    Calculate the cost and gradient for logistic regression
    using regularization (helps with preventing overfitting
    with many features)
    """
    # numpy fmin function only allows flattened arrays instead of
    # matrixes which is stupid so it has to be converted every time
    flattendTheta = sp.asmatrix(flattendTheta)
    (a, b) = flattendTheta.shape
    if a < b:
        theta = flattendTheta.T
    else:
        theta = flattendTheta
    m = sp.shape(y)[0]
    (J, grad) = costFunction(theta, X, y)

    # f is a filter vector that will disregard regularization for theta0
    f = sp.ones((theta.shape[0], 1))
    f[0, 0] = 0
    thetaFiltered = sp.multiply(theta, f)

    J = J + (lmbda/(2.0 * m)) * (thetaFiltered.T.dot(thetaFiltered))
    grad = grad + ((lmbda/m) * thetaFiltered).T

    return (J, grad)

def bondOrientation2sh(atoms,basis,l,neighbs=None,rcut=None,debug=False):
    atoms = array(atoms)
    basis = array(basis)    
    atoms = rectify(atoms,basis)

    if neighbs==None:
        bounds=[[0,basis[0][0]],[0,basis[1][1]],[0,basis[2][2]]]

        if rcut==None:
            rcut = generateRCut(atoms,basis,debug=debug)
            #print "Automatically generating r-cutoff=",rcut

        neighbs = secondShell( neighbors(atoms,bounds,rcut) )

    #sum the spherical harmonic over ever neighbor pair
    a = 4*np.pi / (2*l+1.)
    Ql=list()
    for i,ineighbs in enumerate(neighbs):
        n=len(ineighbs)

        shij = np.vectorize(complex)(zeros(2*l+1)) #spherical harmonic for bond i-j
        for j in ineighbs:
            shij += pairSphereHarms(atoms[i],minImageAtom(atoms[i],atoms[j],basis),l)/n
        shi = a * sum( scipy.real( scipy.multiply(shij,scipy.conj(shij)) ) )
        Ql.append(shi**0.5)
    
    return Ql,rcut

def spectral_radius(net, typed=True, weighted=True):
    '''
    Spectral radius of the graph, defined as the eigenvalue of greatest module.
    
    Parameters
    ----------
    net : :class:`~nngt.Graph` or subclass
        Network to analyze.
    typed : bool, optional (default: True)
        Whether the excitatory/inhibitory type of the connnections should be
        considered.
    weighted : bool, optional (default: True)
        Whether the weights should be taken into account.
    
    Returns
    -------
    the spectral radius as a float.
    '''
    weights = None
    if typed and "type" in net.graph.eproperties.keys():
        weights = net.eproperties["type"].copy()
    if weighted and "weight" in net.graph.eproperties.keys():
        if weights is not None:
            weights = sp.multiply(weights,
                                  net.graph.eproperties["weight"])
        else:
            weights = net.graph.eproperties["weight"].copy()
    mat_adj = adjacency(net.graph,weights)
    eigenval = [0]
    try:
        eigenval = spl.eigs(mat_adj,return_eigenvectors=False)
    except spl.eigen.arpack.ArpackNoConvergence,err:
        eigenval = err.eigenvalues

    def solve_v(self, s8y):
        fl = s8y.flatten()
        self._v = fl / fl.sum()
        v = copy.copy(self._e)
        step = 0

        def write_current_matrix():
            f = open("%s/%s_%03d.v" % (temporary_directory(), self._debug_matrix_file, step), "w")
            x = v.reshape(len(self._p1.modules), len(self._p2.modules))
            for i in xrange(len(self._p1.modules)):
                for j in xrange(len(self._p2.modules)):
                    f.write("%f " % x[i, j])
                f.write("\n")
            f.close()

        while 1:
            if self._debug:
                write_current_matrix()
            new = self.step(v)
            r = v - new
            r = scipy.multiply(r, r)
            s = r.sum()
            if s < 0.0000001 and step >= 10:
                return v
            step += 1
            v = new

 def scale_manual(self, event, val=None):
     a = P4Rm()
     if val is not None:
         P4Rm.ParamDict['DW_multiplication'] = val
     P4Rm.ParamDict['dwp'] = multiply(a.ParamDict['dwp'],
                                      a.ParamDict['DW_multiplication'])
     pub.sendMessage(pubsub_Re_Read_field_paramters_panel, event=event)

 def calc_E( self ):
     filename = os.path.join(self.prefix, 'wf_spectrum_dot_kp8.dat')
     if not os.path.isfile(filename):
        print 'ERROR: file %s not found\n' % (filename)
        sys.exit(1)        
  
     #E=na.genfromtxt(filename,unpack=True)[1][1:]
     E=na.loadtxt(filename,skiprows=1,unpack=True)[1]
     #state=na.genfromtxt(filename,unpack=True)[0][1:]
     state=na.loadtxt(filename,skiprows=1,unpack=True)[0]
     diff_E = []   
     diff_st = []
     for i1 in range(len(E)):
          for i2 in range(len(E)):
               #diff=item2-item1
               diff=abs( E[i2]-E[i1] )
               #if diff > 0 and diff > abs(E[i1]) and diff > abs(E[i2]):
               a=Set([diff])
               b=Set(diff_E)
               if not a.intersection(b):
                  diff_E.append(diff)
                  diff_st.append([state[i1],state[i2]])
     diff_E=sp.multiply(diff_E, 1.0*na.ones(len(diff_E)))
     #print len(diff_E),len(diff_st)
     #print diff_E,diff_st
     #diff_E[0]=sp.true_divide(1240.0*na.ones(len(diff_E[0])),diff_E)
     out_E=list([diff_E,diff_st])
     out_E_T=zip(*out_E)
     out_E_sort=sorted(out_E_T, key=itemgetter(0))
     #out_E.sort()
     #print out_E_sort#, diff_st[0]
     #print diff_E
     return out_E_sort 

 def calc_E( self ):
     load=self.loadDat()
     state_e=load[0]
     state_h=load[1]
     E_e=load[2]
     E_h=load[3]
     diff_E = []   
     diff_st = []
     #we are considering e-h transitions only
     for i1 in range(len(E_e)):
          for i2 in range(len(E_h)):
               #diff=item2-item1
               diff=abs( E_h[i2]-E_e[i1] )
               #if diff > 0 and diff > abs(E[i1]) and diff > abs(E[i2]):
               a=Set([diff])
               b=Set(diff_E)
               if not a.intersection(b):
                  diff_E.append(diff)
                  diff_st.append([state_e[i1],state_h[i2]])
     diff_E=sp.multiply(diff_E,1000.0*na.ones(len(diff_E)))
     #print len(diff_E),len(diff_st)
     #print diff_E,diff_st
     #diff_E[0]=sp.true_divide(1240.0*na.ones(len(diff_E[0])),diff_E)
     out_E=list([diff_E,diff_st])
     out_E_T=zip(*out_E)
     # print 'out_E_T', out_E_T
     out_E_sort=sorted(out_E_T, key=itemgetter(0))
     #out_E.sort()
     print 'sorted', out_E_sort#, diff_st[0]
     #print '\n\n\n'
     print 'out sorted', out_E_sort[0][0]
     #print '\n\n\n'
     return out_E_sort 

 def solve_v(self, s8y):
     fl = s8y.flatten()
     self._v = fl / fl.sum()
     v = copy.copy(self._e)
     step = 0
     def write_current_matrix():
         f = open('%s/%s_%03d.v' % (tempfile.gettempdir(),
                                    self._debug_matrix_file, step), 'w')
         x = v.reshape(len(self._p1.modules),
                       len(self._p2.modules))
         for i in xrange(len(self._p1.modules)):
             for j in xrange(len(self._p2.modules)):
                 f.write('%f ' % x[i,j])
             f.write('\n')
         f.close()
     while 1:
         if self._debug:
             write_current_matrix()
         new = self.step(v)
         r = (v-new)
         r = scipy.multiply(r,r)
         s = r.sum()
         if s < 0.0000001 and step >= 10:
             return v
         step += 1
         v = new

 def _updateSigmas(self, updateSize, lastSample):
     for c in range(self.numberOfCenters):
         self.sigmas[c] *= (1. - updateSize[c])
         dif = self.mus[c] - lastSample
         if self.diagonalOnly:
             self.sigmas[c] += updateSize[c] * multiply(dif, dif)
         else:
             self.sigmas[c] += updateSize[c] * 1.2 * outer(dif, dif)

def rbf(inputs, centroids, weights):
    if len(inputs) > 0:
        icw = np.array([[inputs[i], centroids[i], weights[i]]
                        for i in inputs.keys()])
        sw = np.absolute(np.subtract(icw[:, 0], icw[:, 1]))
        return np.exp(-10 * np.multiply(sw, icw[:, 2]).sum())  # /len(inputs))
    else:
        return 0

 def lineBroadening(self, fromPos, toPos, LB):
     """Applies line broadening of given widh (in Hz) to the FID. Resulting spectrum
     (after fft is called) will be convolution of the original spectrum (fromPos)
     and Lorentzian with full-width-at-half-maximum equal to LB"""
     self.checkToPos(toPos)
     #print("Len Fid Time: {}".format(len(self.fidTime)))
     #print("Len All Fid 0: {}".format(len(self.allFid[fromPos][0])))
     self.allFid[toPos] = sp.multiply(self.allFid[fromPos][:], sp.exp(-self.fidTime[:len(self.allFid[fromPos][0])]*LB*np.pi))

def __kullback_leibler(h1, h2): # 36.3 us
    """
    The actual KL implementation. @see kullback_leibler() for details.
    Expects the histograms to be of type scipy.ndarray.
    """
    result = h1.astype(scipy.float_)
    mask = h1 != 0
    result[mask] = scipy.multiply(h1[mask], scipy.log(h1[mask] / h2[mask]))
    return scipy.sum(result)

def __minowski_low_positive_integer_p(h1, h2, p = 2): # 11..43 us for p = 1..24 \w 100 bins
    """
    A faster implementation of the Minowski distance for positive integer < 25.
    @note do not use this function directly, but the general @link minowski() method.
    @note the passed histograms must be scipy arrays.
    """
    mult = scipy.absolute(h1 - h2)
    dif = mult
    for _ in range(p - 1): dif = scipy.multiply(dif, mult)
    return math.pow(scipy.sum(dif), 1./p)