Pythona zacząłem używać jakieś 5 lat temu. Na początku nie byłem przekonany, ale nie miałem wyboru. Z czasem się przekonałem i zrobiłem kilka projektów małych i dużych. Dodatkowo uczyłem podstaw programowania w pythonie przez 4 lata – do tej pory nie jestem przekonany czy jest to dobry język na początek (trochę za dużo wybacza), ale nie była to moja decyzja.

Wracając do tematu: python ma jedną poważną wadę – jest strasznie wolny. Dziś chcę wam pokazać na przykładzie, jak w prosty sposób można użyć funkcji napisanych w C/C++ do przeprowadzania obliczeń w programie w napisanym w pythonie. Python jest rewelacyjny jeśli chodzi o wczytywanie, zapisywanie i wizualizację danych, z kolei C/C++ jest bardzo wydajny. Połączenie wydaję się idealne i okazuje się całkiem proste.

Dla przykładu weźmy proste zadanie: generujemy listę N liczb losowych z zakresu od 0 do 1. Chcemy obliczyć wszystkie możliwe sumy 3 liczb z tej listy bez powtórzeń i policzyć ile z nich wynosi 1 z pewną ustaloną dokładnością. Oczywiście można ten problem optymalizować, ale dziś nie o tym. Chcemy porównać wydajność czystego pythona i pythona wspieranego przez C++.

W pierwszej kolejności napiszmy sobie funkcję w C++ (plik CALC.cpp):

#include "CALC.h"
#include 
#include 
#include 

void _CALC(double *DATA, int *RESULT, int N,  double treshold)
{
	printf("C++ start\n");

	RESULT[0] = 0;
	for (int a = 0; a < N; a++)
	{
		for (int b = a+1; b < N; b++)
		{
			for (int c = b+1; c < N; c++)
			{
				if (fabs(DATA[a] + DATA[b] + DATA[c] - 1) < treshold)
				{
					RESULT[0]++;
				}
			}
		}
	}

	printf("C++ finish\n");
}

W parze mamy plik nagłówkowy (CALC.h):

void _CALC(double *DATA, int *RESULT, int N, double treshold);

Powyższa funkcja robi dokładnie to co założyliśmy - zlicza ile kombinacji bez powtórzeń daje sumę 1 z dokładnością "treshold". Za chwilę powinno być jasne, czym są argumenty tej funkcji.

Naszym celem jest możliwość wywołania tej funkcji z poziomu pythona w taki sposób (program.py):

import numpy
import CALC
import time
start_time = time.time()

N = 1000
treshold = 1e-7

DATA = numpy.random.uniform(0, 1, size=N)
RESULT = numpy.zeros((1,), dtype=numpy.int)

start_time = time.time()
CALC.CALC(DATA, RESULT, N, treshold)
t1 = (time.time() - start_time)
print("C++:", RESULT[0], t1)


start_time = time.time()
count = 0;
for a in range(0,N):
	for b in range(a+1,N):
		for c in range(b+1,N):
			if abs(DATA[a]+DATA[b]+DATA[c]-1) < treshold:
				count =count + 1
t2 = (time.time() - start_time)
print("Python:", count, t2)

print("Speed:", t2/t1)

Dokładnie chodzi o linię

CALC.CALC(DATA, RESULT, N, treshold)

Potrzebujemy powiązania (interfejsu) pomiędzy pythonem i C++ (CALCmodule.cpp). Jest to kawałek kody, który z naszej funkcji C++ pozwala zrobić moduł dla pythona. Nie pisałem tego sam, złożyłem z kilku tutoriali. Najważniejsza jest funkcja CALC, gdzie definiujemy nasze zmienne i wywołujemy funkcję liczącą:

#include 
#define NPY_NO_DEPRECATED_API NPY_1_10_API_VERSION
#include 
#include "CALC.h"



struct module_state {
    PyObject *error;
};

