W ciągu ostatnich paru tygodni pracowałem nieco nad biblioteką, która pozwalałaby na prowadzenie obliczeń w zakresie geometrii różniczkowej w Ruście. Przez geometrię różniczkową mam tu na myśli głównie rachunek tensorowy, tylko w krzywych przestrzeniach bądź czasoprzestrzeniach. Coś tego typu stworzyłem już swego czasu w C++, postanowiłem jednak spróbować wykorzystać niektóre możliwości Rusta do stworzenia ulepszonej wersji w tym języku.
Co Rust mógłby zrobić lepiej?
Tensory do celów obliczeń najwygodniej zapisywać jako tablice liczb. Problem w tym, że sprowadzenie tensora do liczb wymaga wybrania układu współrzędnych w przestrzeni, w której ten tensor jest zdefiniowany. Różne operacje, np. dodawanie dwóch tensorów, mają sens jedynie wtedy, kiedy tensory biorące w niej udział są zdefiniowane w tym samym układzie współrzędnych. Jedyną możliwością zakodowania tej zasady w C++ było zakodowanie układu współrzędnych jako właściwości obiektu tensora i sprawdzanie w kodzie różnych operatorów, czy wszystko się zgadza. W ten sposób, jeśli zrobimy w kodzie jakiś błąd, program wychwyci go w trakcie wykonania.
Ok, czyli błędy były wykrywalne, co tu można zrobić lepiej? Ano można byłoby zrobić tak, żeby tensory zapisane w różnych układach współrzędnych miały nie tylko różną wartość jednej właściwości, ale były wręcz obiektami różnych typów. W ten sposób o błędzie może nas powiadomić jeszcze kompilator, zanim program w ogóle zostanie przekształcony w postać wykonywalną. W C++ nie byłoby to praktyczne, natomiast system typów Rusta pozwala to zrobić o wiele ciekawiej.
EDIT: Zwrócono mi uwagę, że szablony w C++ też pozwalają na coś takiego. Tym niemniej, zrobienie tego w Ruście było ciekawym eksperymentem :)
Problem z tablicami raz jeszcze
Skoro chcemy statycznie (tzn. jeszcze podczas kompilacji) sprawdzać zgodność tensorów pod względem układu współrzędnych, można by to było nieco rozszerzyć. Dodawanie np. ma sens tylko dla tensorów tego samego rzędu i o tej samej wariancji (tzn. "składzie" indeksów, które mogą być kowariantne lub kontrawariantne) - nie ma sensu choćby dodawanie wektora i kowektora (które mają ten sam rząd równy 1, ale jeden jest kontrawariantny, a drugi - kowariantny). Gdyby do tego również udało się zaprząc system typów, byłoby super, tylko czy to wykonalne?
Okazuje się, że tak, choć nie było to łatwe.
Tensorowe perypetie
Tensory charakteryzowane są w zasadzie przez dwie cechy: wymiar przestrzeni, na której są określone, i rząd / wariancję. Np. wektory i kowektory mają rząd 1, jak już wspomniałem wyżej, z kolei tensory macierzowe (np. transformacje liniowe albo dwuformy) są rzędu 2. Różnice tkwią w wariancji - przekształcenia liniowe mają jeden indeks kontrawariantny i jeden kowariantny, a dwuformy - dwa kowariantne. Wymiar przestrzeni razem z rzędem określają, ile tensor ma współrzędnych. Liczba ta wynosi 
, gdzie 
jest wymiarem przestrzeni, a 
- rzędem tensora. Ponieważ te wielkości będą znane statycznie, chciałoby się również statycznie określić wielkość tablicy współrzędnych. I tu zaczynają się problemy.
Tablice, jak już wspominałem w poprzednim wpisie, nie dopuszczają w Ruście parametryzacji względem długości. Na pierwszy rzut oka musiałbym wobec tego kodować każdy rodzaj tensora oddzielnie, z odpowiadającą mu wielkością tablicy. Na szczęście rozwiązanie tego problemu opisałem również w tamtym wpisie - jest nim struktura GenericArray. Typy, które ta struktura wykorzystuje do określania długości tablicy mają w bibliotece typenum zdefiniowaną całą arytmetykę, więc obliczanie wielkości takich jak nie powinno stanowić problemu. Jest to faktycznie możliwe, ale nie bezproblemowe.
Otóż dokonanie każdej operacji na typach liczbowych 
powoduje niestety utratę wszelkich gwarancji, jakie były założone dla parametrów. Np. jeśli określiłem, że typ reprezentujący wymiar przestrzeni i typ reprezentujący rząd tensora mogą być stosowane jako długość tablicy, ta gwarancja znika po utworzeniu typu będącego wynikiem działania . Muszę zatem określić ją osobno, co powoduje rozrost ograniczeń wypisywanych przy każdym typie i prowadzi do potworków typu czegoś takiego. Nieocenionej pomocy udzielił mi w tej kwestii paholg, autor typenum - dzięki niemu udało mi się wypisać warunki, które pozwoliły na kompilację struktury opisującej tensor.
Po przezwyciężeniu wstępnych przeszkód poszło już z górki. Kolejne operacje na typach tensorów było łatwo zakodować, kiedy już się wiedziało, czego oczekuje kompilator. I tak udało się zaimplementować dodawanie i odejmowanie tensorów, mnożenie ich przez skalary i przez siebie, a także zwężanie (co oznacza wysumowanie elementów tensora po pewnej "przekątnej", np. w przypadku macierzy zwężenie sprowadza się do policzenia śladu).
Kompilator się poddaje
Niestety, po zakodowaniu wszystkiego i napisaniu kilku testów okazało się, że kompilator nie wytrzymał napięcia. Operacja mnożenia tensorów przerosła go i wpuściła w nieskończoną pętlę.
Moje wysiłki zmierzające do identyfikacji i naprawienia problemu spełzły na niczym. Ewidentnie winnym jest mnożenie tensorów (wykomentowanie kodu odpowiedzialnego za tę operację powoduje, że błąd znika), ale jest nie do ustalenia, w jaki sposób powoduje ono problem. W akcie desperacji skompilowałem wersję debug kompilatora i wygenerowałem logi kompilacji. Otrzymałem ponad 500 MB pliku tekstowego, który niestety nic mi nie rozjaśnił. Społeczność reddita poradziła mi, abym spróbował zgłosić problem programistom kompilatora, toteż to zrobiłem. Póki co jednemu z programistów udało się zminimalizować kod wywołujący nieskończoną pętlę, jednak przyczyna wciąż jest nieznana.
Tym niemniej, przy zastosowaniu alternatywnej składni (<Tensor<T,U> as Mul<Tensor<T, V>>::mul(tensor1, tensor2) zamiast tensor1 * tensor2) kompilator sobie radzi i biblioteka jest funkcjonalna. Wczoraj wypuściłem wersję 0.1.
Na koniec: dodatkowe makro w GenericArray
Dodałem również jedno usprawnienie do GenericArray. Otóż do tej pory utworzenie takiej tablicy wiązało się z użyciem metody from_slice, która a) mocno zwiększa objętość kodu, b) sprawdza poprawność długości tablicy dynamicznie. Dzisiaj utworzyłem makro arr!, które nie dość, że skraca kod, to jeszcze pozwala na statyczne sprawdzenie poprawności kodu:
|
1 2 3 4 5 6 7 8 9 |
// stara wersja let array1 = GenericArray::<u32, U3>::from_slice(&[1, 2, 3]); // nowa wersja let array2 = arr![u32; 1, 2, 3]; // kompiluje się, ale wyrzuca błąd w trakcie wykonania let array3 = GenericArray<u32, U2>::from_slice(&[1, 2, 3]); // nie skompiluje się let array4: GenericArray<u32, U2> = arr![u32; 1, 2, 3]; |
Makro zostało wgrane na crates.io w wersji 0.1.1 biblioteki generic-array.
Podsumowanie i dalszy plan
Tak zatem wyglądała moja ostatnia aktywność programistyczna. W najbliższym czasie zamierzam dodać do biblioteki differential-geometry obsługę metryk, metryk odwrotnych i symboli Christoffela, co powinno być wystarczające do zakodowania klona wspomnianego na początku gr-engine. Gdy będę miał to za sobą, napiszę, co z tego wyszło :)