26 czerwca 2016

Rachunek wektorowy i książki

Kilka przydatnych książek, które znalazłem na wikibooks, ucząc się fizyki. Niektórym brakuje rozdziałów, ale i tak jestem pod wrażeniem. Warto się temu przyjrzeć, jako coś dodatkowego do klasycznych podręczników:
Coś z rachunku wektorowego - wirujący wektor, do którego stworzenia potrzebna jest macierz obrotu [1], [2].
Program koduj macierze przekształceń właściwą dla dwóch wymiarów:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

from visual import *
import ImageGrab

display(x=0, y=0, width=400, height=400, userzoom=False,
        center=(0, 0, 1), foreground=(0, 0, 0), background=(1, 1, 1))

ar = arrow(pos=(0, 0, 0), axis=(10, 0, 0), shaftwidth=0.3)

total_angle = 0
angle = 0.1 / (2 * math.pi)
frame = 0
while total_angle <= 2 * math.pi:
    rate(100)

    nx = ar.axis.x * math.cos(angle) - ar.axis.y * math.sin(angle)
    ny = ar.axis.x * math.sin(angle) + ar.axis.y * math.cos(angle)
    nz = 0
    ar.axis = (nx, ny, nz)

    file_name = 'img-' + '{fr:03d}'.format(fr=frame) + '.png'
    frame += 1

    im = ImageGrab.grab((0, 0, 400, 400))
    im.save(file_name)

    total_angle += angle
    print total_angle

exit()

24 czerwca 2016

[python] matplotlib - rysowanie wykresów

Kolejna biblioteka, którą zacząłem się bawić, służąca do rysowania wykresów: matplotlib. Bardzo przyjema dokumentacja, z dużą liczbą przykładów. Dla testu wykreśliłem pochodną i całkę oznaczoną dla funkcji:
$$ f(x) = -2x^{3} - 4x^{2} + 8x + 1 \\
f^\prime(x) = -6x^{2} - 8x + 8 \\
F(x) = \int f(x) dx = -\frac{1}{2}x^{4} - \frac{4}{3}x^{3} + 4x^{2} + x + C \\
F(-3) - F(1) = \int_{-3}^{1} f(x) dx = \left(-\frac{1}{2}x^{4} - \frac{4}{3}x^{3} + 4x^{2} + x + C\right)\Biggr|_{-3}^{1} $$
Kod programu:
#!/usr/bin/env python
# -*- coding: utf-8 -*-


import numpy as np
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt


def main():
    red_patch = mpatches.Patch(color='red', label=r'$f(x)$ - moja funkcja')
    blue_patch = mpatches.Patch(color='blue', label=r'$f^\prime(x)$ - obliczone numerycznie')
    cyan_patch = mpatches.Patch(color='cyan', label=r'$f^\prime(x)$ - obliczone recznie')
    green_patch = mpatches.Patch(color='green', label=r"$\int_{-3}^{1}f(x)$ - obliczone numerycznie")

    plt.legend(handles=[red_patch, blue_patch, cyan_patch, green_patch], loc='lower right')

    x = np.linspace(-4, 2, 10)
    x_dense = np.linspace(-4, 2)

    y_f = func(x)
    plt.plot(x, y_f, 'red', linewidth=2)

    y_deriv = derivative(func, x)
    plt.plot(x, y_deriv, 'blue', linewidth=2)

    y_own_deriv = own_deriv_func(x_dense)
    plt.plot(x_dense, y_own_deriv, 'cyan', linewidth=1)

    summo1 = definite_integral(func, a=-3, b=1)
    summo2 = own_integral_func(a=-3, b=1)
    print 'Całka oznaczona, numerycznie: %f' % summo1
    print 'Całka oznaczona, ręcznie: %f' % summo2

    current_figure = plt.gcf()
    current_figure.savefig('rachunek_rozniczkowy.png')

    plt.show()


def func(x):
    return -2 * (x ** 3) - 4 * (x ** 2) + (8 * x) + 1


def own_deriv_func(x):
    return -6 * (x ** 2) - 8 * x + 8


def own_integral_func(a, b):
    F = lambda x: -(1/2.0) * (x ** 4) - (4/3.0) * (x ** 3) + 4 * (x ** 2) + x
    return F(b) - F(a)


def derivative(fun, x):
    h = 0.2 # dx
    return (fun(x + h) - fun(x)) / h


def definite_integral(fun, a, b):
    axes = plt.gca()
    dx = (b - a) / 20.0

    summo = 0
    x = a
    while x < b:
        axes.add_patch(mpatches.Rectangle(xy=(x, 0), width=dx, height=fun(x), facecolor='green'))
        summo += dx * fun(x)
        x += dx

    return summo


