Several questions regarding LevenbergMarquadt implementation in Scipy
I came across a few oddities while trying to optimize a function using LevenbergMarquadt algorithm in Scipy. I was wondering if anyone had an explanation for them.
A. I have an objective function that accepts two parameters with 2 residuals. As far as I understand, LM for scipy root
using sum(residual^2)
and least_squares
should be the same. However, least_squares
returns an optimized value while root
fails. Are the implementations for the two functions different?
B. When I optimize the function using least_squares
, the second parameter does not change from its initial value. Why does this happen?
Thank you for your help!
See also questions close to this topic

Any way to get barchart (upsidedown) with percent change in python?
I have trade data (import /export of grocery products by worldwide), and I intend to render some plots that indicates worldwide importexport trade flow. To do so, I manipulated my dataframe by using
group_by
,pivot_table
that will be used as input for plotting function. Seems I didn't get correct manipulation on my data, as a result, I didn't get correct plot that I expected for. Can anyone point me out how to make this work? Thanksthis is the original data that I am experimenting with:
based on this original data, I want to show top 10 countries that get beef export from US. How can I get that?
minimal example:
import numpy as np import pandas as pd df = pd.DataFrame(data=[list('EEEIEEIIEI'), ['AR', 'AUC', 'CA', 'CN', 'MX', 'MX', 'AR', 'IT', 'UK', 'RU'], ['ALBANIA', 'PAKISTN', 'UGANDA', 'FRANCE', 'USA', 'RUSSIA', 'COLOMBIA', 'KAZAK', 'KOREA', 'JAPAN'], [20230, 20220, 20120, 20230, 20230, 20220, 20230, 20120, 20130, 20329], list(np.random.randint(10, 100000, 10)), list(np.random.randint(10, 100, 10)), list(np.random.randint(10, 100, 10)), np.random.choice(pd.date_range('1/1/2014', periods=365, freq='D'), 10, replace=False)] ).T df.columns =['ID', 'cty', 'cty_ptn', 'prod_code', 'value','Quantity1', 'Quantity2', 'date'] my_dict={'20230':'Gas', '20220':'Water', '20210': 'Refined', '20120':'Oil', '20239':'Other'} df['prod_label']=df['prod_code'].astype(str).map(my_dict)
my attempt:
I tried the following manipulation for creating bar plot:
import matplotlib.pyplot.plt %matplotlib inline mydf= mydf[(mydf['prod_label']=='Gas') & (mydf['ID']=='E')]
updated attempt:
mydf=mydf.groupby('cty_ptn')['value','Quantity1'].sum().reset_index() mydf=mydf.nlargest(20, 'Quantity1') ## EDA fig=plt.figure() ax=fig.add_axes([0,0,1,1]) ax.scatter(mydf['Quantity1'].values, mydf['value'].values)
my current output
I got this dummy plot that is far from what I expected:
but I didn't get the correct plot here (it is empty, I don't know why). The reason is most likely I got wrong manipulated data which not match for the input for
matplotlib
plot function. Any quick correction?desired plot:
I am trying to achieve the following plot (inspired by trade data report):
how can I get above plot? what sort of manipulation is needed before hitting plot function? Thanks

Combine rows of a single file based on multiple matched columns, while retaining unique content by appending additional columns
I'm sure this code exists and I've read through a LOT of pandas / python documentation, in fact perhaps my answer is contained within Pandas, append column based on unique subset of column values but we can't seem to get it to work as below.
Using the example below, if both the Company and Place match, we want to combine the rest of the columns. IF there is unique data, then we would to retain the data in that column and append it as an additional column.
Here is the visual representation of what we need:

indices of nested lists
I am trying to solve a problem that is a part of my genome alignment project. The problem goes as follows: if given a nested list
y = [[1,2,3],[1,2,3],[3,4,5],[6,5,4],[4,2,5],[4,2,5],[1,2,8],[1,2,3]]
extract indices of unique lists into a nested list again.
For example, the output for the above nested list should be
[[0,1,7],[2],[3],[4,5],[6]]
.This is because list
[1,2,3]
is present in0,1,7th
index positions,[3,4,5]
in 2nd index position and so on.Since, I will be dealing with large lists, what could be the most optimal way of achieving this in python? Thanks in advance.

