def step_generation(self, senders, receivers):
        # x_i(t+1) = (a + u(e^i, x(t)))*x_i(t) / (a + u(x(t), x(t)))
        # a is background (lifetime) birthrate -- set to 0

        s_payoffs = self._data['s_payoffs']
        r_payoffs = self._data['r_payoffs']
        s_fitness = [0.] * len(senders)
        r_fitness = [0.] * len(receivers)

        for (s, sp), (r, (rp, rt)) in itertools.product(enumerate(senders), enumerate(receivers)):
            state_acts = self._interactions[(s, r)]
            s_fitness[s] += math.fsum(s_payoffs[state][act] * rp for state, act in state_acts) / 4.
            r_fitness[r] += math.fsum(r_payoffs[rt][state][act] * sp for state, act in state_acts) / 4.

        avg_s = math.fsum(s_fitness[s] * sp for s, sp in enumerate(senders))
        avg_r = math.fsum(r_fitness[r] * rp for r, (rp, rt) in enumerate(receivers))

        new_senders = [s_fitness[s] * sp / avg_s for s, sp in enumerate(senders)]
        new_receivers = [(r_fitness[r] * rp / avg_r, rt) for r, (rp, rt) in enumerate(receivers)]

        for s, sp in enumerate(new_senders):
            if sp < effective_zero:
                new_senders[s] = 0.
        for r, (rp, rt) in enumerate(new_receivers):
            if rp < effective_zero:
                new_receivers[r] = (0., rt)

        return (tuple(new_senders), tuple(new_receivers))

 def calc_ase(self, mod_t, exp_t, args):
     """
     Calculates the normalized average squared difference of the traces.
     
     :param mod_t: the trace obtained from the model as ``list``
     :param exp_t: the input trace as ``list``
     :param args: optional arguments as ``dictionary``
     
     :return: the normalized average squared difference, where the normalization is done by
         the squared range of the input trace
         
     """
     if (args["cov_m"]!=None):
             return self.calc_ase_cov(mod_t,exp_t,args)
     temp = []
     for n in range(min([len(exp_t), len(mod_t)])):
         try:
             temp.append(pow(exp_t[n] - mod_t[n], 2))
         except OverflowError:
             return 1
         #except TypeError:
         #    return 1
     try:
         if self.option.output_level == "1":
             print "ase"
             print fsum(temp) / len(temp) / (pow(max(exp_t) - min(exp_t), 2))
     except OverflowError:
             return 1
     return fsum(temp) / len(temp) / (pow(max(exp_t) - min(exp_t), 2))

def get_device_local_storage_price(device):
    price = math.fsum(s.get_price() for s in device.storage_set.all())
    if not price and device.model and device.model.type in (
            DeviceType.rack_server.id, DeviceType.blade_server.id,
            DeviceType.virtual_server.id):
        try:
            os = OperatingSystem.objects.get(device=device)
            group = ComponentModelGroup.objects.get(name='OS Detected Storage')
        except (OperatingSystem.DoesNotExist,
                ComponentModelGroup.DoesNotExist):
            pass
        else:
            if not group.per_size:
                return group.price or 0
            else:
                storage = os.storage or 0
                remote_storage_size = math.fsum(
                    m.get_size() for m in device.disksharemount_set.all()
                )
                storage -= remote_storage_size
                if storage > 0:
                    return (storage /
                            (group.size_modifier or 1)) * (group.price or 0)
        if device.model.type != DeviceType.virtual_server.id:
            try:
                group = ComponentModelGroup.objects.get(name='Default Disk')
            except ComponentModelGroup.DoesNotExist:
                pass
            else:
                return group.price
    return price

 def yules_k(self):
     freqs = self.self_wordfreq()
     freq_set = list(set(freqs))
     M1 = math.fsum([freqs.count(f)*f for f in freq_set])
     M2 = math.fsum([(f**2)*freqs.count(f) for f in freq_set])
     K = (10000)*(M2 - M1)/(M1**2)
     return K

    def calc_grad_dif(self, mod_t, exp_t, args):
        """
        Calculates the normalized average squared differences of derivatives of the given traces.
        The gradient is calculated as follows:
        ::

            grad_a=((mod_t[i+1]-mod_t[i-1])/(2*dt))

        where dt is the step between to points in the trace
        
        :param mod_t: the trace obtained from the model as ``list``
        :param exp_t: the input trace as ``list``
        :param args: optional arguments as ``dictionary``
        
        :return: the normalized average squared differences of derivatives where the normalization is done
            by the squared range of the input trace    
        """
        dt = self.reader.data.step
        grad_a = 0
        grad_b = 0
        tmp = []
        for i in range(1, min(len(mod_t), len(exp_t)) - 1):
            grad_a = ((mod_t[i + 1] - mod_t[i - 1]) / (2 * dt))
            grad_b = ((exp_t[i + 1] - exp_t[i - 1]) / (2 * dt))
            tmp.append((grad_a - grad_b) ** 2)
        try:
            if self.option.output_level == "1":
                print "grad dif"
                print fsum(tmp) / len(tmp) / (pow(max(grad_b) - min(grad_b), 2))
        except OverflowError:
                return 1
            
            
        return fsum(tmp) / len(tmp) / (pow(max(grad_b) - min(grad_b), 2))  

