8
\$\begingroup\$

This code (~30lines) detects independence between variables in a dataset, by using a bootstrap of a statistic. This is fully-broadcasted, but it's still long.

Reading it, or making it run, do you see some possibilities for more optimisation?

Right here is a MWE with a sample of data with 100 points. Feel free to comment

import numpy as np

########### Data and prerequisites
# Sample of data with strong dependence structure :
data = np.array([[5.61293390e-01, 9.97100450e-01, 4.23530180e-01, 6.08808896e-01],
       [1.22563280e-01, 1.72015130e-01, 8.71145720e-01, 5.40745844e-01],
       [8.51194500e-02, 1.18289130e-01, 8.90346540e-01, 7.22859351e-01],
       [9.83241310e-01, 9.57282690e-01, 7.22347100e-02, 5.43527399e-02],
       [5.49211550e-01, 3.97858250e-01, 6.86380990e-01, 7.91494336e-01],
       [9.94878920e-01, 6.39160920e-01, 2.01045170e-01, 9.86840712e-01],
       [5.04337540e-01, 5.69995040e-01, 3.99087430e-01, 4.32140476e-01],
       [9.28230540e-01, 9.32143440e-01, 1.02748280e-01, 9.92666867e-01],
       [1.77513660e-01, 1.83466350e-01, 7.99027540e-01, 6.30800256e-01],
       [5.14663640e-01, 6.34361690e-01, 5.33889110e-01, 7.90899958e-01],
       [2.58006640e-01, 2.88319290e-01, 7.09604700e-01, 9.02145588e-01],
       [2.04811730e-01, 1.58717810e-01, 8.41421970e-01, 8.84068574e-01],
       [3.11875950e-01, 2.46353420e-01, 7.58289460e-01, 9.65660849e-01],
       [9.36622730e-01, 8.01263020e-01, 9.83931900e-02, 5.44281251e-01],
       [3.45077880e-01, 2.88884330e-01, 6.99352220e-01, 9.71301027e-01],
       [8.18323020e-01, 8.42968360e-01, 2.68607890e-01, 1.52418342e-01],
       [9.07517590e-01, 7.41841580e-01, 3.94466860e-01, 3.33215046e-01],
       [3.84328210e-01, 2.85716010e-01, 7.73390420e-01, 4.79702183e-01],
       [4.74390660e-01, 7.16142340e-01, 4.52819790e-01, 7.98200958e-01],
       [1.91063150e-01, 3.11021770e-01, 8.17990970e-01, 3.65566296e-02],
       [3.84454330e-01, 5.31087070e-01, 6.46125720e-01, 4.26028784e-01],
       [7.76256310e-01, 9.82152950e-01, 2.24877600e-01, 9.03596071e-01],
       [5.21782070e-01, 6.09994810e-01, 5.49719820e-01, 3.77052128e-01],
       [2.44654800e-01, 2.22705930e-01, 6.92217350e-01, 8.04524950e-01],
       [6.26568380e-01, 6.85160450e-01, 5.05651420e-01, 1.06857236e-01],
       [4.57300600e-02, 3.19478800e-02, 9.56385030e-01, 5.62236853e-01],
       [9.85249100e-02, 9.00401200e-02, 9.15879990e-01, 9.37933793e-01],
       [1.13619610e-01, 1.08790380e-01, 8.86476490e-01, 9.07453097e-01],
       [8.61160760e-01, 9.26073490e-01, 1.51885700e-02, 5.61689264e-01],
       [4.72355650e-01, 8.90940390e-01, 4.79858960e-01, 1.31270153e-01],
       [3.29944340e-01, 3.87579980e-01, 5.10804930e-01, 5.19551698e-01],
       [5.29880000e-03, 7.01797000e-03, 9.93625440e-01, 1.85747216e-01],
       [6.23029300e-01, 4.84573370e-01, 2.68151590e-01, 5.60921564e-03],
       [8.03967200e-01, 7.84008540e-01, 2.37276020e-01, 9.47798098e-01],
       [4.58525700e-01, 6.45049070e-01, 5.76664620e-01, 5.75709041e-01],
       [8.65658770e-01, 8.99023760e-01, 2.40967370e-01, 5.56589158e-01],
       [9.67638170e-01, 8.83972000e-01, 3.44920000e-04, 9.17016810e-01],
       [7.24319710e-01, 8.63911350e-01, 2.64988220e-01, 2.13078474e-01],
       [9.91532810e-01, 9.85368470e-01, 6.58391400e-02, 4.58927178e-04],
       [1.17465420e-01, 1.14261700e-01, 9.05365980e-01, 1.80863318e-01],
       [8.63220080e-01, 7.91506140e-01, 5.18878970e-01, 5.49666344e-01],
       [3.12470000e-02, 3.19607200e-02, 9.72981330e-01, 8.60758375e-01],
       [4.41205810e-01, 4.87292710e-01, 5.39217100e-01, 5.65980037e-01],
       [1.39334500e-02, 1.40076300e-02, 9.84327060e-01, 9.38894626e-01],
       [2.43659030e-01, 2.09662260e-01, 7.81243000e-01, 1.14951150e-01],
       [9.45841720e-01, 9.45075580e-01, 1.37627400e-02, 5.39213927e-01],
       [7.92013050e-01, 6.55037130e-01, 2.08627580e-01, 1.50823215e-01],
       [6.04095200e-02, 4.57398400e-02, 9.53590740e-01, 6.32755639e-01],
       [5.67334500e-01, 2.75674320e-01, 6.47657510e-01, 4.68491101e-01],
       [1.58600060e-01, 1.22128390e-01, 8.47935330e-01, 4.94281577e-01],
       [8.26576000e-03, 4.07989000e-03, 9.95207870e-01, 8.02447365e-02],
       [7.25499790e-01, 7.05574910e-01, 3.94566850e-01, 8.90077195e-01],
       [8.30398180e-01, 7.65006390e-01, 1.04508490e-01, 9.44908637e-01],
       [6.84425700e-01, 8.13177160e-01, 2.55194010e-01, 1.71608600e-01],
       [3.69045830e-01, 4.10810940e-01, 6.39276590e-01, 9.22700243e-01],
       [1.42119190e-01, 1.45086850e-01, 8.87464990e-01, 2.35533293e-01],
       [8.51399930e-01, 8.63513030e-01, 4.28106000e-02, 7.49796027e-01],
       [7.26388760e-01, 9.23435870e-01, 3.21152410e-01, 5.59389176e-01],
       [2.68165680e-01, 2.21699530e-01, 7.13336850e-01, 8.28847266e-01],
       [4.67212100e-02, 6.18397600e-02, 9.08459550e-01, 1.73109978e-01],
       [8.12353540e-01, 6.14787930e-01, 2.36200930e-01, 6.70979632e-01],
       [3.56200600e-01, 2.86300900e-01, 6.87996620e-01, 7.68872468e-01],
       [4.27617260e-01, 4.08906890e-01, 4.65987670e-01, 1.67199623e-01],
       [6.63373240e-01, 9.66214910e-01, 1.39582640e-01, 9.85382902e-01],
       [5.51993350e-01, 4.93202560e-01, 5.63663960e-01, 1.69990831e-01],
       [8.04742160e-01, 7.23388830e-01, 1.97937550e-01, 5.06756753e-01],
       [1.07240370e-01, 1.15115720e-01, 9.07925810e-01, 3.46134208e-01],
       [3.61709450e-01, 2.16649010e-01, 7.91721970e-01, 5.22621049e-01],
       [9.83195600e-01, 9.35189250e-01, 1.09384140e-01, 4.87989100e-01],
       [1.07405620e-01, 1.05033440e-01, 8.76795260e-01, 2.44237928e-01],
       [6.75897130e-01, 6.50329960e-01, 3.04297580e-01, 3.60810270e-01],
       [7.02020600e-02, 4.96392100e-02, 9.33498520e-01, 7.17513612e-01],
       [4.84155500e-01, 6.88098980e-01, 3.46669530e-01, 2.16784063e-01],
       [6.04164790e-01, 7.48494480e-01, 9.49017500e-02, 2.69127829e-03],
       [5.92501140e-01, 7.18188940e-01, 4.79787090e-01, 4.72203718e-01],
       [6.47244640e-01, 9.12962170e-01, 3.94908800e-02, 1.89967176e-02],
       [7.52063710e-01, 8.36582980e-01, 2.56381510e-01, 1.82552057e-01],
       [7.33809600e-01, 5.88942430e-01, 3.17564930e-01, 4.83186793e-02],
       [6.37782580e-01, 7.91589180e-01, 3.08634220e-01, 1.83951279e-01],
       [7.32009020e-01, 9.14051250e-01, 1.80915920e-01, 2.45163585e-01],
       [1.53493780e-01, 1.90967590e-01, 8.19005590e-01, 7.55056039e-01],
       [5.36161820e-01, 5.13641150e-01, 5.01637010e-01, 3.47079632e-01],
       [6.06637230e-01, 6.67565790e-01, 3.33999130e-01, 2.51786198e-01],
       [7.25650010e-01, 8.41152620e-01, 2.36374270e-01, 2.61322095e-01],
       [6.52008490e-01, 8.66015010e-01, 1.90032370e-01, 5.14531432e-01],
       [2.59336300e-02, 3.60464100e-02, 9.42735970e-01, 8.76251330e-01],
       [3.91414850e-01, 3.16164320e-01, 6.36344310e-01, 2.11938819e-01],
       [6.43722130e-01, 5.38235890e-01, 1.13523690e-01, 3.54529909e-01],
       [7.90799970e-01, 7.44277280e-01, 3.24458070e-01, 1.60302427e-01],
       [7.10510700e-02, 8.50407900e-02, 9.08863250e-01, 9.18056054e-02],
       [8.27656880e-01, 7.68024600e-01, 6.35402600e-02, 1.39203186e-01],
       [4.22585470e-01, 4.66851210e-01, 5.36839920e-01, 9.51087042e-02],
       [8.74929100e-02, 9.12235300e-02, 8.91159090e-01, 3.88725280e-02],
       [9.36443830e-01, 8.16299420e-01, 2.90021130e-01, 2.89175878e-01],
       [9.26354550e-01, 9.51074570e-01, 8.49412000e-03, 1.76092602e-01],
       [5.72510800e-02, 3.56584500e-02, 9.67113730e-01, 7.74680782e-01],
       [1.10646890e-01, 1.10709490e-01, 8.85004310e-01, 8.08970193e-01],
       [9.30983140e-01, 9.85525760e-01, 6.47764700e-02, 3.51535913e-01],
       [9.16176650e-01, 8.13787300e-01, 2.80715970e-01, 7.69428516e-01],
       [9.34773620e-01, 8.73895270e-01, 1.23538120e-01, 5.72796569e-01]])

