Cook and Felderman analysis to perform data reduction on heated surfaces

This is a "pain", but I wonder if I can get any real speedup with the follow Python script - it is taking well over 10 minutes to calculate the way it is written now. (And that's only when feeding around 1-2 promille of the data).

I have an with shape (roughly): [260, 300, 300]

Now I am doing a cougar=feldermann analysis on this:

import numpy as np
import math
#COOK AND FELDERMAN
def CandF(dT, time, rho, c, k):
sqrt = np.sqrt
sumdata=np.zeros(len(time))
for n in range(len(time)):
for data in range(1,n+1):
sumdata[n] += dT[data] / (sqrt(time[n]-time[data])+sqrt(time[n]-time[data-1]))

qs = (2*sqrt(rho*c*k))/(sqrt(math.pi)) * sumdata
return qs

def data_reduction_constant_time(dT, timestep, rho, cp, k):
#print (len(dT), timestep)
time = np.linspace(0, len(dT)*timestep, int(len(dT)))
return data_reduction(dT, time, rho, cp, k )

def data_reduction(dT, time, rho, cp, k):
return CandF(dT, time, rho, cp, k)

class data:
def __init__(self):

self.integration = 400
self.cols = 320
self.rows = 260
self.depth = 300

self.ml_delta_temp = 3 * np.random.random_sample((self.cols, self.rows, self.depth))

def calcQML(self):
print("--- Calculating q ---")
print (self.cols, self.rows, self.depth)
v = np.zeros([self.cols,self.rows,self.depth])
T = 590
P = 28 * 100000
R = 287
c = 1005
k = 1.4
q = 20
rho = P/(R*T)
print(v)

print("-------------")
for x in range(self.cols ):
print ("calculating column: ", x)
for y in range(self.rows ):
v_tmp = data_reduction_constant_time(self.ml_delta_temp[x,y,:],
self.integration,
rho, c, k)
v[x,y,:] = v_tmp
print("--- q done ---")
print(v)

if (__name__ == "__main__"):
dat = data()
dat.calcQML()


self.ml_delta_temp is a (260, 300, 300) shaped numpy array. cols, and rows are the size of the array, and integration is the step size (how much time is between each 3rd axis slice).

Now the problem is: on my computer it takes almost a minute to calculate a single row.

How can I speed it up, or as the same question: in what part is this function losing the most time?

I can show the mathematical formula for the Cook and Felderman equation:

$$qCF = \frac{2\sqrt{\rho ck}}{\sqrt{\pi}}\sum_{i=1}^n \frac{T_i-T_{i-1}}{\sqrt{t_n-t_i}+\sqrt{t_n-t_{i+1}}}$$

Berry, S.A. and Horvath, T.J. Discrete Roughness Transition for Hypersonic Flight Vehicles. Journal of Spacecraft and Rockets, 2:216{ 227, March-April 2008. DOI 10.2514/1.30970.

Now this is to get the total heat at time "t_n", I am however interested in the heat flux profile, how it changes over time. So I need to calculate this formula for each n. Also I do know not T, but only delta-T (but as one can see that is only beneficial).

• Something seems wrong with CandF(): sumdata[n] is invariant in the most inner loop, so it is overwritten (n-1) times - ? Could you please describe why you calculate (n-1) values and only keep 1? – user1016274 Mar 21 '15 at 23:28
• OK, got it: should be sumdata[n] += dT[data] / (np.sqrt(.... To save some time you should declare local functions to spare module lookups: sqrt = np.sqrt and then use sqrt() within CandF(). – user1016274 Mar 21 '15 at 23:35
• @user1016274 what do you mean only keep 1? When you multiply a scalar with an array you do generate an element wise multiplication. – paul23 Mar 21 '15 at 23:35
• @user1016274 and here is the cook and felderman formula. – paul23 Mar 21 '15 at 23:49
• This is difficult to answer without even having a runnable program to profile. Can you post the full code and a description of the inputs you're using (even if it's just "each value in delta_temp is between some lower and upper bound, so random data can be generated as input). – Yuushi Mar 22 '15 at 4:16

I suggest the following replacement for function CandF():

def CandF3(dT, time, fac):
sqrt = np.sqrt
sumdata = np.zeros(len(time))
for n in range(len(time)):
tn = time[n]
# precompute sqrt(tn-time[0])
last_sq = sqrt(tn - time[0])
sum_n = 0
for j in range(1, n + 1):
sq_j = sqrt(tn - time[j])
sum_n += dT[j] / (sq_j + last_sq)
last_sq = sq_j
sumdata[n] = sum_n
return sumdata * fac


The most rewarding idea here was that from the 2 square roots per innermost loop one is the second is computed repeatedly and can be saved for the next loop (here: last_sq). Halving the number of np.sqrt() calls nearly halves the runtime for large numbers of time steps (len(time)).
More optimizations are the use of a local variable for a loop-invariant (tn), use of a scalar when indexed variable is not indexed incrementally (sum_n), and precomputing the final constant sqrt(rho*c*k*2.0) in the caller. This seems natural as the contributing constants are all defined there.

Generally, you want to try and avoid looping over things as much as possible when working with data like this. The more you can break things into operations on entire arrays/matrices, the faster it will be, generally thanks to SIMD.

Clearly, most of the time is being spent in the CandF function, so this is where time would be best spent optimizing.

Using array and matrix primitives instead of looping inside this function gives a roughly ~45x speedup on my machine:

def CandF(dT, time, rho, c, k):
sqrt = np.sqrt
sumdata=np.zeros(len(time))
for n in xrange(len(time)):
z = np.repeat(time[n], n)
sumdata[n] = np.sum(dT[1:n+1] / (sqrt(z - time[1:n+1]) + sqrt(z - time[0:n])))

qs = (2*sqrt(rho*c*k))/(sqrt(math.pi)) * sumdata
return qs


Potentially you could squeeze out some more performance by recognizing that the calculation of np.sqrt(z - time[1:n+1]) calculates almost everything you need for sqrt(z - time[0:n])) (only difference is the first value, and the rest of the values need to be shifted by one), but hopefully the above gets you started at least.