if __name__ == '__main__':
    main()
Wynik:
Całka oznaczona, numerycznie: -26.080000
Całka oznaczona, ręcznie: -25.333333

19 czerwca 2016

MathJax - Hello World

Test MathJax do zapisywania formuł matematycznych.
Niech \(f\) będzie funkcją ciągła zdefiniowaną dla \( a \leq x \leq b \). Przedział \( [a, b] \) będzie podzielony \( n \) podprzedziałów o równej długości wynoszącej \( \Delta x = \frac{(b - a)}{n} \), tak że \(x_{0} = a \), a \( x_{n} = b \), będącymi końcami tych podprzedziałów i niech \( x_{1}^*, x_{2}^*,... , x_{n}^* \) będą dowolnymi punktami próbkującymi w tych przedziałach, \( x_{i}^* \in [x_{i-1}, x_{i}] \). Wówczas całkę oznaczoną z funkcji \( f \) w przedziale od \( a \) do \( b \) oznaczamy i definiujemy jako:
$$ \int_{a}^{b} f(x) dx = \lim_{n\to\infty} f(x)\sum_{i=1}^{n} f(x_{i}^*) \Delta x $$

3 maja 2016

{Machine Learning} Recipes

Google Developers na swoim kanale, rozpoczął serię filmów traktujących o nauczaniu maszynowym. Krótka i ciekawa forma, bardzo mi się podoba.



Odrobinę przerobiony przykład "Hello Wrold"
# http://scikit-learn.org/stable/
# https://www.continuum.io/
# sudo apt-get install python-scikits-learn

from sklearn import tree

# Collect training data
SMOOTH = 1
BUMPY = 0
APPLE = 0
ORANGE = 1

features = [[140, SMOOTH], [130, SMOOTH], [150, BUMPY], [170, BUMPY]]
labels = [APPLE, APPLE, ORANGE, ORANGE]

# Train classifier
clf = tree.DecisionTreeClassifier()
clf = clf.fit(features, labels)

# Make predictions
result = clf.predict([[160, 0]])

if result[0] == APPLE:
    print 'APPLE'
elif result[0] == ORANGE:
    print 'ORANGE'
else:
    'Result unknown'
Wynik:
ORANGE

28 marca 2016

vpython, gify i mechanika klasyczna

Od dawna planuje zabrać się za naukę mechaniki klasycznej. To kolejne podejście, tym razem jednak, zacząłem od poszukiwania narzędzi, które pomogą mi w nauce. Zaczęło się od bardzo fajnego artykułu na temat biblioteki vpython, linki poniżej i inne przydatne:
Za cel eksperymentu postawiłem sobie, wygenerowanie gif-a, z animacją ruchu. Liczyłem, że wystarczy do tego ImageMagic, jednak tworzenie gif-ów jest na tyle zagmatwane, że trzeba było skorzystać również z pakietu FFmpeg. Efekt nie jest do końca satysfakcjonujący. Od pewnego czasu, za bardzo skupiam się nad idealną wersją prototypu. Co więcej staram się znaleźć odpowiedzi na wszystkie pytania jakie się pojawią podczas jego wytwarzania, szczególnie te poboczne. Nie taka był idea tego bloga, czas więc na refleksję.

W każdym razie, oto przerobiony przykład "bouncing ball" z niewielkimi modyfikacjami. Po pierwsze, skorzystałem z biblioteki PIL (Python Image Library), aby tworzyć zrzuty ekranu kolejnych klatek animacji. Biblioteka działa jedynie pod systemem Windows, a jej zamiennik pyscreenshot nie zdał rezultatu. Jakoś się z tym pogodziłem. Druga sprawa, to warunek przerwania animacji, czyli moment gdy piłka znów znajdzie się w tej samej pozycji.
from visual import *
import ImageGrab    # from PIL


starting_height = 4
floor = box (pos=(0, 0, 0), length=4, height=0.5, width=4, color=color.blue)
ball = sphere (pos=(0, starting_height, 0), radius=1, color=color.red)
ball.velocity = vector(0, -1, 0)
dt = 0.01

frame = 0
while 1:
    rate (100)
    ball.pos = ball.pos + ball.velocity*dt
    if ball.y < ball.radius:
        ball.velocity.y = abs(ball.velocity.y)
    else:
        ball.velocity.y = ball.velocity.y - 9.8*dt

    file_name = 'img-' + '{fr:03d}'.format(fr=frame)  + '.png'
    frame += 1

    im = ImageGrab.grab((0, 0, 500, 500))  # screen box from (0,0)-(500,500)
    im.save(file_name)                     # save image to disk

    if ball.pos.y > starting_height:
        exit()
