[LNT] Python 3 support: use Python 2's division behavior

[thopre]

Summary

Adapt all divisions to use a call to old_div to maintain Python 2's
division behavior in legacy code. Also import Python 3's division
behavior by default. This was produced by running futurize's stage2
libfuturize.fixes.fix_division_safe.

[hubert.reinterpretcast]

examples/functions.py

33 ↗

(On Diff #220980)

Clearly float division.

44 ↗

(On Diff #220980)

Same.

lnt/external/stats/pstat.py

141

Clearly floor division (len(source) // len(addon)). Also for the rest of the changes in the file.

lnt/external/stats/stats.py

295

Except for division by zero (which is equally bad between floor division and float division), this is float division.

343

There are a few instances in the patch where the change is made where integer division is clearly stated to be the intent, and many cases where it is clear which division operation would have taken place in Python 2 (or cases where the floor division occurs only in degenerate cases and float division is fine anyway).

[thopre]

Decide between float and integer division, getting rid of old_div

[thopre]

examples/functions.py

31 ↗

(On Diff #220980)

math.pi is float as pointed out by Hubert and thus this is float division

42 ↗

(On Diff #220980)

math.pi is float as pointed out by Hubert and thus this is float division

lnt/external/stats/pstat.py

141

len returns an integer as pointed out by Hubert so this is a floor division

146

len returns an integer as pointed out by Hubert so this is a floor division

153

len returns an integer as pointed out by Hubert so this is a floor division

158

len returns an integer as pointed out by Hubert so this is a floor division

lnt/external/stats/stats.py

295

As pointed by Hubert, either sum is 0 and thus both divisions give the same result, or sum is float by virtue of the 1.0/item and thus this should be float division.

343

len returns an int, therefore this is a floor division.

419

Second operand is a float due to being a number to the power 1.5, therefore this is a float division.

429

Denominator is a float due to being a number to the power 1.5, therefore this is a float division.

524

binsize is obtained by float division, therefore this is a float division.

700

math.sqrt returns a float so this is a float division.

710

mean refers to lmean in the same module which returns a float so this is a float division.

lnt/server/db/rules/rule_update_profile_stats.py

38

st_size being the size of f in bytes, it is an integer and thus this computes the floor.

lnt/server/reporting/analysis.py

83 ↗

(On Diff #220980)

This is float division:

the only construction of ComparisonResult (the type of self) that can be the source of pct_delta uses in LNT seems to have samples and prev_samples (and thus ultimately self.delta and value) made of floats (samples or result of calling calc_geomean)
the only uses of pct_delta are by flask (e.g. lnt/server/ui/templates/v4_global_status.html) and max, both of which are happy with float values
pct_delta is set to 0.0 when the condition of the outermost "if" is false

lnt/server/reporting/dailyreport.py

49

day_start_offset is of type datetime.timedelta whose second attribute is an integer, therefore this is computing the floor.

lnt/server/reporting/summaryreport.py

50 ↗

(On Diff #220980)

self.sum elements are initialized by 0. in _initialize and only ever changed by adding a value to it. Therfore this is float division.

68 ↗

(On Diff #220980)

value ** 1.0 will yield a float, so this is clearly float division.

85 ↗

(On Diff #220980)

While uses of NormalizedMean seems have values all come from calling lnt.util.stats.median and thus be integer, this makes more sense to be a float division because:

Mean._append uses the list elements in a float addition as explained above;
the same data that gets appended to NormalizedMean also gets appended to Mean in some cases

This suggests that the intent is to have float division.

505 ↗

(On Diff #220980)

values is constructed by calling getvalue and all aggregation functions have their value of type float. This is therefore a float division

lnt/server/ui/util.py

15–19 ↗

(On Diff #220980)

No use of safediv remains in LNT so this gets deleted in a separate patch.

73 ↗

(On Diff #220980)

This method is only referenced 3 times, all in lnt/server/ui/views.py where in all cases it ends up being called with data which can be either integer or float AFAIK ("column" field of SampleField):

• called once directly to fill the derived_stat attribute of the MatrixDataRequest which is not used anywhere
• as an aggregate function when constructing a ComparisonResult that is being used in get_cell_value macro in lnt/server/ui/templates/v4_matrix.html by printing it with the %.4f specifier

It should therefore be consider float division, which makes sense for a mean.

315 ↗

(On Diff #220980)

Like for ComparisonResult. all uses of pct_delta are fine with float and the initialization of pct_delta above suggest this is the intent.

lnt/server/ui/views.py

981 ↗

(On Diff #220980)

values in samples comes from the column field of SampleField which it appears can be either integer or float. However this should be made a float division unconditionally IMHO. mean is used to mark the baseline line on the graph via the markings attribute of the grid attribute of the plot options passed to JQuery flot's plot function. This seems to be fine as a floating-point value and being a mean seems more sensible to be a floating-point value.

985 ↗

(On Diff #220980)

color_offset is float so this is float division.

1153 ↗

(On Diff #220980)

x iterate over elements of xs which is built from the first element in each tuple composing pts, which in turns is built by calling determine_x_value. Since that method return a float, this is a float division.

1898 ↗

(On Diff #220980)

Both floor and float divisions behave the same so choosing to keep the existing code to minimize the diff.

[thopre]

lnt/external/stats/stats.py

852

As pointed out by the comment, r_den is a float (at least due to the call to math.sqrt) and thus this is a float division.

854

Denominator contains float operand and is thus float, this is a float division.

869

All variables involved are computed with float operations, this is thus a float division.

889

rs is computed with a float operand, it is therefore a float and this is a float division.

891

As mentioned in the comment, t is a float (bein computed from a sqrt call) and this is thus a float division.

922

xmean and ymean are float by virtue of being the result of a call to mean. This is a float division.

924–925

Denominator for the first line has some float operations, this is a float division and thus the line below too.

957–959

All divisions involve calls to sqrt or floating-point operands, these are float divisions.

981

r_num involves a float operand so this is a float division.

982

Float operands => float division.

984–985

float operand => first line uses a float division and thus the second line as well.

1011

Denominator is a sqrt return value, this is a float division

1042–1043

Denominator in first line is a sqrt return value, so it is a float division and therefore the second line as well.

1078–1079

denominator in first line is the return value of a sqrt call and so division on that line is a float division and thus the second line too

1172

sd is a sqrt return value, numerator has float operators => float division

1218

sqrt return value as denominator => float division

1251

se = sqrt return val => float division

1486

qap has float operands => float division

1490

em = float => float division

1493

Likewise.

1497–1499

ap, bp and app are constructed with the d values above which are float => float division

1523

coeff contains float numbers => float division

1546

Float operands in both numerator and denominator => float division

1581

Numerator and denominator formed from float division => float division

1750

As per comment, n is return value from len so this is integer division

1761

gap initialized to integer (see above) => integer division

2085

• 1st outermost if: s is float from the reduce parameter => float division
• 2nd outermost if: size is float => float division
• outermost else:
  • 1st 2nd-level if: s initialized from float and is thus going to be a float at the end
  • 2nd-level else: s initialized from float with some entries replaced by float as well and potentially reshaped => float division

All cases are float division

2125

Denom always a float or float array (outermost else) => float division.

2167

shape contains integer and result of division used for indexing => integer division

2170

Likewise, not sure why the indx assignment is not hoisted out of the if but anyway, this is a float division

2234

Obviously float division

2352

sqrt => float division

2409

denom initialized with something to the power 1.5 => float division

2428

amoment returns a float => float division

2471–2472

Both lines involve float literals => float division

2474–2475

Involve numpy's sqrt => float division

2477

alpha and delta produced by float division => float division

2500

Uses result of numpy's sqrt => float division

2501–2502

numerator has a float literal => float division

2504–2505

Float literals => float division

2507

Numerator involves float division => float division

2508

innermost division: float literal => float division
outermost division: numpy's sqrt => float division

2746

mean returns a float so m is a float and thus this is a float division

2850

sqrt denominator => float division

2862

amean returns a float => float division

2889

mns is the result of amean which is of type float => float division

2978

denominator = numpy's sqrt => float division

3057

total is float => float division

3080–3084

withindf and betwdf are float => withinms and betweenms use float division => rho, t and prob use float division

3100

covar is the result of float division => float division

3118

r_den is the result of sqrt => float division

3120–3121

Division on 1st line involves float literals => float division
Division on 2nd line involves t which involves sqrt => float division

3139

Involves float literals => float division

3141

involves t which involves sqrt => float division

3170

numerator's operands are the result of amean, thus float => float division

3172–3173

1st line: float literals => float division
2nd line: t involves sqrt, thus float => float division

3205–3207

involves sqrt or float iterals => float division

3238–3239

1st line: r_den involves sqrt => float division
2nd line: float literals => float division

3241–3243

1st line: involves float literals => float division
2nd line: involves t which involves sqrt => float division
3rd line: involves float value => float division

3279

x & y are made of floats => r_num is float => float division

3281

float literals => float division

3283–3284

1st line: float literals => float division
2nd line: t involves sqrt and thus float => float division

3288

r_num is float as per above => float division

3315–3316

1st line: involves sqrt => float division
2nd line: involves result from 1st line => float division

3355

Involves numpy's sqrt => float division

3401

step is made of floats => float division

3436

df is float => float division

3439

denom is either 1 or float values from earlier float division and when it is 1 it is substituted by 1.0 below => float division gives same result

3474

f_exp is float => float division

3542

nominator involves float literals => float division

3588

sqrt => float division

3618–3619

se is result from sqrt => float division

3756

z cast to float => float division

3875

float literals => float division

3903

qap involves float literals => qap = float value => float division

3909

involves em which is float => float division

3912

Involves em which is float => float division

3916–3918

All variables involved defined by operations involving d which is the result from ealier float division => float division

3949

coeff is full of float literals => float division

3980

Involves float literals => float division

3984

float literals => float division

4027

Involves numpy's sqrt => float division

4055

Both variables are the result of float division => float division

4104

Both variables are the result of float divisions => float division

4263

n is the result of len => integer division

4274

gap is initialized by an integer division above => integer division

[thopre]

Ping?

[thopre]

Ping?

[thopre]

Ping?

[cmatthews]

Thanks for the comments. LGTM.

[thopre]

Committed as 6d4476fb9eb0f2faa88ff3e791c9827f1b79d492

lnt/external/stats/pstat.py

"""
 
     if not isinstance(source, (list, tuple)):
         source = [source]
     for addon in args:
         if not isinstance(addon, (list, tuple)):
             addon = [addon]
         if len(addon) < len(source):                # is source list longer?
             if len(source) % len(addon) == 0:        # are they integer multiples?
-                repeats = len(source)/len(addon)    # repeat addon n times
+                repeats = len(source) // len(addon)    # repeat addon n times
                 origadd = copy.deepcopy(addon)
                 for i in range(repeats-1):
                     addon = addon + origadd
             else:
-                repeats = len(source)/len(addon)+1  # repeat addon x times,
+                repeats = len(source) // len(addon) + 1  # repeat addon x times,
                 origadd = copy.deepcopy(addon)      #    x is NOT an integer
                 for i in range(repeats-1):
                     addon = addon + origadd
                     addon = addon[0:len(source)]
         elif len(source) < len(addon):                # is addon list longer?
             if len(addon) % len(source) == 0:        # are they integer multiples?
-                repeats = len(addon)/len(source)    # repeat source n times
+                repeats = len(addon) // len(source)    # repeat source n times
                 origsour = copy.deepcopy(source)
                 for i in range(repeats-1):
                     source = source + origsour
             else:
-                repeats = len(addon)/len(source)+1  # repeat source x times,
+                repeats = len(addon) // len(source) + 1  # repeat source x times,
                 origsour = copy.deepcopy(source)    #   x is NOT an integer
                 for i in range(repeats-1):
                     source = source + origsour
                 source = source[0:len(addon)]
 
         source = simpleabut(source,addon)
     return source
 

lnt/external/stats/stats.py

 Calculates the harmonic mean of the values in the passed list.
 That is:  n / (1/x1 + 1/x2 + ... + 1/xn).  Assumes a '1D' list.
 
 Usage:   lharmonicmean(inlist)
 """
     sum = 0
     for item in inlist:
         sum = sum + 1.0/item
     return len(inlist) / sum
 
 
 def lmean (inlist):
     """
 Returns the arithematic mean of the values in the passed list.
 Assumes a '1D' list, but will function on the 1st dim of an array(!).
 
 Usage:   lmean(inlist)

 Returns the 'middle' score of the passed list.  If there is an even
 number of scores, the mean of the 2 middle scores is returned.
 
 Usage:   lmedianscore(inlist)
 """
 
     newlist = sorted(copy.deepcopy(inlist))
     if len(newlist) % 2 == 0:   # if even number of scores, average middle 2
-        index = len(newlist)/2  # integer division correct
+        index = len(newlist) // 2  # integer division correct
         median = float(newlist[index] + newlist[index-1]) /2
     else:
-        index = len(newlist)/2  # int divsion gives mid value when count from 0
+        index = len(newlist) // 2  # int divsion gives mid value when count from 0
         median = newlist[index]
     return median
 
 
 def lmode(inlist):
     """
 Returns a list of the modal (most common) score(s) in the passed
 list.  If there is more than one such score, all are returned.  The

 
 def lskew(inlist):
     """
 Returns the skewness of a distribution, as defined in Numerical
 Recipies (alternate defn in CRC Standard Probability and Statistics, p.6.)
 
 Usage:   lskew(inlist)
 """
     return moment(inlist,3)/pow(moment(inlist,2),1.5)
 
 
 def lkurtosis(inlist):
     """
 Returns the kurtosis of a distribution, as defined in Numerical
 Recipies (alternate defn in CRC Standard Probability and Statistics, p.6.)
 
 Usage:   lkurtosis(inlist)
 """
     return moment(inlist,4)/pow(moment(inlist,2),2.0)
 
 
 def ldescribe(inlist):
     """
 Returns some descriptive statistics of the passed list (assumed to be 1D).
 
 Usage:   ldescribe(inlist)
 Returns: n, mean, standard deviation, skew, kurtosis
if (defaultreallimits is not None):
             upperreallimit = 1.000001 * max(inlist)
         else: # assume both limits given
             lowerreallimit = defaultreallimits[0]
             upperreallimit = defaultreallimits[1]
         binsize = (upperreallimit-lowerreallimit)/float(numbins)
     else:     # no limits given for histogram, both must be calc'd
         estbinwidth=(max(inlist)-min(inlist))/float(numbins) +1e-6 #1=>cover all
         binsize = ((max(inlist)-min(inlist)+estbinwidth))/float(numbins)
         lowerreallimit = min(inlist) - binsize/2 #lower real limit,1st bin
     bins = [0]*(numbins)
     extrapoints = 0
     for num in inlist:
         try:
             if (num-lowerreallimit) < 0:
                 extrapoints = extrapoints + 1
             else:
                 bintoincrement = int((num-lowerreallimit)/float(binsize))
def lsem (inlist):
     """
 Returns the estimated standard error of the mean (sx-bar) of the
 values in the passed list.  sem = stdev / sqrt(n)
 
 Usage:   lsem(inlist)
 """
     sd = stdev(inlist)
     n = len(inlist)
     return sd/math.sqrt(n)
 
 
 def lz (inlist, score):
     """
 Returns the z-score for a given input score, given that score and the
 list from which that score came.  Not appropriate for population calculations.
 
 Usage:   lz(inlist, score)
 """
     z = (score-mean(inlist))/samplestdev(inlist)
     return z
 
 
 def lzs (inlist):
     """
 Returns a list of z-scores, one for each score in the passed list.
 
 Usage:   lzs(inlist)
if len(x) != len(y):
         raise ValueError('Input values not paired in pearsonr.  Aborting.')
     n = len(x)
     x = list(map(float, x))
     y = list(map(float, y))
     xmean = mean(x)
     ymean = mean(y)
     r_num = n*(summult(x,y)) - sum(x)*sum(y)
     r_den = math.sqrt((n*ss(x) - square_of_sums(x))*(n*ss(y)-square_of_sums(y)))
     r = (r_num / r_den)  # denominator already a float
     df = n-2
     t = r*math.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
     prob = betai(0.5*df,0.5,df/float(df+t*t))
     return r, prob
 
 
 def llincc(x,y):
     """
 Calculates Lin's concordance correlation coefficient.
 
 Usage:   alincc(x,y)    where x, y are equal-length arrays
 Returns: Lin's CC
 """
     covar = lcov(x,y)*(len(x)-1)/float(len(x))  # correct denom to n
     xvar = lvar(x)*(len(x)-1)/float(len(x))  # correct denom to n
     yvar = lvar(y)*(len(y)-1)/float(len(y))  # correct denom to n
     lincc = (2 * covar) / ((xvar+yvar) +((amean(x)-amean(y))**2))
     return lincc
 
 
 def lspearmanr(x,y):
     """
 Calculates a Spearman rank-order correlation coefficient.  Taken
 from Heiman's Basic Statistics for the Behav. Sci (1st), p.192.
 
 Usage:   lspearmanr(x,y)      where x and y are equal-length lists
 Returns: Spearman's r, two-tailed p-value
 """
     TINY = 1e-30
     if len(x) != len(y):
         raise ValueError('Input values not paired in spearmanr.  Aborting.')
     n = len(x)
     rankx = rankdata(x)
     ranky = rankdata(y)
     dsq = sumdiffsquared(rankx,ranky)
     rs = 1 - 6*dsq / float(n*(n**2-1))
     t = rs * math.sqrt((n-2) / ((rs+1.0)*(1.0-rs)))
     df = n-2
     probrs = betai(0.5*df,0.5,df/(df+t*t))  # t already a float
 # probability values for rs are from part 2 of the spearman function in
 # Numerical Recipies, p.510.  They are close to tables, but not exact. (?)
     return rs, probrs
 
 
 def lpointbiserialr(x,y):
     """
 Calculates a point-biserial correlation coefficient and the associated
else:   # there are 2 categories, continue
         codemap = pstat.abut(categories,range(2))
         recoded = pstat.recode(data,codemap,0)
         x = pstat.linexand(data,0,categories[0])
         y = pstat.linexand(data,0,categories[1])
         xmean = mean(pstat.colex(x,1))
         ymean = mean(pstat.colex(y,1))
         n = len(data)
         adjust = math.sqrt((len(x)/float(n))*(len(y)/float(n)))
         rpb = (ymean - xmean)/samplestdev(pstat.colex(data,1))*adjust
         df = n-2
         t = rpb*math.sqrt(df/((1.0-rpb+TINY)*(1.0+rpb+TINY)))
         prob = betai(0.5*df,0.5,df/(df+t*t))  # t already a float
         return rpb, prob
 
 
 def lkendalltau(x,y):
     """
 Calculates Kendall's tau ... correlation of ordinal data.  Adapted
 from function kendl1 in Numerical Recipies.  Needs good test-routine.@@@
 
for j in range(len(x)-1):
                     iss = iss + 1
                 else:
                     iss = iss -1
             else:
                 if (a1):
                     n1 = n1 + 1
                 else:
                     n2 = n2 + 1
     tau = iss / math.sqrt(n1*n2)
     svar = (4.0*len(x)+10.0) / (9.0*len(x)*(len(x)-1))
     z = tau / math.sqrt(svar)
     prob = erfcc(abs(z)/1.4142136)
     return tau, prob
 
 
 def llinregress(x,y):
     """
 Calculates a regression line on x,y pairs.  
 
 Usage:   llinregress(x,y)      x,y are equal-length lists of x-y coordinates
 Returns: slope, intercept, r, two-tailed prob, sterr-of-estimate
 """
     TINY = 1.0e-20
     if len(x) != len(y):
         raise ValueError('Input values not paired in linregress.  Aborting.')
     n = len(x)
     x = list(map(float, x))
     y = list(map(float, y))
     xmean = mean(x)
     ymean = mean(y)
     r_num = float(n*(summult(x,y)) - sum(x)*sum(y))
     r_den = math.sqrt((n*ss(x) - square_of_sums(x))*(n*ss(y)-square_of_sums(y)))
     r = r_num / r_den
     z = 0.5*math.log((1.0+r+TINY)/(1.0-r+TINY))
     df = n-2
     t = r*math.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
     prob = betai(0.5*df,0.5,df/(df+t*t))
     slope = r_num / float(n*ss(x) - square_of_sums(x))
     intercept = ymean - slope*xmean
     sterrest = math.sqrt(1-r*r)*samplestdev(y)
     return slope, intercept, r, prob, sterrest
 
 
 ####################################
 #####  INFERENTIAL STATISTICS  #####

 Usage:   lttest_1samp(a,popmean,Name='Sample',printit=0,writemode='a')
 Returns: t-value, two-tailed prob
 """
     x = mean(a)
     v = var(a)
     n = len(a)
     df = n-1
     svar = ((n-1)*v)/float(df)
     t = (x-popmean)/math.sqrt(svar*(1.0/n))
     prob = betai(0.5*df,0.5,float(df)/(df+t*t))
 
     if printit != 0:
         statname = 'Single-sample T-test.'
         outputpairedstats(printit,writemode,
                           'Population','--',popmean,0,0,0,
                           name,n,x,v,min(a),max(a),
                           statname,t,prob)
"""
     x1 = mean(a)
     x2 = mean(b)
     v1 = stdev(a)**2
     v2 = stdev(b)**2
     n1 = len(a)
     n2 = len(b)
     df = n1+n2-2
     svar = ((n1-1)*v1+(n2-1)*v2)/float(df)
     t = (x1-x2)/math.sqrt(svar*(1.0/n1 + 1.0/n2))
     prob = betai(0.5*df,0.5,df/(df+t*t))
 
     if printit != 0:
         statname = 'Independent samples T-test.'
         outputpairedstats(printit,writemode,
                           name1,n1,x1,v1,min(a),max(a),
                           name2,n2,x2,v2,min(b),max(b),
                           statname,t,prob)
     return t,prob
"""
     v2 = var(b)
     n = len(a)
     cov = 0
     for i in range(len(a)):
         cov = cov + (a[i]-x1) * (b[i]-x2)
     df = n-1
     cov = cov / float(df)
     sd = math.sqrt((v1+v2 - 2.0*cov)/float(n))
     t = (x1-x2)/sd
     prob = betai(0.5*df,0.5,df/(df+t*t))
 
     if printit != 0:
         statname = 'Related samples T-test.'
         outputpairedstats(printit,writemode,
                           name1,n,x1,v1,min(a),max(a),
                           name2,n,x2,v2,min(b),max(b),
                           statname,t,prob)
     return t, prob
"""
     u1 = n1*n2 + (n1*(n1+1))/2.0 - sum(rankx)  # calc U for x
     u2 = n1*n2 - u1                            # remainder is U for y
     bigu = max(u1,u2)
     smallu = min(u1,u2)
     T = math.sqrt(tiecorrect(ranked))  # correction factor for tied scores
     if T == 0:
         raise ValueError('All numbers are identical in lmannwhitneyu')
     sd = math.sqrt(T*n1*n2*(n1+n2+1)/12.0)
     z = abs((bigu-n1*n2/2.0) / sd)  # normal approximation for prob calc
     return smallu, 1.0 - zprob(z)
 
 
 def ltiecorrect(rankvals):
     """
 Corrects for ties in Mann Whitney U and Kruskal Wallis H tests.  See
 Siegel, S. (1956) Nonparametric Statistics for the Behavioral Sciences.
 New York: McGraw-Hill.  Code adapted from |Stat rankind.c code.
"""
     n1 = len(x)
     n2 = len(y)
     alldata = x+y
     ranked = rankdata(alldata)
     x = ranked[:n1]
     y = ranked[n1:]
     s = sum(x)
     expected = n1*(n1+n2+1) / 2.0
     z = (s - expected) / math.sqrt(n1*n2*(n1+n2+1)/12.0)
     prob = 2*(1.0 -zprob(abs(z)))
     return z, prob
 
 
 def lwilcoxont(x,y):
     """
 Calculates the Wilcoxon T-test for related samples and returns the
 result.  A non-parametric T-test.
"""
     for i in range(len(absd)):
         if d[i] < 0:
             r_minus = r_minus + absranked[i]
         else:
             r_plus = r_plus + absranked[i]
     wt = min(r_plus, r_minus)
     mn = count * (count+1) * 0.25
     se =  math.sqrt(count*(count+1)*(2.0*count+1.0)/24.0)
     z = math.fabs(wt-mn) / se
     prob = 2*(1.0 -zprob(abs(z)))
     return wt, prob
 
 
 def lkruskalwallish(*args):
     """
 The Kruskal-Wallis H-test is a non-parametric ANOVA for 3 or more
 groups, requiring at least 5 subjects in each group.  This function

 """
     ITMAX = 200
     EPS = 3.0e-7
 
     bm = az = am = 1.0
     qab = a+b
     qap = a+1.0
     qam = a-1.0
     bz = 1.0-qab*x/qap
     for i in range(ITMAX+1):
         em = float(i+1)
         tem = em + em
         d = em*(b-em)*x/((qam+tem)*(a+tem))
         ap = az + d*am
         bp = bz+d*bm
         d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
         app = ap+d*az
         bpp = bp+d*bz
         aold = az
         am = ap/bpp
         bm = bp/bpp
         az = app/bpp
         bz = 1.0
         if (abs(az-aold)<(EPS*abs(az))):
             return az
     print('a or b too big, or ITMAX too small in Betacf.')
 
 
 def lgammln(xx):
     """
 Returns the gamma function of xx.
     Gamma(z) = Integral(0,infinity) of t^(z-1)exp(-t) dt.
 (Adapted from: Numerical Recipies in C.)
 
 Usage:   lgammln(xx)
 """
 
     coeff = [76.18009173, -86.50532033, 24.01409822, -1.231739516,
              0.120858003e-2, -0.536382e-5]
     x = xx - 1.0
     tmp = x + 5.5
     tmp = tmp - (x+0.5)*math.log(tmp)
     ser = 1.0
     for j in range(len(coeff)):
         x = x + 1
         ser = ser + coeff[j]/x
     return -tmp + math.log(2.50662827465*ser)
 
 
 def lbetai(a,b,x):
     """
 Returns the incomplete beta function:
 
     I-sub-x(a,b) = 1/B(a,b)*(Integral(0,x) of t^(a-1)(1-t)^(b-1) dt)
 
 where a,b>0 and B(a,b) = G(a)*G(b)/(G(a+b)) where G(a) is the gamma
 function of a.  The continued fraction formulation is implemented here,
 using the betacf function.  (Adapted from: Numerical Recipies in C.)
 
 Usage:   lbetai(a,b,x)
 """
     if (x<0.0 or x>1.0):
         raise ValueError('Bad x in lbetai')
     if (x==0.0 or x==1.0):
         bt = 0.0
     else:
         bt = math.exp(gammln(a+b)-gammln(a)-gammln(b)+a*math.log(x)+b*
                       math.log(1.0-x))
     if (x<(a+1.0)/(a+b+2.0)):
         return bt*betacf(a,b,x)/float(a)
     else:
         return 1.0-bt*betacf(b,a,1.0-x)/float(b)
 
 
 ####################################
 #######  ANOVA CALCULATIONS  #######
 ####################################
"""
     for list in lists:
         ssbn = ssbn + asquare_of_sums(N.array(list))/float(len(list))
     ssbn = ssbn - (asquare_of_sums(alldata)/float(bign))
     sswn = sstot-ssbn
     dfbn = a-1
     dfwn = bign - a
     msb = ssbn/float(dfbn)
     msw = sswn/float(dfwn)
     f = msb/msw
     prob = fprob(dfbn,dfwn,f)
     return f, prob
 
 
 def lF_value (ER,EF,dfnum,dfden):
     """
 Returns an F-statistic given the following:
         ER  = error associated with the null hypothesis (the Restricted model)

 Shellsort algorithm.  Sorts a 1D-list.
 
 Usage:   lshellsort(inlist)
 Returns: sorted-inlist, sorting-index-vector (for original list)
 """
     n = len(inlist)
     svec = copy.deepcopy(inlist)
     ivec = range(n)
-    gap = n/2   # integer division needed
+    gap = n // 2   # integer division needed
     while gap >0:
         for i in range(gap,n):
             for j in range(i-gap,-1,-gap):
                 while j>=0 and svec[j]>svec[j+gap]:
                     temp        = svec[j]
                     svec[j]     = svec[j+gap]
                     svec[j+gap] = temp
                     itemp       = ivec[j]
                     ivec[j]     = ivec[j+gap]
                     ivec[j+gap] = itemp
-        gap = gap / 2  # integer division needed
+        gap = gap // 2  # integer division needed
 # svec is now sorted inlist, and ivec has the order svec[i] = vec[ivec[i]]
     return svec, ivec
 
 
 def lrankdata(inlist):
     """
 Ranks the data in inlist, dealing with ties appropritely.  Assumes
 a 1D inlist.  Adapted from Gary Perlman's |Stat ranksort.
else: # must be a SEQUENCE of dims to average over
             while incr(idx, loopcap) != -1:
                 s[idx] = asum(1.0/tinarray[idx])
             size = N.multiply.reduce(N.take(inarray.shape,dims))
             if keepdims == 1:
                 shp = list(inarray.shape)
                 for dim in dims:
                     shp[dim] = 1
                 s = N.reshape(s,shp)
     return size / s
 
 
  def amean (inarray,dimension=None,keepdims=0):
     """
 Calculates the arithmatic mean of the values in the passed array.
 That is:  1/n * (x1 + x2 + ... + xn).  Defaults to ALL values in the
 passed array.  Use dimension=None to flatten array first.  REMEMBER: if
 dimension=0, it collapses over dimension 0 ('rows' in a 2D array) only, and
else: # must be a TUPLE of dims to average over
         for dim in dims:
             sum = N.add.reduce(sum,dim)
         denom = N.array(N.multiply.reduce(N.take(inarray.shape,dims)),N.float_)
         if keepdims == 1:
             shp = list(inarray.shape)
             for dim in dims:
                 shp[dim] = 1
             sum = N.reshape(sum,shp)
     return sum/denom
 
 
  def amedian (inarray,numbins=1000):
     """
 Calculates the COMPUTED median value of an array of numbers, given the
 number of bins to use for the histogram (more bins approaches finding the
 precise median value of the array; default number of bins = 1000).  From
 G.W. Heiman's Basic Stats, or CRC Probability & Statistics.

 Usage:   amedianscore(inarray,dimension=None)
 Returns: 'middle' score of the array, or the mean of the 2 middle scores
 """
     if dimension is None:
         inarray = N.ravel(inarray)
         dimension = 0
     inarray = N.sort(inarray,dimension)
     if inarray.shape[dimension] % 2 == 0:   # if even number of elements
         indx = inarray.shape[dimension]/2   # integer division correct
         median = N.asarray(inarray[indx]+inarray[indx-1]) / 2.0
     else:
         indx = inarray.shape[dimension] / 2 # integer division correct
         median = N.take(inarray,[indx],dimension)
         if median.shape == (1,):
             median = median[0]
     return median
 
 
  def amode(a, dimension=None):
     """
"""
      elif limits[0] is None and limits[1] is not None:
          mask = upperfcn(a,limits[1])
      elif limits[0] is not None and limits[1] is None:
          mask = lowerfcn(a,limits[0])
      elif limits[0] is not None and limits[1] is not None:
          mask = lowerfcn(a,limits[0])*upperfcn(a,limits[1])
      s = float(N.add.reduce(N.ravel(a*mask)))
      n = float(N.add.reduce(N.ravel(mask)))
      return s/n
 
 
  def atvar(a,limits=None,inclusive=(1,1)):
      """
 Returns the sample variance of values in an array, (i.e., using N-1),
 ignoring values strictly outside the sequence passed to 'limits'.  
 Note: either limit in the sequence, or the value of limits itself,
 can be set to None.  The inclusive list/tuple determines whether the lower
"""
      elif limits[0] is None and limits[1] is not None:
          mask = upperfcn(a,limits[1])
      elif limits[0] is not None and limits[1] is None:
          mask = lowerfcn(a,limits[0])
      elif limits[0] is not None and limits[1] is not None:
          mask = lowerfcn(a,limits[0])*upperfcn(a,limits[1])
      term1 = N.add.reduce(N.ravel(a*a*mask))
      n = float(N.add.reduce(N.ravel(mask)))
      return sd/math.sqrt(n)
 
 
 #####################################
 ############  AMOMENTS  #############
 #####################################
 
  def amoment(a,moment=1,dimension=None):
     """

 Usage:   askew(a, dimension=None)
 Returns: skew of vals in a along dimension, returning ZERO where all vals equal
 """
     denom = N.power(amoment(a,2,dimension),1.5)
     zero = N.equal(denom,0)
     if isinstance(denom, N.ndarray) and asum(zero) != 0:
         print("Number of zeros in askew: ", asum(zero))
     denom = denom + zero  # prevent divide-by-zero
     return N.where(zero, 0, amoment(a,3,dimension)/denom)
 
 
  def akurtosis(a,dimension=None):
     """
 Returns the kurtosis of a distribution (normal ==> 3.0; >3 means
 heavier in the tails, and usually more peaked).  Use akurtosistest()
 to see if it's close enough.  Dimension can equal None (ravel array
 first), an integer (the dimension over which to operate), or a
 sequence (operate over multiple dimensions).
 
 Usage:   akurtosis(a,dimension=None)
 Returns: kurtosis of values in a along dimension, and ZERO where all vals equal
 """
     denom = N.power(amoment(a,2,dimension),2)
     zero = N.equal(denom,0)
     if isinstance(denom, N.ndarray) and asum(zero) != 0:
         print("Number of zeros in akurtosis: ", asum(zero))
     denom = denom + zero  # prevent divide-by-zero
     return N.where(zero,0,amoment(a,4,dimension)/denom)
 
 
  def adescribe(inarray,dimension=None):
      """
 Returns several descriptive statistics of the passed array.  Dimension
 can equal None (ravel array first), an integer (the dimension over
 which to operate), or a sequence (operate over multiple dimensions).
 

 Usage:   askewtest(a,dimension=None)
 Returns: z-score and 2-tail z-probability
 """
     if dimension is None:
         a = N.ravel(a)
         dimension = 0
     b2 = askew(a,dimension)
     n = float(a.shape[dimension])
     y = b2 * N.sqrt(((n+1)*(n+3)) / (6.0*(n-2)) )
     beta2 = ( 3.0*(n*n+27*n-70)*(n+1)*(n+3) ) / ( (n-2.0)*(n+5)*(n+7)*(n+9) )
     W2 = -1 + N.sqrt(2*(beta2-1))
     delta = 1/N.sqrt(N.log(N.sqrt(W2)))
     alpha = N.sqrt(2/(W2-1))
     y = N.where(y==0,1,y)
     Z = delta*N.log(y/alpha + N.sqrt((y/alpha)**2+1))
     return Z, (1.0-zprob(Z))*2
 
 
  def akurtosistest(a,dimension=None):
     """
 Tests whether a dataset has normal kurtosis (i.e.,
 kurtosis=3(n-1)/(n+1)) Valid only for n>20.  Dimension can equal None
 (ravel array first), an integer (the dimension over which to operate),
 or a sequence (operate over multiple dimensions).
 
 Usage:   akurtosistest(a,dimension=None)
 Returns: z-score and 2-tail z-probability, returns 0 for bad pixels
 """
     if dimension is None:
         a = N.ravel(a)
         dimension = 0
     n = float(a.shape[dimension])
     if n<20:
         print("akurtosistest only valid for n>=20 ... continuing anyway, n=", n)
     b2 = akurtosis(a,dimension)
     E = 3.0*(n-1) /(n+1)
     varb2 = 24.0*n*(n-2)*(n-3) / ((n+1)*(n+1)*(n+3)*(n+5))
     x = (b2-E)/N.sqrt(varb2)
     sqrtbeta1 = 6.0*(n*n-5*n+2)/((n+7)*(n+9)) * N.sqrt((6.0*(n+3)*(n+5))/
                                                        (n*(n-2)*(n-3)))
     A = 6.0 + 8.0/sqrtbeta1 *(2.0/sqrtbeta1 + N.sqrt(1+4.0/(sqrtbeta1**2)))
     term1 = 1 -2/(9.0*A)
     denom = 1 +x*N.sqrt(2/(A-4.0))
     denom = N.where(N.less(denom,0), 99, denom)
     term2 = N.where(N.equal(denom,0), term1, N.power((1-2.0/A)/denom,1/3.0))
     Z = ( term1 - term2 ) / N.sqrt(2/(9.0*A))
     Z = N.where(N.equal(denom,99), 0, Z)
     return Z, (1.0-zprob(Z))*2
 
 
  def anormaltest(a,dimension=None):
     """
 Tests whether skew and/OR kurtosis of dataset differs from normal
 curve.  Can operate over multiple dimensions.  Dimension can equal

 sequence (operate over multiple dimensions).
 
 Usage:   asignaltonoise(instack,dimension=0):
 Returns: array containing the value of (mean/stdev) along dimension,
          or 0 when stdev=0
 """
     m = mean(instack,dimension)
     sd = stdev(instack,dimension)
     return N.where(sd==0,0,m/sd)
 
 
  def acov (x,y, dimension=None,keepdims=0):
     """
 Returns the estimated covariance of the values in the passed
 array (i.e., N-1).  Dimension can equal None (ravel array first), an
 integer (the dimension over which to operate), or a sequence (operate
 over multiple dimensions).  Set keepdims=1 to return an array with the
if dimension is None:
         inarray = N.ravel(inarray)
         dimension = 0
     if isinstance(dimension, list):
         n = 1
         for d in dimension:
             n = n*inarray.shape[d]
     else:
         n = inarray.shape[dimension]
     s = asamplestdev(inarray,dimension,keepdims) / N.sqrt(n-1)
     return s
 
 
  def az (a, score):
     """
 Returns the z-score of a given input score, given thearray from which
 that score came.  Not appropriate for population calculations, nor for
 arrays > 1D.
 
 Usage:   az(a, score)
 """
     z = (score-amean(a)) / asamplestdev(a)
     return z
 
 
  def azs (a):
     """
 Returns a 1D array of z-scores, one for each score in the passed array,
 computed relative to the passed array.
 

 Returns an array of z-scores the shape of scores (e.g., [x,y]), compared to
 array passed to compare (e.g., [time,x,y]).  Assumes collapsing over dim 0
 of the compare array.
 
 Usage:   azs(scores, compare, dimension=0)
 """
     mns = amean(compare,dimension)
     sstd = asamplestdev(compare,0)
     return (scores - mns) / sstd
 
 
 #####################################
 #######  ATRIMMING FUNCTIONS  #######
 #####################################
 
 ## deleted around() as it's in numpy now
 
def acorrelation(X):
     """
 Computes the correlation matrix of a matrix X.  Requires a 2D matrix input.
 
 Usage:   acorrelation(X)
 Returns: correlation matrix of X
 """
     C = acovariance(X)
     V = N.diagonal(C)
     return C / N.sqrt(N.multiply.outer(V,V))
 
 
  def apaired(x,y):
     """
 Interactively determines the type of data in x and y, and then runs the
 appropriated statistic for paired group data.
 
 Usage:   apaired(x,y)     x,y = the two arrays of values to be compared

 
 Usage:  dices(x,y)
 """
     import sets
     x = sets.Set(x)
     y = sets.Set(y)
     common = len(x.intersection(y))
     total = float(len(x) + len(y))
     return 2*common/total
 
 
  def icc(x,y=None,verbose=0):
     """
 Calculates intraclass correlation coefficients using simple, Type I sums of squares.
 If only one variable is passed, assumed it's an Nx2 matrix
 
 Usage:   icc(x,y=None,verbose=0)
 Returns: icc rho, prob ####PROB IS A GUESS BASED ON PEARSON
 """
     TINY = 1.0e-20
     if y:
         all = N.concatenate([x,y],0)
     else:
         all = x+0
         x = all[:,0]
         y = all[:,1]
     totalss = ass(all-mean(all))
     pairmeans = (x+y)/2.
     withinss = ass(x-pairmeans) + ass(y-pairmeans)
     withindf = float(len(x))
     betwdf = float(len(x)-1)
     withinms = withinss / withindf
     betweenms = (totalss-withinss) / betwdf
     rho = (betweenms-withinms)/(withinms+betweenms)
     t = rho*math.sqrt(betwdf/((1.0-rho+TINY)*(1.0+rho+TINY)))
     prob = abetai(0.5*betwdf,0.5,betwdf/(betwdf+t*t),verbose)
     return rho, prob
 
 
  def alincc(x,y):
     """
 Calculates Lin's concordance correlation coefficient.
 
 Usage:   alincc(x,y)    where x, y are equal-length arrays
 Returns: Lin's CC
 """
     x = N.ravel(x)
     y = N.ravel(y)
     covar = acov(x,y)*(len(x)-1)/float(len(x))  # correct denom to n
     xvar = avar(x)*(len(x)-1)/float(len(x))  # correct denom to n
     yvar = avar(y)*(len(y)-1)/float(len(y))  # correct denom to n
     lincc = (2 * covar) / ((xvar+yvar) +((amean(x)-amean(y))**2))
     return lincc
 
 
  def apearsonr(x,y,verbose=1):
     """
 Calculates a Pearson correlation coefficient and returns p.  Taken
 from Heiman's Basic Statistics for the Behav. Sci (2nd), p.195.
 
 Usage:   apearsonr(x,y,verbose=1)      where x,y are equal length arrays
 Returns: Pearson's r, two-tailed p-value
 """
     TINY = 1.0e-20
     n = len(x)
     xmean = amean(x)
     ymean = amean(y)
     r_num = n*(N.add.reduce(x*y)) - N.add.reduce(x)*N.add.reduce(y)
     r_den = math.sqrt((n*ass(x) - asquare_of_sums(x))*(n*ass(y)-asquare_of_sums(y)))
     r = (r_num / r_den)
     df = n-2
     t = r*math.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
     prob = abetai(0.5*df,0.5,df/(df+t*t),verbose)
     return r,prob
 
 
  def aspearmanr(x,y):
     """
 Calculates a Spearman rank-order correlation coefficient.  Taken
 from Heiman's Basic Statistics for the Behav. Sci (1st), p.192.
 
 Usage:   aspearmanr(x,y)      where x,y are equal-length arrays
 Returns: Spearman's r, two-tailed p-value
 """
     TINY = 1e-30
     n = len(x)
     rankx = rankdata(x)
     ranky = rankdata(y)
     dsq = N.add.reduce((rankx-ranky)**2)
     rs = 1 - 6*dsq / float(n*(n**2-1))
     t = rs * math.sqrt((n-2) / ((rs+1.0)*(1.0-rs)))
     df = n-2
     probrs = abetai(0.5*df,0.5,df/(df+t*t))
 # probability values for rs are from part 2 of the spearman function in
 # Numerical Recipies, p.510.  They close to tables, but not exact.(?)
     return rs, probrs
 
 
  def apointbiserialr(x,y):
     """
 Calculates a point-biserial correlation coefficient and the associated
else:   # there are 2 categories, continue
         codemap = pstat.aabut(categories,N.arange(2))
         recoded = pstat.arecode(data,codemap,0)
         x = pstat.alinexand(data,0,categories[0])
         y = pstat.alinexand(data,0,categories[1])
         xmean = amean(pstat.acolex(x,1))
         ymean = amean(pstat.acolex(y,1))
         n = len(data)
         adjust = math.sqrt((len(x)/float(n))*(len(y)/float(n)))
         rpb = (ymean - xmean)/asamplestdev(pstat.acolex(data,1))*adjust
         df = n-2
         t = rpb*math.sqrt(df/((1.0-rpb+TINY)*(1.0+rpb+TINY)))
         prob = abetai(0.5*df,0.5,df/(df+t*t))
         return rpb, prob
 
 
  def akendalltau(x,y):
     """
 Calculates Kendall's tau ... correlation of ordinal data.  Adapted
 from function kendl1 in Numerical Recipies.  Needs good test-cases.@@@
 
for j in range(len(x)-1):
                     iss = iss + 1
                 else:
                     iss = iss -1
             else:
                 if (a1):
                     n1 = n1 + 1
                 else:
                     n2 = n2 + 1
     tau = iss / math.sqrt(n1*n2)
     svar = (4.0*len(x)+10.0) / (9.0*len(x)*(len(x)-1))
     z = tau / math.sqrt(svar)
     prob = erfcc(abs(z)/1.4142136)
     return tau, prob
 
 
  def alinregress(*args):
     """
 Calculates a regression line on two arrays, x and y, corresponding to x,y
 pairs.  If a single 2D array is passed, alinregress finds dim with 2 levels
"""
     else:
         x = args[0]
         y = args[1]
     n = len(x)
     xmean = amean(x)
     ymean = amean(y)
     r_num = n*(N.add.reduce(x*y)) - N.add.reduce(x)*N.add.reduce(y)
     r_den = math.sqrt((n*ass(x) - asquare_of_sums(x))*(n*ass(y)-asquare_of_sums(y)))
     r = r_num / r_den
     z = 0.5*math.log((1.0+r+TINY)/(1.0-r+TINY))
     df = n-2
     t = r*math.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
     prob = abetai(0.5*df,0.5,df/(df+t*t))
     slope = r_num / (float(n)*ass(x) - asquare_of_sums(x))
     intercept = ymean - slope*xmean
     sterrest = math.sqrt(1-r*r)*asamplestdev(y)
     return slope, intercept, r, prob, sterrest, n
 
  def amasslinregress(*args):
     """
 Calculates a regression line on one 1D array (x) and one N-D array (y).
 
"""
     shp = N.ones(len(y.shape))
     shp[0] = len(x)
     x.shape = shp
     print(x.shape, y.shape)
     r_num = n*(N.add.reduce(x*y,0)) - N.add.reduce(x)*N.add.reduce(y,0)
     r_den = N.sqrt((n*ass(x) - asquare_of_sums(x))*(n*ass(y,0)-asquare_of_sums(y,0)))
     zerodivproblem = N.equal(r_den,0)
     r_den = N.where(zerodivproblem,1,r_den)  # avoid zero-division in 1st place
     r = r_num / r_den  # need to do this nicely for matrix division
     r = N.where(zerodivproblem,0.0,r)
     z = 0.5*N.log((1.0+r+TINY)/(1.0-r+TINY))
     df = n-2
     t = r*N.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
     prob = abetai(0.5*df,0.5,df/(df+t*t))
 
     ss = float(n)*ass(x)-asquare_of_sums(x)
     s_den = N.where(ss==0,1,ss)  # avoid zero-division in 1st place
     slope = r_num / s_den
     intercept = ymean - slope*xmean
     sterrest = N.sqrt(1-r*r)*asamplestdev(y,0)
     return slope, intercept, r, prob, sterrest, n
 
 
 #####################################
 #####  AINFERENTIAL STATISTICS  #####
 #####################################

 """
     if not isinstance(a, N.ndarray):
         a = N.array(a)
     x = amean(a)
     v = avar(a)
     n = len(a)
     df = n-1
     svar = ((n-1)*v) / float(df)
     t = (x-popmean)/math.sqrt(svar*(1.0/n))
     prob = abetai(0.5*df,0.5,df/(df+t*t))
 
     if printit != 0:
         statname = 'Single-sample T-test.'
         outputpairedstats(printit,writemode,
                           'Population','--',popmean,0,0,0,
                           name,n,x,v,N.minimum.reduce(N.ravel(a)),
                           N.maximum.reduce(N.ravel(a)),
                           statname,t,prob)
"""
     v1 = avar(a,dimension)
     v2 = avar(b,dimension)
     n1 = a.shape[dimension]
     n2 = b.shape[dimension]
     df = n1+n2-2
     svar = ((n1-1)*v1+(n2-1)*v2) / float(df)
     zerodivproblem = N.equal(svar,0)
     svar = N.where(zerodivproblem,1,svar)  # avoid zero-division in 1st place
     t = (x1-x2)/N.sqrt(svar*(1.0/n1 + 1.0/n2))  # N-D COMPUTATION HERE!!!!!!
     t = N.where(zerodivproblem,1.0,t)     # replace NaN/wrong t-values with 1.0
     probs = abetai(0.5*df,0.5,float(df)/(df+t*t))
 
     if isinstance(t, N.ndarray):
         probs = N.reshape(probs,t.shape)
     if probs.shape == (1,):
         probs = probs[0]
         
Returns: an array of t-values with the shape of pval
     step = N.ones(pval.shape,N.float_)*25
     print("Initial ap2t() prob calc")
     prob = abetai(0.5*df,0.5,float(df)/(df+t*t))
     print('ap2t() iter: ', end=' ')
     for i in range(10):
         print(i, ' ', end=' ')
         t = N.where(pval<prob,t+step,t-step)
         prob = abetai(0.5*df,0.5,float(df)/(df+t*t))
         step = step/2
     print()
     # since this is an ugly hack, we get ugly boundaries
     t = N.where(t>99.9,1000,t)      # hit upper-boundary
     t = t+signs
     return t #, prob, pval
 
 
  def attest_rel (a,b,dimension=None,printit=0,name1='Samp1',name2='Samp2',writemode='a'):
"""
     x1 = amean(a,dimension)
     x2 = amean(b,dimension)
     v1 = avar(a,dimension)
     v2 = avar(b,dimension)
     n = a.shape[dimension]
     df = float(n-1)
     d = (a-b).astype('d')
 
     denom = N.sqrt((n*N.add.reduce(d*d,dimension) - N.add.reduce(d,dimension)**2) /df)
     zerodivproblem = N.equal(denom,0)
     denom = N.where(zerodivproblem,1,denom)  # avoid zero-division in 1st place
     t = N.add.reduce(d,dimension) / denom      # N-D COMPUTATION HERE!!!!!!
     t = N.where(zerodivproblem,1.0,t)     # replace NaN/wrong t-values with 1.0
     probs = abetai(0.5*df,0.5,float(df)/(df+t*t))
     if isinstance(t, N.ndarray):
         probs = N.reshape(probs,t.shape)
     if probs.shape == (1,):
         probs = probs[0]
 
     if printit != 0:

 Usage:   achisquare(f_obs, f_exp=None)   f_obs = array of observed cell freq.
 Returns: chisquare-statistic, associated p-value
 """
 
     k = len(f_obs)
     if f_exp is None:
         f_exp = N.array([sum(f_obs)/float(k)] * len(f_obs),N.float_)
     f_exp = f_exp.astype(N.float_)
     chisq = N.add.reduce((f_obs-f_exp)**2 / f_exp)
     return chisq, achisqprob(chisq, k-1)
 
 
  def aks_2samp (data1,data2):
     """
 Computes the Kolmogorov-Smirnof statistic on 2 samples.  Modified from
 Numerical Recipies in C, page 493.  Returns KS D-value, prob.  Not ufunc-
 like.
"""
     u1 = n1*n2 + (n1*(n1+1))/2.0 - sum(rankx)  # calc U for x
     u2 = n1*n2 - u1                            # remainder is U for y
     bigu = max(u1,u2)
     smallu = min(u1,u2)
     T = math.sqrt(tiecorrect(ranked))  # correction factor for tied scores
     if T == 0:
         raise ValueError('All numbers are identical in amannwhitneyu')
     sd = math.sqrt(T*n1*n2*(n1+n2+1)/12.0)
     z = abs((bigu-n1*n2/2.0) / sd)  # normal approximation for prob calc
     return smallu, 1.0 - azprob(z)
 
 
  def atiecorrect(rankvals):
     """
 Tie-corrector for ties in Mann Whitney U and Kruskal Wallis H tests.
 See Siegel, S. (1956) Nonparametric Statistics for the Behavioral
 Sciences.  New York: McGraw-Hill.  Code adapted from |Stat rankind.c
"""
     n1 = len(x)
     n2 = len(y)
     alldata = N.concatenate((x,y))
     ranked = arankdata(alldata)
     x = ranked[:n1]
     y = ranked[n1:]
     s = sum(x)
     expected = n1*(n1+n2+1) / 2.0
     z = (s - expected) / math.sqrt(n1*n2*(n1+n2+1)/12.0)
     prob = 2*(1.0 - azprob(abs(z)))
     return z, prob
 
 
  def awilcoxont(x,y):
     """
 Calculates the Wilcoxon T-test for related samples and returns the
 result.  A non-parametric T-test.
"""
     for i in range(len(absd)):
         if d[i] < 0:
             r_minus = r_minus + absranked[i]
         else:
             r_plus = r_plus + absranked[i]
     wt = min(r_plus, r_minus)
     mn = count * (count+1) * 0.25
     se =  math.sqrt(count*(count+1)*(2.0*count+1.0)/24.0)
     z = math.fabs(wt-mn) / se
     z = math.fabs(wt-mn) / se
     prob = 2*(1.0 -zprob(abs(z)))
     return wt, prob
 
 
  def akruskalwallish(*args):
     """
 The Kruskal-Wallis H-test is a non-parametric ANOVA for 3 or more
 groups, requiring at least 5 subjects in each group.  This function
#            print z, e, s
             e = N.ones(probs.shape,N.float_)
         else:
             z = 0.5 *N.ones(probs.shape,N.float_)
             e = 1.0 / N.sqrt(N.pi) / N.sqrt(a) * N.ones(probs.shape,N.float_)
         c = 0.0
         mask = N.zeros(probs.shape)
         a_notbig_frozen = -1 *N.ones(probs.shape,N.float_)
         while asum(mask) != totalelements:
             e = e * (a/z.astype(N.float_))
             c = c + e
             z = z + 1.0
 #            print '#2', z, e, c, s, c*y+s2
             newmask = N.greater(z,chisq)
             a_notbig_frozen = N.where(newmask*N.equal(mask,0)*(1-a_big),
                                       c*y+s2, a_notbig_frozen)
             mask = N.clip(newmask+mask,0,1)
         probs = N.where(N.equal(probs,1),1,
def afprob (dfnum, dfden, F):
     """
 Returns the 1-tailed significance level (p-value) of an F statistic
 given the degrees of freedom for the numerator (dfR-dfF) and the degrees
 of freedom for the denominator (dfF).  Can handle multiple dims for F.
 
 Usage:   afprob(dfnum, dfden, F)   where usually dfnum=dfbn, dfden=dfwn
 """
     if isinstance(F, N.ndarray):
         return abetai(0.5*dfden, 0.5*dfnum, dfden/(1.0*dfden+dfnum*F))
     else:
         return abetai(0.5*dfden, 0.5*dfnum, dfden/float(dfden+dfnum*F))
 
 
  def abetacf(a,b,x,verbose=1):
     """
 Evaluates the continued fraction form of the incomplete Beta function,
 betai.  (Adapted from: Numerical Recipies in C.)  Can handle multiple
else:
         arrayflag = 0
         frozen = N.array([-1])
         x = N.array([x])
     mask = N.zeros(x.shape)
     bm = az = am = 1.0
     qab = a+b
     qap = a+1.0
     qam = a-1.0
     bz = 1.0-qab*x/qap
     for i in range(ITMAX+1):
         if N.sum(N.ravel(N.equal(frozen,-1)))==0:
             break
         em = float(i+1)
         tem = em + em
         d = em*(b-em)*x/((qam+tem)*(a+tem))
         ap = az + d*am
         bp = bz+d*bm
         d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
         app = ap+d*az
         bpp = bp+d*bz
         aold = az*1
         am = ap/bpp
         bm = bp/bpp
         az = app/bpp
         bz = 1.0
         newmask = N.less(abs(az-aold),EPS*abs(az))
         frozen = N.where(newmask*N.equal(mask,0), az, frozen)
         mask = N.clip(mask+newmask,0,1)
     noconverge = asum(N.equal(frozen,-1))
     if noconverge != 0 and verbose:
         print('a or b too big, or ITMAX too small in Betacf for ', noconverge, ' elements')
     if arrayflag:
"""
     coeff = [76.18009173, -86.50532033, 24.01409822, -1.231739516,
              0.120858003e-2, -0.536382e-5]
     x = xx - 1.0
     tmp = x + 5.5
     tmp = tmp - (x+0.5)*N.log(tmp)
     ser = 1.0
     for j in range(len(coeff)):
         x = x + 1
         ser = ser + coeff[j]/x
     return -tmp + N.log(2.50662827465*ser)
 
 
  def abetai(a,b,x,verbose=1):
     """
 Returns the incomplete beta function:
 
     I-sub-x(a,b) = 1/B(a,b)*(Integral(0,x) of t^(a-1)(1-t)^(b-1) dt)
"""
 
     bt = N.where(N.equal(x,0)+N.equal(x,1), 0, -1)
     exponents = ( gammln(a+b)-gammln(a)-gammln(b)+a*N.log(x)+b*
                   N.log(1.0-x) )
     # 746 (below) is the MAX POSSIBLE BEFORE OVERFLOW
     exponents = N.where(N.less(exponents,-740),-740,exponents)
     bt = N.exp(exponents)
     if isinstance(x, N.ndarray):
         ans = N.where(N.less(x,(a+1)/(a+b+2.0)),
                       bt*abetacf(a,b,x,verbose)/float(a),
                       1.0-bt*abetacf(b,a,1.0-x,verbose)/float(b))
     else:
         if x<(a+1)/(a+b+2.0):
             ans = bt*abetacf(a,b,x,verbose)/float(a)
         else:
             ans = 1.0-bt*abetacf(b,a,1.0-x,verbose)/float(b)
     return ans
 
 
 #####################################
 #######  AANOVA CALCULATIONS  #######
b = N.dot(N.dot(LA.inv(N.dot(N.transpose(x),x)),  # i.e., b=inv(X'X)X'Y
               data)
     diffs = (data - N.dot(x,b))
     s_sq = 1./(n-len(p)) * N.dot(N.transpose(diffs), diffs)
 
     if len(p) == 2:  # ttest_ind
         c = N.array([1,-1])
         df = n-2
         fact = asum(1.0/asum(x,0))  # i.e., 1/n1 + 1/n2 + 1/n3 ...
         t = N.dot(c,b) / N.sqrt(s_sq*fact)
         probs = abetai(0.5*df,0.5,float(df)/(df+t*t))
         return t, probs
 
 
  def aF_oneway(*args):
     """
 Performs a 1-way ANOVA, returning an F-value and probability given
 any number of groups.  From Heiman, pp.394-7.
"""
     for a in args:
         ssbn = ssbn + asquare_of_sums(N.array(a))/float(len(a))
     ssbn = ssbn - (asquare_of_sums(alldata)/float(bign))
     sswn = sstot-ssbn
     dfbn = na-1
     dfwn = bign - na
     msb = ssbn/float(dfbn)
     msw = sswn/float(dfwn)
     f = msb/msw
     prob = fprob(dfbn,dfwn,f)
     return f, prob
 
 
  def aF_value (ER,EF,dfR,dfF):
     """
 Returns an F-statistic given the following:
         ER  = error associated with the null hypothesis (the Restricted model)

 where ER and EF are matrices from a multivariate F calculation.
 """
      if isinstance(ER, (int, float)):
          ER = N.array([[ER]])
      if isinstance(EF, (int, float)):
          EF = N.array([[EF]])
      n_um = (LA.det(ER) - LA.det(EF)) / float(dfnum)
      d_en = LA.det(EF) / float(dfden)
      return n_um / d_en
 
 
 #####################################
 #######  ASUPPORT FUNCTIONS  ########
 #####################################
 
  def asign(a):
     """

 Shellsort algorithm.  Sorts a 1D-array.
 
 Usage:   ashellsort(inarray)
 Returns: sorted-inarray, sorting-index-vector (for original array)
 """
     n = len(inarray)
     svec = inarray *1.0
     ivec = range(n)
-    gap = n/2   # integer division needed
+    gap = n // 2   # integer division needed
     while gap >0:
         for i in range(gap,n):
             for j in range(i-gap,-1,-gap):
                 while j>=0 and svec[j]>svec[j+gap]:
                     temp        = svec[j]
                     svec[j]     = svec[j+gap]
                     svec[j+gap] = temp
                     itemp       = ivec[j]
                     ivec[j]     = ivec[j+gap]
                     ivec[j+gap] = itemp
-        gap = gap / 2  # integer division needed
+        gap = gap // 2  # integer division needed
 #    svec is now sorted input vector, ivec has the order svec[i] = vec[ivec[i]]
     return svec, ivec
 
 
  def arankdata(inarray):
     """
 Ranks the data in inarray, dealing with ties appropritely.  Assumes
 a 1D inarray.  Adapted from Gary Perlman's |Stat ranksort.

lnt/server/db/rules/rule_update_profile_stats.py

def update_profile_stats(session, ts, run_id):
     blocks = subprocess.check_output("du -s -k %s" % profile_path,
                                      shell=True).split('\t')[0]
     kb = float(blocks)  # 1024 byte blocks.
 
     history.append((dt, kb))
 
     for f in glob.glob('%s/*.lntprof' % profile_path):
         mtime = os.stat(f).st_mtime
-        sz = os.stat(f).st_size / 1000
+        sz = os.stat(f).st_size // 1000
         age.append([mtime, sz])
 
     open(history_path, 'w').write(json.dumps(history))
     open(age_path, 'w').write(json.dumps(age))
 
 
 post_submission_hook = update_profile_stats

lnt/server/reporting/dailyreport.py

def __init__(self, ts, year, month, day, num_prior_days_to_include=3,
         self.result_table = None
         self.nr_tests_table = None
 
     def get_query_parameters_string(self):
         query_params = [
             "{}={}".format(
                 urllib.quote_plus(query_param), urllib.quote_plus(str(value)))
             for query_param, value in (
-                ("day_start", self.day_start_offset.seconds / 3600),
+                ("day_start", self.day_start_offset.seconds // 3600),
                 ("num_days", self.num_prior_days_to_include),
                 ("filter-machine-regex", self.filter_machine_regex_str),)
             if value is not None]
         # Use & instead of plain & to make sure legal XML is
         # produced, so that the unit tests can parse the output using
         # python's built-in XML parser.
         return "&".join(query_params)