SPGD V0.03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SPGD: Search Party Gradient Descent algorithm, a Simple Gradient-Based Parallel Algorithm for Bound-Constrained Optimization\n",
    "\n",
    "Link: https://www.mdpi.com/2227-7390/10/5/800"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This is a preliminary code written for the SPGD paper.\n",
    "The code is not efficient, and has lots of room for improvement.\n",
    "We plan to release an updated version of this code in future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "import threading\n",
    "import time\n",
    "import math, statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #F1Ackley \n",
    "    x,y = x[0],x[1]\n",
    "    '''\n",
    "    Ackley Function\n",
    "    wikipedia: https://en.wikipedia.org/wiki/Ackley_function\n",
    "    global minium at f(x=0, y=0) = 0\n",
    "    bounds: -35<=x,y<=35\n",
    "    '''\n",
    "    return (-20 * np.exp(-0.02 * np.sqrt(0.5 * (x*x + y*y))) -\n",
    "            np.exp(0.5 * (np.cos(2.0*np.pi*x) + np.cos(2*np.pi*y))) + np.e + 20) #Ackley\n",
    "mrnge = [-35,35]\n",
    "optimum = 0\n",
    "dim = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #F2beale\n",
    "    x,y = x[0],x[1]\n",
    "    '''\n",
    "    Beale Function\n",
    "    global minimum: f(x=3, y=0.5) = 0\n",
    "    bounds: -4.5 <= x, y <= 4.5\n",
    "    '''\n",
    "    return ((1.500 - x + x*y)**2 +\n",
    "            (2.250 - x + x*y**2)**2 +\n",
    "            (2.625 - x + x*y**3)**2)  #Beale\n",
    "mrnge = [-5,5]\n",
    "optimum = 0\n",
    "dim = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #F3Egg Holder\n",
    "    x,y = x[0],x[1]\n",
    "    '''\n",
    "    Eggholder Function\n",
    "    global minimum: f(x=512, y=404.2319) = -959.6407\n",
    "    bounds: -512 <= x, y <= 512\n",
    "    '''\n",
    "\n",
    "    return (-(y+47)*np.sin(np.sqrt(abs((x/2.0) + (y+47)))) -\n",
    "            x*np.sin(np.sqrt(abs(x-(y+47)))))\n",
    "mrnge = [-512,512]\n",
    "optimum = -959.6407\n",
    "dim = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(xy):#F4Gold-Stein\n",
    "    '''\n",
    "    Goldstein-Price Function\n",
    "    global minimum: f(0, -1) = 3\n",
    "    bounds: -2 <= x, y <= 2\n",
    "    '''\n",
    "    x, y = xy[0], xy[1]\n",
    "    return ((1 + (x + y + 1)**2 * (19 - 14*x + 3*x**2 - 14*y + 6*x*y + 3*y**2)) *\n",
    "            (30 + (2*x-3*y)**2 * (18 - 32*x + 12*x**2 + 48*y - 36*x*y + 27*y**2)))\n",
    "mrnge = [-2,2]\n",
    "optimum = 3\n",
    "dim = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x):#F5Matyas\n",
    "    x,y = x[0],x[1]\n",
    "    '''\n",
    "    Matyas Function\n",
    "    global minimum: f(x=0, y=0) = 0\n",
    "    bounds: -10 <= x, y <= 10\n",
    "    '''\n",
    "    return 0.26*(x**2 + y**2) - 0.48*x*y #matyas\n",
    "mrnge = [-10,10]\n",
    "optimum = 0\n",
    "dim = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #F6Sahaffer\n",
    "    x,y = x[0],x[1]\n",
    "    '''\n",
    "    Schaffer Function N.4\n",
    "    Reference: https://www.sfu.ca/~ssurjano/schaffer4.html\n",
    "    global minimum:\n",
    "        f(x=0, y= 1.25313) = 0.292579\n",
    "        f(x=0, y=-1.25313) = 0.292579\n",
    "    bounds: -100 <= x, y <= 100\n",
    "    '''\n",
    "    return 0.5 + (np.cos(np.sin(x**2 - y**2))**2 - 0.5)/(1+0.001*(x**2 + y**2))**2 #N.4\n",
    "mrnge = [-100,100]\n",
    "optimum = 0.292579\n",
    "dim = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(xy): #Colville\n",
    "    '''Colville Function\n",
    "    Parameters\n",
    "    ----------\n",
    "        xy : list\n",
    "    Returns\n",
    "    -------\n",
    "        float\n",
    "    Notes\n",
    "    -----\n",
    "    Bounds: The Colville function is a 4-dimensional function usually evaluated\n",
    "    on the hypercube defined by x_i in [-10, 10] for i=1,2,3,4.\n",
    "    Global minimum: f(x)=0 at x=[1,1,1,1]\n",
    "    References\n",
    "    ----------\n",
    "    A.-R. Hedar, “Global Optimization Test Problems”\n",
    "    '''\n",
    "\n",
    "    x1, x2, x3, x4 = xy[0], xy[1], xy[2], xy[3]\n",
    "    a = 100*(x1**2 - x2)**2\n",
    "    b = (x1-1)**2\n",
    "    c = (x3-1)**2\n",
    "    d = 90*(x3**2 - x4)**2\n",
    "    e = 10.1*((x2-1)**2 + (x4-1)**2)\n",
    "    f = 19.8*(x2-1)*(x4-1)\n",
    "    return a + b + c + d + e + f\n",
    "mrnge = [-10,10]\n",
    "optimum = 0\n",
    "dim = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(xy): #griewank\n",
    "    '''Griwank Function\n",
    "    Bounds: x_i in [-600, 600] for all i=1,...,d\n",
    "    Global minimum: f(x)=0 at x=(0,...,0)\n",
    "\n",
    "    '''\n",
    "    a, b, = 0, 1\n",
    "    for i, v in enumerate(xy):\n",
    "        a += v**2 / 4000.0\n",
    "        b *= np.cos(v/np.sqrt(i+1))\n",
    "    return a - b + 1\n",
    "mrnge = [-600,600]\n",
    "optimum = 0\n",
    "dim = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x, m=10): #Michaelwicz\n",
    "    '''Michalewicz Function\n",
    "    Notes\n",
    "    -----\n",
    "    The parameter m defines the steepness of they valleys and ridges; a larger m leads to a more difficult search. The recommended value of m is m = 10.\n",
    "    global minimum for 2D: f(x)=-1.8013 at x*=(2.20,1.57)\n",
    "    bounds: x_i in [0, π] for i=1,..., d\n",
    "    '''\n",
    "    d = len(x)\n",
    "    result = 0\n",
    "    for i in range(d):\n",
    "        result -= np.sin(x[i])*(np.sin((i+1)*x[i]**2/np.pi))**(2*m)\n",
    "    return result\n",
    "mrnge = [0,4]\n",
    "optimum = -4.687658\n",
    "dim = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #Rosenbrock\n",
    "    '''\n",
    "    Rosenbrock Function\n",
    "    wikipedia: https://en.wikipedia.org/wiki/Rosenbrock_function\n",
    "    global minimum:\n",
    "        n=2 -> f(1,1)=0\n",
    "    bounds:\n",
    "        -5 <= x_i <= +10\n",
    "        1 <= i <= n\n",
    "    '''\n",
    "    total = 0\n",
    "    for i in range(len(x)-1):\n",
    "        total += 100*(x[i+1] - x[i]*x[i])**2 + (1-x[i])**2\n",
    "    return total\n",
    "mrnge = [-30,30]\n",
    "optimum = 0\n",
    "dim = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #10D_zakharov\n",
    "    '''Zakharov Function\n",
    "    '''\n",
    "    a, b = 0, 0\n",
    "    for i, val in enumerate(x):\n",
    "        a += val**2\n",
    "        b += 0.5*i*val\n",
    "    return a + b**2 + b**4\n",
    "mrnge = [-5,10]\n",
    "optimum = 0\n",
    "dim = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(xy): #Rotated Ellipsoid\n",
    "    '''Rotated Hyper-Ellipsoid\n",
    "    Bounds: the rotated hyper-ellipsoid is usually evaluated on the hypercube\n",
    "    defined by x_i in [-65.536, 65.536] for all i=1,...,d\n",
    "    Global minimum: f(x)=0 at x=(0,...,0)\n",
    "    References\n",
    "    ----------\n",
    "    Molga, M., & Smutnicki, C. Test functions for optimization needs (2005)\n",
    "    '''\n",
    "    a = 0\n",
    "    for i in range(0, len(xy)): a += sum([xy[j]**2 for j in range(0, i)])\n",
    "    return a\n",
    "mrnge = [-66,66]\n",
    "optimum = 0\n",
    "dim = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(x): #Rastrigin\n",
    "    '''Rastrigin Function\n",
    "    Notes\n",
    "    -----\n",
    "    Bounds: -5.12 <= x_i <= 5.12 for all i=1,...,d\n",
    "    Global minimum: f(x)=0 at x=(0,...,0)\n",
    "    References\n",
    "    ----------\n",
    "    wikipedia: https://en.wikipedia.org/wiki/Rastrigin_function\n",
    "    '''\n",
    "    return len(x)*10.0 +  sum([item*item - 10.0*np.cos(2.0*np.pi*item) for item in x]) #Rastrigin\n",
    "mrnge = [-6,6]\n",
    "optimum = 0\n",
    "dim = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(xy): #Holder Table\n",
    "    '''\n",
    "    Holder Table Function\n",
    "    global minimum:\n",
    "        f(x= 8.05502, y= 9.66459) = -19.2085\n",
    "        f(x=-8.05502, y= 9.66459) = -19.2085\n",
    "        f(x= 8.05502, y=-9.66459) = -19.2085\n",
    "        f(x=-8.05502, y=-9.66459) = -19.2085\n",
    "    bounds: -10 <= x, y <= 10\n",
    "    '''\n",
    "    x, y = xy[0], xy[1]\n",
    "    return -abs(np.sin(x)*np.cos(y)*np.exp(abs(1-(np.sqrt(x**2 + y**2)/np.pi))))\n",
    "mrnge = [-10,10]\n",
    "optimum = -19.2085"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def f(xy): #tripod\n",
    "    '''Tripod Function\n",
    "    Parameters\n",
    "    ----------\n",
    "        xy : list\n",
    "    Returns\n",
    "    -------\n",
    "        float\n",
    "    Notes\n",
    "    -----\n",
    "    Bounds: x_i in [-100, 100] for i=1,2\n",
    "    Global minimum: f(x)=0 at x=[0,-50]\n",
    "    References\n",
    "    ----------\n",
    "    S. Rahnamyan, H. R. Tizhoosh, N. M. M. Salama, “A Novel Population\n",
    "    Initialization Method for Accelerating Evolutionary Algorithms”\n",
    "    Computers and Mathematics with Applications, vol. 53, no. 10,\n",
    "    pp. 1605-1614, 2007.\n",
    "    '''\n",
    "\n",
    "    x1, x2 = xy[0], xy[1]\n",
    "    p_x1 = 1 if x1>=0 else 0\n",
    "    p_x2 = 1 if x2>=0 else 0\n",
    "\n",
    "    a = p_x2 * (1 + p_x1)\n",
    "    b = abs(x1 + 50*p_x2*(1-2*p_x1))\n",
    "    c = abs(x2 + 50*(1-2*p_x2))\n",
    "    return a + b + c\n",
    "mrnge = [-100,100]\n",
    "optimum = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SPGD Begin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def derivative (arr, pos):\n",
    "    h = 0.000000001\n",
    "    temp = arr.copy()\n",
    "    temp[pos] = temp[pos] + h\n",
    "    num = f(temp) - f(arr)\n",
    "    return num/h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def printing (minimalist):\n",
    "    lst=[]\n",
    "    for i in range (len(minimalist)):\n",
    "        lst.append(f(minimalist[i]))\n",
    "    minima_arr.append(min (lst))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def descent(tid, w, a, e):\n",
    "    for i in range(e*Iter_P_Ep, e * Iter_P_Ep + Iter_P_Ep):\n",
    "        for pos in range(dim):\n",
    "            if (w[pos] < rnge[0] or w[pos] > rnge[1]):\n",
    "                w[pos] = random.randint(rnge[0],rnge[1])\n",
    "                \n",
    "        V[tid][i] = w\n",
    "        for pos in range(dim):\n",
    "            w[pos] = w[pos] - a * derivative(w,pos)\n",
    "            \n",
    "             \n",
    "        global CLM\n",
    "        if f( V[tid][i]) < fval[0]:\n",
    "            fval[0] = f( V[tid][i])\n",
    "            CLM = V[tid][i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def check_pre_convrg ():\n",
    "    pre_convg = []       \n",
    "    global Turn, skip\n",
    "    for i in range(ep-5, len(minimalist)):\n",
    "        pre_convg.append(f(minimalist[i]))\n",
    "\n",
    "    if(statistics.stdev(pre_convg) == 0 and Turn == 0):\n",
    "        lrn = np.linspace(0.1, 0.001, num=Tc)\n",
    "        Turn +=1\n",
    "        skip = 3\n",
    "\n",
    "    if(skip > 0):\n",
    "        skip = skip -1\n",
    "\n",
    "    elif(statistics.stdev(pre_convg) == 0 and Turn == 1 and skip ==0):\n",
    "        lrn = np.linspace(0.01, 0.0001, num=Tc)\n",
    "        Turn +=1\n",
    "        skip = 3\n",
    "\n",
    "    elif(statistics.stdev(pre_convg) == 0 and Turn == 2 and skip == 0):\n",
    "        lrn = np.linspace(0.001,0.000001, num = Tc)\n",
    "        Turn +=1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "code_folding": [
     0,
     1
    ]
   },
   "outputs": [],
   "source": [
    "def reduce_search_circle ():\n",
    "    if(ep % 5 == 0 and ep !=0):\n",
    "        dummy = 5\n",
    "        tempmin, tempmax=[],[]\n",
    "        for i in range(ep-dummy, len(minimalist)):\n",
    "            tempmin.append(math.floor(min(minimalist[i])))\n",
    "            tempmax.append(math.ceil(max(minimalist[i])))\n",
    "        rnge[0]= min(tempmin)\n",
    "        rnge[1] = max(tempmax)\n",
    "        if (rnge[0] == rnge[1]):\n",
    "            rnge[0],rnge[1] = mrnge[0],mrnge[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "code_folding": [
     3
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run:  0\n",
      "0 20 0.0\n",
      "Run:  1\n",
      "1 16 0.0\n",
      "Finished in 0.67 second(s)\n"
     ]
    }
   ],
   "source": [
    "minima_arr = []\n",
    "avg_epp = []\n",
    "start = time.perf_counter()\n",
    "for run in range(50):    #Running SPGD 50 times for statistcal analysis\n",
    "    print(\"Run: \", run)\n",
    "    \n",
    "    #Begin Intialization of various variables\n",
    "    rnge = [mrnge[0],mrnge[1]]\n",
    "    fval, CLM = [9999999],np.zeros(dim)\n",
    "    Turn, skip = 0,-1\n",
    "    Episodes, Iter_P_Ep = 50, 20\n",
    "    LocMinima,minimalist, epsilon, Tc =[] ,[], 1, 25 #Threadcount\n",
    "    decay = epsilon/Episodes\n",
    "\n",
    "    V,T = np.zeros((Tc, Episodes * Iter_P_Ep, dim)), [0] * Tc\n",
    "    lrn = np.linspace(1, 0.001, num=Tc)\n",
    "    init_xy = np.random.uniform(rnge[0],rnge[1], size=(Tc, dim))\n",
    "    \n",
    "    #End Intialization\n",
    "    \n",
    "    #Begin SPGD Execution\n",
    "    for ep in range(Episodes):\n",
    "        \n",
    "        for i in range(Tc):\n",
    "            T[i] = threading.Thread(target=descent, args=(i, init_xy[i],lrn[i], ep))    \n",
    "            T[i].start()  \n",
    "        for i in range(Tc):\n",
    "            T[i].join()\n",
    "\n",
    "        Edge_Loc = ep*Iter_P_Ep + (Iter_P_Ep-1)\n",
    "\n",
    "        minimalist.append(CLM) #Remember Local Minima Encountered\n",
    "        \n",
    "        reduce_search_circle()        \n",
    "        if(ep > 5):\n",
    "            check_pre_convrg()\n",
    "            \n",
    "        #Exploitation Vs Exploration\n",
    "        if(random.uniform(0,1) >= epsilon): #Exploitation\n",
    "            for Tid in range(Tc):\n",
    "                if (random.uniform(0,1) >= 0.2):\n",
    "                    init_xy[Tid] = CLM #Assemble at best known location\n",
    "                else: #Chooses to explore 20% of the time, to avoid loc minima.\n",
    "                    for pos in range(dim):\n",
    "                        init_xy[Tid][pos] = np.random.triangular(rnge[0],CLM[pos],rnge[1])\n",
    "\n",
    "        else:\n",
    "            #EXPLORATION\n",
    "            for Tid in range (Tc):\n",
    "                #Continue as and where you are\n",
    "                if (random.uniform(0,1) >= 0.5): #Dont Change\n",
    "                    init_xy[Tid] = V[Tid][Edge_Loc]\n",
    "                else: #Pure Exploration\n",
    "                    for pos in range(dim):\n",
    "                        init_xy[Tid][pos] = np.random.triangular(rnge[0],CLM[pos],rnge[1])\n",
    "                        \n",
    "                        \n",
    "        if(ep >10 ):\n",
    "            convergence = []        \n",
    "            for i in range(ep-10, len(minimalist)):\n",
    "                convergence.append(f(minimalist[i]))\n",
    "            if(statistics.stdev(convergence) == 0):\n",
    "                print(run, ep, f(CLM))\n",
    "                avg_epp.append(ep)\n",
    "                break\n",
    "        epsilon = epsilon - decay\n",
    "    printing(minimalist)\n",
    "finish = time.perf_counter()\n",
    "print(f'Finished in {round(finish-start, 2)} second(s)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.063286299237916e-08 2.0419683007588674e-05 1.0255157935290526e-05 1.4374809212325405e-05\n",
      "0.02\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "49.44"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Summary\n",
    "import statistics\n",
    "print(min(minima_arr), max(minima_arr),statistics.mean(minima_arr) ,statistics.stdev(minima_arr))\n",
    "j=0\n",
    "for i in range(len(minima_arr)):\n",
    "    if abs(minima_arr[i] - optimum )< 0.000001 :\n",
    "        #print(j,\"success\", minima_arr[i])\n",
    "        j+=1\n",
    "print(j/50)\n",
    "dum = avg_epp\n",
    "for _ in range(50 - len(avg_epp)):\n",
    "    dum.append(50)\n",
    "statistics.mean(dum) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}