Skrypt wygenerował około 140 zrzutów ekranu, które trzeba było poddać obróbce, nieoceniony okazał się tutaj ImageMagic. Kluczowym problem jest wysterowanie czasu pomiędzy kolejnymi klatkami.
Wszystkie obrazki zawierały ramki Windows-owego okienka, które trzeba było usunąć:
Przycinanie zdjęcia za pomocą mogrify (z ImageMagic). Pozycje i rozmiar dobrałem eksperymentalnie.
mogrify -crop 413x411+9+30 +repage $(ls -v *.png)
Chociaż ImageMagic potrafi generować gif-y, jednak za żadne skarby nie mogłem skonfigurować właściwego opóźnienia, pomiędzy kolejnymi klatkami. Szperając po internecie udało się ten problem ominąć, wymaga to jednak kroków pośrednich, w postaci utworzenia filmiku z animacją, a następnie przepuszczenia go jeszcze raz przez ffmpeg i convert (ImageMagic). Magiczny zestaw komend:
# Filmik ze zdjęć
ffmpeg -r 100 -f image2 -i img-%3d.png video.avi

# Wyciąganie klatek z filmiku i zlepianie w gif-a
ffmpeg -i video.avi -vf scale=320:-1 -r 10 -f image2pipe -vcodec ppm - | convert -delay 10 -loop 0 - gif:- | convert -layers Optimize - output.gif

# Alternatywna wersja, też działa
ffmpeg -i video.avi -vf scale=320:-1 -r 10 -f image2pipe -vcodec ppm - | convert -delay 10 -loop 0 - output2.gif
W animacji jest jakiś uszczerbek jakości, tak jakby brakowało wszystkich potrzebnych klatek. Być może uda mi się kiedyś rozwiązać tą zagadką i odpowiednio dobrać parametry, wtedy też o nich napiszę. Próbowałem także innych rozwiązań, bezskutecznie próbując wysterować opóźnienie pomiędzy kolejnymi klatkami.
# Gif generowany tylko przez convert (ImageMagic)
convert -dispose none -delay 1.5 $(ls -v *.png) -coalesce output3.gif

# Gif generowany tylko przez ffmpeg
ffmpeg -f image2 -framerate 100 -i img-%3d.png output4.gif

16 marca 2016

Klient OpenVPN - pierwsze kroki

Postawienie dobrze zabezpieczonego połączenia VPN, jest jak się okazuje czynnością dość pracochłonną. Szczególnie warty uwagi jest pierwszy link, który opowiada o dobrym zabezpieczeniu całej usługi:
Moje pierwsze kroki, sprowadziły się do rozwiązania kilku problemów, związanych z systemem operacyjnym (Ubuntu) jak i usługodawcą VPN. Pierwsza próba odpalenia zakończyła się błędem:
Wed Mar 16 10:10:10 2016 ERROR: Cannot ioctl TUNSETIFF tap: Operation not permitted (errno=1)
Wed Mar 16 10:10:10 2016 Exiting due to fatal error
A więc klient, musi w jakiś sposób uzyskać prawa super użytkownika, aby działać na urządzeniu tun/tab, na co jest kilka rozwiązań:
Najprościej jest oczywiście, gdy posiadamy prawa root-a. Możemy wtedy wykonać komendę (gdzie w pliku secret.txt, w kolejnych linijkach zapisany jest login i hasło):
openvpn --config vpn_provider/conf.ovpn --ca vpn_provider/ca.crt \
        --auth-nocache --auth-user-pass secret.txt
Kolejny problem, z którym się spotkałem to niedziałające DNS-y. Czasami usługodawca VPN nie oferuje takiej usługi ze swojej strony, a ISP nie zgadza się na łączenie ze swoimi serwerami DNS, dla osób spoza sieci. DNS-owe zapytania, można przerzucić na router (192.168.1.1), co może nie być najlepszym pomysłem ze względów bezpieczeństwa, można też skorzystać z publicznych serwerów DNS (np. google). Nie wiem jak dokładnie wybierane są serwery DNS-owe z pliku /etc/resolve.conf, przez aplikacje, ale to co zaobserwowałem (przynajmniej w stosunku do komendy ping), warto wpisy o nowych serwerach podać jako pierwsze.
# Dynamic resolv.conf(5) file for glibc resolver(3) generated by resolvconf(8)
#     DO NOT EDIT THIS FILE BY HAND -- YOUR CHANGES WILL BE OVERWRITTEN
nameserver 8.8.8.8
namseerver 8.8.4.4