# Remove the dependence of the last dimension :
data[:,3] = np.random.uniform(size=100)

# use this specific breakpoint :
breakpoint = np.array([0.57425735, 0.52599749, 0.43842457, 0.66851334])

# Dimensions
d = data.shape[1] # Number of variables
M = 2 ** d
n = data.shape[0] # Number of observationsa
N = 999 # Number of bootstrap resamples

# compute prerequisites : Thoose prerequistes dont need optimisation.
binary_repr = np.arange(2 ** d)[:, None] >> np.arange(d)[::-1] & 1 # this is already optimized
min = breakpoint * binary_repr  # (M,d)
max = breakpoint ** (1 - binary_repr)  # (M,d)
lambda_l = np.prod(max - min, axis=1)  # (M,)

########### Algorithm
# small function : the 3 imputs needs to be broadcastable to each other,
# and it returns something that is shaped as the broadcast of the 3 inputs
def is_in(data,min,max): return np.logical_and(data >= min, data < max)

def compute_statistic(A, B, C): return np.sum(B ** 2 / C - 2 * A * B,axis=0)

# Now we get to the algorithm itself.
are_in = is_in(data[None,], min[:, None, ], max[:, None, ])  # (M,n,d)

f_l = np.mean(np.all(are_in, axis=2), axis=1)  # (M,)

