def __call__(self, values, clip=True, out=None):
     values = _prepare(values, clip=clip, out=out)
     np.multiply(values, self.exp, out=values)
     np.add(values, 1., out=values)
     np.log(values, out=values)
     np.true_divide(values, np.log(self.exp + 1.), out=values)
     return values

  def movementCompute(self, displacement, noiseFactor = 0):
    """
    Shift the current active cells by a vector.

    @param displacement (pair of floats)
    A translation vector [di, dj].
    """

    if noiseFactor != 0:
      displacement = copy.deepcopy(displacement)
      xnoise = np.random.normal(0, noiseFactor)
      ynoise = np.random.normal(0, noiseFactor)
      displacement[0] += xnoise
      displacement[1] += ynoise


    # Calculate delta in the module's coordinates.
    phaseDisplacement = (np.matmul(self.rotationMatrix, displacement) *
                         self.phasesPerUnitDistance)

    # Shift the active coordinates.
    np.add(self.activePhases, phaseDisplacement, out=self.activePhases)

    # In Python, (x % 1.0) can return 1.0 because of floating point goofiness.
    # Generally this doesn't cause problems, it's just confusing when you're
    # debugging.
    np.round(self.activePhases, decimals=9, out=self.activePhases)
    np.mod(self.activePhases, 1.0, out=self.activePhases)

    self._computeActiveCells()
    self.phaseDisplacement = phaseDisplacement

def morphological_laplace(input, size = None, footprint = None,
                          structure = None, output = None,
                          mode = "reflect", cval = 0.0, origin = 0):
    """Multi-dimensional morphological laplace.

    Either a size or a footprint, or the structure must be provided. An
    output array can optionally be provided. The origin parameter
    controls the placement of the filter. The mode parameter
    determines how the array borders are handled, where cval is the
    value when mode is equal to 'constant'.
    """
    tmp1 = grey_dilation(input, size, footprint, structure, None, mode,
                         cval, origin)
    if isinstance(output, numpy.ndarray):
        grey_erosion(input, size, footprint, structure, output, mode,
                     cval, origin)
        numpy.add(tmp1, output, output)
        del tmp1
        numpy.subtract(output, input, output)
        return numpy.subtract(output, input, output)
    else:
        tmp2 = grey_erosion(input, size, footprint, structure, None, mode,
                            cval, origin)
        numpy.add(tmp1, tmp2, tmp2)
        del tmp1
        numpy.subtract(tmp2, input, tmp2)
        numpy.subtract(tmp2, input, tmp2)
        return tmp2

def makeFeatureVec(words, model, num_features):

	# Pre-initialize an empty numpy array (for speed)
	featureVec = np.zeros((num_features,),dtype="float32")

	# Count number of words
	nwords = 0.

	# Loop over word by word
	# If in vocabulary, add its feature vector to the total
	for word in words.split():

		if word in model: 
			nwords += 1.
			featureVec = np.add(featureVec,model[word])
		else:
			missingWord = handleMissingWord(word, model, num_features)
			featureVec = np.add(featureVec, missingWord)
			nwords += 1.

	# Divide the result by the number of words to get the average
	featureVec = np.divide(featureVec,nwords)
	
	# If number of words zero
	if nwords == 0:
		featureVec = characterVec(words, model, num_features)
	
	return featureVec

    def test_pi_ops_nat(self):
        idx = PeriodIndex(['2011-01', '2011-02', 'NaT', '2011-04'],
                          freq='M', name='idx')
        expected = PeriodIndex(['2011-03', '2011-04', 'NaT', '2011-06'],
                               freq='M', name='idx')

        self._check(idx, lambda x: x + 2, expected)
        self._check(idx, lambda x: 2 + x, expected)
        self._check(idx, lambda x: np.add(x, 2), expected)

        self._check(idx + 2, lambda x: x - 2, idx)
        self._check(idx + 2, lambda x: np.subtract(x, 2), idx)

        # freq with mult
        idx = PeriodIndex(['2011-01', '2011-02', 'NaT', '2011-04'],
                          freq='2M', name='idx')
        expected = PeriodIndex(['2011-07', '2011-08', 'NaT', '2011-10'],
                               freq='2M', name='idx')

        self._check(idx, lambda x: x + 3, expected)
        self._check(idx, lambda x: 3 + x, expected)
        self._check(idx, lambda x: np.add(x, 3), expected)

        self._check(idx + 3, lambda x: x - 3, idx)
        self._check(idx + 3, lambda x: np.subtract(x, 3), idx)

