gauss.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra # do operacji macierzowych i testowania\n",
    "using RandomNumbers.MersenneTwisters # do MT19937()\n",
    "using Printf\n",
    "using DelimitedFiles\n",
    "\n",
    "epsilon = 10e-8;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test data\n",
    "function Ab_read_from_file(path)    \n",
    "    local D1 = (open(readdlm, path))\n",
    "    local T = []\n",
    "    for i=1:50\n",
    "        T = [T;(M =(A=convert(Array{Float64,2}, D1[(i-1)*100 + i+1:i*100+i, 1:end-1]), b=convert(Array{Float64,1}, D1[(i-1)*100 + i+1:i*100+i, end])), desc =\"Big 100x100 \"*string(i))]\n",
    "    end\n",
    "    return T\n",
    "    end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test na macierzach. \n",
    "# f: funkcja zwracającą wynik do przetestowania z testem\n",
    "# name: nazwa metody liczącej rząd macierzy\n",
    "# rng: generator liczb losowych. Można przekazać np. MT19937(<seed>) aby generować te same zestawy losowych danych\n",
    "# exact_f: funkcja do liczenia dokładnego wyniku, chyba że test zawiera pole `exact`, wtedy to ono jest użyte\n",
    "# equal_f: opcjonalne do porównywania wartości wyliczonej przez f() i exact_f() lub test.exact_f(). Domyślnie: `==`\n",
    "function test(f, f_name, exact_f; rng=nothing, equal_f=(x, y)->x==y)\n",
    "    local _rand = rng==nothing ? rand : (args...)->rand(rng, args...) # custom rand if generator is set\n",
    "      \n",
    "    local Id_100x100          = (b=_rand(-10.0:10.0, 100), A=Matrix{Float64}(I, 100, 100))\n",
    "    local small_3x3           = (A=[1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 10.0], b=[10.0 12.0 13.0]')\n",
    "    local big_100x100         = Ab_read_from_file(\"gauss1.in\")\n",
    "    # tests consist of fieds M=matrix, desc=description of test and exact=wanted\n",
    "    tests = [\n",
    "        (M=Id_100x100,          desc=\"Identity 100x100 matrix with random b\")\n",
    "        (M=small_3x3 ,          desc=\"Small hardcoded 3x3 test\");\n",
    "        big_100x100\n",
    "    ];\n",
    "    \n",
    "    @printf(\">>> Random test ID: [\\e[1;35m%.16f\\e[m]\\n\", _rand())\n",
    "\n",
    "    function _do_test(test)\n",
    "        local t1 = @elapsed exact = haskey(test, :exact_f) ? test.exact_f(test.M) : exact_f(test.M)\n",
    "        local t2 = @elapsed computed = f(test.M) \n",
    "        local success = equal_f(exact, computed)\n",
    "         @printf(\"[%5s] [exact:\\ttime:%f] [%s:\\ttime:%f] <<< %s\\n\", \n",
    "                 success ? \"\\e[1;32mPASS\\e[m\" : \"\\e[1;31mFAIL\\e[m\", t1, f_name, t2, test.desc)\n",
    "    end\n",
    "\n",
    "    @printf(\"> Hardcoded tests:\\n\")\n",
    "    for test in tests\n",
    "        _do_test(test)\n",
    "    end\n",
    "    \n",
    "    flush(stdout)\n",
    "    @printf(\">>> TESTS END\\n\")\n",
    "end\n",
    "\n",
    "function getPrecisionU(T=Float64) return 0.5*2.0^-precision(T) end\n",
    "# zamienia liczby zmiennoprzecinkowe w macierzy na zera jeśli |x|≦|maxError|\n",
    "function filterZeros(matrix, maxError=epsilon) map(x->abs(x)<=abs(maxError) ? 0.0 : x, matrix) end\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solves Ax=b where Aₙₓₙ matrix of coefficients, bₙₓ₁ matrix of results for linear equations\n",
    "function solve!(A, b)\n",
    "    local M = [A b]'\n",
    "    local n = length(M[1, :])\n",
    "    local isZero = x->abs(x)<=epsilon\n",
    "    # Gauss reduce on columns\n",
    "    for col in 1:n\n",
    "        # swap-in column with nonzero value\n",
    "        local allzero = true\n",
    "        for j in col:n \n",
    "            if (!isZero(M[col, j])) \n",
    "                M[:, col], M[:, j], allzero = M[:, j], M[:, col], false\n",
    "                break\n",
    "            end\n",
    "        end\n",
    "        if (allzero) continue end # skip to next step iff row with zeros till the end\n",
    "        \n",
    "        # gauss reduce columns\n",
    "        for next_col in col+1:n\n",
    "            M[:, next_col] -= M[:, col]*M[col, next_col]/M[col, col] # denom !=0 cuz we swapped\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    for col in n:-1:1\n",
    "        if (isZero(M[col, col]))\n",
    "            throw(:UNSOLVABLE)\n",
    "#             if (!isZero(M[end, col])) return NaN*ones(length(b)) end # no solutions fuz we have 0*x = y (!=0)\n",
    "#             return Inf*ones(length(b)) # inf solutions cuz we have 0*x = 0\n",
    "        end\n",
    "        M[:, col] /= M[col, col] # calculate b[col]\n",
    "        for j in col-1:-1:1\n",
    "            M[:, j] -= M[col, j]*M[:, col]\n",
    "        end\n",
    "    end\n",
    "    return M[end, :]\n",
    "    \n",
    "end\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">>> Random test ID: [\u001b[1;35m0.9610688406006067\u001b[m]\n",
      "> Hardcoded tests:\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.000024] [solve!:\ttime:0.077138] <<< Identity 100x100 matrix with random b\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.000013] [solve!:\ttime:0.000007] <<< Small hardcoded 3x3 test\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.029928] [solve!:\ttime:0.014657] <<< Big 100x100 1\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.009936] [solve!:\ttime:0.016688] <<< Big 100x100 2\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.028179] [solve!:\ttime:0.014899] <<< Big 100x100 3\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001478] [solve!:\ttime:0.016826] <<< Big 100x100 4\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.004257] [solve!:\ttime:0.017867] <<< Big 100x100 5\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.008792] [solve!:\ttime:0.014898] <<< Big 100x100 6\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.038288] [solve!:\ttime:0.012812] <<< Big 100x100 7\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.013450] [solve!:\ttime:0.017904] <<< Big 100x100 8\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001066] [solve!:\ttime:0.017286] <<< Big 100x100 9\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.007394] [solve!:\ttime:0.017833] <<< Big 100x100 10\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.031534] [solve!:\ttime:0.012543] <<< Big 100x100 11\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.019626] [solve!:\ttime:0.015567] <<< Big 100x100 12\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.011630] [solve!:\ttime:0.016707] <<< Big 100x100 13\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.021579] [solve!:\ttime:0.018692] <<< Big 100x100 14\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006599] [solve!:\ttime:0.013553] <<< Big 100x100 15\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001660] [solve!:\ttime:0.014863] <<< Big 100x100 16\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001302] [solve!:\ttime:0.016782] <<< Big 100x100 17\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.016191] [solve!:\ttime:0.016815] <<< Big 100x100 18\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.034389] [solve!:\ttime:0.013797] <<< Big 100x100 19\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001158] [solve!:\ttime:0.014430] <<< Big 100x100 20\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.000822] [solve!:\ttime:0.016859] <<< Big 100x100 21\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.000824] [solve!:\ttime:0.017159] <<< Big 100x100 22\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.000770] [solve!:\ttime:0.014640] <<< Big 100x100 23\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006314] [solve!:\ttime:0.014348] <<< Big 100x100 24\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006114] [solve!:\ttime:0.016494] <<< Big 100x100 25\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006066] [solve!:\ttime:0.016940] <<< Big 100x100 26\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006658] [solve!:\ttime:0.016683] <<< Big 100x100 27\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.008010] [solve!:\ttime:0.013922] <<< Big 100x100 28\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001368] [solve!:\ttime:0.015833] <<< Big 100x100 29\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.014416] [solve!:\ttime:0.016858] <<< Big 100x100 30\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006858] [solve!:\ttime:0.017316] <<< Big 100x100 31\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.010688] [solve!:\ttime:0.013077] <<< Big 100x100 32\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.007279] [solve!:\ttime:0.015398] <<< Big 100x100 33\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.015880] [solve!:\ttime:0.016691] <<< Big 100x100 34\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001844] [solve!:\ttime:0.016978] <<< Big 100x100 35\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.007591] [solve!:\ttime:0.015996] <<< Big 100x100 36\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.003007] [solve!:\ttime:0.014691] <<< Big 100x100 37\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.005483] [solve!:\ttime:0.016941] <<< Big 100x100 38\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.005739] [solve!:\ttime:0.017059] <<< Big 100x100 39\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.002455] [solve!:\ttime:0.015538] <<< Big 100x100 40\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.002339] [solve!:\ttime:0.014520] <<< Big 100x100 41\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.000937] [solve!:\ttime:0.016869] <<< Big 100x100 42\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001547] [solve!:\ttime:0.017309] <<< Big 100x100 43\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006841] [solve!:\ttime:0.017265] <<< Big 100x100 44\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.010988] [solve!:\ttime:0.012344] <<< Big 100x100 45\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001839] [solve!:\ttime:0.016139] <<< Big 100x100 46\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.006317] [solve!:\ttime:0.016759] <<< Big 100x100 47\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.001697] [solve!:\ttime:0.016818] <<< Big 100x100 48\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.005746] [solve!:\ttime:0.012656] <<< Big 100x100 49\n",
      "[\u001b[1;32mPASS\u001b[m] [exact:\ttime:0.011035] [solve!:\ttime:0.015238] <<< Big 100x100 50\n",
      ">>> TESTS END\n"
     ]
    }
   ],
   "source": [
    "test(M->try solve!(M.A, M.b) catch; :UNSOLVABLE end, \"solve!\",\n",
    "    M->try \\(M.A, M.b) catch; :UNSOLVABLE end;\n",
    "    equal_f=(v,w)->(v==w==:UNSOLVABLE || abs(norm(v-w))<=epsilon))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.1",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}