Feed variable number of parameters to fit 2D fourier series

Question

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)

Answer 1

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])