def numeric_gemm_var1_flat(A, B, C, mc, kc, nc, mr=1, nr=1):
  M, N = C.shape
  K = A.shape[0]

  mc = min(mc, M)
  kc = min(kc, K)
  nc = min(nc, N)  

  tA = numpy.zeros((mc, kc), dtype = numpy.float)
  tB = numpy.zeros((kc, N), dtype = numpy.float)  

  for k in range(0, K, kc):
    # Pack B into tB
    tB[:,:] = B[k:k+kc:,:]
    
    for i in range(0, M, mc):
      imc = i+mc
      # Pack A into tA
      tA[:,:] = A[i:imc,k:k+kc]
      
      for j in range(0, N): # , nc):
        # Cj += ABj + Cj
        # jnc = j+nc
        ABj = numpy.matrixmultiply(tA, tB[:,j])
        numpy.add(C[i:imc:,j], ABj, C[i:imc:,j])
        
        # Store Caux into memory
  return

def wavefunction(coords, mocoeffs, gbasis, volume):
    """Calculate the magnitude of the wavefunction at every point in a volume.
    
    Attributes:
        coords -- the coordinates of the atoms
        mocoeffs -- mocoeffs for one eigenvalue
        gbasis -- gbasis from a parser object
        volume -- a template Volume object (will not be altered)
    """
    bfs = getbfs(coords, gbasis)
    
    wavefn = copy.copy(volume)
    wavefn.data = numpy.zeros( wavefn.data.shape, "d")

    conversion = convertor(1,"bohr","Angstrom")
    x = numpy.arange(wavefn.origin[0], wavefn.topcorner[0]+wavefn.spacing[0], wavefn.spacing[0]) / conversion
    y = numpy.arange(wavefn.origin[1], wavefn.topcorner[1]+wavefn.spacing[1], wavefn.spacing[1]) / conversion
    z = numpy.arange(wavefn.origin[2], wavefn.topcorner[2]+wavefn.spacing[2], wavefn.spacing[2]) / conversion

    for bs in range(len(bfs)):
        data = numpy.zeros( wavefn.data.shape, "d")
        for i,xval in enumerate(x):
            for j,yval in enumerate(y):
                for k,zval in enumerate(z):
                    data[i, j, k] = bfs[bs].amp(xval,yval,zval)
        numpy.multiply(data, mocoeffs[bs], data)
        numpy.add(wavefn.data, data, wavefn.data)
    
    return wavefn

def average_perceptron(feature_matrix, labels, T):
    theta = np.empty_like(feature_matrix[0])
    theta.fill(0.)
    theta_sum = theta  
    theta_0 = 0.0
    theta_0_sum = theta_0
    ticker = 0
    update_track = 0
    
    while ticker < T:
        
        for i in range(len(feature_matrix)):        

            check_before_label = np.add(np.dot(theta, feature_matrix[i]),theta_0)

            check_mult_label = np.multiply(labels[i], check_before_label)
            if check_mult_label == 0 or check_mult_label < 0:
                update_track += 1                
                (theta, theta_0) = perceptron_single_step_update(feature_matrix[i], labels[i], theta, theta_0)
                theta_sum = np.add(theta, theta_sum)
                theta_0_sum += theta_0

        ticker += 1
        
    theta_average = np.divide(theta_sum, update_track)
    theta_0_average = theta_0_sum/update_track
        
    return (theta_average, theta_0_average)

