Итерация массива в Python с использованием функции brentq

У меня возникают проблемы с повторением каждого элемента массива с помощью функции brentq. q в определенной ниже функции представляет собой файловый массив FITS, и мы используем каждый элемент в этом массиве в качестве входных данных для запуска функции brentq, чтобы найти T.

По сути, моя проблема заключается в том, что я не знаю, где и как реализовать соответствующий цикл for для итерации функции по каждому элементу q.

Любые предложения о том, как решить эту проблему?

def f(T,q,coeff1,coeff2,coeff3):
    return q*const3 - ((exp(const2/T)-1)/(exp(const/T)-1))

a = brentq(f, 10, 435.1, args=(q,4351.041,4262.570,0.206))
print a

newhdu = fits.PrimaryHDU(a)
newhdulist = fits.HDUList([newhdu])
newhdulist.writeto('Temp21DCOT.fits')

Дополнительное объяснение: в основе того, что я пытаюсь сделать, лежит использование brentq для определения значений температуры с использованием значений интенсивности нашего исходного массива (нашего файла FITS).

Уравнение получено из отношения двух длин волн уравнения Планка, поэтому q = B_1/B_2, если мы хотим быть верными физике, где каждый элемент в q является значением интенсивности. brentq решит это аналитически неразрешимое уравнение для T (температуры) для каждого элемента в q и создаст новый массив температур того же размера, что и q. Другими словами, я пытаюсь определить температуру каждого пикселя в файле FITS, используя уравнение Планка.

Примечание. Я повторно опубликовал это, чтобы более эффективно прояснить проблему.


python numpy scipy arrays mathematical-optimization
person Wolfgang    schedule 12.07.2015    source источник
comment
Имена аргументов функции не совпадают с именами переменных в выражении возврата функции f. Вы используете scipy.optimize.brentq?   -  person wwii    schedule 12.07.2015
comment
Да, точно. Я импортировал scipy.optimize как brentq.   -  person Wolfgang    schedule 12.07.2015
comment
Ваша функция принимает 5 аргументов. Вы передаете 4 значения в аргументе brentq args - похоже, f.q в вашем примере равно 1. вы спрашиваете, как выполнить итерацию по массиву и заменить значение в args(1, 4351.041 ... примера значениями из множество?   -  person wwii    schedule 12.07.2015
comment
Да, поэтому мы только итерируем двумерный массив в первый аргумент в arg(q,...), остальные являются константами.   -  person Wolfgang    schedule 13.07.2015

Ответы (1)


arrow_upward
1
arrow_downward

У вас проблемы с итерацией или с эффективностью?

Эта итерация работает для меня:

In [485]: from scipy import optimize

In [486]: def f(T,q,coeff1,coeff2,coeff3):
        return q*coeff3 - ((np.exp(coeff2/T)-1)/(np.exp(coeff1/T)-1))
        # corrected the coeff use

In [487]: q=np.linspace(1,3,10)
# q range chosen to avoid the different signs ValueError

In [488]: A=[optimize.brentq(f, 10, 435.1, args=(i,4351.041,4262.570,0.206),full_output=True) for i in q]

In [489]: A
Out[489]: 
[(55.99858839149886, <scipy.optimize.zeros.RootResults at 0xa861284c>),
 (64.14621536172528, <scipy.optimize.zeros.RootResults at 0xa861286c>),
 (72.98658083834341, <scipy.optimize.zeros.RootResults at 0xa861288c>),
 (82.75638321495505, <scipy.optimize.zeros.RootResults at 0xa86128ac>),
 (93.73016750496367, <scipy.optimize.zeros.RootResults at 0xa86128cc>),
 (106.25045004489637, <scipy.optimize.zeros.RootResults at 0xa86128ec>),
 (120.76612665626851, <scipy.optimize.zeros.RootResults at 0xa861290c>),
 (137.88917389176325, <scipy.optimize.zeros.RootResults at 0xa861292c>),
 (158.4854607193551, <scipy.optimize.zeros.RootResults at 0xa861294c>),
 (183.82941862839408, <scipy.optimize.zeros.RootResults at 0xa861296c>)]

In [490]: [a[1].iterations for a in A]
Out[490]: [8, 9, 10, 10, 10, 10, 10, 9, 8, 10]

В документах brentq f возвращает одно значение для одного набора args. Есть некоторые решатели, такие как ode, которые позволяют определить функцию, которая принимает векторную переменную и возвращает соответствующую векторную производную. Не похоже, что этот искатель корней позволяет это. Таким образом, вы застряли в переборе значений args и решении для каждого случая. Я написал итерацию как понимание списка. Возможны другие форматы итераций (цикл for и т.д.). Мы могли бы даже обернуть этот вызов brentq функцией, которую можно было бы передать через np.vectorize. Но это все еще будет итерация с небольшой экономией времени.

Существуют различные способы обработки многомерного массива. Один простой — flatten ввести, выполнить 1-ю итерацию, а затем изменить форму результата. Например:

In [517]: q1=q.reshape(2,5)

In [518]: q1
Out[518]: 
array([[ 1.        ,  1.22222222,  1.44444444,  1.66666667,  1.88888889],
       [ 2.11111111,  2.33333333,  2.55555556,  2.77777778,  3.        ]])

In [519]: np.array([optimize.brentq(f, 10, 435.1, args=(i,4351.041,4262.570,0.206)) for i in q1.flat]).reshape(q1.shape)
Out[519]: 
array([[  55.99858839,   64.14621536,   72.98658084,   82.75638321,
          93.7301675 ],
       [ 106.25045004,  120.76612666,  137.88917389,  158.48546072,
         183.82941863]])

Я отказался от флага full_output, так как это добавляет сложности.

person hpaulj    schedule 12.07.2015
comment
Будет ли это та же самая процедура для двумерного массива? Потому что это то, что q, к сожалению. - person Wolfgang; 13.07.2015