Feed variable number of parameters to fit 2D fourier series
I am working on an optimization problem using 2D Fourier series, basically I am trying to implement this formula:
and I wanted to have it implemented in my script as a function of the order m, in the sense that provided m I can fit some data with a m-order series and obtain the best coefficients of the series (except for C_00 that I already know)
Anyway I am struggling when it comes to use the curve_fit function from Scipy with a model with a variable number of parameters (eg. if m=1 I have 8 coefficients + 2 angular velocities, if m=2 I have 24 coefficients + 2 angular velocities and so on)
What I have written so far is the following:
c00=50
def fourier(x, v, w1, w2):
f=c00
k=0
for i in range(1, m+1):
f=f+v[k]*np.cos(i*w1*x)+v[k+1]*np.sin(i*w1*x)
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+v[k]*np.cos(i*w1*x+j*w2*x)+v[k+1]*np.sin(i*w1*x+j*w2*x)
k+=2
return f
time = [] #some x data I get from a txt file
energy = [] # some y data I get from a txt file
popt, pcov = curve_fit(fourier, time, energy)
But when I run the script it raises the following error:
File "/home/arcalinux/Desktop/fourier_series.py", line 63, in fourier f=f+v[k]np.cos(iw1x)+v[k+1]np.sin(iw1x) IndexError: invalid index to scalar variable.
I also tried this other function:
def fourier2(args):
f=c00
k=3
for i in range(1, m+1):
f=f+args[k]*np.cos(i*args[1]*args[0])+args[k+1]*np.sin(i*args[1]*args[0])
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+args[k]*np.cos(i*args[0]*x+j*args[1]*x)+
+args[k+1]*np.sin(i*args[0]*x+j*args[1]*x)
k*=2
return f
but when I run curve_fit I have:
ValueError: Unable to determine number of fit parameters.
please notice that the fourier series function works correctly, for example (random numbers):
v=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
m=2
x=2
w1=3
w2=4
print(fourier(x,v,w1,w2))
49.136961626348295
And the same goes for fourier2
How can I solve this problem? Thank you!
the whole code is here:
#importing libraries
import numpy as np
import math
from scipy.optimize import curve_fit
from statistics import mean
import matplotlib.pyplot as plt
m=2 #order of series
number_of_coeff=(2*m+1)**2 - 1
N=number_of_coeff
#build series
coefficients=[]
#first part, j=0
for i in range(1,m+1):
coefficients.append('c'+str(i)+'0')
coefficients.append('s'+str(i)+'0')
#second part
for j in range(1, m+1):
for i in range(-m, 0):
coefficients.append('g'+str(-i)+str(j))
coefficients.append('t'+str(-i)+str(j))
for i in range(0, m+1):
coefficients.append('c'+str(i)+str(j))
coefficients.append('s'+str(i)+str(j))
print(len(coefficients))
v=np.ones(N)
m=2
x=2
w1=3
w2=4
arg=[2,3,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
arg3=[3,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
c00=50
def fourier(x, v, w1, w2):
print(type(v))
print('v is ',v)
f=c00
k=0
print(v[0])
for i in range(1, m+1):
f=f+v[k]*np.cos(i*w1*x)+v[k+1]*np.sin(i*w1*x)
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+v[k]*np.cos(i*w1*x+j*w2*x)+
v[k+1]*np.sin(i*w1*x+j*w2*x)
k+=2
return f
def fourier2(args):
f=c00
k=3
for i in range(1, m+1):
f=f+args[k]*np.cos(i*args[1]*args[0])+
args[k+1]*np.sin(i*args[1]*args[0])
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+args[k]*np.cos(i*args[1]*args[0]+
j*args[2]*args[0])+
args[k+1]*np.sin(i*args[1]*args[0]+
j*args[2]*args[0])
k+=2
#print('number of args =',len(args))
return f
def fourier3(x,args):
f=c00
k=2
for i in range(1, m+1):
f=f+args[k]*np.cos(i*args[0]*x)+
args[k+1]*np.sin(i*args[0]*x)
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+args[k]*np.cos(i*args[0]*x+j*args[1]*x)+
args[k+1]*np.sin(i*args[0]*x+j*args[1]*x)
k+=2
#print('number of args =',len(args))
return f
#print(fourier2(arg))
#print(fourier3(x,arg3))
#print(fourier(2,v,3,4))
#getting data
time = []
energy = []
filepath=''
filename=''
with open(filepath+filename, 'r') as f:
for line in f: # iterate through the file
if not line: continue
t, e = line.split()
time.append(float(t)),
energy.append(float(e)) # add a value
c00=mean(energy)
time=np.array(time)
popt, pcov = curve_fit(fourier, time, energy)
print(popt)
According to scipy documentation, you can provide to the curve_fit function the p0 parameter, which is your initial guess to your function's parameters (here, v, w1 and w2).
If you don't, the function will assume an initial value of 1 for your parameters, which is why you get a scalar instead of an array for v.
So what you could do before calling curve_fit is:
initial_v = ... # Replace with value (an array)
initial_w1 = ... # Replace with value (a scalar)
initial_w2 = ... # Replace with value (a scalar)
popt, pcov = curve_fit(fourier, time, energy, p0=[initial_v, initial_w1, initial_w2])