def _wrapx(input, output, nx):
    """
    Wrap the X format column Boolean array into an ``UInt8`` array.

    Parameters
    ----------
    input
        input Boolean array of shape (`s`, `nx`)

    output
        output ``Uint8`` array of shape (`s`, `nbytes`)

    nx
        number of bits
    """

    output[...] = 0  # reset the output
    nbytes = ((nx - 1) // 8) + 1
    unused = nbytes * 8 - nx
    for i in range(nbytes):
        _min = i * 8
        _max = min((i + 1) * 8, nx)
        for j in range(_min, _max):
            if j != _min:
                np.left_shift(output[..., i], 1, output[..., i])
            np.add(output[..., i], input[..., j], output[..., i])

    # shift the unused bits
    np.left_shift(output[..., i], unused, output[..., i])

    def discretize(self, time_slice_length):
        self.time_slice_length = time_slice_length

        # compute the total number of time-slices
        time_delta = (self.end_date - self.start_date)
        time_delta = time_delta.total_seconds()/60
        self.time_slice_count = int(time_delta // self.time_slice_length) + 1

        # parallelize tweet partitioning using a pool of processes (number of processes = number of cores).
        nb_processes = cpu_count()
        nb_tweets_per_process = self.size // nb_processes
        portions = []
        for i in range(0, self.size, nb_tweets_per_process):
            j = i + nb_tweets_per_process if i + nb_tweets_per_process < self.size else self.size
            portions.append((i, j))
        p = Pool()
        results = p.map(self.discretize_job, portions)
        results.sort(key=lambda x: x[0])

        # insert the time-slices number in the data frame and compute the final frequency matrices
        time_slices = []
        self.tweet_count = np.zeros(self.time_slice_count, dtype=np.int)
        self.global_freq = csr_matrix((len(self.vocabulary), self.time_slice_count), dtype=np.short)
        self.mention_freq = csr_matrix((len(self.vocabulary), self.time_slice_count), dtype=np.short)
        for a_tuple in results:
            time_slices.extend(a_tuple[1])
            self.tweet_count = np.add(self.tweet_count, a_tuple[2])
            self.global_freq = np.add(self.global_freq, a_tuple[3])
            self.mention_freq = np.add(self.mention_freq, a_tuple[4])
        self.df['time_slice'] = np.array(time_slices)

 def __add__(self, other):
     if isinstance(other, Raster):
         result = np.add(self.data, other.data)
     else:
         result = np.add(self.data, other)
     
     return Raster(None, result, self.nodata, self.driver, self.georef, self.proj)

def plot_selfish_cooperative(num_runs):
	import seaborn as sns
	for exp in range(num_runs):
		fig  = plt.figure()
		cooperators = [3,4,5]
		selfish = [0,1,2]

		
		fname = Parameters.dirname + '/%i_assembly_%i.dat' % (exp, 0)
		t, pop = np.loadtxt(fname, unpack = True)

		cooperative_pop = [0.0]*len(pop[:-5])
		selfish_pop = [0.0]*len(pop[:-5])

		for c in cooperators:
			fname = Parameters.dirname + '/%i_assembly_%i.dat' % (exp, c)
			t, pop = np.loadtxt(fname, unpack = True)
			cooperative_pop = np.add(cooperative_pop,pop[:-5] )

		for s in selfish:
			fname = Parameters.dirname + '/%i_assembly_%i.dat' % (exp, s)
			t, pop = np.loadtxt(fname, unpack = True)
			selfish_pop = np.add(selfish_pop,pop[:-5] )

		ax = fig.add_subplot(1,1,1)
		ax.plot(t[:-5], selfish_pop, label = 'selfish', color = 'r')
		ax.plot(t[:-5], cooperative_pop, label = 'cooperative', color = 'g')
		ax.legend(loc = 'upper left')
		plt.xlabel('System Time')
		plt.ylabel('Total Abundance')

		#plt.show()
		plt.savefig(Parameters.dirname + '/%i_cooperative_vs_selfish.png' % exp)
		plt.close()

 def createChord(self,*args):
   if len(args) == 1:
     self.notesCombined = args[0]
   if len(args) == 2:
     self.notesCombined = np.add(args[0],args[1])
   if len(args) == 3:
     self.notesCombined = np.add(args[0],np.add(args[1],args[2]))

	def __init__(self, limitsLow, limitsHigh, status = None):
		# Either merge two clusters (bottom-up)
		# or create a new micro-cluster (top-down)
		if status is "merging":
			# The first two input parameters are two newCluster() objects
			first = limitsLow
			second = limitsHigh
			self.limitsLow = [None]*ndim
			self.limitsHigh = [None]*ndim
			for i in xrange(ndim):
				self.limitsLow[i] = min(first.limitsLow[i], second.limitsLow[i])
				self.limitsHigh[i] = max(first.limitsHigh[i], second.limitsHigh[i])
			self.findKeys()
			self.weight = first.weight + second.weight
			self.Sum = np.add(first.Sum, second.Sum)
			self.sqSum = np.add(first.sqSum, second.sqSum)#self.computeSqSum()
			self.computeSSQ()
			self.computeCoG()
		else:
			# The first two parameters are the edges of the cluster
			self.limitsLow = limitsLow
			self.limitsHigh = limitsHigh
			self.findKeys()
			self.computeWeight()
			self.computeSum()
			self.computeSqSum()
			self.computeSSQ()
			self.computeCoG()

def scale_samples(params, bounds):
    '''
    Rescales samples in 0-to-1 range to arbitrary bounds.

    Arguments:
        bounds - list of lists of dimensions num_params-by-2
        params - numpy array of dimensions num_params-by-N,
        where N is the number of samples
    '''
    # Check bounds are legal (upper bound is greater than lower bound)
    b = np.array(bounds)
    lower_bounds = b[:, 0]
    upper_bounds = b[:, 1]

    if np.any(lower_bounds >= upper_bounds):
        raise ValueError("Bounds are not legal")

    # This scales the samples in-place, by using the optional output
    # argument for the numpy ufunctions
    # The calculation is equivalent to:
    #   sample * (upper_bound - lower_bound) + lower_bound
    np.add(np.multiply(params,
                       (upper_bounds - lower_bounds),
                       out=params),
           lower_bounds,
           out=params)

def plot3d(data):
	windowIndex = 0
	msgrateIndex = 1
	intervalLengthIndex = 3
	dirCPUTOTALIndex=11+6 -2
	DecCPUTotalIndex=23+6 -2
	PunCPUToTalIndex=35+6 -2 
	MaxSize=36+6


	intervals = data[:,intervalLengthIndex]
	dirTimes = data[:,dirCPUTOTALIndex]
	decTime = data[:,DecCPUTotalIndex]
	PunTime = data[:,PunCPUToTalIndex]
	totalTime = np.add(dirTimes,decTime)
	totalTime = np.add(totalTime,PunTime)
	rate = data[:,intervalLengthIndex]

	sizes = data[:,MaxSize]

	fig = plt.figure()
	ax = fig.add_subplot(111, projection='3d')

	ax.scatter(rate,intervals,totalTime)
	ax.set_xlabel(' Message rate (msg/s)')
	ax.set_ylabel('interval')
	ax.set_zlabel('time (s)')
	ax.xaxis.set_scale('log')
	ax.yaxis.set_scale('log')

	plt.show()	

def trainClassifier(postingVec, classVec):
    nWGivenC0 = np.zeros(len(postingVec[0]))
    nWGivenC1 = np.zeros(len(postingVec[0]))
    nWs = np.zeros(len(postingVec[0]))
    numAllWsC0 = 0
    numAllWsC1 = 0
    numC1 = 0
    index = 0
    for post in postingVec:
        nWs = np.add(nWs, np.array(post))
        if (classVec[index] == 0):
            nWGivenC0 = np.add(nWGivenC0, np.array(post))
            numAllWsC0 += sum(post)
        else:
            nWGivenC1 = np.add(nWGivenC1, np.array(post))
            numAllWsC1 += sum(post)
            numC1 += 1
        index += 1

    pWGivenC0 = nWGivenC0/numAllWsC0 # probability of each word, given class C0
    pWGivenC1 = nWGivenC1/numAllWsC1 # probability of each word, given class C1

    pC1 = float(numC1) / len(classVec) # probability of class 1
    pC0 = 1 - pC1 # probability of class 0

    return (pC0,pWGivenC0), (pC1,pWGivenC1)

def open_cfd(cfd, which_window):
    if which_window == "watchlists":
        win = watchlists
    elif which_window == "open_positions":
        win = open_positions
    else:
        raise Exception("Window not defined")

    pos = wait_for_window(win, 1)
    if pos:
        screen = load_prev_screenshot()
        moveclick(*add(pos, (5, 5)))  # To activate Window (if not)
    else:
        print "ERROR: " + which_window + " not found"
        return False

    win_pos = pos
    maxval, pos = get_match(screen[win_pos[1] : win_pos[1] + 200, win_pos[0] : win_pos[0] + 400], cfd_list[cfd], 1)
    pos = add(win_pos, pos)
    if maxval > MATCH_OK:
        moveclick(*add(pos, (5, 5)))
    else:
        print "CFD " + cfd + " not found in Screen"
        return False

    return True

def feature_extraction_Doc2Vec(data_pos, data_neg, model): # use the word2vec under the hood
    vec_size = 300
    data_pos_vec, data_neg_vec = [], []
    for term in data_pos:
        raw_vecs = np.zeros(vec_size)
        vec_num = 0
        for item in term:
            try:
                raw_vecs = np.add(raw_vecs, model[item])
                vec_num += 1
            except:
                pass
        doc_vec = raw_vecs / float(vec_num+0.00001)
        data_pos_vec.append(doc_vec.tolist())
    for term in data_neg:
        raw_vecs = np.zeros(vec_size)
        vec_num = 0
        for item in term:
            try:
                raw_vecs = np.add(raw_vecs, model[item])
                vec_num += 1
            except:
                pass
        doc_vec = raw_vecs / float(vec_num + 0.00001)
        data_neg_vec.append(doc_vec.tolist())
    return data_pos_vec, data_neg_vec

	def __step(self,vis=False):
		## liczenie natezenia w kazdym punkcie
		for to_mod_dipole in self.dlist:
			to_mod_dipole.E1=np.array([0.,0.,0.])
			to_mod_dipole.E2=np.array([0.,0.,0.])
			for get_dipole in self.dlist:
				if not to_mod_dipole == get_dipole:
					to_mod_dipole.E1=np.add(to_mod_dipole.E1,GetElectricField(to_mod_dipole.r1,get_dipole))
					to_mod_dipole.E2=np.add(to_mod_dipole.E2,GetElectricField(to_mod_dipole.r2,get_dipole))
					
		## poddajemy dipole dzialaniom sil 
		for i in range(len(self.dlist)):
			if IsInWorld(self.dlist[i],self.world):		# sprawdzenie czy dipol jest w swiecie
				LetElectricForceWork(self.dlist[i],self.dt)
			else:
				raise 0	# jesli wyszedl poza swiat - koniec sym.

		## wizualizacja stepu
		if vis: 
			print "Time: "+str(self.time)+"  Refreshing..."
			self.visdata.NextFrame(self.dlist)



		## zapis danych do pliku 
		self.fobj.write(str(self.time)+"\n")