为什么我的插值不能在我的函数中正常工作？答案

【问题标题】：Why is my interpolation not working properly in my function?为什么我的插值不能在我的函数中正常工作？
【发布时间】：2018-10-14 06:03:54
【问题描述】：

我有一个相当长的代码来处理光谱，并且在此过程中我需要对一些点进行插值。我曾经在没有任何函数的情况下逐行编写所有这些代码，并且一切正常，但现在我将其转换为两个大函数，以便将来可以更轻松地在其他模型上调用它。下面是我的代码（我在这里的最后一行之后有更多代码，它绘制了一些东西，但这与我的问题无关，因为我已经用一堆 print 行测试了这个并了解到我的问题是在我打电话时出现的我的process 函数中的interpolation 函数。

import re
import numpy as np
import scipy.interpolate

# Required files and lists
filename = 'bpass_spectra.txt' # number of columns = 4
extinctionfile = 'ExtinctionLawPoints.txt' # R_V = 4.0
datalist = []
if filename == 'bpass_spectra.txt':
    filetype = 4
else:
    filetype = 1
if extinctionfile == 'ExtinctionLawPoints.txt':
    R_V = 4.0
else:
    R_V = 1.0 #to be determined

# Constants
h = 4.1357e-15 # Planck's constant [eV s]
c = float(3e8) # speed of light [m/s]

# Inputs
beta = 2.0 # power used in extinction law
R = 1.0 # star formation rate [Msun/yr]
z = 1.0 # redshift
M_gas = 1.0 # mass of gas
M_halo = 2e41 # mass of dark matter halo

# Read spectra file
f = open(filename, 'r')
rawlines = f.readlines()
met = re.findall('Z\s=\s(\d*\.\d+)', rawlines[0])
del rawlines[0]
for i in range(len(rawlines)):
    newlist = rawlines[i].split(' ')
    datalist.append(newlist)

# Read extinction curve data file
rawpoints = open(extinctionfile, 'r').readlines()
def interpolate(R_V, rawpoints, Elist, i):
    pointslist = []
    if R_V == 4.0:
        for i in range(len(rawpoints)):
            newlst = re.split('(?!\S)\s(?=\S)|(?!\S)\s+(?=\S)', rawpoints[i])
            pointslist.append(newlst)
        pointslist = pointslist[3:]
        lambdalist = [float(item[0]) for item in pointslist]
        k_abslist = [float(item[4]) for item in pointslist]
        xvallist = [(c*h)/(lamb*1e-6) for lamb in lambdalist]
        k_interp = scipy.interpolate.interp1d(xvallist, k_abslist)
        return k_interp(Elist[i])

# Processing function
def process(interpolate, filetype, datalist, beta, R, z, M_gas, M_halo, met):
    speclist = []
    if filetype == 4:
        metallicity = float(met[0])
        Elist = [float(item[0]) for item in datalist]
        speclambdalist = [h*c*1e9/E for E in Elist]
        met1list = [float(item[1]) for item in datalist]
        speclist.extend(met1list)
        klist, Tlist = [None]*len(speclist), [None]*len(speclist)
        if metallicity > 0.0052:
            DGRlist = [50.0*np.exp(-2.21)*metallicity]*len(speclist) # dust to gas ratio
        elif metallicity <= 0.0052:
            DGRlist = [((50.0*metallicity)**3.15)*np.exp(-0.96)]*len(speclist)
        for i in range(len(speclist)):
            if Elist[i] <= 4.1357e-3: # frequencies <= 10^12 Hz
                klist[i] = 0.1*(float(Elist[i])/(1000.0*h))**beta # extinction law [cm^2/g]
            elif Elist[i] > 4.1357e-3: # frequencies > 10^12 Hz
                klist[i] = interpolate(R_V, rawpoints, Elist, i) # interpolated function's value at Elist[i]
        print "KLIST (INTERPOLATION) ELEMENTS 0 AND 1000:", klist[0], klist[1000]
        return

print 行的输出是KLIST (INTERPOLATION) ELEMENTS 0 AND 1000: 52167.31734159269 52167.31734159269。

当我在没有函数的情况下运行旧代码时，我会像在此处一样打印klist[0] 和klist[1000]，并为每个代码获取不同的值。在这个新代码中，我从这一行得到两个相同的值。这不应该是这种情况，所以它不能在我的函数内部正确插值（也许它没有在循环中的每个点上正确执行它？）。有没有人有任何见解？在这里发布我的整个代码以及所有使用的文本文件（它们非常大）是不合理的，所以我不希望任何人运行它，而是检查我如何使用和调用我的函数。

编辑：下面是我的代码的原始版本，直到没有函数的插值点（有效）。

import re
import numpy as np
import scipy.interpolate

