Проблема в записи вывода дифференциального уравнения с многопроцессорной обработкой в ​​python - PullRequest
Новая
       21

Проблема в записи вывода дифференциального уравнения с многопроцессорной обработкой в ​​python

2 голосов
/ 14 января 2020

Я совершенно новичок в параллельных вычислениях, фактически, в численных методах. Я пытаюсь решить дифференциальное уравнение, используя python solve_ivp следующей формы: 

y''(x) + (a^2 + x^2)y(x) = 0
y(0)=1
y'(0)=0
x=(0,100)

Я хочу решить для диапазона a и записать файл как a[i] y[i](80).

Исходное уравнение довольно сложное, но по существу структура такая же, как определено выше. Я использовал for l oop, и это занимает много времени для вычислений. Ища в Интернете, я наткнулся на этот прекрасный веб-сайт и нашел этот вопрос и связанный с ним ответ, который очень может решить проблему, с которой я столкнулся.

Я попробовал исходный код, предоставленный в решении; однако полученный вывод неправильно отсортирован. Я имею в виду, что второй столбец не в правильном порядке. 

q1    a1    Y1
q1    a2    Y3
q1    a4    Y4
q1    a3    Y3
q1    a5    Y5
...

Я даже пытался с одним l oop с одним параметром, но та же проблема все еще остается. Ниже приведен мой код с тем же методом многопроцессорной обработки, но с solve_ivp 

import numpy as np
import scipy.integrate
import multiprocessing as mp
from scipy.integrate import solve_ivp


def fun(t, y):
    # replace this function with whatever function you want to work with
    #    (this one is the example function from the scipy docs for odeint)
    theta, omega = y
    dydt = [omega, -a*omega - q*np.sin(theta)]
    return dydt

#definitions of work thread and write thread functions
tspan = np.linspace(0, 10, 201)


def run_thread(input_queue, output_queue):
    # run threads will pull tasks from the input_queue, push results into output_queue
    while True:
        try:
            queueitem = input_queue.get(block = False)
            if len(queueitem) == 3:
                a, q, t = queueitem
                sol = solve_ivp(fun, [tspan[0], tspan[-1]], [1, 0], method='RK45', t_eval=tspan)
                F = 1 + sol.y[0].T[157]
                output_queue.put((q, a, F))
        except Exception as e:
            print(str(e))
            print("Queue exhausted, terminating")
            break

def write_thread(queue):    
    # write thread will pull results from output_queue, write them to outputfile.txt
    f1 = open("outputfile.txt", "w")
    while True:
        try:
            queueitem = queue.get(block = False)
            if queueitem[0] == "TERMINATE":
                f1.close()
                break
            else:
                q, a, F = queueitem                
                print("{}  {} {} \n".format(q, a, F))            
                f1.write("{}  {} {} \n".format(q, a, F))            
        except:
            # necessary since it will throw an error whenever output_queue is empty
            pass

# define time point sequence            
t = np.linspace(0, 10, 201)

# prepare input and output Queues
mpM = mp.Manager()
input_queue = mpM.Queue()
output_queue = mpM.Queue()

# prepare tasks, collect them in input_queue
for q in np.linspace(0.0, 4.0, 100):
    for a in np.linspace(-2.0, 7.0, 100):
        # Your computations as commented here will now happen in run_threads as defined above and created below
        # print('Solving for q = {}, a = {}'.format(q,a))
        # sol1 = scipy.integrate.odeint(fun, [1, 0], t, args=( a, q))[..., 0]
        # print(t[157])
        # F = 1 + sol1[157]    
        input_tupel = (a, q, t)
        input_queue.put(input_tupel)

# create threads
thread_number = mp.cpu_count()
procs_list = [mp.Process(target = run_thread , args = (input_queue, output_queue)) for i in range(thread_number)]         
write_proc = mp.Process(target = write_thread, args = (output_queue,))

# start threads
for proc in procs_list:
    proc.start()
write_proc.start()

# wait for run_threads to finish
for proc in procs_list:
    proc.join()

# terminate write_thread
output_queue.put(("TERMINATE",))
write_proc.join()

Пожалуйста, дайте мне знать, что не так в многопроцессорной обработке, чтобы я мог немного узнать о многопроцессорной обработке в python в процессе. Кроме того, я был бы очень признателен, если бы кто-нибудь дал мне знать о самых элегантных / эффективных способах обработки таких вычислений в python. Спасибо

1 Ответ

1 голос
/ 15 января 2020

То, что вы хотите, является чем-то вроде онлайн сортировки . В этом случае вам известен порядок , в котором вы хотите представить результаты ( т.е. , порядок ввода), поэтому просто накапливайте выходные данные в очереди приоритетов и извлекайте из него элементы, когда они соответствуют следующему ключу, который вы ожидаете.

Тривиальный пример без демонстрации параллелизма: 

import heapq
def sort_follow(pairs,ref):
  """
  Sort (a prefix of) pairs so as to have its first components
  be the elements of (sorted) ref in order.
  Uses memory proportional to the disorder in pairs.
  """
  heap=[]
  pairs=iter(pairs)
  for r in ref:
    while not heap or heap[0][0]!=r:
      heapq.heappush(heap,next(pairs))
    yield heapq.heappop(heap)

Преимущество этого подхода - уменьшенный объем памяти - вероятно не имеет значения для небольших результатов, полученных всего за несколько операций с плавающей запятой, но его достаточно легко применить.

...