#define GETSTATE(m) ((struct module_state*)PyModule_GetState(m))

static PyObject* CALC(PyObject* self, PyObject* args)
{
	PyArrayObject *data_obj;
	PyArrayObject *result_obj;
	double treshold;
	int N;

	if (!PyArg_ParseTuple(args, "OOid", &data_obj, &result_obj, &N, &treshold))
	{
		Py_INCREF(Py_None);
		return Py_None;
	}

	double *DATA = static_cast(PyArray_DATA(data_obj));
	int *RESULT = static_cast(PyArray_DATA(result_obj));

	_CALC(DATA, RESULT, N, treshold);

	Py_INCREF(Py_None);
	return Py_None;
}

static PyMethodDef CALCMethods[] = {
    {"CALC", CALC, METH_VARARGS, "..."},
    {NULL, NULL, 0, NULL}
};

static int CALC_traverse(PyObject *m, visitproc visit, void *arg) {
    Py_VISIT(GETSTATE(m)->error);
    return 0;
}

static int CALC_clear(PyObject *m) {
    Py_CLEAR(GETSTATE(m)->error);
    return 0;
}


static struct PyModuleDef moduledef = {
        PyModuleDef_HEAD_INIT,
        "CALC",
        NULL,
        sizeof(struct module_state),
        CALCMethods,
        NULL,
        CALC_traverse,
        CALC_clear,
        NULL
};

extern "C" PyObject * PyInit_CALC(void)
{
	PyObject *module = PyModule_Create(&moduledef);
	if (module == NULL)
        return NULL;
    struct module_state *st = GETSTATE(module);

    st->error = PyErr_NewException("CALC.Error", NULL, NULL);
    if (st->error == NULL) 
    {
        Py_DECREF(module);
        return NULL;
    }
	import_array();
	Py_INCREF(module);
    return module;

}

Na stronie https://docs.python.org/2/c-api/arg.html znajdują się dokładnie opisy typów danych które możemy przekazywać z pythona do C++ i z powrotem. W powyższym przykładzie "OOid" oznacza dwa obiekty, liczbę całkowitą i zmiennoprzecinkową (double). Obiektem może być lista, z której możemy czytać i do której pisać. Nazwa naszego modułu pojawia się w tym kodzie wielokrotnie - zmieniając trzeba pamiętać o wszystkich miejscach.

Ostatnim krokiem jest kompilacja modułu. Potrzebujemy skryptu do kompilacji (setup.py):

from distutils.core import setup, Extension
import numpy.distutils.misc_util
import os

os.environ["CC"] = "g++"
os.environ["CXX"] = "g++"

module1 = Extension('CALC', sources = ['CALCmodule.cpp', 'CALC.cpp'])

setup (name = 'CALC',
       version = '1.0',
       description = 'Package for calculating number sums',
       ext_modules = [module1],
	   include_dirs=numpy.distutils.misc_util.get_numpy_include_dirs())

Teraz możemy kompilować i uruchamiać:

Funkcja w C++ potrzebowała 0.12 sekundy na policzenie, że jest 19 trójek, które z dokładnością 0.0000001 dają sumę 1. Python potrzebował na analogiczne liczenie 108.42 sekund: 876 razy dłużej! Oczywiście można by ten kod zoptymalizować, ale porównujemy dwa identyczne rozwiązania. Różnica w czasie jest kolosalna. W mojej pracy często spotykam się z obliczeniami numerycznymi, dlatego w pierwszej kolejności przygotowuję prototyp w pythonie i fragment obliczeniowy przepisuję na C++.

Nie chcę wnikać w teorię budowania modułów dla pythona bo się na tym nie znam. Sam najlepiej uczę się na przykładach, więc mam nadzieję, że ten Wam się przyda.

Dajcie znać, czy kiedyś używaliście takich rozwiązań, albo czy przydadzą Wam się w przyszłości.

