<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Optimization of Non-Differentiable Functions Using Differential Evolution with Constraints
This code is provided as supplementary material of the lecture Machine Learning and Optimization in Communications (MLOC).<br>
This code illustrates
Step1:  For illustration, we optimize Rosenbrock's function in 2 dimensions. The function is given by
\begin{equation}
f(\boldsymbol{x}) = (1-x_2)^2 + 100(x_2-x_1^2)^2
\end{equation}
such that $x_1 \geq 0$ and $x_2 \geq 0$.
Step2:  Plot the function as a heat map.
Step3:  Helper function to generate a random initial population of $D$ dimensions $N_P$ elements. The elementsof the population are randomly distributed in the interval $[x_{\min},x_{\max}]$ (in every dimension).
Step4:  Carry out differential evolution similar (not identical, slightly modified) to the scheme DE1 described in [1]
[1] R. Storn and K. Price, "Differential Evolution - A simple and efficient adaptive scheme for global optimization
over continuous spaces", Technical Report TR-95-012, March 1995
Step5:  Generate animation.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker
Explanation: Optimization of Non-Differentiable Functions Using Differential Evolution with Constraints
This code is provided as supplementary material of the lecture Machine Learning and Optimization in Communications (MLOC).<br>
This code illustrates:
* Use of differential evolution to optimize Rosenbrock's function for non-negative inputs
* Modification of cost function to take into account constraints
End of explanation
# Rosenbrocks function
def rosenbrock(x,y):
    return (1-x)**2+100*(y-(x**2))**2
    
    
def fun(x,y):
    value = rosenbrock(x,y)
    
    # CONSTRAINT: x and y should be non-negative. Hence, add a penalty for negative numbers that becomes larger the more negative the number is
    if x < 0:
        value += 100*(-x)
    if y < 0:
        value += 100*(-y)
    return value
# vector-version of the function
vrosenbrock = np.vectorize(rosenbrock)
vfun = np.vectorize(fun)
Explanation: For illustration, we optimize Rosenbrock's function in 2 dimensions. The function is given by
\begin{equation}
f(\boldsymbol{x}) = (1-x_2)^2 + 100(x_2-x_1^2)^2
\end{equation}
such that $x_1 \geq 0$ and $x_2 \geq 0$.
End of explanation
# plot map of Griewanks function
x = np.arange(-2.0, 5.1, 0.1)
y = np.arange(-2.0, 5.1, 0.1)
X, Y = np.meshgrid(x, y)
fZ = vrosenbrock(X,Y)
plt.figure(1,figsize=(10,9))
plt.rcParams.update({'font.size': 14})
temp = np.log(fZ)
temp[temp < -1] = 0
plt.contourf(X,Y,temp,levels=20)
plt.fill_betweenx([-2, 5],[-2,-2],[0,0], alpha=0.5, color='white',linewidth=0)
plt.fill_betweenx([-2, 0],[0,0],[5,5], alpha=0.5, color='white',linewidth=0)
#plt.axhspan(-5, 0, alpha=0.5, color='white')
plt.axis('scaled')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.show()
Explanation: Plot the function as a heat map.
End of explanation
def initial_population(D, NP, xmin, xmax):
    v = np.random.rand(NP, D)*(xmax - xmin) + xmin
    return v
Explanation: Helper function to generate a random initial population of $D$ dimensions $N_P$ elements. The elementsof the population are randomly distributed in the interval $[x_{\min},x_{\max}]$ (in every dimension).
End of explanation
#dimension
D = 2
# population
NP  = 15*D
# twiddling parameter
F = 0.8
# cross-over probability
CR = 0.3
# maximum 1000 iterations
max_iter = 1000
# generate initial population
population = initial_population(D, NP, -2, 5)[:]
# compute initial cost
cost = vfun(population[:,0], population[:,1])
best_index = np.argmin(cost)
best_cost = cost[best_index]
iteration = 0
# keep track of population 
save_population = []
while iteration < max_iter:
    # loop over every element from the population
    for k in range(NP):
        # get 4 random elements
        rp = np.random.permutation(NP)[0:4]
        
        # remove ourselves from the list
        rp = [j for j in rp if j !=  k]
        
        # generate new candidate vector
        v = population[rp[0],:] + F*( population[rp[1],:] - population[rp[2],:] ) 
        
        # take vector from population
        u = np.array(population[k,:])
        # cross-over each coordinate with probability CR with entry from candidate vector v
        idx = np.random.rand(D) < CR
        # cross-over
        u[idx] = v[idx]
        new_cost = fun(u[0], u[1])
        if new_cost < cost[k]:
            # better cost? keep!
            cost[k] = new_cost
            population[k,:] = u
            
            
            if new_cost < best_cost:
                best_cost = new_cost
                best_index = k
    
    save_population.append(np.array(population[:]))
    iteration += 1
    if iteration % 100 == 0:
        print('After iteration %d, best cost %1.4f (obtained for (%1.2f,%1.2f))' % (iteration, best_cost, population[best_index,0], population[best_index,1]))
                
plt.figure(1,figsize=(9,9))
plt.rcParams.update({'font.size': 14})
plt.contourf(X,Y,temp,levels=20)
index = 180
cost = vfun(save_population[index][:,0], save_population[index][:,1])
best_index = np.argmin(cost)
    
plt.fill_betweenx([-2, 5],[-2,-2],[0,0], alpha=0.5, color='white',linewidth=0)
plt.fill_betweenx([-2, 0],[0,0],[5,5], alpha=0.5, color='white',linewidth=0)
    
plt.scatter(save_population[index][:,0], save_population[index][:,1], c='w')
plt.scatter(save_population[index][best_index,0], save_population[index][best_index,1], c='r')
    
plt.xlim((-2,5))
plt.ylim((-2,5))
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.savefig('DE_Rosenbrock.pdf',bbox_inches='tight')
Explanation: Carry out differential evolution similar (not identical, slightly modified) to the scheme DE1 described in [1]
[1] R. Storn and K. Price, "Differential Evolution - A simple and efficient adaptive scheme for global optimization
over continuous spaces", Technical Report TR-95-012, March 1995
End of explanation
%matplotlib notebook
# Generate animation
from matplotlib import animation, rc
from matplotlib.animation import PillowWriter # Disable if you don't want to save any GIFs.
font = {'size'   : 18}
plt.rc('font', **font)
fig, ax = plt.subplots(1, figsize=(10,10))
ax.set_xlim(( -2, 5))
ax.set_ylim(( -2, 5))
ax.axis('scaled')
written = False
def animate(i):
    ax.clear()
    ax.contourf(X,Y,temp,levels=20)
        
    cost = vfun(save_population[i][:,0], save_population[i][:,1])
    best_index = np.argmin(cost)
    
    ax.fill_betweenx([-2, 5],[-2,-2],[0,0], alpha=0.5, color='white',linewidth=0)
    ax.fill_betweenx([-2, 0],[0,0],[5,5], alpha=0.5, color='white',linewidth=0)
    
    ax.scatter(save_population[i][:,0], save_population[i][:,1], c='w')
    ax.scatter(save_population[i][best_index,0], save_population[i][best_index,1], c='r')
    ax.set_xlabel(r'$x_1$',fontsize=18)
    ax.set_ylabel(r'$x_2$',fontsize=18)
    ax.set_xlim(( -2, 5))
    ax.set_ylim(( -2, 5))
    
anim = animation.FuncAnimation(fig, animate, frames=300, interval=80, blit=False)
fig.show()
anim.save('differential_evolution_Rosenbrock.gif', writer=PillowWriter(fps=7))
Explanation: Generate animation.
End of explanation
<END_TASK>