def simpsons(a,b,n):
    Xk = (float) (b-a)/n
    x2 = arange(a+2*Xk, b-Xk, 2*Xk) # 3rd, 5th, 7th... (n-2)th x values
    x4 = arange(a+Xk, b, 2*Xk) # 2nd, 4th, 6th... (n-1)th x values
    sum_of_2s = 2*fsum(f(x2))
    sum_of_4s = 4*fsum(f(x4))
    return Xk/3*(f(a) + sum_of_2s + sum_of_4s + f(b))

    def weightedOverhead(self, numpy=False):
        ns = 'weighted.numpy' if numpy else 'weighted'
        factors = {
            'strConcat': 0.8171743283710745,
            'objectInit': 162.78282024428898,
            'getListItem': 12.018827808252004,
            'functionCall': 48.613475069418904,
            'getDictItem': 0.2148913920265531,
            'methodCall': 75.36797944163118,
            'numpy': 1.0,
            }
        runtimes = {
            'strConcat': self.totalStrConcatRuntime(),
            'objectInit': self.totalObjectCreations() * PYTHON_OBJECT_CREATION_COST,
            'functionCall': self.totalPythonCalls() * PYTHON_CALL_COST,
            'getListItem': self.totalGetitemRuntime(),
            }
        if numpy:
            runtimes['numpy'] = self.totalNumpyRuntime()

        total_runtime = fsum(runtimes.values())
        c_runtimes = {key:(runtime / factors[key]) for key, runtime in runtimes.iteritems()}
        total_c_runtime = fsum(c_runtimes.values())

        self.store(ns, 'total_tt', self.total_tt)
        self.store(ns, 'total_runtime', total_runtime)
        self.store(ns, 'total_c_runtime', total_c_runtime)
        self.store(ns, 'runtimes', runtimes)
        self.store(ns, 'factors', factors)
        self.store(ns, 'c_runtimes', c_runtimes)

        return 1 - (total_c_runtime / total_runtime), (total_runtime / self.total_tt)

def genStats(times):
    N = len(times)
    avgTime = math.fsum(times)/ N
    timesMinusMean = [x - avgTime for x in times]
    timesMMSquared = [math.pow(x,2) for x in timesMinusMean]
    var = math.fsum(timesMMSquared) / N
    return avgTime, var

def _generalised_sum(data, func):
    """_generalised_sum(data, func) -> len(data), sum(func(items of data))

    Return a two-tuple of the length of data and the sum of func() of the
    items of data. If func is None, use just the sum of items of data.
    """
    # Try fast path.
    try:
        count = len(data)
    except TypeError:
        # Slow path for iterables without len.
        # We want to support BIG data streams, so avoid converting to a
        # list. Since we need both a count and a sum, we iterate over the
        # items and emulate math.fsum ourselves.
        ap = add_partial
        partials = []
        count = 0
        if func is None:
            # Note: we could check for func is None inside the loop. That
            # is much slower. We could also say func = lambda x: x, which
            # isn't as bad but still somewhat expensive.
            for count, x in enumerate(data, 1):
                ap(x, partials)
        else:
            for count, x in enumerate(data, 1):
                ap(func(x), partials)
        total = math.fsum(partials)
    else: # Fast path continues.
        if func is None:
            # See comment above.
            total = math.fsum(data)
        else:
            total = math.fsum(func(x) for x in data)
    return count, total