Jeszcze przed obroną doktoratu zostałem poproszony o wygłoszenie wykładu o tematyce mojej pracy doktorskiej w ramach studenckiego koła geofizyki na wydziale fizyki uniwersytetu warszawskiego.

Było to dla mnie ważne z dwóch powodów. Po pierwsze był to pierwszy wykład w ramach koła w tym roku akademickim, a po drugie było to moje pierwsze wystąpienie jako doktora. OK, nie czepiajmy się szczegółów, doktorem zostanę po uchwale rady wydziału, ale to tylko kwestia formalna.

Wszystkich zainteresowanych zapraszam do oglądania. Niestety kamera sportowa zawiodła mnie po raz kolejny – tym razem jasność nie pozwala na odczytanie niczego ze slajdów, więc znów dodałem je w montażu. Druga kamera również mnie zawiodła przerywając nagrywanie po 30 minutach. Mam nadzieję, że kiedyś nauczę się jak porządnie nagrywać.

W poniedziałek 17 października broniłem na Wydziale Fizyki Uniwersytetu Warszawskiego swoją pracę doktorską zatytułowaną „Lokalizacja ognisk trzęsień ziemi metodą propagacji wstecznej z wykorzystaniem implementacji metody fast marching w modelu 3D obszaru Polski”.

Praca doktorska była podsumowaniem 4-letnich studiów doktoranckich, podczas których realizowałem cały temat. Oczywiście najważniejszym celem studiów doktoranckich jest zdobycie stopnia naukowego doktora, ale równie ważne jest wykorzystanie tego czasu na wyrabianie sobie “naukowej marki”. Dowodem wyników pracy w nauce są publikacje w czasopismach naukowych (recenzowanych czasopismach z tak zwanej listy filadelfijskiej). W 4 lata zostałem współautorem 8 takich publikacji.

Drugą połowę czwartego roku studiów doktoranckich poświęciłem na zebranie wszystkich zagadnień, nad którymi wcześniej pracowałem w spójną całość – 180 stron pracy doktorskiej. Moją pracę doktorską możecie pobrać z mojej prywatnej strony: http://marcinpolkowski.com/en/ (po prawej stronie w boksie “Education”).

Należy też zaznaczyć, że jednym z ważnych elementów studiów doktoranckich jest dydaktyka – ja uczyłem głównie programowania w pythonie i C++.

Samo przygotowanie dysertacji (wiki) to dopiero początek drogi do uzyskania stopnia doktora. Trzeba otworzyć przewód doktorski, powołać komisje na dwa (lub trzy jak ktoś nie ma certyfikatu z angielskiego) egzaminy, wyznaczyć dwóch recenzentów i powołać komisję, która przyjmie obronę pracy. Każdy z tych elementów przechodzi przez Radę Wydziału. Gotowa praca trafia najpierw do recenzentów, którzy ją opiniują i dopuszczają do dalszego postępowania – mają na to dwa miesiące (recenzje mojej pracy: recenzja prof. Pietsch, recenzja prof. Guterch). Po otrzymaniu recenzji komisja ds. obrony sprawdza wszystkie dokumenty i wyznacza termin obrony. Moja obrona miała miejsce 17 października 2016.

W celach głównie pamiątkowych postanowiłem dyskretnie zarejestrować całe wydarzenie. Małą kamerę sportową schowałem w piórniku koło komputera, żeby nikogo nie kuła w oczy. Całoś trwała równo dwie godziny. Myślę, że możecie być ciekawi jak taka obrona wygląda (moja prezentacja startuje około 8-9 minuty):

Gdybyście mieli pytania o studia doktorancie, proces uzyskiwania stopnia czy samą moją pracę to chętnie będę odpowiadał – komentujcie poniżej!

Nie zapomnijcie subskrybować mojego kanału na YouTube:

