Итериране на масив в 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 стойностите и решаване за всеки случай. Написах итерацията като разбиране на списък. Възможни са други формати на итерация (за цикъл и т.н.). Може дори да успеем да обвием това brentq извикване във функция, която може да бъде предадена през np.vectorize. Но това все още ще бъде повторение с малки спестявания на време.

Има различни начини за работа с многоизмерен масив. Един прост е да flatten входа, да направите 1d итерация и след това да промените формата на резултата. Например:

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