# Draw random numbers :
random_numbers = np.random.uniform(size=(N, n, d))  # (N,n,d)

# build indices selectors :
masks = (np.identity(d) == 0).reshape(d * d)  # (d*d)
rev_masks = np.logical_not(masks)

# simplify broadcasting by repeating a little the data : # <<<<<<<<<<<< Here problems might start !
maxes = np.repeat(max, d, axis=1)  # (M,d*d)
mines = np.repeat(min, d, axis=1)  # (M,d*d)
zes = np.repeat(data, d, axis=1)  # (n,d*d)
random_numbers = np.repeat(random_numbers, d, axis=2)
are_in_rep = np.repeat(are_in, d, axis=2)  # (M,n,d*d)

# Compute a piece of the result :
f_k = np.mean(np.all(are_in_rep[..., masks].reshape(M, n, d - 1, d), axis=2), axis=1)  # (M,d)
lambda_k = np.prod((maxes - mines)[:, masks].reshape((M, d - 1, d)), axis=1)  # (M,d)
rez_k = np.flip(f_k / lambda_k,axis=1)

# The observed values of the statistic can now be computed :
value = compute_statistic(rez_k, f_l[:,None], lambda_l[:, None])  # (d,)

# Now compute the same statistic value, but bootstrapped :
z_reps = np.repeat(zes[None, :, :], N, axis=0)  # (N,n,d*d)
z_reps[:, :, rev_masks] = random_numbers[:, :, rev_masks]  # (N,n,d*d)
are_in_boot = is_in(z_reps[None,], mines[:, None, None, ], maxes[:, None, None, ]).reshape(M, N, n, d, d)
f_ls = np.mean(np.all(are_in_boot, axis=3), axis=2)  # (M,N,d)
values = compute_statistic(rez_k[:,None],f_ls,lambda_l[:, None, None])  # (N,d)

