--- SOLVED ----

u/Standardw and u/kra_pao provided me with the answers: Passing an array-like initial set of parameters to curve_fit helps, i.e.:

# define start values for scalar parameters (so that len(args) is known)
p = [1] * len(expList)

# fit data to model
popt, pcov = curve_fit(polymodel, xSamples, valSamps,  p0 = p, maxfev = 10000)

However, with this, there is a warning message from curve_fit. This new problem seems to originate in my polynomial model. I will work on that and of course post the solution here.

UPDATE:

Fixed it completely and posted resulting code in my comment below.

Original Post:

Hi r/learnpython!

This subreddit has really helped me a lot and continues to do so when I'm trying to code something in python. Thank you all for your great posts!

Now I turn to you for help.

Context

I'm working on my thesis (mechanical engineering, so I don't know a lot about computers or programming). In a nutshell, I'm trying out different ways to optimize the shape of a mechanical structure and compare multiple optimization techniques. Python is the language to go to for that because it's easily automatable, can directly start my FEM programs with alternative parameters and so on. The code will not be published.

My Specific Problem

I have some sample data that I want to fit to a polynomial model. What I'm too stupid for: I can't figure out how to "create the model" during running the program so that it fits my requirements (i.e. number of dimensions and maximum degree). I tried to boil it down a bit, so as an example I use the ackley function here.

So I have sample data (x, y and function values) of the ackley function and want to fit a model to it. (This level of inaccuaracy is fine - I just need the minimum or can worry about the accuaracy later)

When I use the function also as a model, it works. That's how I got the model in the second picture. When I try to dynamically build the model (sorry, I don't know how to express it better), the code fails:

ValueError: Unable to determine number of fit parameters.

This makes sense: I don't specify the number of *args. I do understand that. But how do I fix it?

The Example Code

Remarks for the functions are below the code block. My problem part is set within "###############" above and below. It's at "# fit data to model".As you can see, I try to define my polynomial model on the run. Doesn't seem elegant to me, I just kind of feel helpless here.

import numpy as np
from itertools import product
from scipy.optimize import curve_fit

def myfunc(x, a=20, b=0.2, c=2*np.pi):
    # ackley function
    d = len(x)
    sumOne = sum([xi**2 for xi in x])
    sumTwo = sum([np.cos(c*xi) for xi in x])
    val = -a * np.exp(-b * np.sqrt(1/d * sumOne)) - np.exp( 1/d * sumTwo) + a + np.exp(1)
    return val


def getpolyexponents(nDim, kDegree):
    """gets the basis function exponents for a polynom"""

    # repeat range {0,...,kDegree} nDim times as rangeSet
    rangeSet = [range(0, kDegree+1) for i in range(0, nDim)]

    # create full product set
    expList = [item for item in product(*rangeSet)]

    # keep only items in list with correct degree,
    # i.e. sum(exponents) <= kDegree            
    expList[:] = [item for item in expList if sum(item) <= kDegree]

    # return the polynom exponents
    return expList

def getbasisfunc(x, exponentSet):
    """gets a single basis function for a polynom with an exponent set"""
    # make list of single elements of basis function, e.g. x^0*y^0, x^0*y^1, ...
    productElements = [x[i]**exponentSet[i] for i in range(0,len(exponentSet))]

    # calculate product of all elements and return it
    val = np.prod(productElements)
    return val

def getbasisfuncarray(x, expList):
    """ builds an array (1-D) of basis functions for a polynom"""
    basisfunclist = [getbasisfunc(x, expList[i]) for i in range(0, len(expList))]
    return np.array(basisfunclist)

def main():    
    # dim d, polynomial degree k
    d = 2
    k = 5
    
    # create n random sample points xSamples in [-5, 5]
    n = 100
    xSamples = np.divide(np.random.randint(-500, 500 + 1, (d,n)), 100)

    # values of samples
    valSamps = myfunc(xSamples)

    # find exponents for polynomial model function
    expList = getpolyexponents(d, k)

    # define model
    def polymodel(x, *args):
        coeffArray = np.array([arg for arg in args])
        basisFuncs = getbasisfuncarray(x, expList)        
        # scalar product / dot product
        poly = np.dot(coeffArray, basisFuncs)
        return poly

################################################
    # fit data to model
    #popt, pcov = curve_fit(polymodel, xSamples, valSamps,  maxfev = 1000) #fails
    popt, pcov = curve_fit(myfunc, xSamples, valSamps,  maxfev = 1000) #works
################################################
    

    # print fitting results
    print("fitting results:")
    print(popt)
    print(pcov)


if __name__ == "__main__": main()

• myfunc: just the ackley function

• getpolyexponents: Since I want to create a polynomial model, I need the exponents for all my variables (x1 and x2 or also called x and y in 2-D). For example: In 2D the polynomial basis functions are:

  • x^0*y^0

  • x^0*y^1

  • x^0*y^2

  • x^1*y^0

  • x^1*y^1

  • x^2*y^0

  • I get those exponents from this function. In the 2-D case there would be those 6 sets of exponents in the returned expList

• getbasisfunc: Uses one exponent set from the expList from above function to create the actual basis function

• getbasisfuncarray: puts all required basis functions in a 1-D Array.

I did not include the plotting part because it's in another module and would just be more code in here.

If there is a better way to fit my data to a polynomial model of d dimensions (d being up to 23 in my case) and degree k (k up to 10 approx.) with python I am happy to learn new stuff!

Thank you already in advance. Any tips are greatly appreciated!

Cheers