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

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

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

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

受付中

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

labshirasawa_
labshirasawa_

総合スコア4

Python 3.x

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

0回答

0評価

0クリップ

47閲覧

投稿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/ツールのバージョンなど)

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

良い質問の評価を上げる

以下のような質問は評価を上げましょう

  • 質問内容が明確
  • 自分も答えを知りたい
  • 質問者以外のユーザにも役立つ

評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

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

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

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

teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

  • プログラミングに関係のない質問
  • やってほしいことだけを記載した丸投げの質問
  • 問題・課題が含まれていない質問
  • 意図的に内容が抹消された質問
  • 過去に投稿した質問と同じ内容の質問
  • 広告と受け取られるような投稿

評価を下げると、トップページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

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

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

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

ただいまの回答率
87.20%

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

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

質問する

関連した質問

同じタグがついた質問を見る

Python 3.x

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