# return the result :
result = np.mean(value < values, axis=0)
print(result)

Edit: I tried to split the code into smaller functions, and run cProfiler on it. The code is below. The problem seems to come from the compute_f_l function which is eating more than half the runtime. But this function is soo trivial that i dont know what to do... The is_in function eats 1/4 of the time, same issue for me.

Here is the profiled and splitted up code :

import numpy as np
import cProfile
cp = cProfile.Profile()
cp.enable()

########### Data and prerequisites
# Sample of data with strong dependence structure :
data = np.array([[5.61293390e-01, 9.97100450e-01, 4.23530180e-01, 6.08808896e-01],
       [1.22563280e-01, 1.72015130e-01, 8.71145720e-01, 5.40745844e-01],
       [8.51194500e-02, 1.18289130e-01, 8.90346540e-01, 7.22859351e-01],
       [9.83241310e-01, 9.57282690e-01, 7.22347100e-02, 5.43527399e-02],
       [5.49211550e-01, 3.97858250e-01, 6.86380990e-01, 7.91494336e-01],
       [9.94878920e-01, 6.39160920e-01, 2.01045170e-01, 9.86840712e-01],
       [5.04337540e-01, 5.69995040e-01, 3.99087430e-01, 4.32140476e-01],
       [9.28230540e-01, 9.32143440e-01, 1.02748280e-01, 9.92666867e-01],
       [1.77513660e-01, 1.83466350e-01, 7.99027540e-01, 6.30800256e-01],
       [5.14663640e-01, 6.34361690e-01, 5.33889110e-01, 7.90899958e-01],
       [2.58006640e-01, 2.88319290e-01, 7.09604700e-01, 9.02145588e-01],
       [2.04811730e-01, 1.58717810e-01, 8.41421970e-01, 8.84068574e-01],
       [3.11875950e-01, 2.46353420e-01, 7.58289460e-01, 9.65660849e-01],
       [9.36622730e-01, 8.01263020e-01, 9.83931900e-02, 5.44281251e-01],
       [3.45077880e-01, 2.88884330e-01, 6.99352220e-01, 9.71301027e-01],
       [8.18323020e-01, 8.42968360e-01, 2.68607890e-01, 1.52418342e-01],
       [9.07517590e-01, 7.41841580e-01, 3.94466860e-01, 3.33215046e-01],
       [3.84328210e-01, 2.85716010e-01, 7.73390420e-01, 4.79702183e-01],
       [4.74390660e-01, 7.16142340e-01, 4.52819790e-01, 7.98200958e-01],
       [1.91063150e-01, 3.11021770e-01, 8.17990970e-01, 3.65566296e-02],
       [3.84454330e-01, 5.31087070e-01, 6.46125720e-01, 4.26028784e-01],
       [7.76256310e-01, 9.82152950e-01, 2.24877600e-01, 9.03596071e-01],
       [5.21782070e-01, 6.09994810e-01, 5.49719820e-01, 3.77052128e-01],
       [2.44654800e-01, 2.22705930e-01, 6.92217350e-01, 8.04524950e-01],
       [6.26568380e-01, 6.85160450e-01, 5.05651420e-01, 1.06857236e-01],
       [4.57300600e-02, 3.19478800e-02, 9.56385030e-01, 5.62236853e-01],
       [9.85249100e-02, 9.00401200e-02, 9.15879990e-01, 9.37933793e-01],
       [1.13619610e-01, 1.08790380e-01, 8.86476490e-01, 9.07453097e-01],
       [8.61160760e-01, 9.26073490e-01, 1.51885700e-02, 5.61689264e-01],
       [4.72355650e-01, 8.90940390e-01, 4.79858960e-01, 1.31270153e-01],
       [3.29944340e-01, 3.87579980e-01, 5.10804930e-01, 5.19551698e-01],
       [5.29880000e-03, 7.01797000e-03, 9.93625440e-01, 1.85747216e-01],
       [6.23029300e-01, 4.84573370e-01, 2.68151590e-01, 5.60921564e-03],
       [8.03967200e-01, 7.84008540e-01, 2.37276020e-01, 9.47798098e-01],
       [4.58525700e-01, 6.45049070e-01, 5.76664620e-01, 5.75709041e-01],
       [8.65658770e-01, 8.99023760e-01, 2.40967370e-01, 5.56589158e-01],
       [9.67638170e-01, 8.83972000e-01, 3.44920000e-04, 9.17016810e-01],
       [7.24319710e-01, 8.63911350e-01, 2.64988220e-01, 2.13078474e-01],
       [9.91532810e-01, 9.85368470e-01, 6.58391400e-02, 4.58927178e-04],
       [1.17465420e-01, 1.14261700e-01, 9.05365980e-01, 1.80863318e-01],
       [8.63220080e-01, 7.91506140e-01, 5.18878970e-01, 5.49666344e-01],
       [3.12470000e-02, 3.19607200e-02, 9.72981330e-01, 8.60758375e-01],
       [4.41205810e-01, 4.87292710e-01, 5.39217100e-01, 5.65980037e-01],
       [1.39334500e-02, 1.40076300e-02, 9.84327060e-01, 9.38894626e-01],
       [2.43659030e-01, 2.09662260e-01, 7.81243000e-01, 1.14951150e-01],
       [9.45841720e-01, 9.45075580e-01, 1.37627400e-02, 5.39213927e-01],
       [7.92013050e-01, 6.55037130e-01, 2.08627580e-01, 1.50823215e-01],
       [6.04095200e-02, 4.57398400e-02, 9.53590740e-01, 6.32755639e-01],
       [5.67334500e-01, 2.75674320e-01, 6.47657510e-01, 4.68491101e-01],
       [1.58600060e-01, 1.22128390e-01, 8.47935330e-01, 4.94281577e-01],
       [8.26576000e-03, 4.07989000e-03, 9.95207870e-01, 8.02447365e-02],
       [7.25499790e-01, 7.05574910e-01, 3.94566850e-01, 8.90077195e-01],
       [8.30398180e-01, 7.65006390e-01, 1.04508490e-01, 9.44908637e-01],
       [6.84425700e-01, 8.13177160e-01, 2.55194010e-01, 1.71608600e-01],
       [3.69045830e-01, 4.10810940e-01, 6.39276590e-01, 9.22700243e-01],
       [1.42119190e-01, 1.45086850e-01, 8.87464990e-01, 2.35533293e-01],
       [8.51399930e-01, 8.63513030e-01, 4.28106000e-02, 7.49796027e-01],
       [7.26388760e-01, 9.23435870e-01, 3.21152410e-01, 5.59389176e-01],
       [2.68165680e-01, 2.21699530e-01, 7.13336850e-01, 8.28847266e-01],
       [4.67212100e-02, 6.18397600e-02, 9.08459550e-01, 1.73109978e-01],
       [8.12353540e-01, 6.14787930e-01, 2.36200930e-01, 6.70979632e-01],
       [3.56200600e-01, 2.86300900e-01, 6.87996620e-01, 7.68872468e-01],
       [4.27617260e-01, 4.08906890e-01, 4.65987670e-01, 1.67199623e-01],
       [6.63373240e-01, 9.66214910e-01, 1.39582640e-01, 9.85382902e-01],
       [5.51993350e-01, 4.93202560e-01, 5.63663960e-01, 1.69990831e-01],
       [8.04742160e-01, 7.23388830e-01, 1.97937550e-01, 5.06756753e-01],
       [1.07240370e-01, 1.15115720e-01, 9.07925810e-01, 3.46134208e-01],
       [3.61709450e-01, 2.16649010e-01, 7.91721970e-01, 5.22621049e-01],
       [9.83195600e-01, 9.35189250e-01, 1.09384140e-01, 4.87989100e-01],
       [1.07405620e-01, 1.05033440e-01, 8.76795260e-01, 2.44237928e-01],
       [6.75897130e-01, 6.50329960e-01, 3.04297580e-01, 3.60810270e-01],
       [7.02020600e-02, 4.96392100e-02, 9.33498520e-01, 7.17513612e-01],
       [4.84155500e-01, 6.88098980e-01, 3.46669530e-01, 2.16784063e-01],
       [6.04164790e-01, 7.48494480e-01, 9.49017500e-02, 2.69127829e-03],
       [5.92501140e-01, 7.18188940e-01, 4.79787090e-01, 4.72203718e-01],
       [6.47244640e-01, 9.12962170e-01, 3.94908800e-02, 1.89967176e-02],
       [7.52063710e-01, 8.36582980e-01, 2.56381510e-01, 1.82552057e-01],
       [7.33809600e-01, 5.88942430e-01, 3.17564930e-01, 4.83186793e-02],
       [6.37782580e-01, 7.91589180e-01, 3.08634220e-01, 1.83951279e-01],
       [7.32009020e-01, 9.14051250e-01, 1.80915920e-01, 2.45163585e-01],
       [1.53493780e-01, 1.90967590e-01, 8.19005590e-01, 7.55056039e-01],
       [5.36161820e-01, 5.13641150e-01, 5.01637010e-01, 3.47079632e-01],
       [6.06637230e-01, 6.67565790e-01, 3.33999130e-01, 2.51786198e-01],
       [7.25650010e-01, 8.41152620e-01, 2.36374270e-01, 2.61322095e-01],
       [6.52008490e-01, 8.66015010e-01, 1.90032370e-01, 5.14531432e-01],
       [2.59336300e-02, 3.60464100e-02, 9.42735970e-01, 8.76251330e-01],
       [3.91414850e-01, 3.16164320e-01, 6.36344310e-01, 2.11938819e-01],
       [6.43722130e-01, 5.38235890e-01, 1.13523690e-01, 3.54529909e-01],
       [7.90799970e-01, 7.44277280e-01, 3.24458070e-01, 1.60302427e-01],
       [7.10510700e-02, 8.50407900e-02, 9.08863250e-01, 9.18056054e-02],
       [8.27656880e-01, 7.68024600e-01, 6.35402600e-02, 1.39203186e-01],
       [4.22585470e-01, 4.66851210e-01, 5.36839920e-01, 9.51087042e-02],
       [8.74929100e-02, 9.12235300e-02, 8.91159090e-01, 3.88725280e-02],
       [9.36443830e-01, 8.16299420e-01, 2.90021130e-01, 2.89175878e-01],
       [9.26354550e-01, 9.51074570e-01, 8.49412000e-03, 1.76092602e-01],
       [5.72510800e-02, 3.56584500e-02, 9.67113730e-01, 7.74680782e-01],
       [1.10646890e-01, 1.10709490e-01, 8.85004310e-01, 8.08970193e-01],
       [9.30983140e-01, 9.85525760e-01, 6.47764700e-02, 3.51535913e-01],
       [9.16176650e-01, 8.13787300e-01, 2.80715970e-01, 7.69428516e-01],
       [9.34773620e-01, 8.73895270e-01, 1.23538120e-01, 5.72796569e-01]])