filename = 'bpass_spectra.txt'
extinctionfile = 'ExtinctionLawPoints.txt' # from R_V = 4.0
pointslist = []
datalist = []
speclist = []

# Constants
h = 4.1357e-15 # Planck's constant [eV s]
c = float(3e8) # speed of light [m/s]

# Read spectra file
f = open(filename, 'r')
rawspectra = f.readlines()
met = re.findall('Z\s=\s(\d*\.\d+)', rawspectra[0])
del rawspectra[0]
for i in range(len(rawspectra)):
    newlist = rawspectra[i].split(' ')
    datalist.append(newlist)

# Read extinction curve data file
rawpoints = open(extinctionfile, 'r').readlines()
for i in range(len(rawpoints)):
    newlst = re.split('(?!\S)\s(?=\S)|(?!\S)\s+(?=\S)', rawpoints[i])
    pointslist.append(newlst)
pointslist = pointslist[3:]
lambdalist = [float(item[0]) for item in pointslist]
k_abslist = [float(item[4]) for item in pointslist]
xvallist = [(c*h)/(lamb*1e-6) for lamb in lambdalist]
k_interp = scipy.interpolate.interp1d(xvallist, k_abslist)

# Create new lists
Elist = [float(item[0]) for item in datalist]
speclambdalist = [h*c*1e9/E for E in Elist]
z1list = [float(item[1]) for item in datalist]
speclist.extend(z1list)
met = met[0]
klist = [None]*len(speclist)
Loutlist = [None]*len(speclist)
Tlist = [None]*len(speclist)

# Define parameters 
b = 2.0 # power used in extinction law (beta)
R = 1.0 # star formation ratw [Msun/yr]
z = 1.0 # redshift
Mgas = 1.0 # mass of gas
Mhalo = 2e41 # mass of dark matter halo

if float(met) > 0.0052:
    DGRlist = [50.0*np.exp(-2.21)*float(met)]*len(speclist)
elif float(met) <= 0.0052:
    DGRlist = [((50.0*float(met))**3.15)*np.exp(-0.96)]*len(speclist)
for i in range(len(speclist)):
    if float(Elist[i]) <= 4.1357e-3: # frequencies <= 10^12 Hz
        klist[i] = 0.1*(float(Elist[i])/(1000.0*h))**b # extinction law [cm^2/g]
    elif float(Elist[i]) > 4.1357e-3: # frequencies > 10^12 Hz
        klist[i] = k_interp(Elist[i]) # interpolated function's value at Elist[i]
print "KLIST (INTERPOLATION) ELEMENTS 0 AND 1000:", klist[0], klist[1000]

print 行的输出是KLIST (INTERPOLATION) ELEMENTS 0 AND 1000 7779.275435560996 58253.589270674354。

【问题讨论】：

可以分享没有这个功能的原版吗？
@wowserx 当然，我会在中编辑它
请按照您创建此帐户时的建议阅读并遵循帮助文档中的发布指南。 Minimal, complete, verifiable example 适用于此。在您发布 MCVE 代码并准确描述问题之前，我们无法有效地帮助您。我们应该能够将您发布的代码粘贴到文本文件中并重现您描述的问题。显示你得到的输出、你期望的输出，以及你对中间结果的跟踪或其他有用的调试尝试。
您能否确认您正在输入interpolate函数（Elist[i] > 4.1357e-3）以及R_V = 4.0在process内？
是的，我已经在我的错误检查过程中这样做了，并且可以确认这两个。

标签： python python-2.7 function interpolation linear-interpolation

【解决方案1】：

您将i 作为参数传递给interpolate，然后还在interpolate 内的循环中使用i。一旦在interpolate 的for i in range(len(rawpoints)) 循环中使用i，它将被设置为某个值：len(rawpoints)-1。然后，interpolate 函数将始终返回相同的值 k_interp(Elist[i])，这等效于 k_interp(Elist[len(rawpoints)-1])。您需要在循环中定义一个新变量（例如for not_i in range(len(rawpoints))），或者为Elist 参数使用不同的变量。考虑对interpolate 进行以下更改：

def interpolate(R_V, rawpoints, Elist, j):
    pointslist = []
    if R_V == 4.0:
        for i in range(len(rawpoints)):
            newlst = re.split('(?!\S)\s(?=\S)|(?!\S)\s+(?=\S)', rawpoints[i])
            pointslist.append(newlst)
        pointslist = pointslist[3:]
        lambdalist = [float(item[0]) for item in pointslist]
        k_abslist = [float(item[4]) for item in pointslist]
        xvallist = [(c*h)/(lamb*1e-6) for lamb in lambdalist]
        k_interp = scipy.interpolate.interp1d(xvallist, k_abslist)
        return k_interp(Elist[j])

【讨论】：