質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Q&A

0回答

480閲覧

モデルフィッティングに関して

labshirasawa_

総合スコア4

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

0グッド

0クリップ

投稿2020/01/17 04:30

前提・実現したいこと

ここに質問の内容を詳しく書いてください。
初めてpythonを使用して、で計測したデータに対してあるモデルでフィッティングするプログラムを作っています。
様々なカーブフィッティングのサイトをオマージュしているのですが,配列とスカラーの相違でエラーが発生中です。
説明不足な点がありますが、pythonに詳しい方手助けお願いします。

発生している問題・エラーメッセージ

TypeError Traceback (most recent call last) <ipython-input-10-341c810d3f60> in <module> 162 b1=([0,0,1.1],[np.inf,370,2.2]) 163 --> 164 para_opt, cov = scipy.optimize.curve_fit(microfacet_EX, inc, me, para_ini) 165 166 #得られたフィッティングパラメータの表示 ~\Anaconda3\envs\research\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs) 750 # Remove full_output from kwargs, otherwise we're passing it in twice. 751 return_full = kwargs.pop('full_output', False) --> 752 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs) 753 popt, pcov, infodict, errmsg, ier = res 754 cost = np.sum(infodict['fvec'] ** 2) ~\Anaconda3\envs\research\lib\site-packages\scipy\optimize\minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag) 381 if not isinstance(args, tuple): 382 args = (args,) --> 383 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n) 384 m = shape[0] 385 ~\Anaconda3\envs\research\lib\site-packages\scipy\optimize\minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape) 24 def _check_func(checker, argname, thefunc, x0, args, numinputs, 25 output_shape=None): ---> 26 res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) 27 if (output_shape is not None) and (shape(res) != output_shape): 28 if (output_shape[0] != 1): ~\Anaconda3\envs\research\lib\site-packages\scipy\optimize\minpack.py in func_wrapped(params) 456 if transform is None: 457 def func_wrapped(params): --> 458 return func(xdata, *params) - ydata 459 elif transform.ndim == 1: 460 def func_wrapped(params): <ipython-input-10-341c810d3f60> in microfacet_EX(inangle, a, d, n2) 151 theta_h = theta_i + theta_r / LA.norm(theta_i+theta_r,1) 152 --> 153 model = D(a, theta_h) * G(theta_i, theta_r, theta_h) / (4 * math.cos(math.radians(theta_i) * math.cos(math.radians(theta_r)))) 154 I = IridescenceTerm(d, np.dot(theta_h,theta_i), n1, n2, eta3, k, wavelength) 155 <ipython-input-10-341c810d3f60> in D(expo, Theta_h) 19 #法線分布項 20 def D(expo,Theta_h): ---> 21 return (expo+2) * (math.cos(math.radians(Theta_h))**expo) / (2*math.pi) 22 23 #幾何学減衰項 TypeError: only size-1 arrays can be converted to Python scalars

該当のソースコード

