symbolic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After many trials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "double_skin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9414781941701408       Mean Squared Error: 0.3618988795335194          2.26412100000000\n",
      "0.9592118867490954       Mean Squared Error: 0.184067494019574          eta_o + 1.4976299\n",
      "0.9800631724197708       Mean Squared Error: 0.11725973700607613          eta_i + eta_o + 1.1657709\n",
      "0.9817337585727615       Mean Squared Error: 0.1055603904824611          eta_i + 1.16053319072701*eta_o + 1.0398949\n",
      "0.9960951534890364       Mean Squared Error: 0.09338918587619967          1.0973018**(eta_o/Di_D) + eta_i + eta_o\n",
      "1.0042132692411005       Mean Squared Error: 0.08964808077071601          eta_i + eta_o + (0.28830558 + eta_o/Di_D)**0.29875413\n",
      "1.0063315779300441       Mean Squared Error: 0.08371387536429944          eta_i + eta_o + 0.775073980830572*(eta_o/(Di_D*fc_fyo))**0.19528726\n",
      "1.0044671728665633       Mean Squared Error: 0.0821596621692013          eta_i + eta_o + (eta_o/(Di_D*fc_fyo*(L_Do + 0.8891143)))**0.1924652\n",
      "0.9946438717091605       Mean Squared Error: 0.0777503146339095          eta_i + eta_o + 2.46042508331603*(eta_o/(Di_D*fc_fyo*fyo))**0.20827866\n",
      "0.9962615963357033       Mean Squared Error: 0.07650050585523047          eta_i + eta_o + 2.87118808723674*(0.018589601584206*fc_fyo + eta_o/(Di_D*fc_fyo*fyo))**0.26466385\n",
      "0.9957780489917618       Mean Squared Error: 0.07675880768188613          eta_i + eta_o + ((eta_o + fc_fyo)/(fc_fyo*(0.0191056168340136*Di_D*fyo + 0.317549512593324)))**0.2788979\n",
      "1.0026741164148825       Mean Squared Error: 0.07415670654312659          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**2 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "0.9934642987567288       Mean Squared Error: 0.07382578948545589          eta_i + eta_o + 1.49779743724763*(0.26038945*Di_D**2 - (eta_o + fc_fyo)/(Di_D*fc_fyo*(2.4148731 - 0.114188455*fyo)))**0.30023974\n",
      "1.0004916287709327       Mean Squared Error: 0.07370264592864345          eta_i + eta_o + 2.99648502989275*(0.0403324563562754*Di_D**3 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "0.9966508486509429       Mean Squared Error: 0.0735849138124225          eta_i + eta_o + 2.99648502989275*(0.0535265644674524*Di_D**4 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "0.9931060686344766       Mean Squared Error: 0.07235784304943822          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**5/eta_i + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "0.9999300249443553       Mean Squared Error: 0.07202512854885683          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**4*(Di_D/eta_i + eta_o) + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "0.9972575372886341       Mean Squared Error: 0.0717578636965337          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**5*(Di_D + eta_o)/eta_i + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "0.9989553145469031       Mean Squared Error: 0.07169917814258617          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**5*(1.21380582748178*Di_D + eta_o)/eta_i + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n",
      "1.0013072156955065       Mean Squared Error: 0.07082807183304216          eta_i + eta_o + (Di_D**5*(eta_o + 0.77493906)/eta_i + (eta_o + fc_fyo)/(fc_fyo*(0.025856011593324*Di_D*fyo - 0.2264328*eta_i**2)))**0.30023974\n",
      "1.0008298592986065       Mean Squared Error: 0.07027056615535476          eta_i + eta_o + (Di_D**5*(Di_D + eta_o)/eta_i + (eta_o + fc_fyo)/(fc_fyo*(0.025856011593324*Di_D*fyo - 0.36121497687792*eta_i**2)))**0.30023974\n",
      "1.0032255718430911       Mean Squared Error: 0.07012380259309721          eta_i + eta_o + (Di_D**5*(Di_D + eta_o + fc_fyo)/eta_i + (eta_o + fc_fyo)/(fc_fyo*(0.025856011593324*Di_D*fyo - 0.36121497687792*eta_i**2)))**0.30023974\n",
      "1.004346088293504       Mean Squared Error: 0.0698131707707705          eta_i + eta_o + (Di_D**5*(eta_o + 0.77493906)/eta_i + (eta_o + fc_fyo)/(fc_fyo*(0.025856011593324*Di_D*fyo - 0.2264328*eta_i*(eta_i + 0.7931713) - 0.2264328*eta_i + 0.06928833490524)))**0.30023974\n"
     ]
    }
   ],
   "source": [
    "from pysr import PySRRegressor\n",
    "import numpy as np\n",
    "PC = pd.read_excel(\"PDS.xlsx\")\n",
    "PC['foAo']=(PC['Do']-PC['to'])*PC['to']*PC['fyo']*pi\n",
    "PC['fiAi']=(PC['Di']-PC['ti'])*PC['ti']*PC['fyi']*pi\n",
    "PC['fcAc']=((PC['Do']-2.0*PC['to'])**2.0-PC['Di']**2.0)*PC['fc']/4.0*pi#error\n",
    "\n",
    "PC['eta_i']=(PC['fiAi'])/(PC['fcAc'])\n",
    "PC['eta_o']=(PC['foAo'])/PC['fcAc']\n",
    "PC['xi_i']=(PC['fcAc'])/(PC['fiAi'])\n",
    "PC['xi_o']=PC['fcAc']/(PC['foAo'])\n",
    "PC['ti_D']=PC['ti']/PC['Di']\n",
    "PC['to_D']=PC['to']/PC['Do']\n",
    "PC['L_Do']=PC['L']/PC['Do']\n",
    "PC['D_ti']=PC['Di']/PC['ti']\n",
    "PC['D_to']=PC['Do']/PC['to']\n",
    "PC['fyi_fc']=(PC['fyi'])/(PC['fc'])\n",
    "PC['fyo_fc']=(PC['fyo'])/(PC['fc'])\n",
    "PC['fc_fyo']=PC['fc']/PC['fyo']\n",
    "PC['Di_D']=PC['Di']/PC['Do']\n",
    "\n",
    "PC['P_PcS']=PC['P']/(PC['fcAc'])*1000.0\n",
    "X=PC[['L_Do','eta_i','eta_o','fyi','fyo','fc_fyo','Di_D']]\n",
    "y=PC['P_PcS']\n",
    "\n",
    "model = PySRRegressor(niterations=100,populations=40,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{i} + \\eta_{o} + 2.99648502989275 \\left(0.0535265644674524 Di_{D}^{4} + \\frac{\\eta_{o} + fc_{fyo}}{Di_{D} fc_{fyo} fyo}\\right)^{0.30023974}$"
      ],
      "text/plain": [
       "eta_i + eta_o + 2.99648502989275*(0.0535265644674524*Di_D**4 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974"
      ]
     },
     "execution_count": 217,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=14\n",
    "model.sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9966508486509429       Mean Squared Error: 0.0735849138124225          eta_i + eta_o + 2.99648502989275*(0.0535265644674524*Di_D**4 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=i)\n",
    "aa=y_pred/y.values\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{i} + \\eta_{o} + 2.99648502989275 \\left(\\frac{0.025856011593324 Di_{D}^{5} \\cdot \\left(1.21380582748178 Di_{D} + \\eta_{o}\\right)}{\\eta_{i}} + \\frac{\\eta_{o} + fc_{fyo}}{Di_{D} fc_{fyo} fyo}\\right)^{0.30023974}$"
      ],
      "text/plain": [
       "eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**5*(1.21380582748178*Di_D + eta_o)/eta_i + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974"
      ]
     },
     "execution_count": 224,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=18\n",
    "model.sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9989553145469031       Mean Squared Error: 0.07169917814258617          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**5*(1.21380582748178*Di_D + eta_o)/eta_i + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=i)\n",
    "aa=y_pred/y.values\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bestttttttttttttttttttttt          no length in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{i} + \\eta_{o} + 2.99648502989275 \\left(0.0535265644674524 Di_{D}^{4} + \\frac{\\eta_{o} + fc_{fyo}}{Di_{D} fc_{fyo} fyo}\\right)^{0.30023974}$"
      ],
      "text/plain": [
       "eta_i + eta_o + 2.99648502989275*(0.0535265644674524*Di_D**4 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974"
      ]
     },
     "execution_count": 229,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=14\n",
    "model.sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9966508486509429       Mean Squared Error: 0.0735849138124225          eta_i + eta_o + 2.99648502989275*(0.0535265644674524*Di_D**4 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=i)\n",
    "aa=y_pred/y.values\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{i} + \\eta_{o} + 1.49779743724763 \\left(0.26038945 Di_{D}^{2} - \\frac{\\eta_{o} + fc_{fyo}}{Di_{D} fc_{fyo} \\left(2.4148731 - 0.114188455 fyo\\right)}\\right)^{0.30023974}$"
      ],
      "text/plain": [
       "eta_i + eta_o + 1.49779743724763*(0.26038945*Di_D**2 - (eta_o + fc_fyo)/(Di_D*fc_fyo*(2.4148731 - 0.114188455*fyo)))**0.30023974"
      ]
     },
     "execution_count": 235,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=12\n",
    "model.sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9934642987567288       Mean Squared Error: 0.07382578948545589          eta_i + eta_o + 1.49779743724763*(0.26038945*Di_D**2 - (eta_o + fc_fyo)/(Di_D*fc_fyo*(2.4148731 - 0.114188455*fyo)))**0.30023974\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=i)\n",
    "aa=y_pred/y.values\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{i} + \\eta_{o} + 2.99648502989275 \\left(0.025856011593324 Di_{D}^{2} + \\frac{\\eta_{o} + fc_{fyo}}{Di_{D} fc_{fyo} fyo}\\right)^{0.30023974}$"
      ],
      "text/plain": [
       "eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**2 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974"
      ]
     },
     "execution_count": 237,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=11\n",
    "model.sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0026741164148825       Mean Squared Error: 0.07415670654312659          eta_i + eta_o + 2.99648502989275*(0.025856011593324*Di_D**2 + (eta_o + fc_fyo)/(Di_D*fc_fyo*fyo))**0.30023974\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=i)\n",
    "aa=y_pred/y.values\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{i} + \\eta_{o} + \\left(\\frac{\\eta_{o} + fc_{fyo}}{fc_{fyo} \\left(0.0191056168340136 Di_{D} fyo + 0.317549512593324\\right)}\\right)^{0.2788979}$"
      ],
      "text/plain": [
       "eta_i + eta_o + ((eta_o + fc_fyo)/(fc_fyo*(0.0191056168340136*Di_D*fyo + 0.317549512593324)))**0.2788979"
      ]
     },
     "execution_count": 239,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=10\n",
    "model.sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9957780489917618       Mean Squared Error: 0.07675880768188613          eta_i + eta_o + ((eta_o + fc_fyo)/(fc_fyo*(0.0191056168340136*Di_D*fyo + 0.317549512593324)))**0.2788979\n"
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=i)\n",
    "aa=y_pred/y.values\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1364: FutureWarning: Index.is_object is deprecated.Use pandas.api.types.is_object_dtype instead\n",
      "  if X.columns.is_object() and X.columns.str.contains(\" \").any():\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling Julia backend...\n",
      "0.9241601015047198       Mean Squared Error: 0.3854855700832014          2.37126970000000\n",
      "0.9518183270344674       Mean Squared Error: 0.19669389472485754          eta_o + 1.5578655\n",
      "0.9936488609784905       Mean Squared Error: 0.14720998156797915          1.6691351*eta_o + 1.073433\n",
      "1.0138385195631137       Mean Squared Error: 0.1174950523445366          1.5289881**eta_i + 1.58925192735733*eta_o\n",
      "0.9808174550070446       Mean Squared Error: 0.10092625263878512          eta_i + eta_o + (Di_D/eta_o)**(-0.34887975)\n",
      "1.0075194408819588       Mean Squared Error: 0.09577021504554631          eta_i + eta_o + (fc_fyo/eta_o)**(fc_fyo - 0.2728227)\n",
      "1.0025772940700524       Mean Squared Error: 0.08789301396166309          eta_i + eta_o + 0.739557453620369/(Di_D*fc_fyo/eta_o)**0.21822666\n",
      "1.000081622978971       Mean Squared Error: 0.08151263478979859          eta_i + eta_o + 2.83165201420837/(Di_D*fc_fyo*fyo/eta_o)**0.24206564\n",
      "1.0025211387343393       Mean Squared Error: 0.07938281134011868          eta_i + eta_o + 2.71012307306078/(Di_D*fc_fyo*fyo/eta_o)**0.3510882 + 0.4139581\n",
      "1.0050033803286584       Mean Squared Error: 0.07793149072660568          eta_i + eta_o + 2.47148372623729/(Di_D*fyo*(fc_fyo - 0.015594775)/eta_o)**0.36451754 + 0.48317894\n",
      "1.0043739456380232       Mean Squared Error: 0.07805543456610362          1.0635232*eta_i + eta_o + 2.3891767554797/(Di_D*fyo*(fc_fyo - 0.01926961)/eta_o)**0.31217787 + 0.34339663\n",
      "1.0110553825425908       Mean Squared Error: 0.07698455182973858          eta_i + eta_o + 2.47148372623729/(Di_D*fyo*(fc_fyo - 0.015594775)/(eta_o + fc_fyo**2))**0.36451754 + 0.48317894\n",
      "1.0061789194027961       Mean Squared Error: 0.07766388731899551          eta_i**1.0595839 + eta_o + (Di_D*(-2.79676345108975*eta_i + 0.052617528*fyo)*(fc_fyo - 0.019685557)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0091984665540088       Mean Squared Error: 0.07718288801330266          eta_i**1.0595839 + eta_o + (Di_D*(fc_fyo - 0.019685557)*(-2.79676345108975*eta_i - 2.79676345108975*fc_fyo + 0.052617528*fyo)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0077723898071038       Mean Squared Error: 0.07724175535491455          eta_i**1.0595839 + eta_o + (Di_D*(fc_fyo - 0.019685557)*(-2.79676345108975*eta_i - 2.79676345108975*fc_fyo + 0.052617528*fyo + 0.193358811102569)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0115616592763002       Mean Squared Error: 0.0757605663081962          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(0.052617528*eta_o + 0.052617528*fyo + 2.0562003 - 2.2983227*(eta_i + fc_fyo)/Di_D)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0086814848469505       Mean Squared Error: 0.07588847803237658          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(0.052617528*fyo + 2.51109848381152 - 2.2983227*(eta_i + fc_fyo)/Di_D)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0099501784183087       Mean Squared Error: 0.07584742254821257          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(0.052617528*fyo + 2.51109848381152 - 2.2983227*(eta_i + 1.1668587*fc_fyo)/Di_D)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.006186129387978       Mean Squared Error: 0.07574631674580205          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(0.052617528*L_Do + 0.052617528*fyo + 3.25003172 - 2.2983227*(1.2068630924853*eta_i + fc_fyo)/Di_D)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.006442721145776       Mean Squared Error: 0.07571236559817152          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(0.052617528*L_Do + 0.052617528*eta_o + 0.052617528*fyo + 3.16099394613294 - 2.2983227*(1.2068630924853*eta_i + fc_fyo)/Di_D)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0069482045869163       Mean Squared Error: 0.07569459993376895          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(0.0267390859292424*L_Do + 0.052617528*eta_o + 0.052617528*fyo + 3.16099394613294 - 2.2983227*(1.2068630924853*eta_i + fc_fyo)/Di_D)/eta_o)**(-0.29722482) + 0.32054222\n",
      "1.0060718143574066       Mean Squared Error: 0.07551246869245769          eta_i + eta_o + (Di_D*(fc_fyo - 0.019685557)*(Di_D + 0.052617528*L_Do + 0.052617528*fyo + 2.43783563523907 - 2.2983227*(1.2068630924853*eta_i + fc_fyo)/Di_D + 0.052617528*eta_o/Di_D**2)/eta_o)**(-0.29722482) + 0.32054222\n"
     ]
    }
   ],
   "source": [
    "#recent\n",
    "from pysr import PySRRegressor\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from math import pi\n",
    "PC = pd.read_excel(\"PDS.xlsx\")\n",
    "PC['foAo']=(PC['Do']-PC['to'])*PC['to']*PC['fyo']*pi\n",
    "PC['fiAi']=(PC['Di']-PC['ti'])*PC['ti']*PC['fyi']*pi\n",
    "PC['fcAc']=((PC['Do']-2.0*PC['to'])**2.0-PC['Di']**2.0)*PC['fc']/4.0*pi#error\n",
    "\n",
    "PC['eta_i']=(PC['fiAi'])/(PC['fcAc'])\n",
    "PC['eta_o']=(PC['foAo'])/PC['fcAc']\n",
    "PC['xi_i']=(PC['fcAc'])/(PC['fiAi'])\n",
    "PC['xi_o']=PC['fcAc']/(PC['foAo'])\n",
    "PC['ti_D']=PC['ti']/PC['Di']\n",
    "PC['to_D']=PC['to']/PC['Do']\n",
    "PC['L_Do']=PC['L']/PC['Do']\n",
    "PC['D_ti']=PC['Di']/PC['ti']\n",
    "PC['D_to']=PC['Do']/PC['to']\n",
    "PC['fyi_fc']=(PC['fyi'])/(PC['fc'])\n",
    "PC['fyo_fc']=(PC['fyo'])/(PC['fc'])\n",
    "PC['fc_fyo']=PC['fc']/PC['fyo']\n",
    "PC['Di_D']=PC['Di']/PC['Do']\n",
    "\n",
    "PC['P_PcS']=PC['P']/(PC['fcAc'])*1000.0\n",
    "X=PC[['L_Do','eta_i','eta_o','fyi','fyo','fc_fyo','Di_D']]\n",
    "y=PC['P_PcS']\n",
    "\n",
    "model = PySRRegressor(niterations=100,populations=40,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Circular section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 800,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8216253460418619       Mean Squared Error: 0.4198046397263072          1.65663020000000\n",
      "0.9443871795769245       Mean Squared Error: 0.14483188957039458          eta + 1.23486\n",
      "1.0063605642118172       Mean Squared Error: 0.08779813786010932          1.5040045*eta + 1.0224893\n",
      "1.0022173367911176       Mean Squared Error: 0.08273599406764118          eta/lam**0.19515437 + 1.0221514\n",
      "1.0012046432939752       Mean Squared Error: 0.07744581617425582          eta + 0.95936775 + 0.00438432*t_D_1/fc_fy\n",
      "1.0007703960208565       Mean Squared Error: 0.0771250484915514          0.75847876**lam + eta + 0.004347205*t_D_1/fc_fy\n",
      "1.0007521169101503       Mean Squared Error: 0.07669168027657838          0.82219374**lam + eta**1.0325893 + 0.0042239283*t_D_1/fc_fy\n",
      "1.0043360230569658       Mean Squared Error: 0.07666439029949486          0.38505226**(lam**2) + eta**1.0437362 + 0.0042239283*t_D_1/fc_fy\n",
      "1.002029610056716       Mean Squared Error: 0.07650233436697544          0.7942968**lam + eta**1.0437362 + 0.0042239283*(0.09525341**eta + t_D_1)/fc_fy\n",
      "1.0002639503405706       Mean Squared Error: 0.07634779193948978          0.67415303**(lam - 0.709392492014766/t_D_1) + eta**1.0437362 + 0.0042239283*t_D_1/fc_fy\n",
      "1.0023588575510896       Mean Squared Error: 0.07632241422290248          0.67415303**(lam - 0.709392492014766/t_D_1) + eta**1.0437362 + 0.0042239283*(fc_fy + t_D_1)/fc_fy\n",
      "1.0017864341548133       Mean Squared Error: 0.07630426410099758          0.67415303**(lam - 0.709392492014766/(eta + t_D_1)) + eta**1.0437362 + 0.0042239283*(fc_fy + t_D_1)/fc_fy\n",
      "1.0015727412999524       Mean Squared Error: 0.07600492763065551          0.7942968**lam + eta**1.0437362 + 0.0042239283*(fc_fy**(1.4160115*lam*(eta + t_D_1) - 0.93048961848115) + t_D_1)/fc_fy\n",
      "1.0036676485104716       Mean Squared Error: 0.07596457015127979          0.7942968**lam + eta**1.0437362 + 0.0042239283*(fc_fy + fc_fy**(1.4160115*lam*(eta + t_D_1) - 0.93048961848115) + t_D_1)/fc_fy\n",
      "1.0014082630963186       Mean Squared Error: 0.0760218181902397          eta - 0.00292870506551424*(eta + t_D_1)/((1.0437362*eta)**t_D_1 + 0.6615228) + 0.9593839 + 0.00920364398411421/lam + 0.0045390404*(t_D_1 - 1.3702877)/fc_fy\n",
      "0.9986165059491032       Mean Squared Error: 0.07586142190455158          eta - 0.00288311757464071*(fc_fy + t_D_1)/((1.0437362*eta)**t_D_1 + 0.6615228) + 0.9593839 + 0.00906038236276849/lam + 0.0044683868318937*(t_D_1 - 1.3702877)/fc_fy\n",
      "0.998665728232018       Mean Squared Error: 0.0758611220478734          eta - 0.00288311757464071*(fc_fy + t_D_1)/((1.0437362*eta)**t_D_1 + 0.6615228) + 0.959483146111117 + 0.00906038236276849/lam + 0.0044683868318937*(t_D_1 - 1.3702877)/fc_fy\n",
      "1.0011321350648186       Mean Squared Error: 0.07597628139389546          eta - 0.00292870506551424*(t_D_1 + 1.7402545)/((1.0437362*eta)**t_D_1 + 0.6615228) + 0.9593839 + 0.00920364398411421/lam + 0.0045390404*(0.10720302**(2*eta) + t_D_1 - 1.3702877)/fc_fy\n",
      "0.9997396928751288       Mean Squared Error: 0.07599143114842062          eta - 0.00292870506551424*(t_D_1 + 1.7402545)/((1.0437362*eta)**t_D_1 + 0.6615228) + 0.956576340819521 + 0.00920364398411421/lam + 0.0045390404*(0.10720302**(2*eta) + t_D_1 - 1.3702877)/fc_fy\n",
      "1.0003799805853373       Mean Squared Error: 0.07579251155680806          eta + 0.9593839 - 0.00292870506551424*(t_D_1 + 1.7402545)/(eta**t_D_1 + 0.6615228) + 0.00920364398411421/lam + 0.0045390404*(0.10720302**(eta + fc_fy*t_D_1 - 0.75960031) + t_D_1 - 1.3702877)/fc_fy\n",
      "1.0008374125252286       Mean Squared Error: 0.07577857940071767          eta - 0.00292870506551424*(t_D_1 + 1.7402545)/((1.0437362*eta)**t_D_1 + 0.6615228) + 0.9593839 + 0.00920364398411421/lam + 0.0045390404*(0.10720302**(eta + fc_fy*t_D_1 - 0.75960031) + t_D_1 - 1.3702877)/fc_fy\n",
      "number:   27\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import math\n",
    "from math import pi\n",
    "PC  = pd.read_excel(\"PC.xlsx\")\n",
    "PC['P_PcS']=PC['P']/((PC['D']-PC['t'])*PC['t']*PC['Fy']+(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0)/pi*1000.0\n",
    "#X=PC[['D', 't', 'L', 'Fy', 'fc']];y=PC['P_PcS']\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "PC['Acfc']=(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0*pi\n",
    "PC['Asfy']=(PC['D']-PC['t'])*PC['t']*PC['Fy']*pi\n",
    "PC['I_s']=(PC['D'])**4.0/64.0*pi-(PC['D']-2.0*PC['t'])**4.0/64.0*pi;PC['I_c'] =(PC['D']-2.0*PC['t'])**4.0/64.0*pi\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "#print(PC['lambda'] )\n",
    "PC['D_t']=PC['D']/PC['t']\n",
    "PC['t_D']=PC['t']/PC['D']\n",
    "\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])#\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0#+PC['Asfy']\n",
    "PC['L_D']=PC['L']/PC['D']\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "PC['fc_fy']=PC['fc']/PC['Fy']\n",
    "\n",
    "PC['D_t_1']=PC['D_t']*PC['Fy']/208000.0\n",
    "PC['t_D_1']=1.0/PC['D_t_1']\n",
    "\n",
    "X=PC[['eta','lam','fc_fy','t_D_1']]\n",
    "y=PC['P_PcS']\n",
    "'''\n",
    "condition = PC['D_t_1'] < 0.07\n",
    "PC_filtered = PC[condition]\n",
    "\n",
    "X = PC_filtered[['eta','lam','fc_fy']]\n",
    "y = PC_filtered['P_PcS']'''\n",
    "dc=[]\n",
    "dc.append(X)\n",
    "dc.append(y)\n",
    "model = PySRRegressor(niterations=100,populations=60,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},parsimony= 0.1,\n",
    "                      constraints={\"^\":(-1,10)},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))\n",
    "modely.append(model)\n",
    "print(\"number:  \",len(modely)-1)\n",
    "XX.append(dc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bestttttttttttttttttttttttttt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 811,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0007703960208565       Mean Squared Error: 0.0771250484915514\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.75847876^{lam} + \\eta + \\frac{0.004347205 t_{D 1}}{fc_{fy}}$"
      ],
      "text/plain": [
       "0.75847876**lam + eta + 0.004347205*t_D_1/fc_fy"
      ]
     },
     "execution_count": 811,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=5\n",
    "j=27\n",
    "y_pred = modely[j].predict(XX[j][0],index=i)\n",
    "aa=(y_pred)/(XX[j][1].values)\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa))\n",
    "modely[j].sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 806,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0002639503405706       Mean Squared Error: 0.07634779193948978\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.67415303^{lam - \\frac{0.709392492014766}{t_{D 1}}} + \\eta^{1.0437362} + \\frac{0.0042239283 t_{D 1}}{fc_{fy}}$"
      ],
      "text/plain": [
       "0.67415303**(lam - 0.709392492014766/t_D_1) + eta**1.0437362 + 0.0042239283*t_D_1/fc_fy"
      ]
     },
     "execution_count": 806,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=9\n",
    "j=27\n",
    "y_pred = modely[j].predict(XX[j][0],index=i)\n",
    "aa=(y_pred)/(XX[j][1].values)\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa))\n",
    "modely[j].sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Circular section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 812,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8216253956380449       Mean Squared Error: 0.4198046397263072          1.65663030000000\n",
      "0.9443871795769245       Mean Squared Error: 0.14483188957039458          eta + 1.23486\n",
      "1.0064287229566484       Mean Squared Error: 0.08780374441629848          1.5043957*eta + 1.0223649\n",
      "0.9996844769197096       Mean Squared Error: 0.08356808004087882          1.0030426**(t_D_1/fc_fy) + eta\n",
      "1.0018095156578086       Mean Squared Error: 0.07915249549375578          eta + 0.9209952 + 0.004829349*t_D_1/fc_fy\n",
      "0.9985246089125102       Mean Squared Error: 0.0774498812378814          eta + 0.9533107 + 0.00439166375011599*t_D_1/fc_fy\n",
      "1.003155656525772       Mean Squared Error: 0.07731605621390675          0.73557305**lam + eta + 0.0044478072481248*t_D_1/fc_fy\n",
      "1.0010361034287507       Mean Squared Error: 0.07700841989000941          0.7014923**lam + eta + 0.024019573 + 0.00420136440961104*t_D_1/fc_fy\n",
      "0.9986326239510949       Mean Squared Error: 0.07693198631023143          eta + 0.9209952 + 0.00649905418222447/lam + 0.00420136440961104*t_D_1/fc_fy\n",
      "1.000032595466487       Mean Squared Error: 0.07620892147328628          eta + 0.9564133 + 0.094253674*(0.18459801*eta + 0.0340764252959601)/lam + 0.00307181814107904*t_D_1/fc_fy\n",
      "number:   28\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import math\n",
    "from math import pi\n",
    "PC  = pd.read_excel(\"PC.xlsx\")\n",
    "PC['P_PcS']=PC['P']/((PC['D']-PC['t'])*PC['t']*PC['Fy']+(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0)/pi*1000.0\n",
    "#X=PC[['D', 't', 'L', 'Fy', 'fc']];y=PC['P_PcS']\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "PC['Acfc']=(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0*pi\n",
    "PC['Asfy']=(PC['D']-PC['t'])*PC['t']*PC['Fy']*pi\n",
    "PC['I_s']=(PC['D'])**4.0/64.0*pi-(PC['D']-2.0*PC['t'])**4.0/64.0*pi;PC['I_c'] =(PC['D']-2.0*PC['t'])**4.0/64.0*pi\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "#print(PC['lambda'] )\n",
    "PC['D_t']=PC['D']/PC['t']\n",
    "PC['t_D']=PC['t']/PC['D']\n",
    "\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])#\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0#+PC['Asfy']\n",
    "PC['L_D']=PC['L']/PC['D']\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "PC['fc_fy']=PC['fc']/PC['Fy']\n",
    "\n",
    "PC['D_t_1']=PC['D_t']*PC['Fy']/208000.0\n",
    "PC['t_D_1']=1.0/PC['D_t_1']\n",
    "\n",
    "X=PC[['eta','lam','fc_fy','t_D_1']]\n",
    "y=PC['P_PcS']\n",
    "'''\n",
    "condition = PC['D_t_1'] < 0.07\n",
    "PC_filtered = PC[condition]\n",
    "\n",
    "X = PC_filtered[['eta','lam','fc_fy']]\n",
    "y = PC_filtered['P_PcS']'''\n",
    "dc=[]\n",
    "dc.append(X)\n",
    "dc.append(y)\n",
    "model = PySRRegressor(niterations=100,populations=60,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},parsimony= 0.5,\n",
    "                      constraints={\"^\":(-1,10)},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))\n",
    "modely.append(model)\n",
    "print(\"number:  \",len(modely)-1)\n",
    "XX.append(dc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BESTTTTTTTTTT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 822,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9986326239510949       Mean Squared Error: 0.07693198631023143\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta + 0.9209952 + \\frac{0.00649905418222447}{lam} + \\frac{0.00420136440961104 t_{D 1}}{fc_{fy}}$"
      ],
      "text/plain": [
       "eta + 0.9209952 + 0.00649905418222447/lam + 0.00420136440961104*t_D_1/fc_fy"
      ]
     },
     "execution_count": 822,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=8\n",
    "j=28\n",
    "y_pred = modely[j].predict(XX[j][0],index=i)\n",
    "aa=(y_pred)/(XX[j][1].values)\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa))\n",
    "modely[j].sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 821,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000032595466487       Mean Squared Error: 0.07620892147328628\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta + 0.9564133 + \\frac{0.094253674 \\cdot \\left(0.18459801 \\eta + 0.0340764252959601\\right)}{lam} + \\frac{0.00307181814107904 t_{D 1}}{fc_{fy}}$"
      ],
      "text/plain": [
       "eta + 0.9564133 + 0.094253674*(0.18459801*eta + 0.0340764252959601)/lam + 0.00307181814107904*t_D_1/fc_fy"
      ]
     },
     "execution_count": 821,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=9\n",
    "j=28\n",
    "y_pred = modely[j].predict(XX[j][0],index=i)\n",
    "aa=(y_pred)/(XX[j][1].values)\n",
    "print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa))\n",
    "modely[j].sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1364: FutureWarning: Index.is_object is deprecated.Use pandas.api.types.is_object_dtype instead\n",
      "  if X.columns.is_object() and X.columns.str.contains(\" \").any():\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8489915312465924       Mean Squared Error: 0.4077891582265128          1.76145720000000\n",
      "0.9573549793584742       Mean Squared Error: 0.14550337962164875          eta + 1.2818646\n",
      "1.004874412352299       Mean Squared Error: 0.08581642365856139          1.5163773*eta + 1.0167091\n",
      "1.0084818692881365       Mean Squared Error: 0.08243652978463031          1.66404895126144*eta**0.8875846 + 0.913833615006899\n",
      "1.0119787091935013       Mean Squared Error: 0.08003052127760675          eta*t_D_1**0.11581023 + fc_fy**(-0.06846441)\n",
      "1.0150325882585436       Mean Squared Error: 0.07842686154308219          eta*t_D_1**0.12227979 + 0.900005088574492*(1/fc_fy)**0.12227979\n",
      "1.0073611090122265       Mean Squared Error: 0.07723402567459944          eta*t_D_1**0.11335925 + fc_fy**(-0.018443204/(fc_fy + lam))\n",
      "1.0072440064875252       Mean Squared Error: 0.07709061373439328          eta*(t_D_1 - 0.84668136)**0.11395208 + fc_fy**(-0.019155843/(fc_fy + lam))\n",
      "1.0079706899726517       Mean Squared Error: 0.07608147506066158          1.03576189527509*eta*(-eta + 0.741818630685891*t_D_1)**0.117653534 + fc_fy**(-0.019112363/(fc_fy + lam))\n",
      "1.0019266906327755       Mean Squared Error: 0.07603622989986278          1.03041170408248*eta*(-eta + 0.775201189902818*t_D_1)**0.117653534 - 0.013386534 + fc_fy**(-0.019112363/(fc_fy + lam))\n",
      "1.0061118836271303       Mean Squared Error: 0.07634267383805193          1.01435196494342*eta*(0.8899984*lam + 0.8899984*t_D_1 - (eta + 4*fc_fy)*(fc_fy + 1.1973821))**0.12227979 + fc_fy**(-0.018132256/(fc_fy + lam))\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'modely' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[1;32mIn [4]\u001b[0m, in \u001b[0;36m<cell line: 51>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     49\u001b[0m     aa\u001b[38;5;241m=\u001b[39my_pred\u001b[38;5;241m/\u001b[39my\u001b[38;5;241m.\u001b[39mvalues\n\u001b[0;32m     50\u001b[0m     \u001b[38;5;28mprint\u001b[39m(np\u001b[38;5;241m.\u001b[39mmean(aa),\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m      Mean Squared Error:\u001b[39m\u001b[38;5;124m\"\u001b[39m,  np\u001b[38;5;241m.\u001b[39mstd(aa)\u001b[38;5;241m/\u001b[39mnp\u001b[38;5;241m.\u001b[39mmean(aa),\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m        \u001b[39m\u001b[38;5;124m\"\u001b[39m,model\u001b[38;5;241m.\u001b[39msympy(i))\n\u001b[1;32m---> 51\u001b[0m \u001b[43mmodely\u001b[49m\u001b[38;5;241m.\u001b[39mappend(model)\n\u001b[0;32m     52\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnumber:  \u001b[39m\u001b[38;5;124m\"\u001b[39m,\u001b[38;5;28mlen\u001b[39m(modely)\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[0;32m     53\u001b[0m XX\u001b[38;5;241m.\u001b[39mappend(dc)\n",
      "\u001b[1;31mNameError\u001b[0m: name 'modely' is not defined"
     ]
    }
   ],
   "source": [
    "#recent\n",
    "import pandas as pd\n",
    "import math\n",
    "from math import pi\n",
    "PC  = pd.read_excel(\"PC.xlsx\")\n",
    "PC['P_PcS']=PC['P']/((PC['D']-PC['t'])*PC['t']*PC['Fy']+(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0)/pi*1000.0\n",
    "#X=PC[['D', 't', 'L', 'Fy', 'fc']];y=PC['P_PcS']\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "PC['Acfc']=(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0*pi\n",
    "PC['Asfy']=(PC['D']-PC['t'])*PC['t']*PC['Fy']*pi\n",
    "PC['I_s']=(PC['D'])**4.0/64.0*pi-(PC['D']-2.0*PC['t'])**4.0/64.0*pi;PC['I_c'] =(PC['D']-2.0*PC['t'])**4.0/64.0*pi\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "#print(PC['lambda'] )\n",
    "PC['D_t']=PC['D']/PC['t']\n",
    "PC['t_D']=PC['t']/PC['D']\n",
    "\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])#\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0#+PC['Asfy']\n",
    "PC['L_D']=PC['L']/PC['D']\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "PC['fc_fy']=PC['fc']/PC['Fy']\n",
    "\n",
    "PC['D_t_1']=PC['D_t']*PC['Fy']/208000.0\n",
    "PC['t_D_1']=1.0/PC['D_t_1']\n",
    "\n",
    "X=PC[['eta','lam','fc_fy','t_D_1']]\n",
    "y=PC['P_PcS']\n",
    "'''\n",
    "condition = PC['D_t_1'] < 0.07\n",
    "PC_filtered = PC[condition]\n",
    "\n",
    "X = PC_filtered[['eta','lam','fc_fy']]\n",
    "y = PC_filtered['P_PcS']'''\n",
    "dc=[]\n",
    "dc.append(X)\n",
    "dc.append(y)\n",
    "model = PySRRegressor(niterations=100,populations=60,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},parsimony= 0.5,\n",
    "                      constraints={\"^\":(-1,10)},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))\n",
    "modely.append(model)\n",
    "print(\"number:  \",len(modely)-1)\n",
    "XX.append(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "rectangular"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 875,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.843454786192332       Mean Squared Error: 0.41727078646391447          1.72742160000000\n",
      "0.9937138604343894       Mean Squared Error: 0.082076394554808          eta + 1.0395808\n",
      "1.0039568901680438       Mean Squared Error: 0.07864283141631527          1.0739423*eta + 0.98694366\n",
      "1.00245109077856       Mean Squared Error: 0.07321405964299886          B_t_i**(-0.19963549) + eta + 0.15244283\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:86: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Failed to evaluate the expression. If you are using a custom operator, make sure to define it in `extra_sympy_mappings`, e.g., `model.set_params(extra_sympy_mappings={'inv': lambda x: 1/x})`, where `lambda x: 1/x` is a valid SymPy function defining the operator. You can then run `model.refresh()` to re-load the expressions.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1936\u001b[0m, in \u001b[0;36mPySRRegressor.predict\u001b[1;34m(self, X, index)\u001b[0m\n\u001b[0;32m   1933\u001b[0m         \u001b[39mreturn\u001b[39;00m np\u001b[39m.\u001b[39mstack(\n\u001b[0;32m   1934\u001b[0m             [eq[\u001b[39m\"\u001b[39m\u001b[39mlambda_format\u001b[39m\u001b[39m\"\u001b[39m](X) \u001b[39mfor\u001b[39;00m eq \u001b[39min\u001b[39;00m best_equation], axis\u001b[39m=\u001b[39m\u001b[39m1\u001b[39m\n\u001b[0;32m   1935\u001b[0m         )\n\u001b[1;32m-> 1936\u001b[0m     \u001b[39mreturn\u001b[39;00m best_equation[\u001b[39m\"\u001b[39;49m\u001b[39mlambda_format\u001b[39;49m\u001b[39m\"\u001b[39;49m](X)\n\u001b[0;32m   1937\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mException\u001b[39;00m \u001b[39mas\u001b[39;00m error:\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\export_numpy.py:36\u001b[0m, in \u001b[0;36mCallableEquation.__call__\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m     35\u001b[0m         X \u001b[39m=\u001b[39m X[:, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_selection]\n\u001b[1;32m---> 36\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_lambda(\u001b[39m*\u001b[39;49mX\u001b[39m.\u001b[39;49mT) \u001b[39m*\u001b[39m np\u001b[39m.\u001b[39mones(expected_shape)\n",
      "File \u001b[1;32m<lambdifygenerated-1851>:2\u001b[0m, in \u001b[0;36m_lambdifygenerated\u001b[1;34m(lam, B_t_i, H_B, eta)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m_lambdifygenerated\u001b[39m(lam, B_t_i, H_B, eta):\n\u001b[1;32m----> 2\u001b[0m     \u001b[39mreturn\u001b[39;00m amax((\u001b[39m0.9757709\u001b[39;49m,B_t_i\u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49m(\u001b[39m-\u001b[39;49m\u001b[39m0.20232658\u001b[39;49m) \u001b[39m+\u001b[39;49m eta), axis\u001b[39m=\u001b[39;49m\u001b[39m0\u001b[39;49m) \u001b[39m+\u001b[39m \u001b[39m0.15244283\u001b[39m\n",
      "File \u001b[1;32m<__array_function__ internals>:5\u001b[0m, in \u001b[0;36mamax\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:2754\u001b[0m, in \u001b[0;36mamax\u001b[1;34m(a, axis, out, keepdims, initial, where)\u001b[0m\n\u001b[0;32m   2641\u001b[0m \u001b[39m\u001b[39m\u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m   2642\u001b[0m \u001b[39mReturn the maximum of an array or maximum along an axis.\u001b[39;00m\n\u001b[0;32m   2643\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m   2752\u001b[0m \u001b[39m5\u001b[39;00m\n\u001b[0;32m   2753\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m-> 2754\u001b[0m \u001b[39mreturn\u001b[39;00m _wrapreduction(a, np\u001b[39m.\u001b[39;49mmaximum, \u001b[39m'\u001b[39;49m\u001b[39mmax\u001b[39;49m\u001b[39m'\u001b[39;49m, axis, \u001b[39mNone\u001b[39;49;00m, out,\n\u001b[0;32m   2755\u001b[0m                       keepdims\u001b[39m=\u001b[39;49mkeepdims, initial\u001b[39m=\u001b[39;49minitial, where\u001b[39m=\u001b[39;49mwhere)\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:86\u001b[0m, in \u001b[0;36m_wrapreduction\u001b[1;34m(obj, ufunc, method, axis, dtype, out, **kwargs)\u001b[0m\n\u001b[0;32m     84\u001b[0m             \u001b[39mreturn\u001b[39;00m reduction(axis\u001b[39m=\u001b[39maxis, out\u001b[39m=\u001b[39mout, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mpasskwargs)\n\u001b[1;32m---> 86\u001b[0m \u001b[39mreturn\u001b[39;00m ufunc\u001b[39m.\u001b[39mreduce(obj, axis, dtype, out, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mpasskwargs)\n",
      "\u001b[1;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32mc:\\Users\\osama\\OneDrive\\Desktop\\work\\svm - Copy.ipynb Cell 443\u001b[0m in \u001b[0;36m<cell line: 34>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/osama/OneDrive/Desktop/work/svm%20-%20Copy.ipynb#Z1241sZmlsZQ%3D%3D?line=33'>34</a>\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(\u001b[39mlen\u001b[39m(model\u001b[39m.\u001b[39mequations_)):\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/osama/OneDrive/Desktop/work/svm%20-%20Copy.ipynb#Z1241sZmlsZQ%3D%3D?line=34'>35</a>\u001b[0m     model\u001b[39m.\u001b[39msympy(i)\n\u001b[1;32m---> <a href='vscode-notebook-cell:/c%3A/Users/osama/OneDrive/Desktop/work/svm%20-%20Copy.ipynb#Z1241sZmlsZQ%3D%3D?line=35'>36</a>\u001b[0m     y_pred \u001b[39m=\u001b[39m model\u001b[39m.\u001b[39;49mpredict(X,index\u001b[39m=\u001b[39;49mi)\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/osama/OneDrive/Desktop/work/svm%20-%20Copy.ipynb#Z1241sZmlsZQ%3D%3D?line=36'>37</a>\u001b[0m     aa\u001b[39m=\u001b[39my_pred\u001b[39m/\u001b[39my\u001b[39m.\u001b[39mvalues\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/osama/OneDrive/Desktop/work/svm%20-%20Copy.ipynb#Z1241sZmlsZQ%3D%3D?line=37'>38</a>\u001b[0m     \u001b[39mprint\u001b[39m(np\u001b[39m.\u001b[39mmean(aa),\u001b[39m\"\u001b[39m\u001b[39m      Mean Squared Error:\u001b[39m\u001b[39m\"\u001b[39m,  np\u001b[39m.\u001b[39mstd(aa)\u001b[39m/\u001b[39mnp\u001b[39m.\u001b[39mmean(aa),\u001b[39m\"\u001b[39m\u001b[39m        \u001b[39m\u001b[39m\"\u001b[39m,model\u001b[39m.\u001b[39msympy(i))\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1938\u001b[0m, in \u001b[0;36mPySRRegressor.predict\u001b[1;34m(self, X, index)\u001b[0m\n\u001b[0;32m   1936\u001b[0m     \u001b[39mreturn\u001b[39;00m best_equation[\u001b[39m\"\u001b[39m\u001b[39mlambda_format\u001b[39m\u001b[39m\"\u001b[39m](X)\n\u001b[0;32m   1937\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mException\u001b[39;00m \u001b[39mas\u001b[39;00m error:\n\u001b[1;32m-> 1938\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[0;32m   1939\u001b[0m         \u001b[39m\"\u001b[39m\u001b[39mFailed to evaluate the expression. \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m   1940\u001b[0m         \u001b[39m\"\u001b[39m\u001b[39mIf you are using a custom operator, make sure to define it in `extra_sympy_mappings`, \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m   1941\u001b[0m         \u001b[39m\"\u001b[39m\u001b[39me.g., `model.set_params(extra_sympy_mappings=\u001b[39m\u001b[39m{\u001b[39m\u001b[39m'\u001b[39m\u001b[39minv\u001b[39m\u001b[39m'\u001b[39m\u001b[39m: lambda x: 1/x})`, where \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m   1942\u001b[0m         \u001b[39m\"\u001b[39m\u001b[39m`lambda x: 1/x` is a valid SymPy function defining the operator. \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m   1943\u001b[0m         \u001b[39m\"\u001b[39m\u001b[39mYou can then run `model.refresh()` to re-load the expressions.\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m   1944\u001b[0m     ) \u001b[39mfrom\u001b[39;00m \u001b[39merror\u001b[39;00m\n",
      "\u001b[1;31mValueError\u001b[0m: Failed to evaluate the expression. If you are using a custom operator, make sure to define it in `extra_sympy_mappings`, e.g., `model.set_params(extra_sympy_mappings={'inv': lambda x: 1/x})`, where `lambda x: 1/x` is a valid SymPy function defining the operator. You can then run `model.refresh()` to re-load the expressions."
     ]
    }
   ],
   "source": [
    "PC  = pd.read_excel(\"PR.xlsx\")#['P', 'EC4', 'FE', 'ref', 'Author ', 'Year', 'L', 'H', 'B', 't', 'Fy', 'fc', 'D/B', 'L/B', 'D/t']\n",
    "PC['Acfc']=(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t'])*(PC['fc'])\n",
    "PC['Asfy']=(PC['H']*PC['B']-(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t']))*PC['Fy']\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])\n",
    "PC['L_B']=PC['L']/PC['B']\n",
    "PC['H_B']=PC['H']/PC['B']\n",
    "PC['B_t']=PC['B']/PC['t']\n",
    "PC['H_t']=PC['H']/PC['t']\n",
    "\n",
    "PC['I_s']=(PC['B'])**3.0*PC['H']/12.0-(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0;    PC['I_c'] =(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0\n",
    "\n",
    "PC['B_t_i']=PC['B']/PC['t']*(PC['Fy']/208000)**0.5\n",
    "\n",
    "\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "X=PC[['lam','B_t_i','H_B','eta']]\n",
    "y=PC['P_PcS']\n",
    "\n",
    "dc=[]\n",
    "dc.append(X)\n",
    "dc.append(y)\n",
    "model = PySRRegressor(niterations=100,populations=60,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},parsimony= 0.5,\n",
    "                      constraints={\"^\":(-1,10)},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\",\"max\",\"min\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))\n",
    "modelyy.append(model)\n",
    "print(\"number:  \",len(modelyy)-1)\n",
    "XX_r.append(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 876,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number:   4\n"
     ]
    }
   ],
   "source": [
    "modelyy.append(model)\n",
    "print(\"number:  \",len(modelyy)-1)\n",
    "XX_r.append(dc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 883,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\max\\left(0.22359633, \\eta\\right) + 0.8592236 + \\frac{0.30603927}{\\max\\left(B_{t i}, 0.7431388 \\max\\left(H_{B}, H_{B}^{\\eta^{2}}\\right)\\right)}$"
      ],
      "text/plain": [
       "Max(0.22359633, eta) + 0.8592236 + 0.30603927/Max(B_t_i, 0.7431388*Max(H_B, H_B**(eta**2)))"
      ]
     },
     "execution_count": 883,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=9\n",
    "j=4\n",
    "modelyy[j].sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 881,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.09016864443202777"
      ]
     },
     "execution_count": 881,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(PC['eta'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bestttttttttttttttttttt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 884,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta + 0.8600225 + \\frac{0.30696008}{\\max\\left(B_{t i}, 0.7451394 \\max\\left(H_{B}, H_{B}^{\\eta^{2}}\\right)\\right)}$"
      ],
      "text/plain": [
       "eta + 0.8600225 + 0.30696008/Max(B_t_i, 0.7451394*Max(H_B, H_B**(eta**2)))"
      ]
     },
     "execution_count": 884,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=8\n",
    "j=4\n",
    "modelyy[j].sympy(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 886,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eta + 0.8600225 + 0.30696008/Max(B_t_i, 0.7451394*Max(H_B, H_B**(eta**2)))\n"
     ]
    }
   ],
   "source": [
    "print(model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1230: UserWarning: Note: Using a large maxsize for the equation search will be exponentially slower and use significant memory. You should consider turning `use_frequency` to False, and perhaps use `warmup_maxsize_by`.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1296: UserWarning: Note: it looks like you are running in Jupyter. The progress bar will be turned off.\n",
      "  warnings.warn(\n",
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1364: FutureWarning: Index.is_object is deprecated.Use pandas.api.types.is_object_dtype instead\n",
      "  if X.columns.is_object() and X.columns.str.contains(\" \").any():\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8527575160339996       Mean Squared Error: 0.4045322065784776          1.75936700000000\n",
      "0.9981491251248705       Mean Squared Error: 0.08302860312019467          eta + 1.0544628\n",
      "1.0063395589701236       Mean Squared Error: 0.07912290649415168          1.0748614*eta + 0.996135\n",
      "1.0029319525806824       Mean Squared Error: 0.07305953441108799          1.15540951714082/B_t_i**0.1642774 + eta\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:86: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Failed to evaluate the expression. If you are using a custom operator, make sure to define it in `extra_sympy_mappings`, e.g., `model.set_params(extra_sympy_mappings={'inv': lambda x: 1/x})`, where `lambda x: 1/x` is a valid SymPy function defining the operator. You can then run `model.refresh()` to re-load the expressions.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1936\u001b[0m, in \u001b[0;36mPySRRegressor.predict\u001b[1;34m(self, X, index)\u001b[0m\n\u001b[0;32m   1933\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39mstack(\n\u001b[0;32m   1934\u001b[0m             [eq[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mlambda_format\u001b[39m\u001b[38;5;124m\"\u001b[39m](X) \u001b[38;5;28;01mfor\u001b[39;00m eq \u001b[38;5;129;01min\u001b[39;00m best_equation], axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m\n\u001b[0;32m   1935\u001b[0m         )\n\u001b[1;32m-> 1936\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mbest_equation\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mlambda_format\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m   1937\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m error:\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\export_numpy.py:36\u001b[0m, in \u001b[0;36mCallableEquation.__call__\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m     35\u001b[0m         X \u001b[38;5;241m=\u001b[39m X[:, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_selection]\n\u001b[1;32m---> 36\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_lambda\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mX\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mT\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mones(expected_shape)\n",
      "File \u001b[1;32m<lambdifygenerated-38>:2\u001b[0m, in \u001b[0;36m_lambdifygenerated\u001b[1;34m(lam, B_t_i, H_B, eta)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_lambdifygenerated\u001b[39m(lam, B_t_i, H_B, eta):\n\u001b[1;32m----> 2\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mamax\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m1.1305108\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m1.15658199986782\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mB_t_i\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m0.1655198\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43meta\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36mamax\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:2793\u001b[0m, in \u001b[0;36mamax\u001b[1;34m(a, axis, out, keepdims, initial, where)\u001b[0m\n\u001b[0;32m   2680\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m   2681\u001b[0m \u001b[38;5;124;03mReturn the maximum of an array or maximum along an axis.\u001b[39;00m\n\u001b[0;32m   2682\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m   2791\u001b[0m \u001b[38;5;124;03m5\u001b[39;00m\n\u001b[0;32m   2792\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m-> 2793\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_wrapreduction\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmaximum\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mmax\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mout\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m   2794\u001b[0m \u001b[43m                      \u001b[49m\u001b[43mkeepdims\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeepdims\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minitial\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43minitial\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwhere\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mwhere\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:86\u001b[0m, in \u001b[0;36m_wrapreduction\u001b[1;34m(obj, ufunc, method, axis, dtype, out, **kwargs)\u001b[0m\n\u001b[0;32m     84\u001b[0m             \u001b[38;5;28;01mreturn\u001b[39;00m reduction(axis\u001b[38;5;241m=\u001b[39maxis, out\u001b[38;5;241m=\u001b[39mout, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mpasskwargs)\n\u001b[1;32m---> 86\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ufunc\u001b[38;5;241m.\u001b[39mreduce(obj, axis, dtype, out, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mpasskwargs)\n",
      "\u001b[1;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Input \u001b[1;32mIn [5]\u001b[0m, in \u001b[0;36m<cell line: 34>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     34\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(\u001b[38;5;28mlen\u001b[39m(model\u001b[38;5;241m.\u001b[39mequations_)):\n\u001b[0;32m     35\u001b[0m     model\u001b[38;5;241m.\u001b[39msympy(i)\n\u001b[1;32m---> 36\u001b[0m     y_pred \u001b[38;5;241m=\u001b[39m \u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredict\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43mindex\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mi\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     37\u001b[0m     aa\u001b[38;5;241m=\u001b[39my_pred\u001b[38;5;241m/\u001b[39my\u001b[38;5;241m.\u001b[39mvalues\n\u001b[0;32m     38\u001b[0m     \u001b[38;5;28mprint\u001b[39m(np\u001b[38;5;241m.\u001b[39mmean(aa),\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m      Mean Squared Error:\u001b[39m\u001b[38;5;124m\"\u001b[39m,  np\u001b[38;5;241m.\u001b[39mstd(aa)\u001b[38;5;241m/\u001b[39mnp\u001b[38;5;241m.\u001b[39mmean(aa),\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m        \u001b[39m\u001b[38;5;124m\"\u001b[39m,model\u001b[38;5;241m.\u001b[39msympy(i))\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1938\u001b[0m, in \u001b[0;36mPySRRegressor.predict\u001b[1;34m(self, X, index)\u001b[0m\n\u001b[0;32m   1936\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m best_equation[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mlambda_format\u001b[39m\u001b[38;5;124m\"\u001b[39m](X)\n\u001b[0;32m   1937\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m error:\n\u001b[1;32m-> 1938\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[0;32m   1939\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFailed to evaluate the expression. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1940\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIf you are using a custom operator, make sure to define it in `extra_sympy_mappings`, \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1941\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124me.g., `model.set_params(extra_sympy_mappings=\u001b[39m\u001b[38;5;124m{\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124minv\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m: lambda x: 1/x})`, where \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1942\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m`lambda x: 1/x` is a valid SymPy function defining the operator. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1943\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mYou can then run `model.refresh()` to re-load the expressions.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1944\u001b[0m     ) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01merror\u001b[39;00m\n",
      "\u001b[1;31mValueError\u001b[0m: Failed to evaluate the expression. If you are using a custom operator, make sure to define it in `extra_sympy_mappings`, e.g., `model.set_params(extra_sympy_mappings={'inv': lambda x: 1/x})`, where `lambda x: 1/x` is a valid SymPy function defining the operator. You can then run `model.refresh()` to re-load the expressions."
     ]
    }
   ],
   "source": [
    "#recent\n",
    "\n",
    "PC  = pd.read_excel(\"PR.xlsx\")#['P', 'EC4', 'FE', 'ref', 'Author ', 'Year', 'L', 'H', 'B', 't', 'Fy', 'fc', 'D/B', 'L/B', 'D/t']\n",
    "PC['Acfc']=(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t'])*(PC['fc'])\n",
    "PC['Asfy']=(PC['H']*PC['B']-(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t']))*PC['Fy']\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])\n",
    "PC['L_B']=PC['L']/PC['B']\n",
    "PC['H_B']=PC['H']/PC['B']\n",
    "PC['B_t']=PC['B']/PC['t']\n",
    "PC['H_t']=PC['H']/PC['t']\n",
    "\n",
    "PC['I_s']=(PC['B'])**3.0*PC['H']/12.0-(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0;    PC['I_c'] =(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0\n",
    "\n",
    "PC['B_t_i']=PC['B']/PC['t']*(PC['Fy']/208000)**0.5\n",
    "\n",
    "\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "X=PC[['lam','B_t_i','H_B','eta']]\n",
    "y=PC['P_PcS']\n",
    "\n",
    "dc=[]\n",
    "dc.append(X)\n",
    "dc.append(y)\n",
    "model = PySRRegressor(niterations=100,populations=60,population_size=60,maxsize=50,\n",
    "                      nested_constraints={\"^\":{\"^\":0,\"/\":1},\"/\":{\"/\":0,\"^\":1}},parsimony= 0.5,\n",
    "                      constraints={\"^\":(-1,10)},\n",
    "                      binary_operators=[\"+\", \"*\",\"^\",\"/\",\"max\",\"min\"],#special(x,y)=x^2.0+y\n",
    "                      denoise=True,loss=\"loss(x, y) = abs(x/y-1)\",\n",
    "                      )\n",
    "model.fit(X,y)\n",
    "for i in range(len(model.equations_)):\n",
    "    model.sympy(i)\n",
    "    y_pred = model.predict(X,index=i)\n",
    "    aa=y_pred/y.values\n",
    "    print(np.mean(aa),\"      Mean Squared Error:\",  np.std(aa)/np.mean(aa),\"        \",model.sympy(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\max\\left(1.1316425, \\max\\left(\\eta, 0.620444171806376 \\left(H_{B} \\eta\\right)^{0.59587467}\\right) + \\left(0.42914733 B_{t i} - 0.0136135618038255 \\eta\\right)^{- \\frac{0.17197762}{H_{B}}}\\right)$"
      ],
      "text/plain": [
       "Max(1.1316425, Max(eta, 0.620444171806376*(H_B*eta)**0.59587467) + (0.42914733*B_t_i - 0.0136135618038255*eta)**(-0.17197762/H_B))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.sympy(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:86: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Failed to evaluate the expression. If you are using a custom operator, make sure to define it in `extra_sympy_mappings`, e.g., `model.set_params(extra_sympy_mappings={'inv': lambda x: 1/x})`, where `lambda x: 1/x` is a valid SymPy function defining the operator. You can then run `model.refresh()` to re-load the expressions.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1936\u001b[0m, in \u001b[0;36mPySRRegressor.predict\u001b[1;34m(self, X, index)\u001b[0m\n\u001b[0;32m   1933\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39mstack(\n\u001b[0;32m   1934\u001b[0m             [eq[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mlambda_format\u001b[39m\u001b[38;5;124m\"\u001b[39m](X) \u001b[38;5;28;01mfor\u001b[39;00m eq \u001b[38;5;129;01min\u001b[39;00m best_equation], axis\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m\n\u001b[0;32m   1935\u001b[0m         )\n\u001b[1;32m-> 1936\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mbest_equation\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mlambda_format\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m   1937\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m error:\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\export_numpy.py:36\u001b[0m, in \u001b[0;36mCallableEquation.__call__\u001b[1;34m(self, X)\u001b[0m\n\u001b[0;32m     35\u001b[0m         X \u001b[38;5;241m=\u001b[39m X[:, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_selection]\n\u001b[1;32m---> 36\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_lambda\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mX\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mT\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mones(expected_shape)\n",
      "File \u001b[1;32m<lambdifygenerated-43>:2\u001b[0m, in \u001b[0;36m_lambdifygenerated\u001b[1;34m(lam, B_t_i, H_B, eta)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_lambdifygenerated\u001b[39m(lam, B_t_i, H_B, eta):\n\u001b[1;32m----> 2\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mamax\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m1.1478395\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43meta\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m0.41471872\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mB_t_i\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m0.175183872971254\u001b[39;49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43mH_B\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32m<__array_function__ internals>:180\u001b[0m, in \u001b[0;36mamax\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:2793\u001b[0m, in \u001b[0;36mamax\u001b[1;34m(a, axis, out, keepdims, initial, where)\u001b[0m\n\u001b[0;32m   2680\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m   2681\u001b[0m \u001b[38;5;124;03mReturn the maximum of an array or maximum along an axis.\u001b[39;00m\n\u001b[0;32m   2682\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m   2791\u001b[0m \u001b[38;5;124;03m5\u001b[39;00m\n\u001b[0;32m   2792\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m-> 2793\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_wrapreduction\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmaximum\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mmax\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mout\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m   2794\u001b[0m \u001b[43m                      \u001b[49m\u001b[43mkeepdims\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mkeepdims\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minitial\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43minitial\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwhere\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mwhere\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:86\u001b[0m, in \u001b[0;36m_wrapreduction\u001b[1;34m(obj, ufunc, method, axis, dtype, out, **kwargs)\u001b[0m\n\u001b[0;32m     84\u001b[0m             \u001b[38;5;28;01mreturn\u001b[39;00m reduction(axis\u001b[38;5;241m=\u001b[39maxis, out\u001b[38;5;241m=\u001b[39mout, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mpasskwargs)\n\u001b[1;32m---> 86\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ufunc\u001b[38;5;241m.\u001b[39mreduce(obj, axis, dtype, out, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mpasskwargs)\n",
      "\u001b[1;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Input \u001b[1;32mIn [13]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[0m y_pred \u001b[38;5;241m=\u001b[39m \u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpredict\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43mindex\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m6\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Users\\osama\\anaconda3\\lib\\site-packages\\pysr\\sr.py:1938\u001b[0m, in \u001b[0;36mPySRRegressor.predict\u001b[1;34m(self, X, index)\u001b[0m\n\u001b[0;32m   1936\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m best_equation[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mlambda_format\u001b[39m\u001b[38;5;124m\"\u001b[39m](X)\n\u001b[0;32m   1937\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m error:\n\u001b[1;32m-> 1938\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[0;32m   1939\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFailed to evaluate the expression. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1940\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIf you are using a custom operator, make sure to define it in `extra_sympy_mappings`, \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1941\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124me.g., `model.set_params(extra_sympy_mappings=\u001b[39m\u001b[38;5;124m{\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124minv\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m: lambda x: 1/x})`, where \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1942\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m`lambda x: 1/x` is a valid SymPy function defining the operator. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1943\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mYou can then run `model.refresh()` to re-load the expressions.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   1944\u001b[0m     ) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01merror\u001b[39;00m\n",
      "\u001b[1;31mValueError\u001b[0m: Failed to evaluate the expression. If you are using a custom operator, make sure to define it in `extra_sympy_mappings`, e.g., `model.set_params(extra_sympy_mappings={'inv': lambda x: 1/x})`, where `lambda x: 1/x` is a valid SymPy function defining the operator. You can then run `model.refresh()` to re-load the expressions."
     ]
    }
   ],
   "source": [
    "y_pred = model.predict(X,index=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimizing the paramters of above equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 896,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated Parameters: [0.84552551 0.29532943 0.83170079]\n"
     ]
    }
   ],
   "source": [
    "PC  = pd.read_excel(\"PR.xlsx\")#['P', 'EC4', 'FE', 'ref', 'Author ', 'Year', 'L', 'H', 'B', 't', 'Fy', 'fc', 'D/B', 'L/B', 'D/t']\n",
    "PC['Acfc']=(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t'])*(PC['fc'])\n",
    "PC['Asfy']=(PC['H']*PC['B']-(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t']))*PC['Fy']\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])\n",
    "PC['L_B']=PC['L']/PC['B']\n",
    "PC['H_B']=PC['H']/PC['B']\n",
    "PC['B_t']=PC['B']/PC['t']\n",
    "PC['H_t']=PC['H']/PC['t']\n",
    "\n",
    "PC['I_s']=(PC['B'])**3.0*PC['H']/12.0-(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0;    PC['I_c'] =(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0\n",
    "\n",
    "PC['B_t_i']=PC['B']/PC['t']*(PC['Fy']/208000)**0.5\n",
    "\n",
    "\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "X=PC[['lam','B_t_i','H_B','eta']]\n",
    "y=PC['P_PcS']\n",
    "\n",
    "X= np.array(X.values)\n",
    "y= np.array(y.values)\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# Define the non-linear model\n",
    "def eqn(parameters, X):\n",
    "    a, b, c = parameters\n",
    "    return X[:,3] + a + b/np.max((X[:,1],c*np.max((X[:,2], X[:,2]**(X[:,3]**2)),axis=0)),axis=0)\n",
    "\n",
    "# Define the likelihood function\n",
    "def fncy(parameters, X, y):\n",
    "    y_pred = eqn(parameters, X)    # Compute the predicted output using the non-linear model\n",
    "    aa=y_pred/y\n",
    "    return np.std(aa)/np.mean(aa)\n",
    "\n",
    "\n",
    "initial_parameters = [0.8600225,0.30696008,0.7451394]# Define initial parameter values\n",
    "result = minimize(fncy, initial_parameters, args=(X, y), method='BFGS')# Optimize the parameters using MLE\n",
    "estimated_parameters = result.x# Retrieve the estimated parameters\n",
    "\n",
    "# Print the estimated parameters\n",
    "print(\"Estimated Parameters:\", estimated_parameters)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated Parameters: [0.85319727 0.31348255 0.82780286] 1.9350013367552292\n",
      "Estimated Parameters: [0.85, 0.31, 0.83] 1.9380742440398122\n",
      "Data exported to SR_rect.xlsx\n"
     ]
    }
   ],
   "source": [
    "#recent\n",
    "PC  = pd.read_excel(\"PR.xlsx\")#['P', 'EC4', 'FE', 'ref', 'Author ', 'Year', 'L', 'H', 'B', 't', 'Fy', 'fc', 'D/B', 'L/B', 'D/t']\n",
    "PC['Acfc']=(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t'])*(PC['fc'])\n",
    "PC['Asfy']=(PC['H']*PC['B']-(PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t']))*PC['Fy']\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])\n",
    "PC['L_B']=PC['L']/PC['B']\n",
    "PC['H_B']=PC['H']/PC['B']\n",
    "PC['B_t']=PC['B']/PC['t']\n",
    "PC['H_t']=PC['H']/PC['t']\n",
    "X2=PC[['H', 'B', 't', 'L', 'Fy', 'fc']]\n",
    "\n",
    "PC['I_s']=(PC['B'])**3.0*PC['H']/12.0-(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0;    PC['I_c'] =(PC['B']-2.0*PC['t'])**3.0*(PC['H']-2.0*PC['t'])/12.0\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0\n",
    "\n",
    "PC['B_t_i']=PC['B']/PC['t']*(PC['Fy']/208000)**0.5\n",
    "\n",
    "\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "X=PC[['lam','B_t_i','H_B','eta']]\n",
    "y=PC['P_PcS']\n",
    "\n",
    "X= np.array(X.values)\n",
    "y= np.array(y.values)\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# Define the non-linear model\n",
    "def eqn(parameters, X):\n",
    "    a, b, c = parameters\n",
    "    return X[:,3] + a + b/np.max((X[:,1],c*np.max((X[:,2], X[:,2]**(X[:,3]**2)),axis=0)),axis=0)\n",
    "\n",
    "# Define the likelihood function\n",
    "def fncy(parameters, X, y):\n",
    "    y_pred = eqn(parameters, X)    # Compute the predicted output using the non-linear model\n",
    "    aa=y_pred/y\n",
    "    return np.sum((y_pred/y-1)**2.0)#np.std(aa)/np.mean(aa)\n",
    "\n",
    "\n",
    "initial_parameters = [0.8600225,0.30696008,0.81394]# Define initial parameter values\n",
    "result = minimize(fncy, initial_parameters, args=(X, y), method='BFGS')# Optimize the parameters using MLE\n",
    "estimated_parameters = result.x# Retrieve the estimated parameters\n",
    "\n",
    "# Print the estimated parameters\n",
    "print(\"Estimated Parameters:\", estimated_parameters,fncy(estimated_parameters,X, y))\n",
    "pp=[0.85,0.31,0.83]\n",
    "print(\"Estimated Parameters:\",pp, fncy(pp,X, y))\n",
    "\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train2, X_test2, y_train, y_test = train_test_split(X2,y, test_size=0.2, random_state=101)\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2, random_state=101)\n",
    "\n",
    "#save data\n",
    "X1_test=np.concatenate((X_test, X_train), axis=0);X2_test=np.concatenate((X_test2, X_train2), axis=0);y1_test=np.concatenate((y_test, y_train), axis=0)\n",
    "y1_pred = eqn(pp, X1_test)\n",
    "data=np.concatenate((X2_test,y1_test.reshape(-1, 1)),axis=1)\n",
    "data2=np.concatenate((data,y1_pred.reshape(-1, 1)),axis=1)\n",
    "df = pd.DataFrame(data2)\n",
    "df.columns = ['H', 'B', 't', 'L', 'Fy', 'fc','P_exp','P_pred'];PC=df\n",
    "PC['P_exp']=PC['P_exp']*((PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t'])*PC['fc'])/1000.0\n",
    "PC['P_pred']=PC['P_pred']*((PC['H']-2.0*PC['t'])*(PC['B']-2.0*PC['t'])*PC['fc'])/1000.0\n",
    "filename = 'SR_rect.xlsx'\n",
    "df.to_excel(filename, index=False)\n",
    "print(\"Data exported to\", filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated Parameters: [1.17215738 0.94915715 1.54977974 2.67568274] 3.801801040368559\n",
      "Estimated Parameters: [1.17, 0.95, 1.55, 2.7] 3.8022721786549902\n",
      "Data exported to SR_circ.xlsx\n"
     ]
    }
   ],
   "source": [
    "#recent\n",
    "import pandas as pd\n",
    "import math\n",
    "from math import pi\n",
    "PC  = pd.read_excel(\"PC.xlsx\")\n",
    "PC['P_PcS']=PC['P']/((PC['D']-PC['t'])*PC['t']*PC['Fy']+(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0)/pi*1000.0\n",
    "X2=PC[['D', 't', 'L', 'Fy', 'fc']]\n",
    "#X=PC[['D', 't', 'L', 'Fy', 'fc']];y=PC['P_PcS']\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "PC['Acfc']=(PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0*pi\n",
    "PC['Asfy']=(PC['D']-PC['t'])*PC['t']*PC['Fy']*pi\n",
    "PC['I_s']=(PC['D'])**4.0/64.0*pi-(PC['D']-2.0*PC['t'])**4.0/64.0*pi;PC['I_c'] =(PC['D']-2.0*PC['t'])**4.0/64.0*pi\n",
    "PC['N_pl'] = PC['Asfy'] + 0.85 * PC['Acfc']\n",
    "PC['N_cr'] = (pi ** 2.0 * ( 200000.0 * PC['I_s'] + 0.6 * 4400.0*PC['fc']**0.50 * PC['I_c'])) / PC['L']** 2\n",
    "#print(A_s*PC['Fy']/A_c/PC['fc'])\n",
    "PC['lam'] =(PC['N_pl'] / PC['N_cr'])**0.5####################     min(0.5, (N_pl / N_cr))\n",
    "#print(PC['lambda'] )\n",
    "PC['D_t']=PC['D']/PC['t']\n",
    "PC['t_D']=PC['t']/PC['D']\n",
    "\n",
    "PC['eta']=PC['Asfy']/(PC['Acfc'])#\n",
    "PC['P_PcS']=PC['P']/(PC['Acfc'])*1000.0#+PC['Asfy']\n",
    "PC['L_D']=PC['L']/PC['D']\n",
    "#X=np.array(PC[['eta','D_t']].values)\n",
    "PC['fc_fy']=PC['fc']/PC['Fy']\n",
    "\n",
    "PC['D_t_1']=PC['D_t']*PC['Fy']/208000.0\n",
    "PC['t_D_1']=1.0/PC['D_t_1']\n",
    "\n",
    "X=PC[['eta','lam','fc_fy','t_D_1']]\n",
    "y=PC['P_PcS']\n",
    "X= np.array(X.values)\n",
    "y= np.array(y.values)\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# Define the non-linear model\n",
    "def eqn(parameters, X):\n",
    "    a, b, c,d = parameters\n",
    "    #print(1/c/X[:,1]**2.0)\n",
    "    return X[:,0]*a + b + 1/c/1000/X[:,1]**2.0+d/1000.0/X[:,2]*X[:,3]\n",
    "\n",
    "# Define the likelihood function\n",
    "def fncy(parameters, X, y):\n",
    "    y_pred = eqn(parameters, X)    # Compute the predicted output using the non-linear model\n",
    "    aa=y_pred/y\n",
    "    return np.sum((y_pred/y-1)**2.0)#np.std(aa)/np.mean(aa)\n",
    "\n",
    "\n",
    "initial_parameters = [2,1,5,10]# Define initial parameter values\n",
    "result = minimize(fncy, initial_parameters, args=(X, y), method='BFGS')# Optimize the parameters using MLE\n",
    "estimated_parameters = result.x# Retrieve the estimated parameters\n",
    "\n",
    "# Print the estimated parameters\n",
    "print(\"Estimated Parameters:\", estimated_parameters,fncy(estimated_parameters,X, y))\n",
    "pp=[1.17,0.95,1.550,2.7]\n",
    "print(\"Estimated Parameters:\",pp, fncy(pp,X, y))\n",
    "\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train2, X_test2, y_train, y_test = train_test_split(X2,y, test_size=0.2, random_state=42)\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2, random_state=42)\n",
    "\n",
    "#save data\n",
    "X1_test=np.concatenate((X_test, X_train), axis=0);X2_test=np.concatenate((X_test2, X_train2), axis=0);y1_test=np.concatenate((y_test, y_train), axis=0)\n",
    "y1_pred = eqn(pp, X1_test)\n",
    "data=np.concatenate((X2_test,y1_test.reshape(-1, 1)),axis=1)\n",
    "data2=np.concatenate((data,y1_pred.reshape(-1, 1)),axis=1)\n",
    "df = pd.DataFrame(data2)\n",
    "df.columns = ['D', 't', 'L', 'Fy', 'fc','P_exp','P_pred'];PC=df\n",
    "PC['P_exp']=PC['P_exp']*((PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0)*pi/1000.0\n",
    "PC['P_pred']=PC['P_pred']*((PC['D']-2.0*PC['t'])**2.0*PC['fc']/4.0)*pi/1000.0\n",
    "filename = 'SR_circ.xlsx'\n",
    "df.to_excel(filename, index=False)\n",
    "print(\"Data exported to\", filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated Parameters: [5.49746705 0.32765549] 0.05695139871720351\n",
      "Estimated Parameters: [5.5, 0.33] 0.05696835745721406\n",
      "Data exported to SR_DS.xlsx\n"
     ]
    }
   ],
   "source": [
    "#recent\n",
    "from pysr import PySRRegressor\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from math import pi\n",
    "PC = pd.read_excel(\"PDS.xlsx\")\n",
    "PC['foAo']=(PC['Do']-PC['to'])*PC['to']*PC['fyo']*pi\n",
    "PC['fiAi']=(PC['Di']-PC['ti'])*PC['ti']*PC['fyi']*pi\n",
    "PC['fcAc']=((PC['Do']-2.0*PC['to'])**2.0-PC['Di']**2.0)*PC['fc']/4.0*pi#error\n",
    "X2=PC[['Do', 'to', 'Di', 'ti', 'fyo', 'fyi', 'fc']]\n",
    "PC['eta_i']=(PC['fiAi'])/(PC['fcAc'])\n",
    "PC['eta_o']=(PC['foAo'])/PC['fcAc']\n",
    "PC['xi_i']=(PC['fcAc'])/(PC['fiAi'])\n",
    "PC['xi_o']=PC['fcAc']/(PC['foAo'])\n",
    "PC['ti_D']=PC['ti']/PC['Di']\n",
    "PC['to_D']=PC['to']/PC['Do']\n",
    "PC['L_Do']=PC['L']/PC['Do']\n",
    "PC['D_ti']=PC['Di']/PC['ti']\n",
    "PC['D_to']=PC['Do']/PC['to']\n",
    "PC['fyi_fc']=(PC['fyi'])/(PC['fc'])\n",
    "PC['fyo_fc']=(PC['fyo'])/(PC['fc'])\n",
    "PC['fc_fyo']=PC['fc']/PC['fyo']\n",
    "PC['Di_D']=PC['Di']/PC['Do']\n",
    "\n",
    "PC['P_PcS']=PC['P']/(PC['fcAc'])*1000.0\n",
    "X=PC[['L_Do','eta_i','eta_o','fyi','fyo','fc_fyo','Di_D','fc']]\n",
    "y=PC['P_PcS']\n",
    "X= np.array(X.values)#[:195,:]\n",
    "y= np.array(y.values)#[:195]\n",
    "#print(np.maximum(0.2, np.minimum(0.85, X[:, 6])))\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# Define the non-linear model\n",
    "def eqn(parameters, X):\n",
    "    c = np.array(parameters)\n",
    "    X[:,6]=np.maximum(0.2, np.minimum(0.85, X[:, 6]))\n",
    "    #print(1/c/X[:,1]**2.0)\n",
    "    return X[:,1]+X[:,2] + (X[:,6]**2.0+210/c/X[:,7]/X[:,6]*(X[:,2]+X[:,5]))**0.3\n",
    "\n",
    "\n",
    "def fncy1(parameters, X, y):\n",
    "    y_pred = eqn(parameters, X)    # Compute the predicted output using the non-linear model\n",
    "    aa=y_pred/y\n",
    "    return np.mean(abs(y_pred/y-1))#np.std(aa)/np.mean(aa)\n",
    "initial_parameters = [5]# Define initial parameter values\n",
    "result = minimize(fncy1, initial_parameters, args=(X, y), method='BFGS')# Optimize the parameters using MLE\n",
    "estimated_parameters = result.x# Retrieve the estimated parameters\n",
    "\n",
    "# Print the estimated parameters\n",
    "print(\"Estimated Parameters:\", estimated_parameters,fncy1(estimated_parameters,X, y))\n",
    "pp=[5]\n",
    "print(\"Estimated Parameters:\",pp, fncy1(estimated_parameters,X, y))\n",
    "#print(\"Estimated Parameters:\",pp, (fncy1(pp,X, y)))\n",
    "\n",
    "#X,y = np.array(X.values),np.array(y.values)\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train2, X_test2, y_train, y_test = train_test_split(X2,y, test_size=0.2, random_state=42)\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.2, random_state=42)\n",
    "#save data\n",
    "X1_test=np.concatenate((X_test, X_train), axis=0);X2_test=np.concatenate((X_test2, X_train2), axis=0);y1_test=np.concatenate((y_test, y_train), axis=0)\n",
    "y1_pred = eqn(pp, X1_test)\n",
    "data=np.concatenate((X2_test,y1_test.reshape(-1, 1)),axis=1)\n",
    "data2=np.concatenate((data,y1_pred.reshape(-1, 1)),axis=1)\n",
    "df = pd.DataFrame(data2)\n",
    "df.columns = ['Do', 'to', 'Di', 'ti', 'fyo', 'fyi', 'fc','P_exp','P_pred'];PC=df\n",
    "PC['P_exp']=PC['P_exp']*(((PC['Do']-2.0*PC['to'])**2.0-PC['Di']**2.0)*PC['fc']/4.0*pi)/1000.0\n",
    "PC['P_pred']=PC['P_pred']*(((PC['Do']-2.0*PC['to'])**2.0-PC['Di']**2.0)*PC['fc']/4.0*pi)/1000.0\n",
    "filename = 'SR_DS.xlsx'\n",
    "df.to_excel(filename, index=False)\n",
    "print(\"Data exported to\", filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}