python のoptimize 使用時のエラーの原因と解決策が知りたい

前提・実現したいこと

python を使ってデータをフィッティングしているときに生じました。

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

Improper input: N=20 must not exceed M=1

該当のソースコード

python
1from __future__ import print_function
2if hasattr(__builtins__, 'raw_input'):
3    input = raw_input
4from scipy import *
5from scipy.interpolate import interp1d
6import matplotlib.pyplot as plt
7import csv
8import numpy as np
9import random
10from matplotlib import pyplot as plt
11import scipy.optimize as optimize
12%matplotlib inline 
13from sklearn.metrics import r2_score
14from sklearn.metrics import mean_squared_error
15from sklearn.metrics import mean_absolute_error
16from scipy import integrate
17import PyCrystalField as Py
18import scipy.linalg as LA
19import random
20
21
22
23data = np.genfromtxt("Tb_1p51p5_0p5.txt", unpack=True)
24E_data = data[0]
25I_data = data[1]
26er = data[2]
27## len(E_data) = len(I_data) = len(er) = 61
28
29
30def spectra(E, B20, B22, B2m2, B40, B42, B4m2, B44, B4m4, B60, B62, B6m2,
31            B64, B6m4, B66, B6m6, lam, A, a, b, ei):
32    Bs = {'B20':B20/f[0]*alpha, 'B22':B22/f[1]*alpha, 'B2-2':B2m2/f[1]*alpha, 
33          'B40':B40/f[2]*beta, 'B42':B42/f[3]*beta, 'B4-2':B4m2/f[3]*beta, 'B44':B44/f[4]*beta, 
34          'B4-4':B4m4/f[4]*beta, 'B60':B60/f[5]*gamma, 'B62':B62/f[6]*gamma, 'B6-2':B6m2/f[6]*gamma, 
35          'B64':B64/f[7]*gamma, 'B6-4':B6m4/f[7]*gamma, 'B66':B66/f[8]*gamma, 'B6-6':B6m6/f[8]*gamma}
36  
37    Original = Py.CFLevels.Bdict(ion, Bs)
38    jx = Original.opttran.Jy
39    jy = Original.opttran.Jy*1j
40    jz = Original.opttran.Jz
41    CFLevels = Py.CFLevels.Bdict(ion, Bs)
42    h = lam*Mfe
43    H_Zeeman = gJ*mub*h*(n[0]*jx + n[1]*jy + n[2]*jz)
44    Ei = {}
45    Ei['Ei'] = ei
46    resolution = fix_consts(Resfun, Ei)
47    spec = np.zeros((1,len(E)))
48    for site in sites:
49        R = np.dot(LA.expm(1j*site[0]*jz), np.dot(LA.expm(1j*site[1]*jy), LA.expm(1j*site[2]*jz)))
50        CFLevels.H = np.dot(R, np.dot(CFLevels.H, np.conjugate(R.T)))
51        CFLevels.H += H_Zeeman
52        spec += CFLevels.neutronSpectrum2(E, 5, ei, resolution, 'Ei')
53    spec = A*spec/len(sites)
54    spec += a*E + b
55    return spec
56
57constrain = {'B20':Bpara[0], 'B22':Bpara[1], 'B2m2':Bpara[2], 'B40':Bpara[3], 
58             'B42':Bpara[4], 'B4m2':Bpara[5], 'B44':Bpara[6], 'B4m4':Bpara[7], 
59             'B60':Bpara[8], 'B62':Bpara[9], 'B6m2':Bpara[10], 'B64':Bpara[11], 
60             'B6m4':Bpara[12], 'B66':Bpara[13], 'B6m6':Bpara[14], 'lam':1.9}
61shape = {'A':0.1, 'a':0., 'b':0.1, 'ei':16.12}
62
63def fix_consts(func, const_dict):
64    # 引数名の取得
65    arg_names = func.__code__.co_varnames[:func.__code__.co_argcount]
66    # 新しい関数(引数を定数化した関数)の引数名リスト
67    var_names = [n for n in arg_names if n not in const_dict.keys()]
68    # 新しい関数内で元の関数を呼ぶ際に渡す引数リスト
69    args = [str(const_dict.get(name, name)) for name in arg_names]
70    # 新しい関数の生成
71    fixed_const_func = eval(
72        f'[lambda {",".join(var_names)}: f({",".join(args)})'
73        f' for f in [func]][0]', {},
74        {'var_names': var_names, 'args': args, 'func': func})
75    return fixed_const_func
76
77
78P = {}
79C = {}
80
81RMSE = 2
82for roop in np.linspace(1,3,3):
83    for i, v in constrain.items():
84        P[i] = v
85        A = dict(constrain.items() ^ P.items())
86        C = random.sample(A.items(),2)
87        P[C[0][0]] = C[0][1]
88        P[C[1][0]] = C[1][1]
89        A = dict(constrain.items() ^ P.items())
90        spec = fix_consts(spectra, A)
91        spec = fix_consts(spec, shape)
92        p0 = [v, C[0][1], C[1][1]]
93        cf = optimize.curve_fit(f=spec, xdata=E_data, ydata=I_data, p0=p0, sigma = er)
94        error = rmse(cf[0], E_data, I_data)
95        print(error)
96        if error < RMSE:
97            constrain[i] = cf[0][0]
98            constrain[C[0][0]] = cf[0][1]
99            constrain[C[1][0]] = cf[1][1]
100            RMSE = error
101        P.clear()
102        A.clear()
103        C.clear()
104    print('roop' + str(roop))
105