# Remove the dependence of the last dimension :
data[:,3] = np.random.uniform(size=100)

# use this specific breakpoint :
breakpoint = np.array([0.57425735, 0.52599749, 0.43842457, 0.66851334])

# Dimensions
d = data.shape[1] # Number of variables
M = 2 ** d
n = data.shape[0] # Number of observationsa
N = 999 # Number of bootstrap resamples

# compute prerequisites : Thoose prerequistes dont need optimisation.
binary_repr = np.arange(2 ** d)[:, None] >> np.arange(d)[::-1] & 1 # this is already optimized
min = breakpoint * binary_repr  # (M,d)
max = breakpoint ** (1 - binary_repr)  # (M,d)
lambda_l = np.prod(max - min, axis=1)  # (M,)

########### Algorithm
# small function : the 3 imputs needs to be broadcastable to each other,
# and it returns something that is shaped as the broadcast of the 3 inputs


def compute_statistic(A, B, C): return np.sum(B ** 2 / C - 2 * A * B,axis=0)

def is_in(data,min,max): return np.logical_and(data >= min, data < max)

def compute_f_l(are_in,axes = (2,1)):
    return np.mean(np.all(are_in, axis=axes[0]), axis=axes[1])

def compute_repeats(max,min,data,random_numbers,are_in):
    maxes = np.repeat(max, d, axis=1)  # (M,d*d)
    mines = np.repeat(min, d, axis=1)  # (M,d*d)
    zes = np.repeat(data, d, axis=1)  # (n,d*d)
    random_numbers = np.repeat(random_numbers, d, axis=2)
    are_in_rep = np.repeat(are_in, d, axis=2)  # (M,n,d*d)

    return maxes,mines,zes,random_numbers,are_in_rep

