1

I have recorded experimental temperatures at five locations from the surface of a solid. At every time step, I want to fit these readings to a theoretical curve defined by my function: Temp_Function_JLT(X,h).

X is a multi-dimensional array that includes the x_coordinates as well as time, initial temperature and material properties (all independent variables). "h" is the heat transfer coefficient, which for the purpose of this exercise I'm trying to optimize (leaving the physics aside for a moment.)

This is the definition of my temperature function:

 import pandas as pd
 import numpy as np
 import matplotlib.pyplot as plt
 import pickle
 import scipy.optimize as opt
 from scipy.special import erfc


def Temp_Function_JLT(X,ht):
    # Work around the fact that only one independent variable can be passed to optimize.curve_fit
    x,t,T0,q,alpha,rho,c,k = X 

    term_a = q/ht
    term_b = erfc(x/np.sqrt(4*alpha*t))
    term_c = np.exp(((ht*x)/(np.sqrt(alpha)*np.sqrt(k*rho*c)))+((ht**2)/(k*rho*c)))
    term_d = erfc((ht*np.sqrt(t))/(np.sqrt(k*rho*c)) + (x/np.sqrt(4*alpha*t)))
    Temperature = (term_a * (term_b - term_c * term_d)) + T0 - 273

return Temperature

The function works. I can run it with some initial parameters and obtain sensible values. More importantly for this question, if I call it with the following data:

t = 1
x_test = np.linspace(0.004,0.02,5) # TC locations
time_test = range(1,180,30)
T0_test = 25 + 273
q_test = 20000
h_test = 10

I will obtain a numpy array as a solution of shape (1,) which gives an answer to np.ndim of 1 (This has been mentioned in the following previous questions:

Least Linear Squares: scipy.optimize.curve_fit() throws "Result from function call is not a proper array of floats."

Fitting a vector function with curve_fit in Scipy

Fitting a 2D Gaussian function using scipy.optimize.curve_fit - ValueError and minpack.error

The problem arises when I call opt.curve_fit(). indepth_temperatures is a list that contains each test as an array. I iterate over it (to iterate over each test) and then I perform the fit on each row (each time step), according to the following code:

for i,test in enumerate(indepth_temperatures):

    # Iterate over every row
    for j,row in enumerate(test):
        # Define tuple that contains all independent variables
        X = (TC_depth,
             times[i][j],
             T0_temperatures[i] + 273,
             20000,
             pmma_alpha,
             pmma_rho,
             pmma_c,
             pmma_k)
        print(Temp_Function_JLT(X,h0))
        print(row)
        print('---')
        # Call function to optimize curve fit on h
        popt, pcov = opt.curve_fit(Temp_Function_JLT,X,row,h0)
        print(popt)

For the first iteration, I obtain the following result:

[23.2034 23.2034 23.2034 23.2034 23.2034]   # comes from print(Temp_Function_JLT(X,h0))
[23.937 22.619 22.59 24.884 21.987000000000002]  # comes from print(row)

Followed by this error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
<ipython-input-67-9c4545fd257b> in <module>()
     22         print('---')
     23         # Call function to optimize curve fit on h
---> 24         popt, pcov = opt.curve_fit(Temp_Function_JLT,X,row,h0)
     25         print(popt)

~\AppData\Local\Continuum\anaconda2\envs\py36\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    749         # Remove full_output from kwargs, otherwise we're passing it in twice.
    750         return_full = kwargs.pop('full_output', False)
--> 751         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    752         popt, pcov, infodict, errmsg, ier = res
    753         cost = np.sum(infodict['fvec'] ** 2)

~\AppData\Local\Continuum\anaconda2\envs\py36\lib\site-packages\scipy\optimize\minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    392         with _MINPACK_LOCK:
    393             retval = _minpack._lmdif(func, x0, args, full_output, ftol, xtol,
--> 394                                      gtol, maxfev, epsfcn, factor, diag)
    395     else:
    396         if col_deriv:

error: Result from function call is not a proper array of floats.

I have tried returning from my function np.ravel(Temperature) or Temperature.flatten() with no luck. The error remains, and I can't figure out why it's there. As I mentioned, I have checked the dimensions of the return of my function and it is a 1D array.

Any help will be greatly appreciated!

UPDATE: I realized it was hard to replicate this code, so this is a simplified version:

Temp_Function_JLT(X,h0): stays the same. 

pmma_rho = 1200 # kg/m3
pmma_c = 1500 # J/kgK
pmma_k = 0.16 # W/mK
pmma_alpha = pmma_k/(pmma_rho*pmma_c)
x_test = np.linspace(0.004,0.02,5) # TC locations
t = 1
T0_test = 25 + 273
q_test = 20000
h_test = 10

X = (x_test,t,T0_test,q_test,pmma_alpha,pmma_rho,pmma_c,pmma_k)
y_data = [23.937 22.619 22.59 24.884 21.987000000000002]

opt.curve_fit(Temp_Function_JLT, X, y_data, h_test)
SimonS
  • 36
  • 6
  • 1
    I can't run your code because you didn't specify what `indepth_temperatures` is. Also, there are so many other variables which are unknown: `h0`, and all others starting with `pmma___` – Sheldore Sep 10 '18 at 10:32
  • Thanks for the answer! pmma_rho = 1200 pmma_c = 1500 pmma_k = 0.16 pmma_alpha = pmma_k/(pmma_rho * pmma_c) h0 = 10 indepth_temperatures is a list, that contains the temperature readings per test (20 tests in total). Each tests contains five temperature readings (my x_data) per time step (around 1200 time steps). That is why I printed in my results the first row, which is the "y_data" for the curve fit in that first time step. In other words, my curve fit uses: Temp_Function_JLT as function, row (printed in my result) as y_data, and X as x_data. I'll simplify my post. – SimonS Sep 10 '18 at 10:43
  • Mmmm right. So, for some reason, my row (y_data) is a numpy array with dtype = object. I think this is the source of the problem. If I print row, this is what I get: array([23.937, 22.619, 22.59, 24.884, 21.987000000000002], dtype=object) – SimonS Sep 10 '18 at 11:21
  • 1
    Sorry, `indepth_temperatures` is still not defined for me although you explained what it contains. In the interest of time of the people trying to help you, it would be good if you define all variables and copy paste the code above in a new python file and execute it. **Only** if you don't get any undefined variables error, post the code. Otherwise we will keep asking you again and again that this variable is missing and that variable is missing – Sheldore Sep 10 '18 at 11:33

1 Answers1

1

I realized what was wrong with my code. Even though my y_data (row) was defined as a 1-D numpy array, its data type was object. I don't yet understand why this was the cause, but by forcing the data type with np.astype(np.float), opt.curve_fit worked.

SimonS
  • 36
  • 6