JS how to optimize function that finds the text value in gematria
Hello guys i wrote this function code:
function gematria(text){ var newText = text.toLowerCase(); // handle capital letters var startValue=0; //initialize value to zero for(let i=0 ; i<newText.length; i++){ // loop through every letter on the text let currentCharValue = newText.charCodeAt(i)  96; // get the distance from the 'a' startValue +=currentCharValue; } return startValue } gematria("someText");
it works well but is there any way to do it more short? or better? how good is that code for doing the task?

Python performance issue and name check confusion
I have a question about some of my python code and apparently could use some clarification on the functionality of the name check in python ( if name='main').
The way I understand it is that you put the name check in the main python module that calls all of the others to begin your loop/call your programs functions. When you run this 'main' program, it should execute everything. If you call the 'main' module, some lines won't run. In the modules below, it isn't executing my Logger.writeAndPrintLine() lines, but if I take out the name check it does.. what am I not understanding here?
Also, I recently redid this "program" and consolidated all of the functions into just a few modules. When it was created, functions were spread out amongst almost 10 different modules... but it ran much faster. Is it better to keep the functions separate to their own modules? Is there something inherently wrong with my code? I think the problem might be related to all the config file calls, as the variables used to be all stored in their modules, not in a config file. Any input is appreciated.
main.py:
from Logger import Logger from s3nd import sleepLoop #from receiveListener import receiveListener Logger.writeAndPrintLine('Beginning program execution...', 1) if __name__ == "__main__": Logger.writeAndPrintLine('Starting send listener loop...', 1) sleepLoop() s3nd.py: from Logger import Logger import pyodbc from configparser import ConfigParser from datetime import datetime from re import search import requests import time config = ConfigParser() config.read('c0nfig.ini') location = config.get('KENECT', 'location') apiToken = config.get('KENECT', 'apiToken') apiKey = config.get('KENECT', 'apiKey') dbUser = config.get('DB', 'dbUser') dbPassword = config.get('DB', 'dbPassword') dbHost = config.get('DB', 'dbHost') pullId = config.get('QUERY', 'pullId') pullNewest = config.get('QUERY', 'pullNewest') updateFin = config.get('QUERY', 'updateFin') last_sent = '' def sendLoop(): while 0 == 0: sendListener() time.sleep(2.0) def sendListener(): global last_sent cnxn = pyodbc.connect('UID='+dbUser+';PWD='+dbPassword+';DSN='+dbHost) cursor = cnxn.cursor() cursor.execute(pullId) results = cursor.fetchone() if results != last_sent: sendSMS() else: cnxn.close() def sendSMS(): msg, destPhone, msgId, caseNum, sender = pullData() msg = msg destPhone = destPhone sender = sender caseNum = caseNum msgId = str(msgId) sendCT = 'multipart/formdata; boundary=' success='"status":"SENT"' headers = { 'xapitoken': apiToken, 'xapikey': apiKey, 'contenttype': sendCT } formData = { 'locationId': location, 'messageBody': msg, 'contactPhone': destPhone } Logger.writeAndPrintLine('[REQ: ' + msgId + '] ' + sender + ' to ' + destPhone, 1) r=requests.post('https://api.kenect.com/v1/conversations/messages', headers=headers, params=formData).text if search(success, r): Logger.writeAndPrintLine('[RES: ' + msgId + '] ' + r, 1) Logger.writeAndPrintLine('SMS successful', 1) updateSent(msgId) else: Logger.writeAndPrintLine('SMS failure; Kenect error.', 3) def pullData(): cnxn = pyodbc.connect('UID='+dbUser+';PWD='+dbPassword+';DSN='+dbHost) cursor = cnxn.cursor() outbound = cursor.execute(pullNewest) results = cursor.fetchone() collection = {} for index, outbound in enumerate(results): key_name = "col{0}".format(index) collection[key_name] = outbound msg = collection['col1'] destPhone = collection['col2'] msgId = collection['col3'] caseNum = collection['col4'] sender = collection['col10'] cnxn.close() return (msg, destPhone, msgId, caseNum, sender) def updateSent(msgId): global last_sent cnxn = pyodbc.connect('UID='+dbUser+';PWD='+dbPassword+';DSN='+dbHost) cursor = cnxn.cursor() cursor.execute(updateFin + msgId) cnxn.commit() cursor.execute(pullId) last_sent = cursor.fetchone() cnxn.close() sendLoop()

Speeding up a pandas iteration looking testing a condition on subsequent elements
Given a pandas dataframe with three columns (C1,C2,C3) and a Series of positive numbers of equal length (coeff) I am computing the fourth column C4, as follows
def event(data, coeff, rate_low=2, rate_high=2): bot_col_name = 'C4' data[bot_col_name] = 1 I = data.index for k in range(len(I)1): i = I[k] next_val = data.at[ I[k+1], 'C1'] c = coeff.at[i] low_bound = next_val  rate_low*c high_bound = next_val + rate_high*c for j in range(k+1, len(data)): if data.at[ I[j], 'C2'] < low_bound: data.at[i, bot_col_name] = 0 break if data.at[ I[j], 'C3'] >= high_bound: data.at[i, bot_col_name] = 1 break return data
In other words, given a row, we compute a certain upper and lower bound and then set the respective row element depending on whether we first hit the upper bound under C2 or lower bound on C3.
Unfortunately, the above solution is quite inefficient. I've tried optimizing it by computing the values backwards and trying to cache the results (sometimes one can infer the value of C4 from 'past' values) but unfortunately its not significantly better.
In my experience, the best way to gain maximum performance is to try to express as much as possible within the pandas framework.
Is there any meaningful way that one could optimize the above code?

Forecasting / Extrapolating with Scipy Curve_Fit
I need to predict / forecast / extrapolate values past
20010115
using Scipy's curve_fit operation. How do I predict out past20010115
and into20010120
?import pandas as pd import numpy as np from datetime import timedelta from scipy.optimize import curve_fit def hyperbolic_equation(t, qi, b, di): return qi/((1.0+b*di*t)**(1.0/b)) df1 = pd.DataFrame({ 'date': ['20010101','20010102','20010103', '20010104', '20010105', '20010106','20010107','20010108', '20010109', '20010110', '20010111','20010112','20010113', '20010114', '20010115'], 'cumsum_days': [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15], 'prod': [800, 900, 1200, 700, 600, 550, 500, 650, 625, 600, 550, 525, 500, 400, 350]}) df1['date'] = pd.to_datetime(df1['date']) qi = max(df1['prod']) #Hyperbolic curve fit the data to get best fit equation popt_hyp, pcov_hyp = curve_fit(hyperbolic_equation, df1['cumsum_days'], df1['prod'],bounds=(0, [qi,1,20])) #Adding in predicted values back into df1 df1.loc[:,'Hyperbolic_Predicted'] = hyperbolic_equation(df1['cumsum_days'], *popt_hyp)
Here I create a future date df (test set)
df1['future_date'] = df1['date'] ftr = (df1['future_date'] + pd.Timedelta(5, unit='days')).to_frame() #Constructs empty columns for ftr dataframe for col in df1.columns: if col not in ftr.columns: ftr[col] = None #Subset future dataframe to predict on (test set) ftr = ftr[(ftr['future_date'] > max(df1['date']))] ftr['cumsum_days'] = [16,17,18,19,20]
This snippet will join the future data set with the original dataset (if needed)
df1 = pd.concat([df1, ftr], ignore_index=True) print(df1) Hyperbolic_Predicted cumsum_days date future_date prod 0 931.054472 1 20010101 20010101 800 ... 14 409.462743 15 20010115 20010115 350 15 NaN 16 NaT 20010116 None 16 NaN 17 NaT 20010117 None 17 NaN 18 NaT 20010118 None 18 NaN 19 NaT 20010119 None 19 NaN 20 NaT 20010120 None
Once I rerun the
curve_fit
operation, I get an error. How do I predict out past20010115
and into20010120
?popt_hyp, pcov_hyp = curve_fit(hyperbolic_equation, df1['cumsum_days'], df1['prod'],bounds=(0, [qi,1,20]))
Error:
TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Behaviour of affine transform on 3D image with nonuniform resolution with Scipy
I'm looking to apply an affine transformation, defined in homogeneous coordinates on images of different resolutions, but I encounter an issue when one ax is of different resolution of the others.
Normally, as only the translation part of the affine is dependent of the resolution, I normalize the translation part by the resolution and apply the corresponding affine on the image, using scipy.ndimage.affine_transform.
If the resolution of the image is the same for all axes, it works perfectly, you can see below images of the same transformation (being a scale+translation, or rotation+translation, see code below) being applied to an image, its downsampled version (meaning at a lower resolution). Images match (almost) perfectly, differences in voxel values are mainly caused as far as I know by interpolation errors.
But you can see that the shapes overlay between the downsampled transformed image and the transformed (downsampled for comparison) image
Scale affine transformation applied on the same image, at two different (uniform) resolutions
Rotation affine transformation applied on the same image, at two different (uniform) resolutions
Unfortunately, if one of the image axis has a different resolution than the other (see code below), it works well with affine transform with null nondiagonal terms (like translation, or scaling) but the result of the transformation gives a completely wrong result.
Rotation affine transformation applied on the same image, at two different (nonuniform) resolutions
Here you can see a minimal working example of the code:
import numpy as np import nibabel as nib from scipy.ndimage import zoom from scipy.ndimage import affine_transform import matplotlib.pyplot as plt ################################ #### LOAD ANY 3D IMAGE HERE #### ################################ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TO BE DEFINED BY USER orig_img = any 3D grayscale image ndim = orig_img.ndim ################################ ##### DEFINE RESOLUTIONS ####### #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TO BE DEFINED BY USER # Comment/uncomment to choose the resolutions (in mm) of the images # ORIG_RESOLUTION = [1., 1., 1.] # TARGET_RESOLUTION = [2., 2., 2.] ORIG_RESOLUTION = [1., 0.5, 1.] TARGET_RESOLUTION = [2., 2., 2.] ##################################### ##### DEFINE AFFINE TRANSFORM ####### affine_scale_translation = np.array([[2.0, 0.0, 0.0, 150.], [0.0, 0.8, 0.0, 0. ], [0.0, 0.0, 1.0, 0. ], [0.0, 0.0, 0.0, 1.0]]) a = np.sqrt(2)/2. affine_rotation_translation = np.array([[a , a , 0.0, 50.], [a , a , 0.0, 100. ], [0.0, 0.0, 1.0, 0.0 ], [0.0, 0.0, 0.0, 1.0]]) # #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TO BE DEFINED BY USER # Comment/uncomment to choose the transformation to be applied # affine_tf, name_affine = affine_scale_translation, "Tf scale" affine_tf, name_affine = affine_rotation_translation, "Tf rotation" ###################################################### ######## DOWNSAMPLE IMAGE TO LOWER RESOLUTION ######## ###################################################### downsample_img = zoom(orig_img, zoom=np.array(ORIG_RESOLUTION)/np.array(TARGET_RESOLUTION), prefilter=False, order=1) ############################################################################## ######## APPLY AFFINE TRANSFORMATION TO ORIGINAL AND DOWNSAMPLE IMAGE ######## ############################################################################## affine_st_full_res, affine_st_low_res = affine_tf.copy(), affine_tf.copy() # Inverse transform as affine_transform apply the tf from the target space to the original space affine_st_full_res, affine_st_low_res = np.linalg.inv(affine_st_full_res), np.linalg.inv(affine_st_low_res) # Normalize translation part (normally expressed in millimeters) for the resolution affine_st_full_res[:ndim, ndim] = affine_st_full_res[:ndim, ndim] / ORIG_RESOLUTION affine_st_low_res[:ndim, ndim] = affine_st_low_res[:ndim, ndim] / TARGET_RESOLUTION # Apply transforms on images of different resolutions orig_tf_img = affine_transform(orig_img, affine_st_full_res, prefilter=False, order=1) downsample_tf_img = affine_transform(downsample_img, affine_st_low_res, prefilter=False, order=1) # Downsample result at full resolution to be compared to result on downsample image downsample_orig_tf_img = zoom(orig_tf_img, zoom=np.array( ORIG_RESOLUTION)/np.array(TARGET_RESOLUTION), prefilter=False, order=1) # print(orig_img.shape) # print(downsample_img.shape) # print(orig_tf_img.shape) # print(downsample_orig_tf_img.shape) ############################### ######## VISUALISATION ######## ############################### # We'll visualize in 2D the slice at the middle of the z (third) axis of the image, in both resolution mid_z_slice_full, mid_z_slice_low = orig_img.shape[2]//2, downsample_img.shape[2]//2 fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3) ax1.imshow(orig_img[:, :, mid_z_slice_full], cmap='gray') ax1.axis('off') ax1.set_title('1/ Origin image, at full res: {}'.format(ORIG_RESOLUTION)) ax2.imshow(downsample_img[:, :, mid_z_slice_low], cmap='gray') ax2.axis('off') ax2.set_title('2/ Downsampled image, at low res: {}'.format(TARGET_RESOLUTION)) ax3.imshow(downsample_tf_img[:, :, mid_z_slice_low], cmap='gray') ax3.axis('off') ax3.set_title('3/ Transformed downsampled image') ax4.imshow(orig_tf_img[:, :, mid_z_slice_full], cmap='gray') ax4.axis('off') ax4.set_title('4/ Transformed original image') ax5.imshow(downsample_tf_img[:, :, mid_z_slice_low], cmap='gray') ax5.axis('off') ax5.set_title('5/ Downsampled transformed image') error = ax6.imshow(np.abs(downsample_tf_img[:, :, mid_z_slice_low] \ downsample_orig_tf_img[:, :, mid_z_slice_low]), cmap='hot') ax6.axis('off') ax6.set_title('Error map between 3/ and 5/') fig.colorbar(error) plt.suptitle('Result for {} applied on {} and {} resolution'.format(name_affine, ORIG_RESOLUTION, TARGET_RESOLUTION)) plt.tight_layout() plt.show()

How to use cross_validation with curve_fit?
I use
scipy.optimize.curve_fit
to check models on my data.Normally, I would have a function of several variables, here y_data is a function of 2 variables in x_data:
import numpy as np from scipy.optimize import curve_fit x_data = np.array( [np.arange(10), np.arange(10)*1+5] ) print( x_data )
[[ 0 1 2 3 4 5 6 7 8 9]
[ 5 4 3 2 1 0 1 2 3 4]]
y_data = np.arange(10)  5 print( y_data )
[5 4 3 2 1 0 1 2 3 4]
To fit the function, I just define it, and use
curve_fit
:def lincomb( X, a, b ): x1 = X[0] x2 = X[1] return a*x1*x2 + b popt, pcov = curve_fit( lincomb, x_data, y_data ) print( popt )
[0.17857143 1.57142857]
The latter are the optimized values of the coefficients a and b in the function.
Now, I would like to use cross validation from
sklearn
to do the same fit. For this, I packed my function into a class to be used as an estimator, like this:from sklearn.model_selection import cross_validate class LinComb: def __init__( self, a=None, b=None ): self.a = a self.b = b def _lincomb_background(self, X, a, b): x1 = X[0] x2 = X[1] return a*x1*x2 + b def predict( self, X ): return self._lincomb_background( X, self.a, self.b ) def fit( self, X, y ): from scipy.optimize import curve_fit popt, pcov = curve_fit( self._lincomb_background, X, y ) self.a = popt[0] self.b = popt[1] return self def get_params( self, deep=False ): return { 'a':self.a, 'b':self.b } def set_params( self, **parameters ): for parameter, value in parameters.intems(): setattr( self, parameter, value ) return self
When I then call the cross validation, I get an error in dimensions:
cross_validate( LinComb(), x_data, y_data, cv=5, scoring='neg_mean_squared_error' )  ValueError Traceback (most recent call last) <ipythoninput60e0ff8bb83213> in <module> > 1 cross_validate( LinComb(), x_data, y_data, cv=5, scoring='neg_mean_squared_error' ) /usr/local/lib/python3.7/distpackages/sklearn/model_selection/_validation.py in cross_validate(estimator, X, y, groups, scoring, cv, n_jobs, verbose, fit_params, pre_dispatch, return_train_score, return_estimator, error_score) 215 216 """ > 217 X, y, groups = indexable(X, y, groups) 218 219 cv = check_cv(cv, y, classifier=is_classifier(estimator)) /usr/local/lib/python3.7/distpackages/sklearn/utils/validation.py in indexable(*iterables) 228 else: 229 result.append(np.array(X)) > 230 check_consistent_length(*result) 231 return result 232 /usr/local/lib/python3.7/distpackages/sklearn/utils/validation.py in check_consistent_length(*arrays) 203 if len(uniques) > 1: 204 raise ValueError("Found input variables with inconsistent numbers of" > 205 " samples: %r" % [int(l) for l in lengths]) 206 207 ValueError: Found input variables with inconsistent numbers of samples: [2, 10]
Indeed, the dimensions of x_data and y_data are:
print( x_data.shape, y_data.shape )
(2, 10) (10,)
Still I do not understand why the dimensions work in the first, simple case, and how to achieve this cleanly.
Am I missing something?

Finding params using leastsq method for levenberg marquardt
I am struggling to plot the data for my leastsq model and equation as I keep getting shape errors. I'm just trying to guess the parameters that fit the data but I can't do that if I can't see the graphs.
Would really appreciate any help been stuck on this for days! :(
Here is what I have so far:
from scipy.integrate import odeint Import numpy as np Import matplotlib.pyplot as plt #data didata1 = np.loadtxt("diauxic1.txt") time = didata1[:,0] pop2 = didata1[:,1] # model equations def diauxic_ode(x,t,params): r1,r2,k = params y,S1,S2 = x derivs = [r1*S1*y+(k/(k+S1))*r2*S2*y, r1*S1*y, (k/(k+S1))*r2*S2*y] return derivs # runs a simulation and returns the population size def diauxic_run(pars,t): r1,r2,k,y0,S10,S20 = pars ode_params=[r1,r2,k] ode_starts=[y0,S10,S20] out = odeint(diauxic_ode, ode_starts, t, args=(ode_params,)) return out[:,0] # residual function def diauxic_resid(pars,t,data): r1,r2,k,y0,S10,S20 = pars ode_params=[r1,r2,k] ode_starts=[y0,S10,S20] out = odeint(diauxic_ode, ode_starts, t, args=(ode_params,)) return diauxic_run(pars,t)data p0 =[1,1,1,1,1,1] lsq_out = leastsq(diauxic_resid, p0, args=(time,pop2)) plt.plot(time,pop2,'o',time,diauxic_resid(p0,time,lsq_out[0])) plt.show()

LevenbergMarquardt algorithm with GPU support
For the shallow neural net, the LM algorithm does amazingly well.
However, only MatLab and pyrenn (Python package) seems to have a robust implementation of it. A problem with both of these implementation is that they do not have GPU support for it. I also tried neupy (a python package) but it is not robust and fails when you try to train for longer epoch or large dataset. Do you guys know of a good LM python package for NN that can be trained using GPU?