Tutorial_26a_Py_C_Call.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running Python or C Code from Julia\n",
    "\n",
    "This notebook provides a basic introduction to how to run Python, using the `PyCall.jl` package and also C code (requires no package) from Julia. Please see the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) homepage for instructions for how to either use an existing Python installation, or let the package make one.\n",
    " \n",
    "You need Python's `statsmodels` package installed to run the code below. If you have let PyCall install Python for you, use [Conda.jl](https://github.com/JuliaPy/Conda.jl) to add packages: `import Conda; Conda.add(\"statsmodels\")`.\n",
    "\n",
    "An alternative package (not used here) for running Python is [PythonCall.jl](https://github.com/cjdoris/PythonCall.jl), which seems to be gaining popularity.\n",
    "\n",
    "Another notebooks discusses how to run R code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "using Printf, DelimitedFiles\n",
    "include(\"src/printmat.jl\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OLS coeffs according to Julia\n",
      "    -0.504\n",
      "     1.341\n",
      "\n"
     ]
    }
   ],
   "source": [
    "x = readdlm(\"Data/MyData.csv\",',',skipstart=1)  #reading the csv file\n",
    "\n",
    "(Rme,Rf,R) = (x[:,2],x[:,3],x[:,4])  #creating variables from columns of x\n",
    "y  = R - Rf                          #do R .- Rf if R has several columns\n",
    "\n",
    "c = ones(length(Rme))\n",
    "x = [c Rme]\n",
    "\n",
    "b = x\\y\n",
    "println(\"OLS coeffs according to Julia\")\n",
    "printmat(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python\n",
    "\n",
    "In the next cells we *(a)* load the PyCall.jl package and activates the (Python) package `statsmodels`; *(b)* call some functions (eg. `OLS()`) from statsmodels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "using PyCall\n",
    "sm = pyimport(\"statsmodels.api\");     #activate this package and call it `sm`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyObject <class 'statsmodels.iolib.summary.Summary'>\n",
      "\"\"\"\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.519\n",
      "Model:                            OLS   Adj. R-squared:                  0.518\n",
      "Method:                 Least Squares   F-statistic:                     416.2\n",
      "Date:                Mon, 25 Nov 2024   Prob (F-statistic):           2.72e-63\n",
      "Time:                        09:42:56   Log-Likelihood:                -1241.7\n",
      "No. Observations:                 388   AIC:                             2487.\n",
      "Df Residuals:                     386   BIC:                             2495.\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -0.5042      0.305     -1.654      0.099      -1.103       0.095\n",
      "x1             1.3410      0.066     20.401      0.000       1.212       1.470\n",
      "==============================================================================\n",
      "Omnibus:                      259.682   Durbin-Watson:                   1.870\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             5249.944\n",
      "Skew:                           2.482   Prob(JB):                         0.00\n",
      "Kurtosis:                      20.323   Cond. No.                         4.68\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "\"\"\"\n"
     ]
    }
   ],
   "source": [
    "resultsP = sm.OLS(y, x).fit()        #can use Python functions directly\n",
    "\n",
    "println(resultsP.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[:HC0_se, :HC1_se, :HC2_se, :HC3_se, :_HCCM, :__class__, :__delattr__, :__dict__, :__dir__, :__doc__, :__eq__, :__format__, :__ge__, :__getattribute__, :__getstate__, :__gt__, :__hash__, :__init__, :__init_subclass__, :__le__, :__lt__, :__module__, :__ne__, :__new__, :__reduce__, :__reduce_ex__, :__repr__, :__setattr__, :__sizeof__, :__str__, :__subclasshook__, :__weakref__, :_abat_diagonal, :_cache, :_data_attr, :_data_in_cache, :_get_robustcov_results, :_get_wald_nonlinear, :_is_nested, :_transform_predict_exog, :_use_t, :_wexog_singular_values, :aic, :bic, :bse, :centered_tss, :compare_f_test, :compare_lm_test, :compare_lr_test, :condition_number, :conf_int, :conf_int_el, :cov_HC0, :cov_HC1, :cov_HC2, :cov_HC3, :cov_kwds, :cov_params, :cov_type, :df_model, :df_resid, :diagn, :eigenvals, :el_test, :ess, :f_pvalue, :f_test, :fittedvalues, :fvalue, :get_influence, :get_prediction, :get_robustcov_results, :info_criteria, :initialize, :k_constant, :llf, :load, :model, :mse_model, :mse_resid, :mse_total, :nobs, :normalized_cov_params, :outlier_test, :params, :predict, :pvalues, :remove_data, :resid, :resid_pearson, :rsquared, :rsquared_adj, :save, :scale, :ssr, :summary, :summary2, :t_test, :t_test_pairwise, :tvalues, :uncentered_tss, :use_t, :wald_test, :wald_test_terms, :wresid]\n"
     ]
    }
   ],
   "source": [
    "println(keys(resultsP))              #print all keys (field names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mComparing the estimates in Julia and Python:\u001b[22m\u001b[39m\n",
      "     Julia    Python\n",
      "    -0.504    -0.504\n",
      "     1.341     1.341\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b_P = resultsP.params                #the numerical results are now a Julia vector\n",
    "\n",
    "printblue(\"Comparing the estimates in Julia and Python:\")\n",
    "printmat([b b_P];colNames=[\"Julia\",\"Python\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " -0.5041626034967046\n",
       "  1.3410486453848383"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#we can run blocks of code like this, notice: $x and $y\n",
    "py\"\"\"\n",
    "import numpy as np\n",
    "xx = np.matmul(np.matrix.transpose($x),$x)\n",
    "xy = np.matmul(np.matrix.transpose($x),$y)\n",
    "b_p = np.linalg.solve(xx,xy)\n",
    "\"\"\"\n",
    "\n",
    "py\"b_p\"               #to print the result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# C\n",
    "\n",
    "This section shows some simple examples of how to call a C function. The functions are in the file `My_C_Stuff.c` (printed in the next cell). The first function `c_dot` defines a dot product between two vectors and the second function `c_ols` a simple linear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#include <stddef.h>\n",
      "\n",
      "// calculate the inner (dot) product of vectors Y and Y, returns the result (Sxy)\n",
      "double c_dot(size_t n, double *Y, double *X) {\n",
      "    double Sxy = 0.0;\n",
      "    for (size_t i = 0; i < n; ++i) {\n",
      "        Sxy += X[i]*Y[i];\n",
      "    }\n",
      "    return Sxy;\n",
      "}\n",
      "\n",
      "// calculate a simple regression, Y = a + b*X + u, puts (a,b) in vector ab, returns nothing\n",
      "void c_ols(size_t n, double *Y, double *X, double *ab) {\n",
      "    double Sx = 0.0, Sy = 0.0, Sxx = 0.0, Sxy = 0.0;\n",
      "    for (size_t i = 0; i < n; ++i) {\n",
      "        Sx  += X[i];\n",
      "        Sy  += Y[i];\n",
      "        Sxx += X[i]*X[i];\n",
      "        Sxy += X[i]*Y[i];\n",
      "    }\n",
      "    ab[1] = (Sxy-Sx*Sy/n)/(Sxx-Sx*Sx/n);   //slope\n",
      "    ab[0] = (Sy - ab[1]*Sx)/n;             //intercept\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "println(read(\"Data/My_C_Stuff.c\",String))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compile to a dynamlic library (dll on windows), I use gcc (for x86_64) from [mingw-64](http://mingw-w64.org)\n",
    "and run the following in the mingw terminal\n",
    "```\n",
    "gcc -shared -fPIC My_C_Stuff.c -o My_C_Stuff.dll\n",
    "```\n",
    "\n",
    "To call the C functions, place the dll file in the current folder and then run the following cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mylibc = \"My_C_Stuff.dll\"\n",
    "x2     = x[:,2];               #get a vector with the regressor values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Function which Returns a Number\n",
    "\n",
    "In the next example, we a function `c_dot` in `My_C_Stuff.dll`. The function calculates the inner product of two vectors.\n",
    "\n",
    "The details are:\n",
    "1. `mylibc.c_dot` is the library.function\n",
    "2. `length(y)::Csize_t` is the first input and its type (an integer indicating the number of elements in `y`)\n",
    "3. `y::Ptr{Float64}` is the second input (a pointer to an array of Floats) and similarly for `x2`\n",
    "4. `Float64` is the type of the output\n",
    "\n",
    "(We could potentially wrap this in a Julia function that checks for the right input types and outputs the result.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The inner product of x2 and y in Julia and C:   11071.648           11071.648\n"
     ]
    }
   ],
   "source": [
    "z = @ccall mylibc.c_dot(length(y)::Csize_t, y::Ptr{Float64}, x2::Ptr{Float64})::Float64\n",
    "\n",
    "printlnPs(\"The inner product of x2 and y in Julia and C:  \",x2'y,\" \",z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Function which Returns a Vector\n",
    "\n",
    "The details are as above, except that \n",
    "1. `mylibc.c_ols` is the library.function\n",
    "2. `Cvoid` is the type of the output, which here indicates that the function does not have an output. Rather, the function modifies the vector `b_c` by putting the OLS results there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comparing the estimates in Julia and C\n",
      "     Julia         C\n",
      "    -0.504    -0.504\n",
      "     1.341     1.341\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b_c = zeros(2)          #where C will store the regression results\n",
    "\n",
    "@ccall mylibc.c_ols(length(y)::Csize_t, y::Ptr{Float64}, x2::Ptr{Float64}, b_c::Ptr{Float64})::Cvoid\n",
    "\n",
    "println(\"Comparing the estimates in Julia and C\")\n",
    "printmat([b b_c];colNames=[\"Julia\",\"C\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.11.1",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}