def compute_first_piece(are_in_rep,masks,maxes,mines):
    f_k = np.mean(np.all(are_in_rep[..., masks].reshape(M, n, d - 1, d), axis=2), axis=1)  # (M,d)
    lambda_k = np.prod((maxes - mines)[:, masks].reshape((M, d - 1, d)), axis=1)  # (M,d)
    rez_k = np.flip(f_k / lambda_k, axis=1)

    return rez_k

def compute_z_reps(zes,rev_masks,random_numbers):
    z_reps = np.repeat(zes, N, axis=0)  # (N,n,d*d)
    z_reps[:, :, rev_masks] = random_numbers[:, :, rev_masks]  # (N,n,d*d)
    return z_reps

def compute_bootstrap(z_reps,mines,maxes,rez_k,lambda_l):

    are_in_boot = is_in(z_reps, mines, maxes).reshape(M, N, n, d, d)
    f_ls = compute_f_l(are_in_boot,axes=(3,2)) # (M,N,d)
    values = compute_statistic(rez_k,f_ls,lambda_l)  # (N,d)
    return values

for i in np.arange(20): # just repeat 20 times to augment values of times.

    # Now we get to the algorithm itself.
    are_in = is_in(data[None,], min[:, None, ], max[:, None, ])  # (M,n,d)
    f_l = compute_f_l(are_in)  # (M,)

    # Draw random numbers :
    random_numbers = np.random.uniform(size=(N, n, d))  # (N,n,d)

    # build indices selectors :
    masks = (np.identity(d) == 0).reshape(d * d)  # (d*d)
    rev_masks = np.logical_not(masks)

    # simplify broadcasting by repeating a little the data : # <<<<<<<<<<<< Here problems might start !
    maxes,mines,zes,random_numbers,are_in_rep = compute_repeats(max,min,data,random_numbers,are_in)

    # Compute a piece of the result :
    rez_k = compute_first_piece(are_in_rep,masks,maxes,mines)

    # The observed values of the statistic can now be computed :
    value = compute_statistic(rez_k, f_l[:,None], lambda_l[:, None])  # (d,)

    z_reps = compute_z_reps(zes[None, :, :],
                            rev_masks,
                            random_numbers)

    # Now compute the same statistic value, but bootstrapped :
    values = compute_bootstrap(z_reps[None,],
                               mines[:, None, None, ],
                               maxes[:, None, None, ],
                               rez_k[:,None],
                               lambda_l[:, None, None])

    # return the result :
    result = np.mean(value < values, axis=0)
    print(result)


cp.disable()
cp.print_stats()

Edit : Was originally posted on SO but after comments we thought it was better suited for code review.

\$\endgroup\$
1
  • \$\begingroup\$ You can try "Numba" but it's rare it will do any good since everything seems to be vectorized here. There's although few places where "Numba" stands out since it optimizes your python code. Otherwise tensorflow and gpu. Not an answer since I don't know if that will work. \$\endgroup\$ Jul 10, 2020 at 15:45

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.