Przez ostatnie 4 lata prowadziłem na Wydziale Fizyki ćwiczenie do kilku kursów programowania w pythonie i c++. Kilka razy spotkałem się z pytaniem o liczby zmiennoprzecinkowe. Poruszamy na zajęciach kwestię zapisu binarnego liczb całkowitych i wtedy pojawia się pytanie o liczby rzeczywiste. Do tej pory odpowiadałem bardzo skrótowo, że “jest to skomplakowne”, bardzo sprytne, ale nie pozbawione wad. W przypadku liczb całkowitych sprawa jest prosta. W danej ilości bitów możemy przechować dany zakres liczb całkowitych. Z liczbami zmiennoprzecinkowymi problem polega na tym, że w każdym przedziale jest ich nieskończenie wiele, a my mamy tylko 32 albo 64 bity. Więc jak to właściwie jest? O samej formie zapisu można oczywiście przeczytać wszystko w Wikipedii – ja chcę zrobić eksperyment. Na pierwszy ogień liczby typu float, czyli 32 bitowe (1 bit na znak, 8 bitów na wykładnik i 23 na mantysę). 32 bity daje jakieś 4 294 967 296 kombinacji – to na tyle mało, że mogę sprawdzić wszystkie. W ramach testu sprawdzam ile liczb mieści się w 4 przedziałach: 0-1, 1-2, 1000-1001 i 1000000-1000001:

int main()
{
	float f;
	uint i;	 
	std::vector range = {0, 1, 1000, 1000000};
	std::vector count = {0, 0, 0, 0};
	for(i = 0; i=range[a] && f < range[a]+1)
				count[a]++;
		}
	}
	for(int a = 0; a < range.size(); a++)
	{
		std::cout << range[a] << "-" << range[a]+1 << ":\t" << count[a] << std::endl;
	}
}

Wynik mamy po niecałych 2 minutach. Policzone ilości liczb są bardzo ciekawe:

Z liczbami w pojedynczej precyzji sprawa była prosta - udało się sprawdzić wszystkie kombinacje w niecałe 2 minuty. Przy podwójnej precyzji mamy 64 bity, więc sprawdzenie wszystkich zajęłoby 16000 lat. Szkoda czasu. Możemy natomiast pobawić się takim kodem:

int main()
{
	union
	{
		uint64_t input;
		uint64_t   output;
	} data;

	uint64_t x = 5;
	data.input = x;
	std::bitset   bits(data.output);

	for(uint64_t a = 0; a<=2048; a++)
	{
		for(uint64_t b = 0; b<=pow(2,52)-1; b+=pow(2,52)-1)
		{
			x = b+(a<<52);
			data.input = x;
			bits = data.output;
			double f;
			memcpy(&f, &x, sizeof(f));
			std::cout << std::setw(4)<< a<< '\t'<< std::setw(16)<<  b << '\t'<< std::setw(12) << f << '\t' << bits  << std::endl;
		}
		std::cout << std::endl;
	}
}

Powyższy kod wypisze zakres liczb dla każdego z 2^11 wykładników. Mantysa ma 52 bity (~4.5*10^15), więc od 1 do 2 mamy właśnie 2^52 liczb, Tyle samo od 2 do 4, 4 do 8 itd:

Na powyższym przykładzie pierwsza kolumna to wykładnik, druga to mantysa, trzecia to wartość liczby, a czwarta to reprezentacja binarna.

Nie chcę tutaj wnikać w szczegóły samego formatu. Zależy mi na pokazaniu specyfiki takiego pamiętania liczb. Możecie modyfikować i uruchamiać powyższe kody, a dowiecie się więcej o liczbach zmiennoprzecinkowych. Dodatkowo powyższe przykłady powinny dobrze obrazować kolosalną różnicę pomiędzy pojedynczą i podwójną precyzją.

Musimy pamiętać, że operując na bardzo dużych liczbach, gdzie dokładność zapisu zmiennoprzecinkowego jest mała, możemy wprowadzać błędy numeryczne - w szczególności gdy łączymy wielkie i małe liczby.

