<jupyter_start><jupyter_text>MTH 9842 HW1 Group L: <NAME>, <NAME>, <NAME>, <NAME><jupyter_code>import numpy as np
from scipy.optimize import least_squares,rosen, differential_evolution
import pandas as pd
from math import *
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Problem 1<jupyter_code>data=pd.read_csv('Assignment1_data.csv',header=None,names=['T','value'])
data<jupyter_output><empty_output><jupyter_text>Define function $g(x, T)=x_1+x_2\exp(x_3 T)+x_4\exp(x_5 T), x=(x_1,x_2,x_3,x_4,x_5) \in \mathbb{R}^5$<jupyter_code>def g(x, T):
    '''
    input: x, of shape (5,); independent variables
           T, vector
    return a vector of len(T)
    the code relies on numpy broadcasting across vector T
    '''
    assert x.shape==(5,)
    result=x[0]+x[1]*np.exp(x[2]*T)+x[3]*np.exp(x[4]*T)
    assert result.shape==T.shape
    return result
# rosenbrock function, not used, for testing purpose
def fun_rosenbrock_2(x):
    return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
# computes vector of residuals
def g_residuals(x, *args, **kwargs):
    # assert x.shape==(5,)
    T=data['T']
    v=data['value']
    resid=np.zeros(len(T))
    resid=g(x, T)-v
    return resid
# solve unconstrained non-linear least square problem using Levenberg-Marquardt algorithm

# use all zeros as initial guess of independent variables
x_ini=np.zeros(5)

res_1=least_squares(g_residuals, x_ini, method='lm')
# reason for algorithm termination
print("Success", res_1.success);
print(res_1.message)
# solutions found
res_1.x
# value of the cost function at the solution
res_1.cost
# vector of residuals at the solution
res_1.fun
# number of function evaluations done
res_1.nfev<jupyter_output><empty_output><jupyter_text>Problem 2<jupyter_code>def rosen_der(x):
    '''
    computes gradient of rosenbrock function
    input: x(2,)
    output: array of length 2 representing the gradient
    '''
    return np.array([2*(x[0]-1)-400*(x[1]-x[0]**2)*x[0],
                     200*(x[1]-x[0]**2)])

def rosenbrock(x):
    '''
    rosenbrock function
    input: x(2,)
    '''
    return (x[0]-1)**2+100*(x[1]-x[0]**2)**2

# keep stats during descent
l_x=[]  # list of x
l_cmin=[]  # list function values
l_resid=[]  # list of residuals
l_ba=[]  # list of learning rates used
l_der=[]  # list of gradients

# step size/learning rate to try
alpha=[0.001,0.002,0.003,0.004,0.005,0.1,0.5]

def rosenbrock_opt(x,alpha,tol1=1e-12,tol2=1e-6):
    '''
    implement gradient descent algorithm on rosenbrock function
    input: x(2,) initial guesses
           alpha() list of learning rates, the algorithm will choose one by means of a line search
    output: 
    '''
    it=0  # number of iteration
    #global l_x, l_cmin, l_resid, l_ba, l_der
    l_x.clear()
    l_cmin.clear()
    l_resid.clear()
    l_ba.clear()
    l_der.clear()
    while True:
        it+=1
        crosen=rosenbrock(x)
        der=rosen_der(x)
        # select alpha
        ba=alpha[0]
        cmin=rosenbrock(x-der*ba)
        for i in range(1,len(alpha)):
            val=rosenbrock(x-der*alpha[i])
            if val<cmin:
                ba=alpha[i]
                cmin=val
        # keep trails of the algorithm
        l_der.append(der)
        l_ba.append(ba)
        l_cmin.append(cmin)
        l_x.append(x)
        # update x estimates
        x=x-der*ba
        if abs(der).all()<tol1 and abs(cmin-crosen)<tol2:
        #if abs(cmin-crosen)<tol2:
            return x,cmin,'tol exit',it
        elif it>10000:
            return x,cmin,'iter exit',it