def collapse_duplicates(raw_data):
	# Create dictionary of lists of duplicates
	dup_data = raw_data.get_array()
	set_sp = {}
	set_ab = {}
	set_co = {}
	set_sz = {}
	set_plasmids = {}
	for sp,ab,co in dup_data:
		if 'taxid' in sp: # retain useful information
			name = sp.rpartition('|')[0] # last segment is usually the original chromosome etc name
		else:
			name = sp.partition('_gi|')[0].partition('|')[0].partition('_gca')[0] #the prepended strain name
		set_sp.setdefault(name,[]).append(sp)
		set_ab.setdefault(name,[]).append(ab)
		set_co.setdefault(name,[]).append(co)
		set_sz.setdefault(name,[]).append(co)

	assert(set_ab.keys() == set_co.keys() == set_sp.keys())

	# New, clean dataset for data without duplicates
	undupe = Dataset()

	# Note: we include plasmids in the count total solely because i100 was simulated to include 1x plasmid coverage.
	for k,v in set_sp.items():
		if len(v) == 1: # just add record directly if it has no duplicates
			undupe.add_record(k,set_ab[k][0],set_co[k][0],set_sz[k][0])
		else: # sum counts and average abundances
			undupe.add_record(k,math.fsum(set_ab[k])/len(v),math.fsum(set_co[k]),math.fsum(set_sz[k]))

	print "Number of entries after combining duplicates: {0}".format(len(undupe.species))

	return undupe

def stats(results, optimal=None):
    sum, notfound, worst = first_stats(results)
    avg = sum / len(results)
    varianza = fsum([(x - avg) ** 2 for x in results]) / len(results)
    scarto = fpformat.fix(sqrt(varianza), 2)
    valori = set(results)
    frequenze = dict(zip(valori, [results.count(v) for v in valori]))
    sorted_frequenze = sorted(frequenze, key=frequenze.get, reverse=True)
    sorted_frequenze = sorted_frequenze[:10]
    if optimal:
        opt_sum, opt_nf, opt_worst = first_stats(optimal)
        opt_avg = opt_sum / len(optimal)
        opt_scarto = fpformat.fix(sqrt(fsum([(x - opt_avg) ** 2 for x in optimal]) / len(optimal)), 2)
        ratio_avg = avg / opt_avg
        ratio_worst = worst / opt_worst
        ratio_scarto = fpformat.fix((float(scarto) / float(opt_scarto)), 2)

    print "-------------------------------------------------"
    print "Statistiche:\t\t\t\tOffline\tOnline\tRapporto"
    print "Numero di test eseguiti:\t\t " + str(len(results)) + "\t" + str(len(optimal))
#    print "Carburante esaurito:\t\t\t " + str(notfound)
    print "Caso peggiore:\t\t\t\t " + str(worst) + "\t" + str(opt_worst) + "\t" + str(ratio_worst)
    print "Media aritmetica dei risultati:\t\t " + str(avg) + "\t" + str(opt_avg) + "\t" + str(ratio_avg)
    print "Scarto quadratico medio:\t\t " + str(scarto) + "\t" + str(opt_scarto) + "\t" + str(ratio_scarto)
    print "I dieci risultati piu' riscontrati:"
    print "Costo:\tOttenuto:\tSotto la media?"
    for el in sorted_frequenze:
        sotto = "media"
        if el < avg:
            sotto = "si"
        elif el > avg:
            sotto = "no"
        print str(el) + "\t" + str(frequenze[el]) + "\t\t" + sotto