import numpy as np from numpy import linalg as LA import matplotlib.pyplot as plt #import seaborn as sns import scipy.optimize import math data = np.loadtxt("fitting_data.csv", delimiter = ",")#計測データの読み取り incident = -data[0]#計測データから入射角を取り出す inc = incident[:5]#初め6列 print(inc) me = data[3,:5]; print(me) plt.plot(inc, me, 'o', label='Raw data') def microfacet_EX(inangle, a, d, n2): #法線分布項 def D(expo,Theta_h): return (expo+2) * (math.cos(math.radians(Theta_h))**expo) / (2*math.pi) #幾何学減衰項 def G(theta_I, theta_R, theta_H): if theta_R > theta_I: theta_rh = theta_R-theta_H else: theta_rh = theta_R+theta_H num_tmp = np.array([0]*3) num_tmp[0] = 1 num_tmp[1] = 2 * math.cos(math.radians(theta_H)) * math.cos(math.radians(theta_R)) / math.cos(math.radians(theta_rh)) num_tmp[2] = 2 * math.cos(math.radians(theta_H)) * math.cos(math.radians(theta_I)) / math.cos(math.radians(theta_rh)) return min(num_tmp) #構造色に関する部分 def IridescenceTerm(h, ct1, eta1, eta2, eta3, kappa3, wavelength): def fresnelPhaseExact(cost, eta_1, eta_2, kappa2): sinThetaSqr = 1.0 - cost**2 A = (eta_2**2)*(1.0-kappa2**2) - (eta1**2)*sinThetaSqr B = math.sqrt((A**2) + (2*(eta2**2)*kappa2)**2) U = math.sqrt((A+B)/2.0) V = math.sqrt((B-A)/2.0) phiS = math.atan2(2*eta1*V*cost, (U**2)+(V**2)-(eta1*cost)**2) phiP = math.atan2(2*eta1*(eta2**2)*cost * (2*kappa2*U - (1.0-(kappa2**2) * V)), ((eta2**2)*(1.0+(kappa2**2))*cost)**2 - (eta1**2)*((U**2)+(V**2))) return phiP, phiS def fresnelConductorExact(cosThetaI, eta, k): cosThetaI2 = cosThetaI*cosThetaI sinThetaI2 = 1-cosThetaI2 sinThetaI4 = sinThetaI2*sinThetaI2 temp1 = eta*eta - k*k - sinThetaI2 a2pb = temp1*temp1 + 4*k*k*eta*eta if a2pb>0.0: a2 = a2pb else: a2 = 0.0 a2pb2 = math.sqrt(a2) A = 0.5 * (a2pb2+temp1) if A>0.0: AA = A else: AA = 0.0 a = math.sqrt(AA) term1 = a2pb2 + cosThetaI2 term2 = 2*a*cosThetaI Rs2 = (term1 - term2) / (term1 + term2) term3 = a2pb2*cosThetaI2 + sinThetaI4 term4 = term2*sinThetaI2 Rp2 = Rs2 * (term3 - term4) / (term3 + term4) return Rp2, Rs2 R12p=T121p=R23p=R12s=T121s=R23s=ct2 = 0.0 scale = eta1 / eta2 cosThetaTSqr = 1 - (1-ct1**2) * (scale**2) if cosThetaTSqr <= 0.0: R12s = 1.0 R12p = 1.0 T121p = 0.0 T121s = 0.0 else: ct2 = math.sqrt(cosThetaTSqr) result_film = fresnelConductorExact(ct1, eta2/eta1, 0.0) R12p = result_film[0] R12s = result_film[1] #ベースの反射項 result_base = fresnelConductorExact(ct2, eta3/eta2, kappa3/eta2) R23p = result_base[0] R23s = result_base[1] #透過項 T121p = 1.0 - R12p T121s = 1.0 - R12s #光路差 OPD = 2.0 * eta2 * h * ct2 OPDphi = 2.0 * math.pi * OPD / wavelength #位相シフト shif_1to2 = fresnelPhaseExact(ct1, 1.0, eta2, 0.0,) shif_2to3 = fresnelPhaseExact(ct2, eta2, eta3, kappa3,) phi21p = math.pi - shif_1to2[0] phi21s = math.pi - shif_1to2[1] phi23p = shif_2to3[0] phi23s = shif_2to3[1] r123p = math.sqrt(R12p*R23p) r123s = math.sqrt(R12s*R23s) Rs = ((T121p**2)*R23p) / (1.0 - R12p*R23p) cosP = math.cos(OPDphi + phi23p + phi21p) irid = (r123p * cosP - r123p**2) / (1.0 - 2.0*r123p*cosP + r123p**2) I = R12p + Rs + 2.0*(Rs - T121p) * irid Rs = ((T121s**2)*R23s) / (1.0 - R12s*R23s) cosP = math.cos(OPDphi + phi23s + phi21s) irid = (r123s * cosP - r123s**2) / (1.0 - 2.0*r123s*cosP + r123s**2) I += R12s + Rs + 2.0*(Rs - T121s) * irid return I data = np.loadtxt("fitting_data.csv", delimiter = ",")#計測データの読み取り view = data[1,1] Eta3 = np.loadtxt("Fe_eta.csv", delimiter = ",", usecols=1) eta3 = Eta3[1] K = np.loadtxt("Fe_k.csv", delimiter = ",", usecols=1) k = K[1] wavelength = 390 n1 = 1.000277;#空気の屈折率(波長依存性はほぼないとみる) theta_i = inangle theta_r = view theta_h = theta_i + theta_r / LA.norm(theta_i+theta_r,1) model = D(a, theta_h) * G(theta_i, theta_r, theta_h) / (4 * math.cos(math.radians(theta_i) * math.cos(math.radians(theta_r)))) I = IridescenceTerm(d, np.dot(theta_h,theta_i), n1, n2, eta3, k, wavelength) PI = I * model return PI # initial guess and boundary para_ini = (0.001, 100, 1.2) b1=([0,0,1.1],[np.inf,370,2.2]) para_opt, cov = scipy.optimize.curve_fit(microfacet_EX, inc, me, para_ini) #得られたフィッティングパラメータの表示 print ("Fitted a0 =", str(para_opt[0])) print ("Fitted d0 =", str(para_opt[1])) print ("Fitted n_20' =", str(para_opt[2]))

試したこと

関数microfacet_EXのtheta_iの前にfor文を書きtheta_iを一つずつ計算した方がいいのか検討中

補足情報(FW/ツールのバージョンなど)

ここにより詳細な情報を記載してください。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問