# gradient descent with initial guess (0,0)
rosenbrock_opt(np.array([0.,0.]),alpha,tol1=1,tol2=1e-12)
# plot function value during each iteration
obj=pd.Series(l_cmin)
obj.plot()
# check that the learning rates used are reasonable
np.histogram(l_ba,bins=alpha)
# initial guess (1.2, 1.2)
rosenbrock_opt(np.array([1.2,1.2]),alpha,tol2=1e-8)
# plot function value during each iteration
obj=pd.Series(l_cmin)
obj.plot()
# run gradient descent with initial guess x,y vary from -3 to +3.75
xcov=[]
ycov=[]
xexp=[]
yexp=[]
for i in range(25):
    for j in range(25):
        x=-3+0.3*i
        y=-3+0.3*j
        array,funval,status,_=rosenbrock_opt(np.array([x,y]),alpha)
        if status=='tol exit':
            xcov.append(x)
            ycov.append(y)
        elif status=='iter exit':
            xexp.append(x)
            yexp.append(y)
            
plt.scatter(xexp,yexp)
plt.scatter(xcov,ycov,facecolor='red')<jupyter_output><empty_output><jupyter_text>Conclusion: gradient descent algorithm applied to Rosenbrock function generally suffers from slow convergence once it hits the valley. Finding the valley is trivial, but once inside the valley finding global minimum is very slow.  Problem 3<jupyter_code>def ackley(x):
    res=20+exp(1)-20*exp(-0.2*sqrt(0.5*(x[0]**2+x[1]**2)))-exp(0.5*(cos(2*pi*x[0])+cos(2*pi*x[1])))
    return res
bounds1 = [(-1,3)]*2
differential_evolution(rosenbrock,bounds=bounds1,strategy='rand1bin',popsize=10)
bounds2 = [(-3,3)]*2
differential_evolution(ackley,bounds=bounds2,strategy='rand1bin')
Flist=[0.1+0.2*i for i in range(10)]
CRlist=[0.1*i for i in range(1,10)]
#rosen success x
rsx=[]
rsy=[]
#rosen unsuccess x
rux=[]
ruy=[]
asx=[]
asy=[]
aux=[]
auy=[]
for F in Flist:
    for CR in CRlist:
        rres=differential_evolution(rosenbrock,bounds=bounds1,strategy='rand1bin',popsize=10,mutation=F,recombination=CR)
        ares=differential_evolution(ackley,bounds=bounds2,strategy='rand1bin',popsize=10,mutation=F,recombination=CR)
        if rres['success']:
            rsx.append(F)
            rsy.append(CR)
        else:
            rux.append(F)
            ruy.append(CR)
        if ares['success']:
            asx.append(F)
            asy.append(CR)
        else:
            aux.append(F)
            auy.append(CR)
plt.scatter(rsx,rsy,c='b')
plt.scatter(rux,ruy,c='r')
plt.xlabel('F')
plt.ylabel('CR')
plt.scatter(asx,asy,c='b')
plt.scatter(aux,auy,c='r')
plt.xlabel('F')
plt.ylabel('CR')
#rosen success pop
rsp=[]
rup=[]
asp=[]
aup=[]
for n in range(4,22):
    rres=differential_evolution(rosenbrock,bounds=bounds1,strategy='rand1bin',popsize=n,mutation=0.5,recombination=0.1)
    ares=differential_evolution(ackley,bounds=bounds2,strategy='rand1bin',popsize=n,mutation=0.5,recombination=0.1)
    if rres['success']:
            rsp.append(n)
    else:
            rup.append(n)
    if ares['success']:
            asp.append(n)
    else:
            aup.append(n)
rsp
asp<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1495312740126018, 'fraction_numerical': 0.033656062701705854, 'mean_word_length': 3.5194444444444444, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '300939', 'license': 'MIT', 'n_tokens_mistral': 2591, 'n_tokens_neox': 2486, 'n_words': 508, 'path': 'HW1/9842HW1_Group L.ipynb', 'repo_name': 'jjusuf/MTH-9842', 'size': '6528'}