def render_summary(stats):
    """
    Render summary of an event stream run.
    :param stats: Dictionary('clock':list()<float>, 'rss':list()<int>)
    :return: Void.
    """
    print('\nSummary profile from stream execution:')
    print('Samples: %i' % len(stats['clock']))
    if -1 in stats['clock']:
        print('(ERRORS DETECTED: Removing timing samples from aborted invocations.)')
        stats['clock'] = [x for x in stats['clock'] if x > 0]
        print('New sample size: %i' % len(stats['clock']))
    median = sorted(stats['clock'])[math.trunc(len(stats['clock']) / 2)]
    print(stats['clock'])
    mean =  math.fsum(stats['clock'])/len(stats['clock'])
    print('Clock time:\n'
          '\tMin: %ims, Max: %ims, Median: %ims, Median Billing Bucket: %ims, Rounded Standard Deviation: %sms' % (
              min(stats['clock']),
              max(stats['clock']),
              median,
              billing_bucket(median),
              math.trunc(math.ceil(math.sqrt(math.fsum((x-mean)**2 for x in stats['clock'])/(len(stats['clock'])-1))))
          )) if len(stats['clock']) > 0 else print("No valid timing samples!")
    print('Peak resident set size (memory):\n'
          '\tMin: %s, Max: %s' % (
              size(min(stats['rss'])),
              size(max(stats['rss']))
          ))

def inTheGrandSchemeOfThings(nrNodesMin, nrNodesMax, algorithm):
    listOfAverage = [[0.0], []]
    listOfMaximum = [[0.0], []]
    listOfMedian = [[0.0], []]

    for nrNodes in range(nrNodesMin, nrNodesMax, 2):
        listOfNodes = setupDesign(nrNodes, algorithm)
        maxLevel = len(listOfNodes) - 1
        while len(listOfAverage) <= maxLevel and len(listOfMaximum) <= maxLevel and len(listOfMedian) <= maxLevel:
            listOfAverage.append([])
            listOfMaximum.append([])
            listOfMedian.append([])
        for level in range(1, maxLevel + 1):
            average = Measurements.averageDistance(listOfNodes, level)
            listOfAverage[level].append(average / level)
            maximum = Measurements.maximumDistance(listOfNodes, level)
            listOfMaximum[level].append(maximum / level)
            median = Measurements.distanceMedian(listOfNodes, level)
            listOfMedian[level].append(median / level)
    for i in range(1, len(listOfAverage)):
        listOfAverage[i] = math.fsum(listOfAverage[i]) / len(listOfAverage[i])
        listOfMaximum[i] = math.fsum(listOfMaximum[i]) / len(listOfMaximum[i])
        listOfMedian[i] = math.fsum(listOfMedian[i]) / len(listOfMedian[i])
    del listOfAverage[0]
    del listOfMaximum[0]
    del listOfMedian[0]
    plt.plot(range(1, 1 + len(listOfAverage)), listOfAverage, label = 'Average')
    plt.plot(range(1, 1 + len(listOfMaximum)), listOfMaximum, label = 'Maximum')
    plt.plot(range(1, 1 + len(listOfMedian)), listOfMedian, label = 'Median')
    plt.legend(loc=3)
    plt.axis([0, maxLevel, 0, 1.1])
    plt.show()

def compute_volume(mesh):
    if "tetra" in mesh.cells:
        vol = math.fsum(
            get_simplex_volumes(*prune_nodes(mesh.points, mesh.cells["tetra"]))
        )
    elif "triangle" in mesh.cells or "quad" in mesh.cells:
        vol = 0.0
        if "triangle" in mesh.cells:
            # triangles
            vol += math.fsum(
                get_triangle_volumes(*prune_nodes(mesh.points, mesh.cells["triangle"]))
            )
        if "quad" in mesh.cells:
            # quad: treat as two triangles
            quads = mesh.cells["quad"].T
            split_cells = numpy.column_stack(
                [[quads[0], quads[1], quads[2]], [quads[0], quads[2], quads[3]]]
            ).T
            vol += math.fsum(
                get_triangle_volumes(*prune_nodes(mesh.points, split_cells))
            )
    else:
        assert "line" in mesh.cells
        segs = numpy.diff(mesh.points[mesh.cells["line"]], axis=1).squeeze()
        vol = numpy.sum(numpy.sqrt(numpy.einsum("...j, ...j", segs, segs)))

    return vol

