Instantly share code, notes, and snippets.

kvedala/c-algorithms.ipynb

Last active July 21, 2020 22:54
Show Gist options
• Save kvedala/27f1b0b6502af935f6917673ec43bcd7 to your computer and use it in GitHub Desktop.
Compile and run C Algorithms
Display the source blob
Display the rendered blob
Raw
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 { "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "name": "Compile and run C Algorithms", "provenance": [], "collapsed_sections": [], "include_colab_link": true }, "kernelspec": { "name": "python3", "display_name": "Python 3" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "id": "b29wElh1z2Ye", "colab_type": "text" }, "source": [ "# Initializer\n", "Documentation available here: https://kvedala.github.io/C" ] }, { "cell_type": "code", "metadata": { "id": "WJLeiiCj6UZs", "colab_type": "code", "colab": {} }, "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from matplotlib import widgets, animation\n", "from IPython.display import display\n", "import ipywidgets\n", "from plotly import express as px\n", "from plotly import graph_objects as go\n", "from matplotlib.animation import FuncAnimation\n", "rand = np.random.default_rng()" ], "execution_count": 0, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "gALUN9Dbz4iK", "colab_type": "text" }, "source": [ "# Setup environment\n", "1. Clone the [GitHub repository](https://github.com/kvedala/C.git) \n", "2. Compile and build the code" ] }, { "cell_type": "code", "metadata": { "id": "f6IiAkSWVNgI", "colab_type": "code", "colab": {} }, "source": [ "%%bash\n", "apt -qq update\n", "apt -qq install ninja-build cmake\n", "git clone https://github.com/kvedala/C.git\n", "mkdir C/build\n", "cd C/build\n", "cmake -G Ninja ..\n", "cmake --build ." ], "execution_count": 0, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "m8oQJFVe0SkS", "colab_type": "text" }, "source": [ "# [Durand Kerner Roots](https://kvedala.github.io/C/da/d38/durand__kerner__roots_8c.html)\n", "## Setup variables\n", "For Durand-Kerner algorithm, this creates a polynomial of degree N. The coefficients are randomly distributed in the interval [a,b]. The " ] }, { "cell_type": "code", "metadata": { "id": "AbQQ_sMSUzAq", "colab_type": "code", "outputId": "8b6cdfa0-2c9b-4697-c368-95116b4f7d42", "colab": { "base_uri": "https://localhost:8080/", "height": 51 } }, "source": [ "a = -2\n", "b = 2\n", "N = 10\n", "coeffs = rand.random(size=(N,)) * (b - a) + a \n", "\n", "print(\"The polynomial to find roots is: \")\n", "coeffs_str = \"\"\n", "for i, c in enumerate(coeffs):\n", " coeffs_str += \"%.4g \" % c\n", " print(\"%.4g x^%d\" %(c, N-i-1), end=\" + \")\n", "print(\"\\b\\b= 0\")" ], "execution_count": 0, "outputs": [ { "output_type": "stream", "text": [ "The polynomial to find roots is: \n", "1.208 x^9 + 1.19 x^8 + 1.699 x^7 + -0.727 x^6 + -0.8313 x^5 + 1.535 x^4 + 1.814 x^3 + -1.277 x^2 + 1.899 x^1 + 1.384 x^0 + \b\b= 0\n" ], "name": "stdout" } ] }, { "cell_type": "markdown", "metadata": { "id": "aDrGH8hs3uLh", "colab_type": "text" }, "source": [ "## Execute the compiled C-function" ] }, { "cell_type": "code", "metadata": { "id": "a0oTYsz_Y633", "colab_type": "code", "outputId": "a1398b85-350c-472a-baec-675309871b96", "colab": { "base_uri": "https://localhost:8080/", "height": 272 } }, "source": [ "!./build/numerical_methods/durand_kerner_roots \$coeffs_str" ], "execution_count": 0, "outputs": [ { "output_type": "stream", "text": [ "Computing the roots for:\n", "\t(1.208) x^9 + (1.19) x^8 + (1.699) x^7 + (-0.727) x^6 + (-0.8313) x^5 + (1.535) x^4 + (1.814) x^3 + (-1.277) x^2 + (1.899) x^1 + (1.384) x^0 = 0\n", "\n", "Iterations: 149\n", "\t 0.902+0.4628j\n", "\t 0.902-0.4628j\n", "\t 0.36-0.8036j\n", "\t-0.5329 +1.444j\n", "\t-0.9732 -0.524j\n", "\t-0.5329 -1.444j\n", "\t-0.9732 +0.524j\n", "\t 0.36+0.8036j\n", "\t-0.4969+2.816e-17j\n", "absolute average change: 1.065e-13\n", "Time taken: 0.006181 sec\n" ], "name": "stdout" } ] }, { "cell_type": "markdown", "metadata": { "id": "4sMFnaiZ3vy1", "colab_type": "text" }, "source": [ "## Read stored log file\n", "Convert the log file into a pandas data frame and convert complex number values from string format to complex numbers. Save the real and imaginary components of the complex root approximations in separate columns for convenience." ] }, { "cell_type": "code", "metadata": { "id": "O_l3Ei0U6tKA", "colab_type": "code", "colab": {} }, "source": [ "d = pd.read_csv('durand_kerner.log.csv')\n", "root_columns = [c for c in d.columns if c not in ['iter#', 'avg. correction']]" ], "execution_count": 0, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "rK7KGwRg7ct8", "colab_type": "code", "outputId": "672956b8-54bf-4eec-d319-42645e0fac71", "colab": { "base_uri": "https://localhost:8080/", "height": 598 } }, "source": [ "to_complex = lambda s: complex(s.replace(' ', ''))\n", "\n", "for col in root_columns:\n", " d.loc[:,col] = d.loc[:,col].apply(to_complex)\n", " d.loc[:,col + '_r'] = d.loc[:,col].apply(np.real).astype(float)\n", " d.loc[:,col + '_i'] = d.loc[:,col].apply(np.imag).astype(float)\n", " \n", "display(d.head())\n", "display(d.tail())" ], "execution_count": 0, "outputs": [ { "output_type": "display_data", "data": { "text/html": [ "
\n", "\n", "\n", "
iter#root_0root_1root_2root_3root_4root_5root_6root_7root_8avg. correctionroot_0_rroot_0_iroot_1_rroot_1_iroot_2_rroot_2_iroot_3_rroot_3_iroot_4_rroot_4_iroot_5_rroot_5_iroot_6_rroot_6_iroot_7_rroot_7_iroot_8_rroot_8_i
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
001.270000e+09+1.814000e+09j9.655000e+08+1.588000e+09j8.856000e+08+1.458000e+09j1.163000e+09+1.927000e+08j1.318000e+09+1.574000e+09j1.793000e+09+2.053000e+09j9.696000e+08+7.220000e+08j1.593000e+09+1.745000e+09j1.651000e+09+6.870000e+08jNaN1.270000e+091.814000e+099.655000e+081.588000e+098.856000e+081.458000e+091.163000e+09192700000.01.318000e+091.574000e+091.793000e+092.053000e+09969600000.0722000000.01.593000e+091.745000e+091.651000e+09687000000.0
11-5.346000e+13-3.378000e+13j9.796000e+08+1.679000e+09j8.667000e+08+1.449000e+09j1.163000e+09+1.926000e+08j1.281000e+09+1.343000e+09j1.705000e+09+2.012000e+09j9.699000e+08+7.220000e+08j1.824000e+09+1.839000e+09j1.653000e+09+6.835000e+08j6.324000e+13-5.346000e+13-3.378000e+139.796000e+081.679000e+098.667000e+081.449000e+091.163000e+09192600000.01.281000e+091.343000e+091.705000e+092.012000e+09969900000.0722000000.01.824000e+091.839000e+091.653000e+09683500000.0
22-1.044000e+10-9.918000e+09j-1.149000e+11+1.375000e+11j7.891000e+08+1.481000e+09j1.163000e+09+1.931000e+08j1.124000e+09+8.713000e+08j-1.465000e+08+1.580000e+09j9.640000e+08+7.077000e+08j1.967000e+09+1.665000e+09j1.554000e+09+6.766000e+08j6.322000e+13-1.044000e+10-9.918000e+09-1.149000e+111.375000e+117.891000e+081.481000e+091.163000e+09193100000.01.124000e+098.713000e+08-1.465000e+081.580000e+09964000000.0707700000.01.967000e+091.665000e+091.554000e+09676600000.0
33-9.990000e+09-1.028000e+10j2.037000e+09+3.472000e+09j-2.201000e+09+2.112000e+09j1.164000e+09+2.223000e+08j5.561000e+08+2.046000e+09j-9.227000e+07+1.530000e+09j9.263000e+08+6.980000e+08j3.237000e+09+6.670000e+09j1.523000e+09+5.961000e+08j1.779000e+11-9.990000e+09-1.028000e+102.037000e+093.472000e+09-2.201000e+092.112000e+091.164000e+09222300000.05.561000e+082.046000e+09-9.227000e+071.530000e+09926300000.0698000000.03.237000e+096.670000e+091.523000e+09596100000.0
44-5.612000e+09-8.478000e+09j3.001000e+09+1.080000e+10j-2.167000e+09+2.194000e+09j1.165000e+09+2.236000e+08j5.949000e+08+1.913000e+09j-9.958000e+07+1.525000e+09j9.236000e+08+6.988000e+08j1.237000e+10+1.801000e+10j1.520000e+09+5.875000e+08j1.456000e+10-5.612000e+09-8.478000e+093.001000e+091.080000e+10-2.167000e+092.194000e+091.165000e+09223600000.05.949000e+081.913000e+09-9.958000e+071.525000e+09923600000.0698800000.01.237000e+101.801000e+101.520000e+09587500000.0
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "" ], "text/plain": [ " iter# root_0 ... root_8_r root_8_i\n", "0 0 1.270000e+09+1.814000e+09j ... 1.651000e+09 687000000.0\n", "1 1 -5.346000e+13-3.378000e+13j ... 1.653000e+09 683500000.0\n", "2 2 -1.044000e+10-9.918000e+09j ... 1.554000e+09 676600000.0\n", "3 3 -9.990000e+09-1.028000e+10j ... 1.523000e+09 596100000.0\n", "4 4 -5.612000e+09-8.478000e+09j ... 1.520000e+09 587500000.0\n", "\n", "[5 rows x 29 columns]" ] }, "metadata": { "tags": [] } }, { "output_type": "display_data", "data": { "text/html": [ "
\n", "\n", "\n", "
iter#root_0root_1root_2root_3root_4root_5root_6root_7root_8avg. correctionroot_0_rroot_0_iroot_1_rroot_1_iroot_2_rroot_2_iroot_3_rroot_3_iroot_4_rroot_4_iroot_5_rroot_5_iroot_6_rroot_6_iroot_7_rroot_7_iroot_8_rroot_8_i
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
1451450.919300+0.458100j0.897400-0.462900j0.352400-0.784900j-0.534100+1.444000j-0.976800-0.525100j-0.532800-1.443000j-0.955400+0.524600j0.366100+0.777600j-0.493200-0.006509j2.847000e-010.91930.45810.8974-0.46290.3524-0.7849-0.53411.444-0.9768-0.5251-0.5328-1.443-0.95540.52460.36610.7776-0.4932-6.509000e-03
1461460.901200+0.463400j0.902100-0.462700j0.360400-0.804100j-0.532900+1.444000j-0.973200-0.523800j-0.532900-1.444000j-0.973000+0.523600j0.360200+0.803700j-0.496900+0.000004j2.677000e-020.90120.46340.9021-0.46270.3604-0.8041-0.53291.444-0.9732-0.5238-0.5329-1.444-0.97300.52360.36020.8037-0.49694.303000e-06
1471470.902000+0.462800j0.902000-0.462800j0.360000-0.803600j-0.532900+1.444000j-0.973200-0.524000j-0.532900-1.444000j-0.973200+0.524000j0.360000+0.803600j-0.496900-0.000000j9.285000e-040.90200.46280.9020-0.46280.3600-0.8036-0.53291.444-0.9732-0.5240-0.5329-1.444-0.97320.52400.36000.8036-0.4969-2.211000e-12
1481480.902000+0.462800j0.902000-0.462800j0.360000-0.803600j-0.532900+1.444000j-0.973200-0.524000j-0.532900-1.444000j-0.973200+0.524000j0.360000+0.803600j-0.496900+0.000000j4.987000e-070.90200.46280.9020-0.46280.3600-0.8036-0.53291.444-0.9732-0.5240-0.5329-1.444-0.97320.52400.36000.8036-0.49698.900000e-17
1491490.902000+0.462800j0.902000-0.462800j0.360000-0.803600j-0.532900+1.444000j-0.973200-0.524000j-0.532900-1.444000j-0.973200+0.524000j0.360000+0.803600j-0.496900+0.000000j1.065000e-130.90200.46280.9020-0.46280.3600-0.8036-0.53291.444-0.9732-0.5240-0.5329-1.444-0.97320.52400.36000.8036-0.49692.816000e-17