Kto dowiedział się czegoś nowego klika łapkę w górę!

Dostałem dziś mailem pytanie dotyczące wcześniejszych wpisów dotyczących budowy prostego czujnika sejsmicznego (Sejsmograf domowej roboty – pierwsza przymiarka oraz Sejsmograf domowej roboty – krok drugi):

Każdy taki mail jest dla mnie bardzo ważny, bo oznacza, że treść, którą tworzę jest wartościowa. Dziś postanowiłem odpowiedzieć publicznie na blogu – myślę, że odpowiedź może się komuś jeszcze przydać.

Cewka zastosowana w sejsmometrze to zwykły zwój drutu mający dwa końce. W wyniku zmiany pola magnetycznego spowodowanej ruchem magnesu w cewce indukuje się napięcie, czyli różnica potencjału pomiędzy obydwoma końcami drutu. Wartość tego napięcia jest proporcjonalna do tego, jak szybko zmienia się pole magnetyczne – im szybciej rusza się magnes tym większe napięcie. Do pomiaru stałego napięcia można oczywiście użyć multimetru, natomiast w naszym przypadku zależy nam na pomiarze zmian napięcia w czasie. Pomiaru napięcia w funkcji czasu możemy dokonać za pomocą oscyloskopu. W przypadku, gdy nie mamy pod ręką oscyloskopu można ratować się kartą dźwiękową komputera z wejściem mikrofonowym. Sygnał z mikrofonu jest napięciem zmiennym w czasie dokładnie tak jak w przypadku czujnika sejsmicznego. Niektóre karty dźwiękowe w komputerach stacjonarnych mają również wejścia zwane line-in, które działają bardzo podobnie.

Oczywiście karta dźwiękowa nie będzie doskonałym oscyloskopem (i tym bardziej rejestratorem sejsmicznym) ponieważ:

1. Ma ograniczony zakres napięć wejściowych – to bardzo ważne, bo przekroczenie maksymalnego dopuszczalnego napięcia może kartę uszkodzić. Dokładne wartości należy sprawdzić w dokumentacji danego modelu karty.
2. Pracuje w domyślnym dla dźwięku próbowaniu (44.1, 48kHz) jest to zdecydowanie za dużo, ale lepiej za dużo niż za mało.
3. Ma spore szumy, co jest istotne przy rejestracji słabych sygnałów.
4. Ma zazwyczaj niską rodzielczość.

Do samego podłączenia będziemy potrzebowali wtyk jack 3.5mm (mono lub stereo). Do masy podpinamy jeden koniec drutu cewki, a do plusa (w przypadku wtyku stereo do kanału L lub P) drugi. Zamiany miejscami kabli od cewki będzie miała wpływ na kierunek pomiaru: dodatnie napięcie przy ruch w górę lub w dół – na początek zabawy nie ma to żadnego znaczenia, przy głębszej interpretacji warto to wiedzieć i kontrolować. Czasami łatwiej będzie kupić kabel zakończony złączem jack (np. przedłużacz do słuchawek, lub kabel do podłączania odtwarzacza do radia), przeciąć i połączyć kabelki. W przypadku kabla stereo trzeba uważać, żeby nie połączyć cewki do kanałów L i P, bo wtedy nie zadziała – jeden koniec musi być do masy. Ja stosowałem właśnie taki przerobiony kabel:

Na komputrze można użyć dowolnego oprogramowania do rejestracji dźwięku lub specjalnego programu udającego oscyloskop: Zelscope lub Soundcard Oscilloscope. Drugi z nich widać na moich starych zdjęciach (to zdjęcie jest 2006, dlatego w warsztacie miałem monitor CRT):

Pozdrawiam i jak zawsze czekam na Wasze maile!

Zapraszam w środę o 20 na nowy wpis – w tym tygodniu będzie o liczbach zmiennoprzecinkowych.