def get_init_nation_data(nation_set, all_games):
    nation_data = {}
    for nation in nation_set:
        nation_data[nation] = {"games":[],"GF":[],"GA":[],"weight":[], \
            "Off":0, "Def":0, "error":1}
    total_goals = []
    total_weight = []
    for game in all_games:
        team1 = game["team1"]
        team2 = game["team2"]
        nation_data[team1]["games"].append(game)
        nation_data[team1]["GF"].append(game["goals1"]*game["weight"])
        nation_data[team1]["GA"].append(game["goals2"]*game["weight"])
        nation_data[team1]["weight"].append(game["weight"])
        nation_data[team2]["games"].append(game)
        nation_data[team2]["GF"].append(game["goals2"]*game["weight"])
        nation_data[team2]["GA"].append(game["goals1"]*game["weight"])
        nation_data[team2]["weight"].append(game["weight"])
        total_goals.append((game["goals1"]+game["goals2"]) * game["weight"])
        total_weight.append(game["weight"])
    
    for nation in nation_set:
        # initial estimates of Offense/Defense are GF/GA divided by weight
        nation_data[nation]["GF"] = math.fsum(nation_data[nation]["GF"])
        nation_data[nation]["GA"] = math.fsum(nation_data[nation]["GA"])
        nation_data[nation]["Offense"] = nation_data[nation]["GF"]
        nation_data[nation]["weight"]=math.fsum(nation_data[nation]["weight"])
        nation_data[nation]["Defense"] = nation_data[nation]["GA"]
        nation_data[nation]["Offense"] /= nation_data[nation]["weight"]
        nation_data[nation]["Defense"] /= nation_data[nation]["weight"]
        
    total_goals = math.fsum(total_goals)
    total_weight = math.fsum(total_weight)
    return nation_data, total_goals/total_weight/2 # PER TEAM

def calcDistance(*coords):
    """    Calculate distance according to :
    ((sum of all distances^-6)/number of distances)^-1/6
    or (sum of all distances^-6)^-1/6

    calcDistance.method should be set before use
    """
    result = None

    try:
        distance_list = (sqrt(fsum(sub(*coord) ** 2 for coord in izip(*atoms))) for atoms in product(*coords))
        sum6 = fsum(pow(distance, -6) for distance in distance_list)
        if calcDistance.method == 'ave6':
            number_of_distances = reduce(mul, (len(coord) for coord in coords))
        elif calcDistance.method == 'sum6':
            number_of_distances = 1
        result = pow(sum6/number_of_distances, -1./6)
    except(ValueError, TypeError):
        errors.add_error_message("Problem using coordinates : " +
                                 str(coords) + "\n" +
                                 " and distances list : " +
                                 str([distance for distance in distance_list]) + "\n")
    except AttributeError:
        sys.stderr.write("Please set calcDistance.method before using calcDistance()\n")

    return result

    def P_Nbpi(self, obs, Nbpi):
        top=PoissonApprox2(obs.Nmi,obs.Ci*Nbpi)
        #print(top, obs.Nmi, obs.Ci, Nbpi)

        if obs.Nmi>5:
            bottom=1.0/(obs.Ci)
        else:
            poissonrange=(obs.Npi*3+obs.Nmi*3+30)
            #memoize the bottom part
            if obs.Nmi in self.memoizedict2:
                if obs.Ci*Nbpi in self.memoizedict2[obs.Nmi]:
                    bottom=self.memoizedict2[obs.Nmi][obs.Ci*Nbpi]
                else:
                    self.memoizedict2[obs.Nmi][obs.Ci*Nbpi]=math.fsum(\
                        [PoissonApprox2(obs.Nmi,obs.Ci*Nbpi2)\
                        for Nbpi2 in range(poissonrange)])
                    bottom=self.memoizedict2[obs.Nmi][obs.Ci*Nbpi]
            else:
                self.memoizedict2[obs.Nmi]={}
                self.memoizedict2[obs.Nmi][obs.Ci*Nbpi]=math.fsum(\
                    [PoissonApprox2(obs.Nmi,obs.Ci*Nbpi2)\
                    for Nbpi2 in range(poissonrange)])
                bottom=self.memoizedict2[obs.Nmi][obs.Ci*Nbpi]
        if bottom==0.0:
            print("%s/%s" % (top, bottom))
            print("Nbpi: %s" % Nbpi)
            print("Nmi: %s" % obs.Nmi)
            print("Ci: %s" % obs.Ci)
        prob=top/bottom

        return prob

	def test_Poisson_sum_times(self):
		print("Testing Poisson sum times")
		k_max=400000
		lmbda=100000
		print("k_max=%s, lmbda=%s" % (k_max,lmbda))
		gt1=time.clock()
		s=math.fsum([bayes.PoissonApprox2(x,lmbda) for x in range(k_max)])
		gt2=time.clock()
		print("Time: %s seconds" % (gt2-gt1,))
		k_max=400000
		lmbda=200000
		print("k_max=%s, lmbda=%s" % (k_max,lmbda))
		gt1=time.clock()
		s=math.fsum([bayes.PoissonApprox2(x,lmbda) for x in range(k_max)])
		gt2=time.clock()
		print("Time: %s seconds" % (gt2-gt1,))
		k_max=600000
		lmbda=200000
		print("k_max=%s, lmbda=%s" % (k_max,lmbda))
		gt1=time.clock()
		s=math.fsum([bayes.PoissonApprox2(x,lmbda) for x in range(k_max)])
		gt2=time.clock()
		print("Time: %s seconds" % (gt2-gt1,))
		k_max=800000
		lmbda=200000
		print("k_max=%s, lmbda=%s" % (k_max,lmbda))
		gt1=time.clock()
		s=math.fsum([bayes.PoissonApprox2(x,lmbda) for x in range(k_max)])
		gt2=time.clock()
		print("Time: %s seconds" % (gt2-gt1,))

