Orthogonal sampling

Question

I have been staring at what I have produced for almost 5 hours and I still cannot see how and where to improve my implementation. I am implementing an Orthogonal sampling method according to the description of Orthogonal Sampling.

import numpy as np
from numba import jit
import random
def Orthogonal(n):
   assert(np.sqrt(n) % 1 == 0),"Please insert an even number of samples"
   step = int(np.sqrt(n))

   row_constraints = [ False for i in range(n) ]
   column_constraints = [ False for i in range(n) ]
   subbox_constraints = [ False for i in range(n)]

   result = []
   append = result.append

   def map_coord_to_box(x:int, y:int, step:int)->int:
       horz = int(np.floor(x/step))
       vert = int(np.floor(y/step))

       return vert * step + horz

   def remap_value(value:float, olb:int, oub:int, nlb:int, nub:int)->int:
# https://stackoverflow.com/questions/1969240/mapping-a-range-of-values-to-another
       delta_old = abs(oub - olb)
       delta_new = abs(nub - nlb)

       ratio = (value - olb) / delta_old

       return ratio * delta_new + nlb

   while(all(subbox_constraints) == False):
         value = np.random.uniform(low=-2, high=2, size=2)
         x = int(np.floor(remap_value(value[0], -2, 2, 0, n)))
         y = int(np.floor(remap_value(value[1], -2, 2, 0, n)))
         if not (row_constraints[y] or column_constraints[x] or subbox_constraints[map_coord_to_box(x, y, step)]): #check_constraints(row_constraints, column_constraints, subbox_constraints, x, y, step): 
              append(tuple(value))
              row_constraints[y] = True
              column_constraints[x] = True
              subbox_constraints[map_coord_to_box(x, y, step)] = True    

   return result

The problem is obvious when generating 100 samples it takes on average 300 ms, and I need something faster as I need to generate at least 10.000 samples. So I have not sat still. I tried to use jit for the sub-functions but it does not make it faster, but slower. I am aware that these function calls in python have a higher overhead. And so far on my own I thought using these functions are a way to approach the sampling method I want to implement. I have also asked a friend and he came up with a different approach which is on average a factor 100 faster than the above code. So he only prunes every row and columns and after randomly choosing those and stores the indices in to a list which later fills randomly.

def orthogonal_l(n):
    bs = int(np.sqrt(n))
    result = [0 for i in range(n)]
    columns = [[i for i in range(bs)] for j in range(bs)]
    c = 0
    for h in range(bs):
        rows = [i for i in range(bs)]
        for z in range(bs):
            w = random.choice(rows)
            x = random.choice(columns[w])
            columns[w].remove(x)
            rows.remove(w)
            x += w * bs
            result[x] = c 
            c += 1
   return result, bs

How can I make use of pruning with my own code and is it wise to do so? If not, how can I improve the code, if so, where?

Answer 1

So after talking with other students and some drawings I realised that I have three choices: use the definition of my friend which is faster than my original function, use my original function OR create an efficient datastructure to deal with the look-up complexity.

So here goes the fresh function which in comparison with my old function is a factor 100 time faster:

def another_Orthogonal(ns):
   assert(np.sqrt(ns) % 1 == 0),"Please insert an even number of samples"
   n = int(np.sqrt(ns))
   # Making a datastructure of a dict with coordinate tuples of a bigger grid with subcoordinate of sub-grid points
   blocks = {(i,j):[(a,b) for a in range(n) for b in range(n)] for i in range(n) for j in range(n)}
   points = []#np.empty((n,2))
   append = points.append # tips of python to fasten up append call
   for block in blocks:
       point = random.choice(blocks[block])
       lst_row = [(k1, b) for (k1, b), v in blocks.items() if k1 == block[0]]
       lst_col = [(a, k1) for (a, k1), v in blocks.items() if k1 == block[1]]

       for col in lst_col:
           blocks[col] = [a for a in blocks[col] if a[1] != point[1]]

       for row in lst_row:
           blocks[row] = [a for a in blocks[row] if a[0] != point[0]]
      #Adjust the points to fit the grid they fall in  
      point = (point[0] + n * block[0], point[1] + n * block[1])
      append(point)

    return points

So the function does return only the coordinates of where the randomly chosen point is set and where no other points should be set by looking at at the row or column coordinates of the block and eliminating those option as choice.

The random value themselves on these points are calculated in another function that scales the values of the interval (0, numberOfpoints) to a desired interval where one could sample from.

The scaling function is the following :

def scale_points(points):
    p = another(points)
    maximum = points 
    scaling =[ 1/maximum * i for i in range(len(p))]
    min_ = -2.0
    max_ = 2.0
    result = np.zeros((points,2))
    anti_res = np.zeros((points,2)) # this is for antithetic variables

    for idx, scale in enumerate(scaling):

        x =  min_ + np.random.uniform(p[idx][0]/maximum, p[idx][0]/maximum +1/maximum ) *4  # 4 is just max - min which is in my case 4
        y =  min_ + np.random.uniform(p[idx][1]/maximum, p[idx][1]/maximum + 1/maximum ) * 4
        result[idx, :] = [x,y]
        anti_res[idx,:] = [x*-1.0, y*-1.0] # antithetic variables

    return result, anti_res