Совершенствование простой реализации простой весенней сети

Question

Мне нужна простая пружинная система, написанная в numpy. Система будет определена как простая сеть knots, связанная links. Меня не интересует оценка системы с течением времени, но вместо этого я хочу перейти из начального состояния, изменить переменную (как правило, переместить knot в новую позицию) и решить систему до тех пор, пока она не достигнет состояния (последний приложенная сила ниже заданного порога). Узлы не имеют массы, нет силы тяжести, все силы выведены из текущей длины/длины длины каждой линии. И единственная «специальная» переменная состоит в том, что каждый узел может поставить ставку как «привязанный» (не перемещается).

Итак, я написал этот простой решатель ниже и включил простой пример. Перейти к самому концу моего вопроса.

import numpy as np 
from numpy.core.umath_tests import inner1d 

np.set_printoptions(precision=4) 
np.set_printoptions(suppress=True) 
np.set_printoptions(linewidth =150) 
np.set_printoptions(threshold=10) 


def solver(kPos, kAnchor, link0, link1, w0, cycles=1000, precision=0.001, dampening=0.1, debug=False): 
    """ 
    kPos  : vector array - knot position 
    kAnchor : float array - knot's anchor state, 0 = moves freely, 1 = anchored (not moving) 
    link0  : int array - array of links connecting each knot. each index corresponds to a knot 
    link1  : int array - array of links connecting each knot. each index corresponds to a knot 
    w0   : float array - initial link length 
    cycles  : int   - eval stops when n cycles reached 
    precision : float  - eval stops when highest applied force is below this value 
    dampening : float  - keeps system stable during each iteration 
    """ 

    kPos  = np.asarray(kPos) 
    pos   = np.array(kPos) # copy of kPos 
    kAnchor  = 1-np.clip(np.asarray(kAnchor).astype(float),0,1)[:,None] 
    link0  = np.asarray(link0).astype(int) 
    link1  = np.asarray(link1).astype(int) 
    w0   = np.asarray(w0).astype(float) 

    F = np.zeros(pos.shape) 
    i = 0 

    for i in xrange(cycles): 

     # Init force applied per knot 
     F = np.zeros(pos.shape) 

     # Calculate forces 
     AB = pos[link1] - pos[link0] # get link vectors between knots 
     w1 = np.sqrt(inner1d(AB,AB)) # get link lengths 
     AB/=w1[:,None] # normalize link vectors 
     f = (w1 - w0) # calculate force vectors 
     f = f[:,None] * AB 

     # Apply force vectors on each knot 
     np.add.at(F, link0, f) 
     np.subtract.at(F, link1, f) 

     # Update point positions  
     pos += F * dampening * kAnchor 

     # If the maximum force applied is below our precision criteria, exit 
     if np.amax(F) < precision: 
      break 

    # Debug info 
    if debug: 
     print 'Iterations: %s'%i 
     print 'Max Force: %s'%np.amax(F) 

    return pos 

Вот несколько тестовых данных, чтобы показать, как это работает. В этом случае я использую сетку, но в действительности это может быть любой тип сети, как строка с большим количеством узлов, или беспорядок полигонов ...:

import cProfile 

# Create a 5x5 3D knot grid 
z = np.linspace(-0.5, 0.5, 5) 
x = np.linspace(-0.5, 0.5, 5)[::-1] 
x,z = np.meshgrid(x,z) 
kPos = np.array([np.array(thing) for thing in zip(x.flatten(), z.flatten())]) 
kPos = np.insert(kPos, 1, 0, axis=1) 
''' 
array([[-0.5 , 0. , 0.5 ], 
     [-0.25, 0. , 0.5 ], 
     [ 0. , 0. , 0.5 ], 
     ..., 
     [ 0. , 0. , -0.5 ], 
     [ 0.25, 0. , -0.5 ], 
     [ 0.5 , 0. , -0.5 ]]) 
''' 


# Define the links connecting each knots 
link0 = [0,1,2,3,5,6,7,8,10,11,12,13,15,16,17,18,20,21,22,23,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] 
link1 = [1,2,3,4,6,7,8,9,11,12,13,14,16,17,18,19,21,22,23,24,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24] 
AB = kPos[link0]-kPos[link1] 
w0 = np.sqrt(inner1d(AB,AB)) # this is a square grid, each link's initial length will be 0.25 

# Set the anchor states 
kAnchor = np.zeros(len(kPos)) # All knots will be free floating 
kAnchor[12] = 1 # Middle knot will be anchored 

Это то, что сетка выглядит как:

Если запустить свой код, используя эти данные, ничего не произойдет, так как ссылки не толкая или растяжение:

print np.allclose(kPos,solver(kPos, kAnchor, link0, link1, w0, debug=True)) 
# Returns True 
# Iterations: 0 
# Max Force: 0.0 

Теперь давайте двигаться, что средний привязанный узел вверх немного и решить систему:

