fork download
copy 
  1. # One-dimensional linear interpolation. (2.01)
  2. # @note The product of trying to figure out the behavior of numpy.interp
  3. # after reading the documentation. I think this captures the gist of it.
  4.  
  5. import math
  6. import bisect
  7. import numbers
  8.  
  9. # Utility.
  10.  
  11. def is_scalar(x):
  12.     return isinstance(x, numbers.Number)
  13.  
  14. def is_least_float(x):
  15.     return is_scalar(x) and not isinstance(x, numbers.Rational)
  16.  
  17. def to_least_float(x):
  18.     return x if is_least_float(x) else float(x)
  19.  
  20. def lerp(a, b, t):
  21.     return (1 - t) * a + t * b
  22.  
  23. def invlerp(x, a, b):
  24.     if math.isclose(a, b):
  25.         return 0.5
  26.     return (x - a) / (b - a)
  27.  
  28. def modinvlerp(x, a, b, m):
  29.     d = (b - a) % m
  30.     if math.isclose(0, d):
  31.         return 0.5
  32.     return (x - a) % m / d
  33.  
  34. def argsort(seq):
  35.     return sorted(range(len(seq)), key=seq.__getitem__)
  36.  
  37. # Interpolate.
  38.  
  39. def interp(x, xp, fp, left=None, right=None, period=None):
  40.     """
  41.     Linear interpolation for monotonically increasing sample points.
  42.     """
  43.     n = len(xp)
  44.  
  45.     if n != len(fp):
  46.         raise ValueError('xp and fp are not of the same length')
  47.     if n == 0:
  48.         return
  49.  
  50.     xp = list(map(to_least_float, xp))
  51.     fp = list(map(to_least_float, fp))
  52.  
  53.     if period is None:
  54.         return _interp_bounded(n, x, xp, fp, left, right)
  55.     return _interp_periodic(n, x, xp, fp, period)
  56.  
  57. def _interp_bounded(n, x, xp, fp, left, right):
  58.     if left is None:
  59.         left = fp[0]
  60.     else:
  61.         left = to_least_float(left)
  62.     if right is None:
  63.         right = fp[-1]
  64.     else:
  65.         right = to_least_float(right)
  66.  
  67.     first, last = xp[0], xp[-1]
  68.  
  69.     for y in x:
  70.         if y < first:
  71.             yield left
  72.         elif y > last:
  73.             yield right
  74.         else:
  75.             i = bisect.bisect(xp, y)
  76.             j = i-1
  77.             k = min(i, n-1)
  78.             yield lerp(fp[j], fp[k], invlerp(y, xp[j], xp[k]))
  79.  
  80. def _interp_periodic(n, x, xp, fp, period):
  81.     if period == 0:
  82.         raise ValueError('period must be a non-zero value')
  83.  
  84.     # Normalize periodic boundaries; reorder xp and fp.
  85.  
  86.     xp = [x % period for x in xp]
  87.     xp_sorted_index = argsort(xp)
  88.     xp = [xp[i] for i in xp_sorted_index]
  89.     fp = [fp[i] for i in xp_sorted_index]
  90.  
  91.     for y in x:
  92.         i = bisect.bisect(xp, y % period)
  93.         j = (i-1) % n
  94.         k = i % n
  95.         yield lerp(fp[j], fp[k], modinvlerp(y, xp[j], xp[k], period))
  96.  
  97. # Main.
  98.  
  99. import numpy as np
  100.  
  101. x  = [0, 1, 1.5, 2.73, 3, 3.14]
  102. xp = [1, 2, 3]
  103. fp = [3, 2, 0]
  104.  
  105. print(list(np.interp([2.5], xp, fp)))
  106. print(list(interp([2.5], xp, fp)))
  107. print(list(np.interp(x, xp, fp)))
  108. print(list(interp(x, xp, fp)))
  109.  
  110. UNDEF = 99
  111. print(list(np.interp(x, xp, fp, left=-UNDEF, right=UNDEF)))
  112. print(list(interp(x, xp, fp, left=-UNDEF, right=UNDEF)))
  113.  
  114. x  = [-180, -170, -185, 185, -10, -5, 0, 365]
  115. xp = [190, -190, 350, -350]
  116. fp = [5, 10, 3, 4]
  117.  
  118. print(list(np.interp(x, xp, fp, period=360)))
  119. print(list(interp(x, xp, fp, period=360)))
  120.  
  121. x  = [1.5, 4.0]
  122. xp = [2, 3, 5]
  123. fp = [1.0j, 0, 2+3j]
  124.  
  125. print(list(np.interp(x, xp, fp)))
  126. print(list(interp(x, xp, fp)))
Success #stdin #stdout 0.81s 41608KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[1.0]
[1.0]
[3.0, 3.0, 2.5, 0.54, 0.0, 0.0]
[3.0, 3.0, 2.5, 0.54, 0.0, 0.0]
[-99.0, 3.0, 2.5, 0.54, 0.0, 99.0]
[-99.0, 3.0, 2.5, 0.54, 0.0, 99.0]
[7.5, 5.0, 8.75, 6.25, 3.0, 3.25, 3.5, 3.75]
[7.5, 5.0, 8.75, 6.25, 3.0, 3.25, 3.5, 3.75]
[1j, (1+1.5j)]
[1j, (1+1.5j)]
https://ideone.com/VWT1ap
language:
Python 3 (python  3.12)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!