title: Python - Quick Upscaling of Array with Numpy, No Image Libary Allowed
body: <p><strong>Note on duplicate message:</strong></p>
  
  <p>Similar themes, not exactly a duplicate. Esp. since the loop is still the fastest method. Thanks.</p>
  
  <p><strong>Goal:</strong> </p>
  
  <p>Upscale an array from [small,small] to [big,big] by a factor quickly, don't use an image library. Very simple scaling, one small value will become several big values, after it is normalized for the several big values it becomes. In other words, this is "flux conserving" from an astronomical wording - a value of 16 from the small array spread into a big array's 4 values (factor of 2) would be 4 4's so the amount of the value has been retained.</p>
  
  <p><strong>Problem:</strong></p>
  
  <p>I've got some working codes to do the upscaling, but they don't work very fast compared to downscaling. Upscaling is actually easier than downscaling (which requires many sums, in this basic case) - upscaling just requires already-known data to be put in big chunks of a preallocated array.</p>
  
  <p>For a working example, a [2,2] array of [16,24;8,16]:</p>
  
  <blockquote>
    <p>16 , 24</p>
    
    <p>8 , 16</p>
  </blockquote>
  
  <p>Multiplied by a factor of 2 for a [4,4] array would have the values:</p>
  
  <blockquote>
    <p>4 , 4 , 6 , 6</p>
    
    <p>4 , 4 , 6 , 6</p>
    
    <p>2 , 2 , 4 , 4</p>
    
    <p>2 , 2 , 4 , 4</p>
  </blockquote>
  
  <p>The fastest implementation is a for loop accelerated by numba's jit &amp; prange. I'd like to better leverage Numpy's pre-compiled functions to get this job done. I'll also entertain Scipy stuff - but not its resizing functions.</p>
  
  <p>It seems like a perfect problem for strong matrix manipulation functions, but I just haven't managed to make it happen quickly.</p>
  
  <p>Additionally, the single-line numpy call is <strong>way</strong> funky, so don't be surprized. But it's what it took to get it to align correctly.</p>
  
  <p><strong>Code examples:</strong></p>
  
  <p><strong>Check more optimized calls below</strong> Be warned, the case I have here makes a 20480x20480 float64 array that can take up a fair bit of memory - but can show off if a method is too memory intensive (as matrices can be).</p>
  
  <p>Environment: Python 3, Windows, i5-4960K @ 4.5 GHz. Time to run for loop code is ~18.9 sec, time to run numpy code is ~52.5 sec on the shown examples.</p>
  
  <p>% MAIN: To run these</p>
  
  <pre><code>import timeit
  
  timeitSetup = ''' 
  from Regridder1 import Regridder1
  import numpy as np
  
  factor = 10;
  
  inArrayX = np.float64(np.arange(0,2048,1));
  inArrayY = np.float64(np.arange(0,2048,1));
  [inArray, _] = np.meshgrid(inArrayX,inArrayY);
  ''';
  
  print("Time to run 1: {}".format( timeit.timeit(setup=timeitSetup,stmt="Regridder1(inArray, factor,)", number = 10) ));
  
  timeitSetup = ''' 
  from Regridder2 import Regridder2
  import numpy as np
  
  factor = 10;
  
  inArrayX = np.float64(np.arange(0,2048,1));
  inArrayY = np.float64(np.arange(0,2048,1));
  [inArray, _] = np.meshgrid(inArrayX,inArrayY);
  ''';
  
  print("Time to run 2: {}".format( timeit.timeit(setup=timeitSetup,stmt="Regridder2(inArray, factor,)", number = 10) ));
  </code></pre>
  
  <p>% FUN: Regridder 1 - for loop</p>
  
  <pre><code>import numpy as np
  from numba import prange, jit
  
  @jit(nogil=True)
  def Regridder1(inArray,factor):
      inSize = np.shape(inArray);
      outSize = [np.int64(np.round(inSize[0] * factor)), np.int64(np.round(inSize[1] * factor))];
  
      outBlockSize = factor*factor; #the block size where 1 inArray pixel is spread across # outArray pixels
      outArray = np.zeros(outSize); #preallcoate
      outBlocks = inArray/outBlockSize; #precalc the resized blocks to go faster
      for i in prange(0,inSize[0]):
          for j in prange(0,inSize[1]):
              outArray[i*factor:(i*factor+factor),j*factor:(j*factor+factor)] = outBlocks[i,j]; #puts normalized value in a bunch of places
  
      return outArray;
  </code></pre>
  
  <p>% FUN: Regridder 2 - numpy</p>
  
  <pre><code>import numpy as np
  
  def Regridder2(inArray,factor):
      inSize = np.shape(inArray);
      outSize = [np.int64(np.round(inSize[0] * factor)), np.int64(np.round(inSize[1] * factor))];
  
      outBlockSize = factor*factor; #the block size where 1 inArray pixel is spread across # outArray pixels
  
      outArray = inArray.repeat(factor).reshape(inSize[0],factor*inSize[1]).T.repeat(factor).reshape(inSize[0]*factor,inSize[1]*factor).T/outBlockSize;
  
      return outArray;
  </code></pre>
  
  <p>Would greatly appreciate insight into speeding this up. Hopefully code is good, formulated it in the text box.</p>
  
  <p><strong>Current best solution:</strong></p>
  
  <p>On my comp, the numba's jit for loop implementation (Regridder1) with jit applied to only what needs it can run the timeit test at 18.0 sec, while the numpy only implementation (Regridder2) runs the timeit test at 18.5 sec. The bonus is that on the first call, the numpy only implementation doesn't need to wait for jit to compile the code. Jit's cache=True lets it not compile on subsequent runs. The other calls (nogil, nopython, prange) don't seem to help but also don't seem to hurt. Maybe in future numba updates they'll do better or something. </p>
  
  <p>For simplicity and portability, Regridder2 is the best option. It's nearly as fast, and doesn't need numba installed (which for my Anaconda install required me to go install it) - so it'll help portability.</p>
  
  <p>% FUN: Regridder 1 - for loop</p>
  
  <pre><code>import numpy as np
  
  def Regridder1(inArray,factor):
      inSize = np.shape(inArray);
      outSize = [np.int64(np.round(inSize[0] * factor)), np.int64(np.round(inSize[1] * factor))];
  
      outBlockSize = factor*factor #the block size where 1 inArray pixel is spread across # outArray pixels
      outArray = np.empty(outSize) #preallcoate
      outBlocks = inArray/outBlockSize #precalc the resized blocks to go faster
      factor = np.int64(factor) #convert to an integer to be safe (in case it's a 1.0 float)
  
      outArray = RegridderUpscale(inSize, factor, outArray, outBlocks) #call a function that has just the loop
  
      return outArray;
  #END def Regridder1
  
  from numba import jit, prange
  @jit(nogil=True, nopython=True, cache=True) #nopython=True, nogil=True, parallel=True, cache=True
  def RegridderUpscale(inSize, factor, outArray, outBlocks ):
      for i in prange(0,inSize[0]):
          for j in prange(0,inSize[1]):
              outArray[i*factor:(i*factor+factor),j*factor:(j*factor+factor)] = outBlocks[i,j];
          #END for j
      #END for i
      #scales the original data up, note for other languages you need i*factor+factor-1 because slicing
      return outArray; #return success
  #END def RegridderUpscale
  </code></pre>
  
  <p>% FUN: Regridder 2 - numpy based on @ZisIsNotZis's answer</p>
  
  <pre><code>import numpy as np
  
  def Regridder2(inArray,factor):
      inSize = np.shape(inArray);
      #outSize = [np.int64(np.round(inSize[0] * factor)), np.int64(np.round(inSize[1] * factor))]; #whoops
  
      outBlockSize = factor*factor; #the block size where 1 inArray pixel is spread across # outArray pixels
  
      outArray = np.broadcast_to( inArray[:,None,:,None]/outBlockSize, (inSize[0], factor, inSize[1], factor)).reshape(np.int64(factor*inSize[0]), np.int64(factor*inSize[1])); #single line call that gets the job done
  
      return outArray;
  #END def Regridder2
  </code></pre>