# Move the center knot up a little 
kPos[12] = np.array([0,0.3,0]) 

# eval the system 
new = solver(kPos, kAnchor, link0, link1, w0, debug=True) # positions will have moved 
#Iterations: 102 
#Max Force: 0.000976603249133 

# Rerun with cProfile to see how fast it runs 
cProfile.run('solver(kPos, kAnchor, link0, link1, w0)') 
# 520 function calls in 0.008 seconds 

А вот что сетка выглядит как после того, как тянули этот единственный привязанного узел:

Вопрос:

Мои фактические случаи использования немного сложнее, чем этот пример, и решить слишком медленно для моего вкуса: (100-200 узлов с сетью где-то между 200-300 ссылками, решает через несколько секунд).

Как я могу запустить функцию моего решателя быстрее? Я бы рассмотрел Cython, но у меня есть нулевой опыт работы с C. Любая помощь будет принята с благодарностью.

Answer 1

Ваш метод при беглом взгляде представляется явным методом релаксации. Рассчитайте остаточную силу на каждом узле, примените коэффициент этой силы как смещение, повторите до сближения. Это повторение до тех пор, пока конвергенция не займет времени. Чем больше у вас очков, тем дольше будет каждая итерация, но вам также потребуется больше итераций для ограничений на одном конце сетки для распространения на другой.

Вы считали неявный метод? Напишите уравнение для остаточной силы на каждом не ограниченном узле, соберите их в большую матрицу и решите на одном шаге. Информация теперь распространяется по всей проблеме за один шаг. В качестве дополнительного преимущества матрица, которую вы создаете, должна быть редкой, для которой scipy имеет модуль.

Wikipedia: explicit and implicit methods

---

РЕДАКТИРОВАТЬ Пример неявного согласования метода (примерно) вашей проблемы. Это решение линейно, поэтому оно не учитывает влияние рассчитанного смещения на силу. Чтобы вычислить это, вам понадобится повторить (или использовать нелинейные методы). Надеюсь, поможет.

#!/usr/bin/python3 

import matplotlib.pyplot as pp 
from mpl_toolkits.mplot3d import Axes3D 
import numpy as np 
import scipy as sp 
import scipy.sparse 
import scipy.sparse.linalg 

#------------------------------------------------------------------------------# 

# Generate a grid of knots 
nX = 10 
nY = 10 
x = np.linspace(-0.5, 0.5, nX) 
y = np.linspace(-0.5, 0.5, nY) 
x, y = np.meshgrid(x, y) 
knots = list(zip(x.flatten(), y.flatten())) 

# Create links between the knots 
links = [] 
# Horizontal links 
for i in range(0, nY): 
    for j in range(0, nX - 1): 
     links.append((i*nX + j, i*nX + j + 1)) 
# Vertical links 
for i in range(0, nY - 1): 
    for j in range(0, nX): 
     links.append((i*nX + j, (i + 1)*nX + j)) 

# Create constraints. This dict takes a knot index as a key and returns the 
# fixed z-displacement associated with that knot. 
constraints = { 
    0   : 0.0, 
    nX - 1  : 0.0, 
    nX*(nY - 1): 0.0, 
    nX*nY - 1 : 1.0, 
    2*nX + 4 : 1.0, 
    } 

#------------------------------------------------------------------------------# 

# Matrix i-coordinate, j-coordinate and value 
Ai = [] 
Aj = [] 
Ax = [] 

# Right hand side array 
B = np.zeros(len(knots)) 

# Loop over the links 
for link in links: 

    # Link geometry 
    displacement = np.array([ knots[1][i] - knots[0][i] for i in range(2) ]) 
    distance = np.sqrt(displacement.dot(displacement)) 

    # For each node 
    for i in range(2): 

     # If it is not a constraint, add the force associated with the link to 
     # the equation of the knot 
     if link[i] not in constraints: 

      Ai.append(link[i]) 
      Aj.append(link[i]) 
      Ax.append(-1/distance) 

      Ai.append(link[i]) 
      Aj.append(link[not i]) 
      Ax.append(+1/distance) 

     # If it is a constraint add a diagonal and a value 
     else: 

      Ai.append(link[i]) 
      Aj.append(link[i]) 
      Ax.append(+1.0) 
      B[link[i]] += constraints[link[i]] 

# Create the matrix and solve 
A = sp.sparse.coo_matrix((Ax, (Ai, Aj))).tocsr() 
X = sp.sparse.linalg.lsqr(A, B)[0] 

#------------------------------------------------------------------------------# 

# Plot the links 
fg = pp.figure() 
ax = fg.add_subplot(111, projection='3d') 

for link in links: 
    x = [ knots[i][0] for i in link ] 
    y = [ knots[i][1] for i in link ] 
    z = [ X[i] for i in link ] 
    ax.plot(x, y, z) 

pp.show()