修改表面代码以求解 4 个维度而不是 3 个 [编辑]
Modify surface code to solve for 4 dimensions instead of 3 [edited]
我通过一些简洁的代码发现了这个很好的问题，通过一些调整，将 3D 多项式曲面拟合到 space.
中的一组点上
Python 3D polynomial surface fit, order dependent
下面是我的版本
最终，我意识到随着时间的推移我需要拟合一个曲面，即我需要求解一个 4 维曲面，我一直在努力解决这个问题。
我想出了一个非常复杂且计算量大的解决方法。我为每个时间间隔创建一个表面。然后我创建一个点网格并找到每个表面上每个点的 Z 值。所以现在我有一堆 x,y 点，每个点都有一个 z 值列表，这些值需要从一个间隔平稳地流动到下一个间隔。所以我们对 z 值进行回归。现在 z 流是平滑的，我根据 x、y 点以及相关时间间隔的平滑 Z 值为每个时间间隔重新拟合表面。
听起来就是这样。笨拙且次优。生成的曲面流动得更顺畅并且仍然表现良好，但必须有一种方法可以去掉中间人并直接在 fitSurface 函数中求解第 4 维。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import itertools
# Start black magic
def xy_powers(order):
powers = itertools.product(range(order + 1), range(order + 1))
return [tup for tup in powers if sum(tup) <= order]
def fitSurface(x, y, z, order):
ncols = (order + 1)**2
G = np.zeros((x.size, ncols))
ij = xy_powers(order)
for k, (i,j) in enumerate(ij):
G[:,k] = x**i * y**j
m, _, _, _ = np.linalg.lstsq(G, z, rcond=None)
return m
def getZValuesForXYInputs(surface, order, x, y):
order = int(np.sqrt(len(surface))) - 1
ij = xy_powers(order)
z = np.zeros_like(x)
for a, (i,j) in zip(surface, ij):
z += a * x**i * y**j
return z
# End black magic
def showRender_3D(x_raw, y_raw, z_raw, xx, yy, zz):
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x_raw, y_raw, z_raw, color='red', zorder=0)
ax.plot_surface(xx, yy, zz, zorder=10, alpha=0.4)
ax.set_xlabel('X data')
ax.set_ylabel('Y data')
ax.set_zlabel('Z data')
plt.show()
def main(order):
# Make generic data
numdata = 100
x = np.random.random(numdata)
y = np.random.random(numdata)
z = x**2 + y**2 + 3*x**3 + y + np.random.random(numdata)
t = np.random.randint(1, 4, numdata) # Break the data into
# Fit the surface
m = fitSurface(x, y, z, order)
# Sample the surface at regular points so we can more easily plot the surface
nx, ny = 40, 40
xx, yy = np.meshgrid(np.linspace(x.min(), x.max(), nx),
np.linspace(y.min(), y.max(), ny))
zz = getZValuesForXYInputs(m, order, xx, yy)
# Plot it
showRender_3D(x, y, z, xx, yy, zz)
orderForSurfaceFit = 3
main(orderForSurfaceFit)
好吧，我想我已经拨通了这个。我不会详细说明如何，只是说一旦你对代码进行了足够的研究，魔法就不会消失，但模式确实会出现。我只是扩展了这些模式，它看起来很有效。
最终结果
无可否认，这是如此低的分辨率，看起来它没有从 C=1 变为 C=2，但确实如此。加载它，你会看到。 gif 现在应该更清楚地显示流程了。
首先是证明背后的方法论。我找到了一个时髦的表面方程并添加了第三个输入变量 C（in-effect 创建 4D 表面），然后使用原始 3D fitter/renderer 作为可信参考点研究具有固定 C 值的表面形状。
当 C 为 1 时，你从地狱得到半管。一个稍微不平衡的向下倾斜的半管。
当 C 为 2 时，您得到的结果大致相同，但不平衡性更加夸张。
当 C 为 3 时，你会得到一个非常不同的形状。就像上面夸张的半管被切成两半，反转，然后粘在一起。
当你 运行 下面的代码时，你会得到一个带有滑块的 3D 渲染，它允许你从 1 到 3 流动 C 值。1、2 和 3 处的值看起来像实心与参考相匹配。我还向数据添加了一个随机化器，以查看它在从不完美的噪声数据近似表面时的表现如何，我也喜欢我在那里看到的结果。
支持以下问题的代码和想法。
Python 3D polynomial surface fit, order dependent
import itertools
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.widgets import Slider
class Surface4D:
def __init__(self, order, a, b, c, z):
# Setting initial attributes
self.order = order
self.a = a
self.b = b
self.c = c
self.z = z
# Setting surface attributes
self.surface = self._fit_surface()
self.aa = None
self.bb = None
self._sample_surface_grid()
# Setting graph attributes
self.surface_render = None
self.axis_3d = None
# Start black magic math
def _abc_powers(self):
powers = itertools.product(range(self.order + 1), range(self.order + 1), range(self.order + 1))
return [tup for tup in powers if sum(tup) <= self.order]
def _fit_surface(self):
ncols = (self.order + 1)**3
G = np.zeros((self.a.size, ncols))
ijk = self._abc_powers()
for idx, (i,j,k) in enumerate(ijk):
G[:,idx] = self.a**i * self.b**j * self.c**k
m, _, _, _ = np.linalg.lstsq(G, self.z, rcond=None)
return m
def get_z_values(self, a, b, c):
ijk = self._abc_powers()
z = np.zeros_like(a)
for s, (i,j,k) in zip(self.surface, ijk):
z += s * a**i * b**j * c**k
return z
# End black magic math
def render_4d_flow(self):
# Set up the layout of the graph
fig = plt.figure()
self.axis_3d = Axes3D(fig, rect=[0.1,0.2,0.8,0.7])
slider_ax = fig.add_axes([0.1,0.1,0.8,0.05])
self.axis_3d.set_xlabel('X data')
self.axis_3d.set_ylabel('Y data')
self.axis_3d.set_zlabel('Z data')
# Plot the point cloud and initial surface
self.axis_3d.scatter(self.a, self.b, self.z, color='red', zorder=0)
zz = self.get_z_values(self.aa, self.bb, 1)
self.surface_render = self.axis_3d.plot_surface(self.aa, self.bb, zz, zorder=10, alpha=0.4, color="b")
# Setup the slider behavior
slider_start_step = self.c.min()
slider_max_steps = self.c.max()
slider = Slider(slider_ax, 'time', slider_start_step, slider_max_steps , valinit=slider_start_step)
slider.on_changed(self._update)
plt.show()
input("Once youre done, hit any enter to continue.")
def _update(self, val):
self.surface_render.remove()
zz = self.get_z_values(self.aa, self.bb, val)
self.surface_render = self.axis_3d.plot_surface(self.aa, self.bb, zz, zorder=10, alpha=0.4, color="b")
def _sample_surface_grid(self):
na, nb = 40, 40
aa, bb = np.meshgrid(np.linspace(self.a.min(), self.a.max(), na),
np.linspace(self.b.min(), self.b.max(), nb))
self.aa = aa
self.bb = bb
def noisify_array(one_dim_array, randomness_multiplier):
listOfNewValues = []
range = abs(one_dim_array.min()-one_dim_array.max())
for item in one_dim_array:
# What percentage are we shifting the point dimension by
shift = np.random.randint(0, 101)
shiftPercent = shift/100
shiftPercent = shiftPercent * randomness_multiplier
# Is that shift positive or negative
shiftSignIndex = np.random.randint(0, 2)
shifSignOption = [-1, 1]
shiftSign = shifSignOption[shiftSignIndex]
# Shift it
newNoisyPosition = item + (range * (shiftPercent * shiftSign))
listOfNewValues.append(newNoisyPosition)
return np.array(listOfNewValues)
def generate_data():
# Define our range for each dimension
x = np.linspace(-6, 6, 20)
y = np.linspace(-6, 6, 20)
w = [1, 2, 3]
# Populate each dimension
a,b,c,z = [],[],[],[]
for X in x:
for Y in y:
for W in w:
a.append(X)
b.append(Y)
c.append(W)
z.append((1 * X ** 4) + (2 * Y ** 3) + (X * Y ** W) + (4 * X) + (5 * Y))
# Convert them to arrays
a = np.array(a)
b = np.array(b)
c = np.array(c)
z = np.array(z)
return [a, b, c, z]
def main(order):
# Make the data
a,b,c,z = generate_data()
# Show the pure data and best fit
surface_pure_data = Surface4D(order, a, b, c, z)
surface_pure_data.render_4d_flow()
# Add some noise to the data
a = noisify_array(a, 0.10)
b = noisify_array(b, 0.10)
c = noisify_array(c, 0.10)
z = noisify_array(z, 0.10)
# Show the noisy data and best fit
surface_noisy_data = Surface4D(order, a, b, c, z)
surface_noisy_data.render_4d_flow()
# ----------------------------------------------------------------
orderForSurfaceFit = 5
main(orderForSurfaceFit)
[编辑 1：我已经开始将此代码合并到我的实际项目中，并且我发现了一些调整以使事情变得合理。还有一个范围问题，运行time 需要在它仍在 render_4d_flow 函数的范围内时暂停，以便滑块工作。]
[编辑 2：显示从 c=2 到 c=3 的流程的高分辨率 gif]
我通过一些简洁的代码发现了这个很好的问题，通过一些调整，将 3D 多项式曲面拟合到 space.
中的一组点上Python 3D polynomial surface fit, order dependent
下面是我的版本
最终，我意识到随着时间的推移我需要拟合一个曲面，即我需要求解一个 4 维曲面，我一直在努力解决这个问题。
我想出了一个非常复杂且计算量大的解决方法。我为每个时间间隔创建一个表面。然后我创建一个点网格并找到每个表面上每个点的 Z 值。所以现在我有一堆 x,y 点，每个点都有一个 z 值列表，这些值需要从一个间隔平稳地流动到下一个间隔。所以我们对 z 值进行回归。现在 z 流是平滑的，我根据 x、y 点以及相关时间间隔的平滑 Z 值为每个时间间隔重新拟合表面。
听起来就是这样。笨拙且次优。生成的曲面流动得更顺畅并且仍然表现良好，但必须有一种方法可以去掉中间人并直接在 fitSurface 函数中求解第 4 维。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import itertools
# Start black magic
def xy_powers(order):
powers = itertools.product(range(order + 1), range(order + 1))
return [tup for tup in powers if sum(tup) <= order]
def fitSurface(x, y, z, order):
ncols = (order + 1)**2
G = np.zeros((x.size, ncols))
ij = xy_powers(order)
for k, (i,j) in enumerate(ij):
G[:,k] = x**i * y**j
m, _, _, _ = np.linalg.lstsq(G, z, rcond=None)
return m
def getZValuesForXYInputs(surface, order, x, y):
order = int(np.sqrt(len(surface))) - 1
ij = xy_powers(order)
z = np.zeros_like(x)
for a, (i,j) in zip(surface, ij):
z += a * x**i * y**j
return z
# End black magic
def showRender_3D(x_raw, y_raw, z_raw, xx, yy, zz):
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x_raw, y_raw, z_raw, color='red', zorder=0)
ax.plot_surface(xx, yy, zz, zorder=10, alpha=0.4)
ax.set_xlabel('X data')
ax.set_ylabel('Y data')
ax.set_zlabel('Z data')
plt.show()
def main(order):
# Make generic data
numdata = 100
x = np.random.random(numdata)
y = np.random.random(numdata)
z = x**2 + y**2 + 3*x**3 + y + np.random.random(numdata)
t = np.random.randint(1, 4, numdata) # Break the data into
# Fit the surface
m = fitSurface(x, y, z, order)
# Sample the surface at regular points so we can more easily plot the surface
nx, ny = 40, 40
xx, yy = np.meshgrid(np.linspace(x.min(), x.max(), nx),
np.linspace(y.min(), y.max(), ny))
zz = getZValuesForXYInputs(m, order, xx, yy)
# Plot it
showRender_3D(x, y, z, xx, yy, zz)
orderForSurfaceFit = 3
main(orderForSurfaceFit)
好吧，我想我已经拨通了这个。我不会详细说明如何，只是说一旦你对代码进行了足够的研究，魔法就不会消失，但模式确实会出现。我只是扩展了这些模式，它看起来很有效。
最终结果
无可否认，这是如此低的分辨率，看起来它没有从 C=1 变为 C=2，但确实如此。加载它，你会看到。 gif 现在应该更清楚地显示流程了。
首先是证明背后的方法论。我找到了一个时髦的表面方程并添加了第三个输入变量 C（in-effect 创建 4D 表面），然后使用原始 3D fitter/renderer 作为可信参考点研究具有固定 C 值的表面形状。
当 C 为 1 时，你从地狱得到半管。一个稍微不平衡的向下倾斜的半管。
当 C 为 2 时，您得到的结果大致相同，但不平衡性更加夸张。
当 C 为 3 时，你会得到一个非常不同的形状。就像上面夸张的半管被切成两半，反转，然后粘在一起。
当你 运行 下面的代码时，你会得到一个带有滑块的 3D 渲染，它允许你从 1 到 3 流动 C 值。1、2 和 3 处的值看起来像实心与参考相匹配。我还向数据添加了一个随机化器，以查看它在从不完美的噪声数据近似表面时的表现如何，我也喜欢我在那里看到的结果。
支持以下问题的代码和想法。
Python 3D polynomial surface fit, order dependent
import itertools
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.widgets import Slider
class Surface4D:
def __init__(self, order, a, b, c, z):
# Setting initial attributes
self.order = order
self.a = a
self.b = b
self.c = c
self.z = z
# Setting surface attributes
self.surface = self._fit_surface()
self.aa = None
self.bb = None
self._sample_surface_grid()
# Setting graph attributes
self.surface_render = None
self.axis_3d = None
# Start black magic math
def _abc_powers(self):
powers = itertools.product(range(self.order + 1), range(self.order + 1), range(self.order + 1))
return [tup for tup in powers if sum(tup) <= self.order]
def _fit_surface(self):
ncols = (self.order + 1)**3
G = np.zeros((self.a.size, ncols))
ijk = self._abc_powers()
for idx, (i,j,k) in enumerate(ijk):
G[:,idx] = self.a**i * self.b**j * self.c**k
m, _, _, _ = np.linalg.lstsq(G, self.z, rcond=None)
return m
def get_z_values(self, a, b, c):
ijk = self._abc_powers()
z = np.zeros_like(a)
for s, (i,j,k) in zip(self.surface, ijk):
z += s * a**i * b**j * c**k
return z
# End black magic math
def render_4d_flow(self):
# Set up the layout of the graph
fig = plt.figure()
self.axis_3d = Axes3D(fig, rect=[0.1,0.2,0.8,0.7])
slider_ax = fig.add_axes([0.1,0.1,0.8,0.05])
self.axis_3d.set_xlabel('X data')
self.axis_3d.set_ylabel('Y data')
self.axis_3d.set_zlabel('Z data')
# Plot the point cloud and initial surface
self.axis_3d.scatter(self.a, self.b, self.z, color='red', zorder=0)
zz = self.get_z_values(self.aa, self.bb, 1)
self.surface_render = self.axis_3d.plot_surface(self.aa, self.bb, zz, zorder=10, alpha=0.4, color="b")
# Setup the slider behavior
slider_start_step = self.c.min()
slider_max_steps = self.c.max()
slider = Slider(slider_ax, 'time', slider_start_step, slider_max_steps , valinit=slider_start_step)
slider.on_changed(self._update)
plt.show()
input("Once youre done, hit any enter to continue.")
def _update(self, val):
self.surface_render.remove()
zz = self.get_z_values(self.aa, self.bb, val)
self.surface_render = self.axis_3d.plot_surface(self.aa, self.bb, zz, zorder=10, alpha=0.4, color="b")
def _sample_surface_grid(self):
na, nb = 40, 40
aa, bb = np.meshgrid(np.linspace(self.a.min(), self.a.max(), na),
np.linspace(self.b.min(), self.b.max(), nb))
self.aa = aa
self.bb = bb
def noisify_array(one_dim_array, randomness_multiplier):
listOfNewValues = []
range = abs(one_dim_array.min()-one_dim_array.max())
for item in one_dim_array:
# What percentage are we shifting the point dimension by
shift = np.random.randint(0, 101)
shiftPercent = shift/100
shiftPercent = shiftPercent * randomness_multiplier
# Is that shift positive or negative
shiftSignIndex = np.random.randint(0, 2)
shifSignOption = [-1, 1]
shiftSign = shifSignOption[shiftSignIndex]
# Shift it
newNoisyPosition = item + (range * (shiftPercent * shiftSign))
listOfNewValues.append(newNoisyPosition)
return np.array(listOfNewValues)
def generate_data():
# Define our range for each dimension
x = np.linspace(-6, 6, 20)
y = np.linspace(-6, 6, 20)
w = [1, 2, 3]
# Populate each dimension
a,b,c,z = [],[],[],[]
for X in x:
for Y in y:
for W in w:
a.append(X)
b.append(Y)
c.append(W)
z.append((1 * X ** 4) + (2 * Y ** 3) + (X * Y ** W) + (4 * X) + (5 * Y))
# Convert them to arrays
a = np.array(a)
b = np.array(b)
c = np.array(c)
z = np.array(z)
return [a, b, c, z]
def main(order):
# Make the data
a,b,c,z = generate_data()
# Show the pure data and best fit
surface_pure_data = Surface4D(order, a, b, c, z)
surface_pure_data.render_4d_flow()
# Add some noise to the data
a = noisify_array(a, 0.10)
b = noisify_array(b, 0.10)
c = noisify_array(c, 0.10)
z = noisify_array(z, 0.10)
# Show the noisy data and best fit
surface_noisy_data = Surface4D(order, a, b, c, z)
surface_noisy_data.render_4d_flow()
# ----------------------------------------------------------------
orderForSurfaceFit = 5
main(orderForSurfaceFit)
[编辑 1：我已经开始将此代码合并到我的实际项目中，并且我发现了一些调整以使事情变得合理。还有一个范围问题，运行time 需要在它仍在 render_4d_flow 函数的范围内时暂停，以便滑块工作。]
[编辑 2：显示从 c=2 到 c=3 的流程的高分辨率 gif]