fork(1) download 
  1. # One-dimensional linear interpolation. (2.00)
  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 + b * t
  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.     t = (b - a) % m
  30.     if math.isclose(0, t):
  31.         return 0.5
  32.     return (x - a) % m / t
  33.  
  34. # Interpolate.
  35.  
  36. def interp(x, xp, fp, left=None, right=None, period=None):
  37.     """
  38.     Linear interpolation for monotonically increasing sample points.
  39.     """
  40.     n = len(xp)
  41.  
  42.     if n != len(fp):
  43.         raise ValueError('xp and fp are not of the same length')
  44.     if n == 0:
  45.         return
  46.  
  47.     xp = list(map(to_least_float, xp))
  48.     fp = list(map(to_least_float, fp))
  49.  
  50.     if period is None:
  51.         return _interp_bounded(n, x, xp, fp, left, right)
  52.     return _interp_periodic(n, x, xp, fp, period)
  53.  
  54. def _interp_bounded(n, x, xp, fp, left, right):
  55.     if left is None:
  56.         left = fp[0]
  57.     else:
  58.         left = to_least_float(left)
  59.     if right is None:
  60.         right = fp[-1]
  61.     else:
  62.         right = to_least_float(right)
  63.  
  64.     first, last = xp[0], xp[-1]
  65.  
  66.     for y in x:
  67.         if y < first:
  68.             yield left
  69.         elif y > last:
  70.             yield right
  71.         else:
  72.             i = bisect.bisect(xp, y)
  73.             j = i-1
  74.             k = min(i, n-1)
  75.             yield lerp(fp[j], fp[k], invlerp(y, xp[j], xp[k]))
  76.  
  77. def _interp_periodic(n, x, xp, fp, period):
  78.     if period == 0:
  79.         raise ValueError('period must be a non-zero value')
  80.  
  81.     # Normalize periodic boundaries; reorder xp and fp.
  82.  
  83.     xp = [x % period for x in xp]
  84.     xp_sort_index = sorted(range(n), key=lambda i: xp[i])
  85.     xp = [xp[i] for i in xp_sort_index]
  86.     fp = [fp[i] for i in xp_sort_index]
  87.  
  88.     for y in x:
  89.         i = bisect.bisect(xp, y % period)
  90.         j = (i-1) % n
  91.         k = i % n
  92.         yield lerp(fp[j], fp[k], modinvlerp(y, xp[j], xp[k], period))
  93.  
  94. # Main.
  95.  
  96. import numpy as np
  97.  
  98. x  = [0, 1, 1.5, 2.73, 3, 3.14]
  99. xp = [1, 2, 3]
  100. fp = [3, 2, 0]
  101.  
  102. print(list(np.interp([2.5], xp, fp)))
  103. print(list(interp([2.5], xp, fp)))
  104. print(list(np.interp(x, xp, fp)))
  105. print(list(interp(x, xp, fp)))
  106.  
  107. UNDEF = 99
  108. print(list(np.interp(x, xp, fp, left=-UNDEF, right=UNDEF)))
  109. print(list(interp(x, xp, fp, left=-UNDEF, right=UNDEF)))
  110.  
  111. x  = [-180, -170, -185, 185, -10, -5, 0, 365]
  112. xp = [190, -190, 350, -350]
  113. fp = [5, 10, 3, 4]
  114.  
  115. print(list(np.interp(x, xp, fp, period=360)))
  116. print(list(interp(x, xp, fp, period=360)))
  117.  
  118. x  = [1.5, 4.0]
  119. xp = [2, 3, 5]
  120. fp = [1.0j, 0, 2+3j]
  121.  
  122. print(list(np.interp(x, xp, fp)))
  123. print(list(interp(x, xp, fp)))
Success #stdin #stdout 0.79s 42284KB
comments (?)
stdin 
Standard input is empty
stdout 
[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/9mnbtx
language:
Python 3 (python  3.12)
created:
8 hours 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!