Bayesian_Optmization_V2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modification adfter meeting \n",
    "\n",
    "1. delete theta, phi, mu, only keep U, g, r. In addition, add error range for flow rate (fr)\n",
    "2. rewrite class function for more general application\n",
    "3. less initial datapoints, increase iterations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import minimize \n",
    "import sys\n",
    "\n",
    "if not sys.warnoptions:\n",
    "    import warnings\n",
    "    warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# newtonian list add theta, phi, mu\n",
    "property_list_newtonian = [30,1,0.045,0.066,1210]\n",
    "# non_newtonian property_list = [0.53,0.25,0.065,1143.5]\n",
    "x_f = 0.0005\n",
    "x_d = 0.101\n",
    "h_w = 0.1\n",
    "quantile = 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(x_input,ink_type, ink_value,judge = x_f, qt = quantile,Xd = x_d):\n",
    "    \"\"\"\n",
    "    objective function to max rewards \n",
    "    \n",
    "    design this function with quantile weight, inspired by quantle loss regression idea and rewards matrix used in the RL notebook\n",
    "\n",
    "    for slot die coating, two edge cases are: leaking, right: breaking-up\n",
    "\n",
    "    leakings are zero-tolerance defects. goal is to encourage more prediction points on right side\n",
    "\n",
    "    R(x_u,x_f) = (x_f - x_u)*tao if x_u - half point > 0;\n",
    "\n",
    "               = (x_u - 0)*(1-tao) if x_u - half point < 0 \n",
    "    \n",
    "    interpretation : if tao = 0.1, then R will be larger for points located on left side, which potentially increase chance to make x_u on left part of middle point\n",
    "    goal is to max rewards \n",
    "    \"\"\"\n",
    "    # get output from input, input is a matrix with 6 columns\n",
    "    # get ink properties:\n",
    "    mu,n,m,sigma,density = ink_value\n",
    "    y = []\n",
    "    x_location = []\n",
    "    if ink_type == \"Non_Newtonian\":\n",
    "        for x_i in x_input: \n",
    "        # searching values and doing calculations\n",
    "            U, g, r, fr = x_i  \n",
    "            # delete theta, phi, mu. only keep U, g, r. In addition, add error range for fr\n",
    "            thickness = fr/(U*h_w) \n",
    "\n",
    "            a = (U*(n+1)*(2*n+1)/n)**n*(g**(-n-1))\n",
    "            b_s = (U*(g-2*thickness)*(n+1)*(2*n+1)/n)**n * (g**(-2*n-1))\n",
    "            b_l = -(U*(2*thickness-g)*(n+1)*(2*n+1)/n)**n * (g**(-2*n-1))\n",
    "            modified_CA = ((U*m*(U/g)**(n-1))/sigma)**(2/3)      \n",
    "            p_ambient = density*U**2/2\n",
    "\n",
    "            if thickness <= g/2:\n",
    "                xu_estimate = judge - (p_ambient - 1.34*modified_CA*sigma/thickness - (Xd-judge)*m*b_s - sigma*(math.cos(theta)+math.cos(phi)))/a*m\n",
    "\n",
    "            elif thickness > g/2:\n",
    "                xu_estimate = judge - (p_ambient - 1.34*modified_CA*sigma/thickness - (Xd-judge)*m*b_l - sigma*(math.cos(theta)+math.cos(phi)))/a*m\n",
    "            \n",
    "            x_location.append(xu_estimate)\n",
    "            if xu_estimate <= 0:\n",
    "                y.append(xu_estimate*2)\n",
    "            elif xu_estimate >= judge/2 and xu_estimate < judge:\n",
    "                y.append((judge - xu_estimate)*qt)\n",
    "            elif xu_estimate < judge/2 and xu_estimate > 0:\n",
    "                y.append(xu_estimate*(1-qt))\n",
    "            else:\n",
    "                y.append(-xu_estimate*2)\n",
    "        \n",
    "    elif ink_type == \"Newtonian\":\n",
    "        for x_i in x_input: \n",
    "        # searching values and doing calculations\n",
    "            U, g, r, fr = x_i  \n",
    "            # delete theta, phi, mu. only keep U g r. In addition, add error range for fr\n",
    "            thickness = fr/(U*h_w) \n",
    "\n",
    "            a = (U*(n+1)*(2*n+1)/n)**n*(g**(-n-1))\n",
    "            b_s = (U*(g-2*thickness)*(n+1)*(2*n+1)/n)**n * (g**(-2*n-1))\n",
    "            b_l = -(U*(2*thickness-g)*(n+1)*(2*n+1)/n)**n * (g**(-2*n-1))\n",
    "            modified_CA = ((U*m*(U/g)**(n-1))/sigma)**(2/3)      \n",
    "            p_ambient = density*U**2/2\n",
    "\n",
    "            if thickness <= g/2:\n",
    "                xu_estimate = judge - g**2*(p_ambient - 1.34*modified_CA*sigma/thickness - (Xd-judge)*m*b_s - sigma*(math.cos(theta)+math.cos(phi))/g)/(6*mu*U)\n",
    "                     \n",
    "            elif thickness > g/2:\n",
    "                xu_estimate = judge - g**2*(p_ambient - 1.34*modified_CA*sigma/thickness - (Xd-judge)*m*b_l - sigma*(math.cos(theta)+math.cos(phi))/g)/(6*mu*U)\n",
    "            \n",
    "            x_location.append(xu_estimate)\n",
    "            if xu_estimate <= 0:\n",
    "                y.append(xu_estimate*2)\n",
    "            elif xu_estimate >= judge/2 and xu_estimate < judge:\n",
    "                y.append((judge - xu_estimate)*qt)\n",
    "            elif xu_estimate < judge/2 and xu_estimate > 0:\n",
    "                y.append(xu_estimate*(1-qt))\n",
    "            else:\n",
    "                y.append(-xu_estimate*2)\n",
    "    return  y  #,x_location uncommend this if you want to check x_u location\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_data(size = 100):\n",
    "#     theta = np.random.choice(np.arange(120,140),self.batch_size)  # 10 扩大范围 看看有没有奇怪的suggestion\n",
    "#     phi = np.random.choice(np.arange(60,80),self.batch_size)  # 100 需要有一个相互依靠的关系\n",
    "    # U, g, r, fr = x_i\n",
    "    U = np.random.uniform(low = 0.05, high = 0.5, size = size)\n",
    "    g = np.random.choice(np.arange(250,500)*1e-6, size = size)\n",
    "    r = np.random.uniform(low = 1.5, high = 4, size = size)\n",
    "    fr = np.random.uniform(low = 0.05, high = 0.2, size = size)\n",
    "    input_matrix = np.stack([U, g, r, fr], axis=1)\n",
    "    return input_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Bayesian_opt():\n",
    "    def __init__(self, target_func,ink_type, ink_value, x_init, y_init, n_iter, batch_size, input_list):\n",
    "        # initial values\n",
    "        self.x_init = x_init\n",
    "        self.y_init = y_init\n",
    "        self.n_iter = n_iter\n",
    "        self.target_func = target_func\n",
    "        self.batch_size = batch_size  \n",
    "        self.input_list = input_list\n",
    "        self.ink_type = ink_type\n",
    "        self.ink_value = ink_value\n",
    "        self.best_samples = pd.DataFrame(columns = (self.input_list + [\"y\",\"ei\"]))\n",
    "        self.gauss_pr = GaussianProcessRegressor()\n",
    "        self.distances_ = []\n",
    "    # surrogate function using EI method, as a starting point  \n",
    "    def _get_expected_improvement(self, x_new): \n",
    "        # do more research on different methods and figure out why choose EI\n",
    "        # Using estimate from Gaussian surrogate instead of actual function for a new trial data point to avoid cost \n",
    "        mean_y_new, sigma_y_new = self.gauss_pr.predict(np.array([x_new]), return_std=True)\n",
    "        sigma_y_new = sigma_y_new.reshape(-1,1)\n",
    "        if sigma_y_new < 0.0:\n",
    "            return 0.0\n",
    "        \n",
    "        # Using estimates from Gaussian surrogate instead of actual function for entire prior distribution to avoid cost\n",
    "        # add a very small to avoid divide by zero sigma \n",
    "        mean_y = self.gauss_pr.predict(self.x_init)\n",
    "        max_mean_y = np.max(mean_y)\n",
    "        z = (mean_y_new - max_mean_y) / (sigma_y_new+1e-9)\n",
    "        exp_imp = (mean_y_new - max_mean_y) * norm.cdf(z) + sigma_y_new * norm.pdf(z)\n",
    "        \n",
    "        return exp_imp\n",
    "    \n",
    "    def acquisition_function(self,x):\n",
    "        # max (function) = min (-function)\n",
    "        return -self._get_expected_improvement(x)\n",
    "    \n",
    "    ###### modify this\n",
    "    def random_select_inputs(self):\n",
    "#         theta = np.random.choice(np.arange(120,140),self.batch_size)  # 10 扩大范围 看看有没有奇怪的suggestion\n",
    "#         phi = np.random.choice(np.arange(60,80),self.batch_size)  # 100 需要有一个相互依靠的关系\n",
    "        # start with less points and larger iterations\n",
    "        U = np.random.uniform(low = 0.05, high = 1.5, size = self.batch_size)\n",
    "        g = np.random.choice(np.arange(250,500)*1e-6, size = self.batch_size)\n",
    "        r = np.random.uniform(low = 1.5, high = 4, size = self.batch_size)\n",
    "        fr = np.random.uniform(low = 0.05, high = 1, size = self.batch_size)\n",
    "        random_matrix = np.stack([U, g, r, fr], axis=1)\n",
    "        return random_matrix\n",
    "\n",
    "    def _get_next_probable_point(self):\n",
    "        \n",
    "        min_ei = float(sys.maxsize)\n",
    "        x_optimal = None \n",
    "        ###### fix bounds \n",
    "        bounds = [(1e-6, 1.6), (1e-6, 500*1e-6), (1e-6, 4.1),(1e-6, 1.1)]\n",
    "        # Trial with an array of input_simulated_data\n",
    "        random_matrix = self.random_select_inputs()\n",
    "        for x_start in random_matrix:\n",
    "            response = minimize(fun=self.acquisition_function, x0=x_start, method='L-BFGS-B',bounds = bounds)\n",
    "            if response.fun < min_ei:\n",
    "                min_ei = response.fun\n",
    "                x_optimal = response.x\n",
    "\n",
    "        return x_optimal, min_ei  \n",
    "        \n",
    "    # wait for debug for rest of two functions     \n",
    "    def _extend_prior_with_posterior_data(self, x,y):\n",
    "        \n",
    "        self.x_init = np.append(self.x_init, np.array([x]), axis = 0)\n",
    "        self.y_init = np.append(self.y_init, np.array([y]), axis = 0)\n",
    "  \n",
    "    def optimize(self):\n",
    "        y_max_ind = np.argmax(self.y_init)\n",
    "        y_max = self.y_init[y_max_ind]\n",
    "        optimal_x = self.x_init[y_max_ind]\n",
    "        optimal_ei = None\n",
    "        for i in range(self.n_iter):\n",
    "            self.gauss_pr.fit(self.x_init, self.y_init)\n",
    "            x_next, ei = self._get_next_probable_point()\n",
    "            y_next = self.target_func(np.array([x_next]),self.ink_type,self.ink_value)[0]\n",
    "            self._extend_prior_with_posterior_data(x_next,y_next)\n",
    "            # \n",
    "            if y_next > y_max:\n",
    "                y_max = y_next\n",
    "                optimal_x = x_next\n",
    "                optimal_ei = ei\n",
    "\n",
    "            if i == 0:\n",
    "                 prev_x = x_next\n",
    "            else:\n",
    "                self.distances_.append(np.linalg.norm(prev_x - x_next))\n",
    "                prev_x = x_next\n",
    "            \n",
    "            self.best_samples = self.best_samples.append({\"Web_speed(m/s)\":optimal_x[0],\"coating_height(m)\":optimal_x[1],\"ratio\":optimal_x[2],\n",
    "                                                           \"Flow_rate(m/s)\":optimal_x[3],\"y\": y_max, \"ei\": optimal_ei},ignore_index=True)\n",
    "        \n",
    "        return optimal_x, y_max \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.random.choice(np.arange(10,140),size =1)\n",
    "phi = np.random.choice(np.arange(60,100),size =1)\n",
    "init_X = initial_data()\n",
    "init_Y = objective(init_X,ink_type = \"Newtonian\", ink_value = property_list_newtonian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.5552154610249296"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(init_Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[107] [62]\n"
     ]
    }
   ],
   "source": [
    "print(theta,phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([2.53573884e-02, 5.00000000e-04, 1.00000000e-06, 1.00000000e-06]),\n",
       " 3.7137806739023966e-05)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# bayesian for newtonian\n",
    "bopt_newtonian = Bayesian_opt(target_func=objective, \n",
    "                              ink_type = \"Newtonian\",\n",
    "                              ink_value = property_list_newtonian,\n",
    "                              x_init=init_X,\n",
    "                              y_init=init_Y, \n",
    "                              n_iter=300,\n",
    "                              batch_size=30,\n",
    "                              input_list =[\"Web_speed(m/s)\",\"coating_height(um)\",\"ratio\",\"Flow_rate(m/s)\"])\n",
    "bopt_newtonian.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([3.7137806563719154e-05], [0.0004257243868725617])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "objective(x_input = [[2.53573884e-02, 5.00000000e-04, 1.00000000e-06, 1.00000000e-06]],\n",
    "          ink_type = \"Newtonian\", ink_value =property_list_newtonian)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}