Za niecałe dwa tygodnie będę bronił doktorat i skończy się moja przygoda ze studiowaniem. Zaproponowano mi pozostanie na uczelni w charakterze pracownika naukowego, ale mimo tego – trochę z ciekawości – postanowiłem wysłać swoje CV do kilku firm z branży IT w Polsce i w USA. Praca w USA to dla mnie pewien cel na przyszłość. Na razie rozmowy o ew. przeprowadzce zostały odsunięte w czasie, żeby choć trochę odchować potomstwo (a druga córka w drodze).

Google

Pierwsza przygoda z poszukiwaniem pracy związana była z polskim oddziałem Google. Wysłałem CV i po ok. tygodniu odezwał się rekruter z Londynu. Pierwsza rozmowa z rekruterem trwała ok. pół godziny. Po kilku dniach zaproponowano mi termin rozmowy technicznej. Pozwolono mi wybrać język programowania (C++, java lub Python, którego wybrałem). Rozmowa odbyła się przez telefon (wyświetlił mi się numer ze Szwajcarii) i trwała 40 minut. Swoją odpowiedź tworzyłem w specjalnie do tego celu przygotowanym dokumencie google, który w momencie zakończenia rozmowy stał się tylko do odczytu i chwilę później zniknął. Pierwsze zadanie to prosta łamigłówka (nie trzeba było nic kodować):

Po powyższym dostałem zadanie z właściwego programowania i muszę przyznać, że zajęło mi chwile zrozumienie o co chodziło autorowi. Teraz jak patrze z perspektywy czasu treść jest oczywista a rozwiązanie banalne. Zadanie brzmiało:

Po około tygodniu otrzymałem telefon, że proces rekrutacji nie będzie kontynuowany.

Tesla Motors

Druga próba była związana ze stanowiskiem inżyniera oprogramowania map i nawigacji w Tesla Motors. Tutaj cały kontakt był przez e-mail a czas pomiędzy wiadomościami liczony był w tygodniach. Po 3 mailach zostałem poproszony o rozwiązanie zadania programistycznego (Tesla Coding Challenge). Zadanie zostało mi udostępnione do pobrania. Miałem tydzień na zabranianie się za rozwiązywanie, ale od pobrania do odesłania zadania czas miał nie przekroczyć 4h (rekomendowane 3h). Rozwiązanie miało być w C++. Zadanie brzmiało:

Tak dział mój program:

Po trzech tygodniach dostałem maila, że doszli do wniosku, że potrzebują kogoś od zaraz, a mnie trzeba by “sprowadzić”.

HSBC

Wysłałem CV również na stanowisko specjalisty ds. oceny ryzyka w banku – było to jedno z niewielu stanowisk, gdzie oczekiwano znajomości matlaba. Praca niestety w Krakowie, ale zostałem zaproszony na rozmowę przez telekonferencję do siedziby w Warszawie. Zostałem poproszony o próbki mojego kodu w matlabie, a że nie mam nic, co nie byłoby tajemnicą mojego dotychczasowego pracodawcy poproszono mnie o rozwiązanie zadania:

Nie wiem czy dostali moje rozwiązanie, bo wrzuciłem w załącznik duży plik ze słownikiem sjp.pl, ale już się nie odezwali.

Nie publikuję na razie rozwiązań – pobawcie się sami, kiedyś pokażę swoje. Tak na prawdę najbardziej szkoda mi było pracy w google. Ani do USA, ani (tym bardziej) do Krakowa w obecnej sytuacji nie zdecydowałbym się jechać całą rodziną. Google był na miejscu. Niestety była to moja pierwsza rekrutacja – nie miałem pojęcia czego się spodziewać i poległem na banalnym zadaniu. Zachęcali mnie, żeby spróbować w przyszłości i pewnie się jeszcze skuszę. Na razie mam nową pracę.