{ "cells": [ { "cell_type": "markdown", "id": "485b4afe", "metadata": {}, "source": [ "# Module03: Introduksjon til numpy\n", "\n", "I denne modulen skal vi introdusere Python modulen som kalles *numpy*. Denne modulen er spesielt laget for at python skal bli kunne brukes for å gjøre numeriske beregninger, noe som står sentralt i fysikken. \n", "\n", "Den sentrale typen som numpy er basert rundt er den såkalte *numpy.ndarrays*. Denne typen er en array type, dvs. vektorer, matriser og flerdimensjonale arrayer, som blir brukt til å lagre numeriske data som blir brukt i numeriske beregninger. Denne typen kan *bare* lagre numeriske data av de innebygde typene integers, floats eller complex. I motsetning til hva som er tilfelle for de innebygde Python typende har numpy støtte for *både* single og double presisjons flyttall (kalt float og float64); dette er tilsvarende som for de vanlige kompilerte programmeringsspråkene some C, C++, eller Fortran (samt julia). \n", "\n", "Som en første introduksjon til numpy, vil vi bruk deler av de mange gode [online numpy resursene](https://numpy.org/learn/) som finnes. \n", "\n", "Spesielt vil vi anbefale: \n", "\n", "* [NumPy Quickstart](https://numpy.org/doc/stable/user/quickstart.html) \n", "* [NumPy for beginners](https://numpy.org/doc/stable/user/absolute_beginners.html)\n", "* [NumPy fundamentals](https://numpy.org/doc/stable/user/basics.html)\n", "* [NumPy Reference Manual](https://numpy.org/doc/stable/reference/index.html#reference) \n", " \n", "\n", "Det er *viktig* at du mestrer de fundamentale aspektene rundt numpy da det er en viktig, muligens, den viktigste, byggesteinen for å gjøre numeriske beregninger i python (på en effektiv og forholdvis rask måte). Det finnes mange ulike funksjoner og metoder som er del av numpy. Det er derfor ikke trolig at du på kort sikt vil ha full oversikt over dem alle bare ved å følge forelesningene. Jeg oppdager fortsatt nye aspekter ved denne modulen. Dog, det viktigste er å komme igang, og så lærer man mens man jobber med numpy. I starten er det viktigst å mestre \n", "\n", "* opprettelsen av numpy arrays\n", "* indeksering av elementene i arrayen\n", "* oppdatering av enkelt elementer i arrayen\n", "* reshaping \n", "* kjenne til og å kunne bruke sentrale array funksjoner (f.eks. sum, eig,...)\n", "* ...\n", "\n", "\n", "Vi starter nå forsiktig med å ta for oss deler av innholdet i quickstart guiden. Denne, og de tre andre linkene ovenfor, anbefaler jeg *sterkt* at du går gjennom på egen hånd og prøver å kjøre eksemplene som blir gitt der. Det som blir presentert nedenfor gir bare en liten smakebit på numpy, og vil ikke være tilstrekkelig for å mestre de sentrale elementene av denne ytterst viktige modulen i python for numeriske beregninger. " ] }, { "cell_type": "code", "execution_count": null, "id": "1b57c444", "metadata": {}, "outputs": [], "source": [ "import numpy as np # import statement\n", "\n", "print(np.__version__) # Hvilken versjon av numpy kjører vi\n", "\n", "print( np.__doc__ ) # Hva inneholder numpy modulen" ] }, { "cell_type": "code", "execution_count": null, "id": "827fbcdb", "metadata": {}, "outputs": [], "source": [ "# Alternativt for innholdet i en modul: dir( )\n", "dir( np )\n", "\n", "# Eller, hold muspekeren over np." ] }, { "cell_type": "markdown", "id": "7d0096aa", "metadata": {}, "source": [ "### Deklarere en numpy array" ] }, { "cell_type": "code", "execution_count": null, "id": "9abcef14", "metadata": {}, "outputs": [], "source": [ "a = np.array( [1, 2, 3, 4, 5] ) # Integer \n", "#a = np.array( [1.1,2.2,3.3,4.4,5.5] ) # Default: double presisjon\n", "#a = np.array( [1.1,2.2,3.3,4.4,5.5], dtype=np.float32 ) # Single presisjon; do not use float (uten np.) herS\n", "\n", "print(a)\n", "\n", "print( type( a ) )\n", "print( a.dtype ) # typen til elementene" ] }, { "cell_type": "markdown", "id": "ba94d643", "metadata": {}, "source": [ "Summen og produktet av elementene i en np.array er:" ] }, { "cell_type": "code", "execution_count": null, "id": "61000bab", "metadata": {}, "outputs": [], "source": [ "a_sum = a.sum() # summen av elementene\n", "a_prod = a.prod() # produktet av elementene\n", "\n", "print(\" Summen av elementene i numpy array a er : \", a_sum )\n", "print(\" Produktet av elementene i numpy array a er : \", a_prod) \n" ] }, { "cell_type": "markdown", "id": "6662e831", "metadata": {}, "source": [ "np.array generert fra en liste.......(der var forsåvidt hva vi gjorde ovenfor)" ] }, { "cell_type": "code", "execution_count": null, "id": "d14b5eb6", "metadata": {}, "outputs": [], "source": [ "list = [1, 2, 3, 4, 5]\n", "aa =np.array( list )\n", "\n", "print(aa)\n", "\n", "if np.linalg.norm( a-aa) < 5e-14: # Normen | aa -a | \n", " print(\"Arrayene er LIKE!\") \n", "else: \n", " print(\"Arrayene er ULIKE!\") \n", "# end if" ] }, { "cell_type": "markdown", "id": "20c99aa7", "metadata": {}, "source": [ "eller via en tuple" ] }, { "cell_type": "code", "execution_count": null, "id": "c188532a", "metadata": {}, "outputs": [], "source": [ "# np.array via en liste\n", "a = np.array( [1, 2, 3, 4, 5] )\n", "\n", "# np.array via en tuple\n", "b = np.array( (1, 2, 3, 4, 5 ) )\n", "\n", "print( type(b) ) # b er en numpy.ndarray \n", "\n", "a is b # a og b er ikke det samme objektet\n", "\n", "if np.linalg.norm(b-a) < 5e-14: # Normen | aa -a | \n", " print(\"Arrayene er LIKE!\") \n", "else: \n", " print(\"Arrayene er ULIKE!\") \n", "## end if" ] }, { "cell_type": "code", "execution_count": null, "id": "5578143b", "metadata": {}, "outputs": [], "source": [ "# Du kan gjøre det samme på denne måten!\n", "c = np.array( range(1,6))\n", "print(c)" ] }, { "cell_type": "code", "execution_count": null, "id": "967a7f63", "metadata": {}, "outputs": [], "source": [ "# Merk deg\n", "if np.linalg.norm( aa-c) < 5e-14 : # Normen | aa -c | \n", " print(\"Arrayene er LIKE!\") \n", "else: \n", " print(\"Arrayene er ULIKE!\") \n", "# end if\n", "\n", "\n", "# Merk deg\n", "# ----------------------------\n", "# Hva blir skrevet ut her? \n", "#print( aa is c )\n", "#\n", "#c = aa.copy()\n", "#print( aa is c )\n", "#\n", "#c = aa\n", "#print( aa is c )" ] }, { "cell_type": "markdown", "id": "e730f2ff", "metadata": {}, "source": [ "Mao kan man lett initiere en numpy array via (a) en liste; (b) en tuple; eller (c) en range! \n", "\n", "\n", "I numpy er arrays default av type til inputen til np.array funksjonen. Dog kan man overstyre dette med den valgfrie parameteren *dtype*" ] }, { "cell_type": "code", "execution_count": null, "id": "28963720", "metadata": {}, "outputs": [], "source": [ "v=np.array( range(5) )\n", "print( v.dtype )\n", "\n", "v=np.array( range(5), dtype=np.float64 )\n", "print( v.dtype )\n", "\n", "v=np.array( [0,1,2,3,4+0j] )\n", "print( v.dtype )\n", "\n", "v=np.array( range(5), dtype=np.complex64 )\n", "print( v.dtype )\n" ] }, { "cell_type": "markdown", "id": "aec929f4", "metadata": {}, "source": [ "Hva med matriser eller multi-dimensjonalle (ofte kalt n-dimensjonalle) arrays? Hvordan kan vi definere dem? " ] }, { "cell_type": "code", "execution_count": null, "id": "565411d5", "metadata": {}, "outputs": [], "source": [ "# Metode 1: Som flere vektorer.....\n", "A1 = np.array( [[1,2,3],[4,5,6]] ) # Ev. A1 = np.array( [(1,2,3),(4,5,6)] )\n", "\n", "print(\"A1 = \\n\",A1 )" ] }, { "cell_type": "code", "execution_count": null, "id": "ef0af0a9", "metadata": {}, "outputs": [], "source": [ "# Metode 2: bruk av reshape....\n", "tmp = np.array( [1,2,3,4,5,6] ) \n", "A2 = tmp.reshape( (2,3) )\n", "\n", "print(\"A2 =\\n\",A2 )\n", "\n", "# Ev som en kombinert funksjon\n", "A2_0 = np.array(range(1,7)).reshape( (2,3))\n", "print(\"A2_0 =\\n\",A2_0 )" ] }, { "cell_type": "markdown", "id": "37d56e35", "metadata": {}, "source": [ "Hva skjedde her....?" ] }, { "cell_type": "code", "execution_count": null, "id": "daf42e0c", "metadata": {}, "outputs": [], "source": [ "help( np.reshape ) # for hjelp\n", "\n", "#A2_1 = tmp.reshape( (2,3) )\n", "#print(\"A2_1 =\\n\",A2_1 )\n", "\n", "#A2_2 = tmp.reshape( (2,3), order='F')\n", "#print(\"A2_2 =\\n\",A2_2 )\n" ] }, { "cell_type": "code", "execution_count": null, "id": "17727153", "metadata": {}, "outputs": [], "source": [ "# Numpy's range funksjon; Denne genererer en np.array\n", "v = np.arange(5)\n", "print(v, type(v) )\n", "\n", "# \n", "v = np.arange(1,6)\n", "print(v)\n", "\n", "print( np.arange(1,10,2) )\n" ] }, { "cell_type": "code", "execution_count": null, "id": "66edf9c9", "metadata": {}, "outputs": [], "source": [ "# np.arange fungerer også for flyttall\n", "print( np.arange(5.) )\n", "\n", "print( np.arange(0., 1., 0.1) )" ] }, { "cell_type": "markdown", "id": "b3f554af", "metadata": {}, "source": [ "Det finnes også en nær beslektet numpy funksjon til arange som heter linspace som er nyttig. " ] }, { "cell_type": "code", "execution_count": null, "id": "2793f475", "metadata": {}, "outputs": [], "source": [ "help(np.linspace)" ] }, { "cell_type": "code", "execution_count": null, "id": "d8aca003", "metadata": {}, "outputs": [], "source": [ "v1=np.arange(0.,1.,0.1)\n", "v2=np.linspace(0.,1.)\n", "v3=np.linspace(0.,1, 11)\n", "v4=np.linspace(0.,1, 10, endpoint=False)\n", "\n", "print(\"v1 =\", v1)\n", "print(\"v2 =\", v2)\n", "print(\"v3 =\", v3)\n", "print(\"v4 =\", v4)" ] }, { "cell_type": "markdown", "id": "854d7443", "metadata": {}, "source": [ "Merk deg forskjellen i default oppførsel mellom *np.arange* og *np.linspace* når det gjelder hvordan endepunktet blir behandlet. \n", "\n", "\n", "np.linspace funksjonen brukes mye f.eks. skal plotte data. La oss vise hvordan dette kan gjøres for $\\cos(2\\pi x)$ og $\\sin(2\\pi x)$ for $x\\in[0,1]$. " ] }, { "cell_type": "code", "execution_count": null, "id": "89f6ba44", "metadata": {}, "outputs": [], "source": [ "N=100 # antall punkter\n", "x = np.linspace(0.,1.,N)\n", "Cos = np.cos( 2* np.pi * x )\n", "Sin = np.sin( 2* np.pi * x )\n", "\n", "import matplotlib.pyplot as plt\n", "plt.plot(x,Cos)\n", "plt.plot(x,Sin)\n", "\n", "plt.xlabel(\"x\")\n", "plt.ylabel(r\"$\\cos(2\\pi x)$ and $\\sin(2\\pi x)$\") " ] }, { "cell_type": "markdown", "id": "7197931a", "metadata": {}, "source": [ "### Ulike nyttige kommandoer i numpy" ] }, { "cell_type": "code", "execution_count": null, "id": "cff2f24b", "metadata": {}, "outputs": [], "source": [ "A = np.array(range(1,7)).reshape((2,3)) # Vår matrise\n", "\n", "print( A.ndim ) # antall dimensjoner (her en matrise)\n", "print( A.shape ) # Dimensjonen \n", "print( A.size ) \n", "print( A.shape[0] ) # Antall rekker i matrisen (index 0)\n", "print( A.shape[1] ) # Antall kolonner i matrisen (index 1)" ] }, { "cell_type": "code", "execution_count": null, "id": "42bbede1", "metadata": {}, "outputs": [], "source": [ "# Null matrisen\n", "np.zeros( (2,3) ) # default float64\n", "\n", "#np.zeros( (2,3), dtype=np.int32 )\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9e5164db", "metadata": {}, "outputs": [], "source": [ "# ener matriser\n", "np.ones( (2,3) )\n", "\n", "# fra en maske\n", "B = np.ones_like(A)\n", "print(B)" ] }, { "cell_type": "code", "execution_count": null, "id": "3f61f3bc", "metadata": {}, "outputs": [], "source": [ "# Identitets matrisen (matrisen må være kvadratisk)\n", "np.eye( 3 )" ] }, { "cell_type": "markdown", "id": "4127e6c1", "metadata": {}, "source": [ "Det er mange, mange flere nyttige kommandoer i numpy, men for en beskrivelse av dem henviser vi til ekstern dokumentasjon. Når det er sagt, påpeker vi at det er *vell* verdt å kjenne til samt å kunne bruke disse. Derfor bør du jobbe deg gjennom eksemplene som er gitt i de eksterne linke gitt overfor (eller andre steder). " ] }, { "cell_type": "markdown", "id": "5cb73293", "metadata": {}, "source": [ "#### Verktoriell oppførsel av numpy.ndarray objekter" ] }, { "cell_type": "markdown", "id": "93b490da", "metadata": {}, "source": [ "La oss starte med å genere noen random tall:" ] }, { "cell_type": "code", "execution_count": null, "id": "76905723", "metadata": {}, "outputs": [], "source": [ "help( np.random.random )" ] }, { "cell_type": "code", "execution_count": null, "id": "f01068e5", "metadata": {}, "outputs": [], "source": [ "# Fyller matrisen med tilfeldige (random) tall mellom 0 og 1.\n", "a = np.random.random( (2,3) )\n", "print(\"a =\\n\",a)" ] }, { "cell_type": "code", "execution_count": null, "id": "45bafdd4", "metadata": {}, "outputs": [], "source": [ "# En ny random matrise\n", "b = np.random.random( (2,3) )\n", "print(\"b =\\n\",b)" ] }, { "cell_type": "markdown", "id": "928bad6b", "metadata": {}, "source": [ "Vi kan også sette et såkalt *seed* for den tilfeldige tall generatoren. Dette seeded gjør at vi får generert *samme* sekvens av tall. La oss se på følgende eksempel: " ] }, { "cell_type": "code", "execution_count": null, "id": "8e868aaf", "metadata": {}, "outputs": [], "source": [ "np.random.seed(42)\n", "A = np.random.random( (2,3) )\n", "print(\"A =\\n\",A)\n", "\n", "np.random.seed(42)\n", "B = np.random.random( (2,3) )\n", "print(\"B =\\n\",B)" ] }, { "cell_type": "markdown", "id": "bc77e70e", "metadata": {}, "source": [ "La oss ta sinus av alle elementene i $a$, dvs, S_{ij}=$\\sin(a_{ij})$. Man sier at sin-funksjonen er *verktoriell*.\n", "\n", "Mange funksjoner i numpy har denne egenskapen." ] }, { "cell_type": "code", "execution_count": null, "id": "fd52542e", "metadata": {}, "outputs": [], "source": [ "sina = np.sin(a) # La oss ta sinus av alle elementene i a\n", "print(\"sin(a) = \\n\", sina) " ] }, { "cell_type": "code", "execution_count": null, "id": "ded70550", "metadata": {}, "outputs": [], "source": [ "# Et annet eksempel\n", "A = np.arange(6).reshape((2,3)) \n", "print(\"A = \\n\", A)\n", "\n", "Acumsum = A.cumsum()\n", "print(\"Acumsum = \\n\", Acumsum)\n", "\n", "# Hva skjedde her?\n", "#help(np.cumsum)" ] }, { "cell_type": "code", "execution_count": null, "id": "60e3860a", "metadata": {}, "outputs": [], "source": [ "# La oss ta cumulative sum over kolonner\n", "Acumsum = A.cumsum(axis=1)\n", "print(\"Acumsum = \\n\", Acumsum)" ] }, { "cell_type": "markdown", "id": "876b58af", "metadata": {}, "source": [ "La oss se hvordan vi man bruke numpy til å løse et lineært ligningsystem $Ax=b$. Dette vil vi gjøre på følgende måte\n", "\n", " - Generer en kvadratisk ($N \\times N$) matrise $A$ med random elementer\n", " - Tilsvarende la oss generere en løsnings vektor $y$ som også er tilfeldig\n", " - Høresiden generer vi via matrisemultilikasjon som $b=Ay$ \n", " - Ved hjelp av numpy skal vi så løse systemet Ax=b slik at og undersøke om vi får rett løsning, dvs $x=y$.\n", " " ] }, { "cell_type": "code", "execution_count": 282, "id": "9fa9e677", "metadata": {}, "outputs": [], "source": [ "# Set opp likningsystemet\n", "N=25\n", "A=np.random.random( (N,N) )\n", "y=np.random.random( N ) \n", "b=np.matmul(A,y)\n", "\n", "# print(\"A = \\n\", A)\n", "# print(\"b = \\n\", b)" ] }, { "cell_type": "markdown", "id": "3cf2b517", "metadata": {}, "source": [ "Merk: funksjonen np.random.random() gir uniformt fordelte tall over $[0,1)$. For normal fordelte tall (Gaussisk) må man bruke *np.random.normal*." ] }, { "cell_type": "code", "execution_count": null, "id": "6a79b760", "metadata": {}, "outputs": [], "source": [ "# Hva er rutinen for å løse likningsystemet\n", "help( np.linalg )\n", "\n", "#help( np.linalg.solve )\n", "#help( np.linalg.matmul )" ] }, { "cell_type": "code", "execution_count": null, "id": "e49f2771", "metadata": {}, "outputs": [], "source": [ "# Løs likningsystemet\n", "x = np.linalg.solve(A,b)\n", "\n", "error = np.linalg.norm(y-x)\n", "print( \"error norm = \", error)\n", "\n", "# samelikning av svarene\n", "M=6\n", "for i in range(M):\n", " print(i, \" x_i= \", x[i] ) \n", "# end for" ] }, { "cell_type": "markdown", "id": "ad8595ee", "metadata": {}, "source": [ "**Oppgave**: Genere A og y med normal fordelte tall!" ] }, { "cell_type": "markdown", "id": "f2e50fb3", "metadata": {}, "source": [ "Sentrale rutiner i np.linalg\n", "\n", "- norm : normen av en vektor eller matrise\n", "- matmul : matrise multiplikasjon (ev. matrise vektor multiplikasjon)\n", "- inv : inverse av en matrise\n", "- solve : Løser et linært likningsystem $A\\vec{x}=\\vec{b}$\n", "- eig : Finner egenverdien til en matrise\n", "- det : Determinanten\n", "etc\n", "\n", "for mer detaljer bruk *help(np.linalg)* eller *help(np.linalg.funksjons_navn)*.\n", "\n", "\n", "La oss nå se på hvordan noen av disse rutine kan brukes.... \n", "\n", "Vi starter med å generere en random matrise, med uniformt eller normal fordelte tilfeldige tall. " ] }, { "cell_type": "code", "execution_count": null, "id": "5d2e6890", "metadata": {}, "outputs": [], "source": [ "# Oppsett av likningssystemet\n", "import numpy as np\n", "N = 3 # System dimensjon\n", "RHS = 1 # Antall høyresider\n", "\n", "# Uniformt fordelte tall\n", "A = np.random.random( (N,N) ) \n", "# A_copy = A.copy() # en kopi til bruk senere\n", "print(\"A =\\n\",A)\n", "xx = np.random.random( (N,RHS) )\n", "print(\"xx =\\n\",xx)\n", "\n", "bb =np.zeros((N,RHS))\n", "for i in range(RHS):\n", " bb[:,i] = np.matmul(A,xx[:,i])\n", "# end for \n", "print(\"bb =\\n\",bb)" ] }, { "cell_type": "markdown", "id": "2f476bd8", "metadata": {}, "source": [ "Men $\\vec{bb}$ er det samme som:" ] }, { "cell_type": "code", "execution_count": null, "id": "4602e5e7", "metadata": {}, "outputs": [], "source": [ "b = np.matmul(A,xx)\n", "# Alternativ for matrise multiplikasjon\n", "\n", "# b = A @ xx\n", "print(\"b =\\n\",b)" ] }, { "cell_type": "markdown", "id": "7cbfd32c", "metadata": {}, "source": [ "Men er løsning riktig?" ] }, { "cell_type": "code", "execution_count": 3, "id": "ef21b263", "metadata": {}, "outputs": [], "source": [ "# Vi definerer en funksjon\n", "def feil_funksjon(x,y):\n", " feil = np.linalg.norm(x-y)\n", " return feil\n", "# end feil_funksjon" ] }, { "cell_type": "markdown", "id": "a069da13", "metadata": {}, "source": [ "La oss regne ut løsningen $A \\vec{x} = \\vec{b}$ samt normen til feilen:" ] }, { "cell_type": "code", "execution_count": null, "id": "53473105", "metadata": {}, "outputs": [], "source": [ "x = np.linalg.solve(A,b)\n", "print(\"Feilen er :\", feil_funksjon(x,xx) )" ] }, { "cell_type": "markdown", "id": "bab65af8", "metadata": {}, "source": [ "Dersom antall høyresider $RHS>1$, hva skjer da? Hva gjør *np.linalg.solve* da? \n", "\n", "For å finne det ut, la oss gå tilbake å sette en verdi for $RHS$ som er større enn $1$. " ] }, { "cell_type": "code", "execution_count": 20, "id": "efd8175c", "metadata": {}, "outputs": [], "source": [ "import sys\n", "\n", "# Vi definerer en funksjon\n", "def feil_funksjon_2(x,y):\n", " ndim =x.ndim\n", " if ndim == 1:\n", " feil2 = np.linalg.norm(x-y)\n", " elif ndim == 2:\n", " M = x.shape[1]\n", " feil2 = np.zeros(M) \n", " for i in range(M):\n", " feil2[i] = np.linalg.norm(x[:,i] - y[:,i])\n", " # end for\n", " else: \n", " print(\"Arrayene har ikke konsistent dimensjon!\")\n", " sys.exit()\n", " # end if \n", " return feil2\n", "# end feil_funksjon_2" ] }, { "cell_type": "code", "execution_count": null, "id": "f9f2755f", "metadata": {}, "outputs": [], "source": [ "x = np.linalg.solve(A,b)\n", "print(\"Feilen er :\", feil_funksjon_2(x,xx) )" ] }, { "cell_type": "markdown", "id": "24e9ac1f", "metadata": {}, "source": [ "En enklere løsning er å lese grundig manualsiden til np.linalg.norm" ] }, { "cell_type": "code", "execution_count": null, "id": "827c1806", "metadata": {}, "outputs": [], "source": [ "help(np.linalg.norm)" ] }, { "cell_type": "code", "execution_count": null, "id": "5c0a988f", "metadata": {}, "outputs": [], "source": [ "print(\"norm(xx[:,0]) :\", np.linalg.norm(xx[:,0]))\n", "print(\"norm(xx) :\", np.linalg.norm(xx))\n", "print(\"norm(xx,axis=0) :\", np.linalg.norm(xx, axis=0))" ] }, { "cell_type": "markdown", "id": "111f9b3e", "metadata": {}, "source": [ "La oss definere en funksjon som har en valgbart (opsjonelt) argument axis:" ] }, { "cell_type": "code", "execution_count": 7, "id": "0b7328f6", "metadata": {}, "outputs": [], "source": [ "def feil_funksjon_3(x,y, axis=None):\n", " if axis is None: \n", " feil3 = np.linalg.norm(x-y)\n", " else:\n", " feil3 = np.linalg.norm(x-y, axis=axis)\n", " # end if \n", " return feil3\n", "# end feil_funksjon_3" ] }, { "cell_type": "markdown", "id": "85d182ad", "metadata": {}, "source": [ "Mer hvordan kan du lett kan lage en *docstring* for denne funksjonen ved å gi \"\"\" + return:" ] }, { "cell_type": "code", "execution_count": 9, "id": "cc653c63", "metadata": {}, "outputs": [], "source": [ "def feil_funksjon_3(x,y, axis=None):\n", " \n", " # Gi \"\"\" + return \n", " if axis is None: \n", " feil3 = np.linalg.norm(x-y)\n", " else:\n", " feil3 = np.linalg.norm(x-y, axis=axis)\n", " # end if \n", " return feil3\n", "# end feil_funksjon_3" ] }, { "cell_type": "code", "execution_count": null, "id": "17424d0e", "metadata": {}, "outputs": [], "source": [ "x = np.linalg.solve(A,b)\n", "print(\"shape(x) :\", x.shape)\n", "print(\"shape(xx) :\", xx.shape)\n", "print(\"Feilen er :\", feil_funksjon_3(x,xx,axis=0) )\n", "print(\"Feilen er (men hva betyr dette):\", feil_funksjon_3(x,xx) )\n", "\n", "# Enklere syntax\n", "print(\"Feilen er (enklere syntax):\", np.linalg.norm(x-xx,axis=0) )" ] }, { "cell_type": "markdown", "id": "e8236bd5", "metadata": {}, "source": [ "Moralen her er, les DOKUMENTASJONENE!!!" ] }, { "cell_type": "markdown", "id": "a2110885", "metadata": {}, "source": [ "**OPPGAVE**: gjennomfør den samme beregningen, men ved å bruke normal fordelte tall for matrisene! " ] }, { "cell_type": "markdown", "id": "269bb85a", "metadata": {}, "source": [ "Man kan også løse likningssystemet $A\\vec{x}=\\vec{b}$ ved å regne den inverse matrisen til $A$ og matrise multiplisere $\\vec{b}$ med denne, dvs.\n", "\\begin{align*}\n", " \\vec{x}\n", " &=\n", " A^{-1} \\vec{b}.\n", "\\end{align*}\n", "\n", "La oss nå teste dette:" ] }, { "cell_type": "code", "execution_count": null, "id": "306c176f", "metadata": {}, "outputs": [], "source": [ "x = np.linalg.inv(A) @ b \n", "# x = np.matmul( np.linalg.inv(A), b ) # Alternativ\n", "print(\"shape(x) :\", x.shape )\n", "print(\"Feilen er (men hva betyr dette):\", np.linalg.norm(x-xx,axis=0))" ] }, { "cell_type": "markdown", "id": "34524f0c", "metadata": {}, "source": [ "### Egenverdier\n", "\n", "En egenverdi likning er definert slik at $A\\vec{v} = \\lambda \\vec{v}$ hvor $\\vec{v}$ ikke er en null vektor. Spørsmålet er hva disse er $\\vec{v}$ og hva er $\\lambda$. Dette er et standard problem i lineær algebra, så vi vill ikke gå i detalj her. " ] }, { "cell_type": "code", "execution_count": null, "id": "abbadd36", "metadata": {}, "outputs": [], "source": [ "help(np.linalg.eigvals)\n", "help(np.linalg.eig)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1a8f34db", "metadata": {}, "outputs": [], "source": [ "# eigenvalues\n", "w = np.linalg.eigvals(A)\n", "print(\"Eigenvalues (lambda) :\", w)" ] }, { "cell_type": "code", "execution_count": null, "id": "571f087f", "metadata": {}, "outputs": [], "source": [ "# eigenverdier og tilhørende eigenvektorer\n", "w,v = np.linalg.eig(A)\n", "print(\"Eigenvalues (lambda) :\\n\", w)\n", "print(\"Eigenvektorer :\\n\", v)" ] }, { "cell_type": "markdown", "id": "ecce1700", "metadata": {}, "source": [ "La oss teste resultatet $A\\vec{v} = \\lambda \\vec{v}$ " ] }, { "cell_type": "code", "execution_count": null, "id": "5485f0a7", "metadata": {}, "outputs": [], "source": [ "for i in range(w.size):\n", " vektor = A @ v[:,i] - w[i] * v[:,i]\n", " error = np.linalg.norm( vektor )\n", " print( \"vektro : \", vektor)\n", " print( \"error (\",i,\") :\", error)\n", "# end for " ] }, { "cell_type": "markdown", "id": "7ba64483", "metadata": {}, "source": [ "#### Elementer i en numpy.ndarray" ] }, { "cell_type": "code", "execution_count": null, "id": "a0ccffa0", "metadata": {}, "outputs": [], "source": [ "v = np.array( np.arange(6) )\n", "A = (np.array( np.arange(9) )).reshape( (3,3) )\n", "\n", "print( v[0] ) # først element\n", "print( v[1:3] ) # andre til fjerde element\n", "#print( v[:-2] ) # Alle elementer foruten de to siste\n", "#print( v[-2:] ) # Alle elementer foruten de TRE første\n", "\n", "#print( A[1,2] ) # etc...\n", "#print( A[:,2] ) # 3. kolonne" ] }, { "cell_type": "code", "execution_count": null, "id": "6582a759", "metadata": {}, "outputs": [], "source": [ "print( \"A = \\n\",A)" ] }, { "cell_type": "markdown", "id": "40c582d5", "metadata": {}, "source": [ "### Determinater" ] }, { "cell_type": "code", "execution_count": null, "id": "95be6f78", "metadata": {}, "outputs": [], "source": [ "help(np.linalg.det)" ] }, { "cell_type": "code", "execution_count": null, "id": "18bf9391", "metadata": {}, "outputs": [], "source": [ "AA=np.array( [[ 1, 2],\n", " [3,4]])\n", "print( \" AA = \\n\", AA)" ] }, { "cell_type": "code", "execution_count": null, "id": "4c03bbd1", "metadata": {}, "outputs": [], "source": [ "# determinanten til denne matrisen er \n", "AAdet = np.linalg.det(AA)\n", "print(\" |AA| = \", AAdet)" ] }, { "cell_type": "code", "execution_count": null, "id": "7da44a04", "metadata": {}, "outputs": [], "source": [ "print( \" A = \\n\", A, \"\\n\")\n", "Adet = np.linalg.det(A)\n", "print(\" |A| = \", Adet)" ] }, { "cell_type": "markdown", "id": "4c5aa90b", "metadata": {}, "source": [ "#### Kopiering av numpy.ndarray" ] }, { "cell_type": "code", "execution_count": null, "id": "36991a95", "metadata": {}, "outputs": [], "source": [ "v = np.array( np.arange(6) )\n", "\n", "# Intet NYTT Object lages\n", "a = v\n", "\n", "# Nytt object lages (deep copy)\n", "b = v.copy()\n", "c = v[:]FY1008_M04.pdf\n", "\n", "#print( id(v),id(a), id(b), id(c) )\n" ] }, { "cell_type": "markdown", "id": "1444a567", "metadata": {}, "source": [ "Tilsvarende som vi så tidligere, så må vi være forsiktige når vi kopierer np.arrays." ] }, { "cell_type": "markdown", "id": "d7d1005f", "metadata": {}, "source": [ "### Operasjoner på numpy.ndarray" ] }, { "cell_type": "code", "execution_count": null, "id": "facb31e1", "metadata": {}, "outputs": [], "source": [ "# definerer to np.arrays\n", "a = np.array( [1,2,3] )\n", "b = np.array( [1,4,6] )\n", "\n", "c = a + b # sum (element vis)\n", "print (\"c = \", c ) # samme for suptraksjon\n", "\n", "d = a * b # multiplikasjon (element vis)\n", "print (\"d = \", d ) # samme for divisjon (mer resultat er float)\n", "\n", "dot = np.dot(a,b)\n", "print(\"dot = \", dot) # dot-produkt eller skalar produkt\n", "\n", "cross = np.cross(a,b) # cross produktet\n", "print(\"cross = \", cross )\n" ] }, { "cell_type": "markdown", "id": "da946b85", "metadata": {}, "source": [ "Man kan også enkelt lese data from en file i ASCII format." ] }, { "cell_type": "code", "execution_count": null, "id": "923877dd", "metadata": {}, "outputs": [], "source": [ "filename=\"https://web.phys.ntnu.no/~ingves/Teaching/FY1008/Downloads/Data/square.dat\"\n", "\n", "square=np.loadtxt(filename)\n", "\n", "print(\"Data som var lest er : \\n\", square )" ] }, { "cell_type": "code", "execution_count": null, "id": "22afc553", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.plot(square[:,0], square[:,1])\n", "plt.title(\"Graf av innholdet i data filen!\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "6a08f4c5", "metadata": {}, "source": [ "For å lagre data som ASCII finnes det en tilsvarende routine np.savetxt som kan være nyttig." ] }, { "cell_type": "code", "execution_count": null, "id": "23706d32", "metadata": {}, "outputs": [], "source": [ "np.savetxt( \"/tmp/kast.dat\", square )\n", "# \n", "\n", "# Finnes filen? (%ls er en \"magic-command\")\n", "%ls -l /tmp/kast.dat\n", "\n", "# Hva inneholder denne filen?\n", "#%cat /tmp/kast.dat\n", "\n", "# Prøv også dette; Den lagrede filen er komprimert (og binær) \n", "np.savetxt( \"/tmp/kast.dat.gz\", square )" ] }, { "cell_type": "code", "execution_count": null, "id": "244baad3", "metadata": {}, "outputs": [], "source": [ "# Man kan også lagre data vha np.save, men formatet er et binært internt Python format (*.npy *.npz). \n", "help(np.load)" ] }, { "cell_type": "markdown", "id": "985df208", "metadata": {}, "source": [ "## Effektivitet av numpy.ndarray objekter" ] }, { "cell_type": "markdown", "id": "e682fcb5", "metadata": {}, "source": [ "Vi har argumentert for at numpy er en ytterst god module når man skal gjøre numeriske beregninger i Python. Men, hva betyr dette, og hvor mye bedre ytelse vil man typisk kunne få dersom man bruker numpy.ndarrays for sine beregninger. Dette vil vi kort prøve å svare på nedenfor med noen enkle eksempler. " ] }, { "cell_type": "code", "execution_count": null, "id": "f6f370a9", "metadata": {}, "outputs": [], "source": [ "# Beregning med lister\n", "import time\n", "import math\n", "\n", "# Antall elementer\n", "N = 100000\n", "\n", "# Sett starttid\n", "time0 = time.time()\n", "\n", "# Gjør beregningen\n", "Sin = [] # Starter med en tom liste\n", "X2 = []\n", "for x in range(N):\n", " Sin.append( math.sin(x) )\n", " #X2.append(x*x) \n", "# end for\n", "\n", "# Sett sluttid\n", "time1 = time.time()\n", "\n", "time_list = time1-time0\n", "\n", "print(\"Beregningen med lister tok (s) : \",time_list)" ] }, { "cell_type": "code", "execution_count": null, "id": "4d58eb86", "metadata": {}, "outputs": [], "source": [ "# Sett starttid\n", "time0 = time.time()\n", "\n", "# Beregning\n", "x=np.arange(N)\n", "Sin= np.sin(x)\n", "#X2 = x*x\n", "\n", "# Sett sluttid\n", "time1 = time.time()\n", "\n", "time_numpy = time1-time0\n", "\n", "print(\"Beregningen med numpy tok (s) : \",time_numpy)\n", "\n", "print(\"Speedup (time_list/ time_numpy) : \", time_list/time_numpy)" ] }, { "cell_type": "markdown", "id": "c163f7c0", "metadata": {}, "source": [ "Merk deg at disse eksemplene er bare illustrasjoner. For andre beregninger kan forskjellen være både større og mindre! Videre vil tiden en beregning tar kunne variere fra kjøring til kjøring!\n", "\n", "\n", "Her er en annen sammelikning (inspirert av https://www.geeksforgeeks.org/python/python-lists-vs-numpy-arrays/):" ] }, { "cell_type": "code", "execution_count": null, "id": "71efb697", "metadata": {}, "outputs": [], "source": [ "\n", "# size of arrays and lists\n", "size = 1000000\n", "\n", "# declaring lists\n", "list1 = range(size)\n", "list2 = range(size)\n", "\n", "# declaring arrays\n", "array1 = np.arange(size)\n", "array2 = np.arange(size)\n", "\n", "#\n", "# Using lists:\n", "# --------------------------\n", "#\n", "# capturing time before the multiplication of Python lists\n", "initialTime = time.time()\n", "\n", "# multiplying elements of both the lists and stored in another list\n", "resultantList = [(a * b) for a, b in zip(list1, list2)]\n", "\n", "# calculating execution time\n", "dt_list = time.time() - initialTime\n", "print(\"Time taken by Lists to perform multiplication:\",\n", " (dt_list),\"seconds\")\n", "\n", "#\n", "# Using numpy\n", "# --------------------------\n", "#\n", "# capturing time before the multiplication of Numpy arrays\n", "initialTime = time.time()\n", "\n", "# multiplying elements of both the Numpy arrays and stored in another Numpy array\n", "resultantArray = array1 * array2\n", "\n", "# calculating execution time\n", "dt_numpy = time.time() - initialTime\n", "print(\"Time taken by NumPy Arrays to perform multiplication:\",\n", " (dt_numpy), \"seconds\")\n", "\n", "# Speedup\n", "print(\"Speedup : \", dt_list/dt_numpy)" ] }, { "cell_type": "markdown", "id": "65251b71", "metadata": {}, "source": [ "### Eksempel med bruk av numpy.array objekter" ] }, { "cell_type": "markdown", "id": "d7d2423b", "metadata": {}, "source": [ "La oss se på problemet fra Øving 3, hvor en ball starter fra bakkenivå $z=z_0=0$ m med hastigheten $v_0=5$ m/s rett opp i vertikal retning. Den blir bare utsatt for tyngde kreften og vi ble bedt om å finne fremtidig posisjon (og hastighet). Hvordan kan vi løse denne oppgaven ved å bruke numpy.ndarray objekter (istedet for f.eks. lister)? " ] }, { "cell_type": "code", "execution_count": 294, "id": "92ee5327", "metadata": {}, "outputs": [], "source": [ "import scipy as sp\n", "from scipy import constants # for å få konstanten g \n", "\n", "g = constants.g\n", "\n", "def euler(u,v,dt): # Hva med doc streng her?\n", " unew = u + dt*v\n", " vnew = v - dt * g \n", " return unew, vnew\n", "# end euler" ] }, { "cell_type": "code", "execution_count": null, "id": "95d41c1f", "metadata": {}, "outputs": [], "source": [ "import time\n", "\n", "dt = 0.000005 # tidssteg (Prøv med e.g. dt = 0.000005 )\n", "t0 = 0 # start tidspunkt\n", "u0 = 0 # start betingelser\n", "v0 = 5\n", "\n", "U = u0 * np.ones((1,)) # for lagring av resultater\n", "V = v0 * np.ones((1,))\n", "tid = t0 * np.ones((1,))\n", "\n", "u = u0\n", "v = v0 \n", "t = t0\n", "\n", "# Løser diff likningen\n", "start = time.time()\n", "while u>=0:\n", " u,v = euler(u,v,dt)\n", " t += dt\n", " # Dette er fristende og det fungerer MEN er ikke effektivt\n", " tid = np.append(tid, t)\n", " U = np.append(U,u)\n", " V = np.append(V,v)\n", "# end while\n", "stopp = time.time()\n", "tid_cpu = stopp-start\n", "print(\"Tidsforbruk (s): \", tid_cpu)\n", "\n", "\n", "# Plotter resultatene (bør vi definere dette i en funksjon?)\n", "import matplotlib.pyplot as plt\n", "\n", "plt.plot(tid,U)\n", "plt.xlabel(\"t [s]\")\n", "plt.ylabel(\"posisjon [m]\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "76c1d5f5", "metadata": {}, "source": [ "Problemet er at *np.append()* funksjonen generer et nytt *numpy.ndarray* objekt hver gang. Siden denne funksjonen er inne i en loop, er ikke dette effektivt da det blir mye kopieringer og minne allokering." ] }, { "cell_type": "code", "execution_count": null, "id": "2b003549", "metadata": {}, "outputs": [], "source": [ "# For dokumentasjon\n", "help( np.append )" ] }, { "cell_type": "markdown", "id": "6bd82313", "metadata": {}, "source": [ "En alternativ metode er å allokere alle np.arrays før man går inn i loopen, men man kjenner ikke hvor stor disse arrayene skal være. En mulig løsning som er mer effektiv blir presentert nedenfor, selv om noen vil hevde at den ikke er spesielt elegant." ] }, { "cell_type": "code", "execution_count": null, "id": "4927bbcd", "metadata": {}, "outputs": [], "source": [ "\n", "N = 500 # ev. N = int(1.3/dt)\n", "U = np.zeros(N) # for lagring av resultater\n", "V = np.zeros(N)\n", "tid = np.zeros(N)\n", "\n", "u = u0; U[0] = u0\n", "v = v0; V[0] = v0 \n", "t = t0; tid[0] = t0\n", "n = 0\n", "\n", "# Løser diff likningen\n", "start = time.time()\n", "while u>=0:\n", " u,v = euler(u,v,dt)\n", " n += 1\n", " t = n * dt\n", " # Må vi utvide numpy arrayene?\n", " if n>N-1:\n", " null = np.zeros((N,) )\n", " tid = np.append( tid, null )\n", " U = np.append( U, null )\n", " V = np.append( V, null )\n", " N = N+N\n", " # end if\n", " # Her lagrer vi resultatene.....\n", " tid[n] = t\n", " U[n] = u\n", " V[n] = v \n", "# end while\n", "\n", "# Oppdaterer arrays (ved å fjerne hva man ikke har brukt) \n", "U = U[0:n-1]\n", "V = V[0:n-1]\n", "tid = tid[0:n-1]\n", "\n", "stopp = time.time()\n", "\n", "# Rapporter tidsbruk\n", "print(\"Tidsforbruk (s): \", stopp-start)\n", "print(\"Speedup :\", tid_cpu/(stopp-start))\n", "\n", "\n", "# Plotter resultatene \n", "import matplotlib.pyplot as plt\n", "\n", "plt.plot(tid,U)\n", "plt.xlabel(\"t [s]\")\n", "plt.ylabel(\"posisjon [m]\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "45aa6dc9", "metadata": {}, "source": [ "Når vi nå har posisjon og hastighet kan vi lett regne ut kinetisk og potensiell energi for ballen. Dette gjør vi å definere en funksjon. [Merk deg hvordan doc strengen kan lages i code]" ] }, { "cell_type": "code", "execution_count": null, "id": "48164538", "metadata": {}, "outputs": [], "source": [ "def energi_beregning(U,V,m=1):\n", " \"\"\"\n", " Mekanisk energi for en ball som beveger seg vertikalt i tyngdefeltet. \n", " Som default regner man ut mekanisk energi for en enhets masse (m=1kg), \n", " eller om man vil, energi per masse enhet. \n", "\n", " Args:\n", " U (ndarray): posisjon som funksjon av tiden\n", " V (ndarray): hastighet som funksjon av tiden\n", " m (int, optional): ballens masse. Defaults to 1.\n", "\n", " Returns:\n", " Etot (ndarray): Total mekanisk energi som funksjon av tiden\n", " Ekin (ndarray): Kinetisk energi som funksjon av tiden\n", " Epot (ndarray): Potensiell energi som funksjon av tiden \n", " \"\"\"\n", " # Kinetisk energi\n", " Ekin= 0.5 * m * V**2\n", " # Potensiell energi\n", " Epot = m * g * U\n", " # Total energi\n", " Etot = Epot + Ekin\n", " #\n", " return Etot, Ekin, Epot" ] }, { "cell_type": "markdown", "id": "60d69ff7", "metadata": {}, "source": [ "Vi bruker nå denne funksjonen for å beregne energien for ballen og plotter resultatene." ] }, { "cell_type": "code", "execution_count": null, "id": "6f1bd510", "metadata": {}, "outputs": [], "source": [ "# Beregner energien\n", "Etot, Ekin, Epot = energi_beregning(U,V)\n", "\n", "# Plotter energien\n", "plt.plot(tid,Etot,label='total')\n", "plt.plot(tid,Ekin,label='kinetisk')\n", "plt.plot(tid,Epot,label='potensiell')\n", "\n", "plt.legend()\n", "plt.xlabel(\"tid (s)\")\n", "plt.ylabel(\"energi/masse (J/kg)\")\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "1ef00afb", "metadata": {}, "source": [ "### En siste oversikt over NumPy\n", "\n", "\n", "En oppsummering av NumPy kan du finne [her](https://colab.research.google.com/github/QuantEcon/lecture-py-notebooks/blob/master/numpy.ipynb). Det kan være lurt å gå gjennom dette materialet for å oppfriske noe av de sentrale egenskapene med NumPy!\n", "\n", "Bruk NumPy hvor du kan, da dette typiske r raskt for numeriske beregninger. \n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.2" } }, "nbformat": 4, "nbformat_minor": 5 }