def h2jpopx(t,v,maxj,maxv):
    j = np.arange(maxj + 1, dtype = float)  #Array of j values
    g_j = 2.0*j + 1.0
    E_j = h2ejx(v,maxj)                  #cm^-1

    nj = g_j*np.exp(-(E_j*h*c)/(kb*t)) #Herzberg pg. 124

    #Need to properly normalize now (follow Herzberg pgs. 123-125)
    #For most astrophysically relevant temperatures, a max j of 100
    #should be adequate for getting the Q sum to converge

    highj = 25   #TOO MUCH HIGHER AND THE NUMBERS GET SO SMALL THEY ROLLOVER
    E_i = h2ejx(v,highj)
    i = np.arange(highj + 1, dtype = float)

    Qr = (2.0*i + 1.0)*np.exp(-(E_i*h*c)/(kb*t))
    njn = nj/math.fsum(Qr)

    #Get vibrational population (only relevant for large temperatures)
    E_v = np.zeros(maxv+1)

    for m in range(0,maxv+1):
        E_v[m] = h2ejx(m,0)      #cm^-1
    Qv = np.exp(-(E_v*h*c)/(kb*t))
    nvn = Qv/math.fsum(Qv)
    njn = njn*nvn[v]
    return njn

def kl_measures(sender_pop, receiver_pop, n=2):
    msgs = list(itertools.product(range(n), range(n)))
    states = list(itertools.product(range(n), range(n)))
    state_probs = [1. / float(len(states))] * len(states)

    all_cprobs_msg_given_state = collections.defaultdict(list)
    all_cprobs_state_given_msg = collections.defaultdict(list)
    information_contents = collections.defaultdict(list)
    for i, msg in enumerate(msgs):
        cprobs_msg_given_state = []
        for j, state in enumerate(states):
            pr = 0.
            for (sender, sender_prob) in sender_pop:
                if simulation.sender_matrix(sender)[j][i] == 1:
                    pr += sender_prob
            cprobs_msg_given_state.append(pr)
        all_cprobs_msg_given_state[i] = cprobs_msg_given_state

        for j, state in enumerate(states):
            if math.fsum(cprobs_msg_given_state) > 0.:
                prob_state_given_msg = ((state_probs[j] * cprobs_msg_given_state[j]) /
                                        math.fsum(state_probs[k] * cprobs_msg_given_state[k]
                                                    for k in xrange(len(states))))
            else:
                prob_state_given_msg = float('inf')

            all_cprobs_state_given_msg[j].append(prob_state_given_msg)

            if prob_state_given_msg > 0. and not math.isinf(prob_state_given_msg):
                information_contents[i].append(math.log(prob_state_given_msg / state_probs[j]))
            else:
                information_contents[i].append(- float('inf'))

    return (information_contents, all_cprobs_state_given_msg, all_cprobs_msg_given_state)