C++Guns – RoboBlog blogging the bot

Habt ihr euch schon immer einmal gefragt, welche Bedingung in einer if() denn letztlich wahr wurde? Hier ist ein einfacher Weg das zu protokollieren.

inline bool test(const bool b, bool& p) {
  p = b;
  return b;
};

void func() {
  double x = 7;

// I thinks this is the best
    bool p1 = x>5;
    bool p2 = x<10;
    if(p1 and p2) {

// Or choose from this
//  bool p1 = false, p2 = false;
//  if(p=x>5 and p2=x<10) {
//  if(test(x>5, p1) and test(x<10, p2)) {
    if(p1) {
      cout << "p1 was true\n";
    }
    if(p2) {
      cout << "p2 was true\n";
    }
  }
}

Thanks to Juan Jimenez.

Diesen Post möchte ich mit einem kürzlich aufgetretenen Problem in Fortran beginnen.
In Fortran gibt es das forall keyword. Damit ist es möglich Schleifen zu bauen, wie die do Schleife, nur dass es möglich ist, forall Schleifen parallel abzuarbeiten.
Eine Vorraussetzung dafür ist folgende: Wenn data der Vector ist, auf den parallel gearbeitet werden soll, dann müssen die Operationen auf data(i) und data(j) mit i != j unabhängig voneinander sein.
Das ist, für das Beispiel, welches ich hier vorstellen möchte, der Fall.

Beispiel: Jeder Punkt im Array soll transformiert, und in einem neuen Array abgespeichert werden. Hier ist die erste native Implementation:

program LLL
  use naaaah
  implicit none
  integer :: i
  type(Point2D), allocatable :: data(:), data2(:)
  real :: x
  
  allocate(data(100))
  allocate(data2(100))
  
  forall(i=1:size(data))
    x = data(i)%x * cos(data(i)%y)
    data2(i)%x = x
    data2(i)%y = data(i)%y
  end forall
end program

Und die Compiler Fehlermeldung

x = data(i)%x * cos(data(i)%y)
Warning: The FORALL with index ‘i’ is not used on the left side of the assignment at (1) and so might cause multiple assignment to this object

Das ist vollkommen richtig, WENN die forall Schleife parallel ausgeführt wird, so wird parallel auf x geschrieben. Aber x ist nur ein scalar. Daher der Fehler. Abhilfe würde schaffen, x als Array zu declarieren. Aber hier gehen bei HPC sofort alle roten Alarmglocken an und die Notbremse wird gezogen.
Okay.
Untersuchen wir jetzt das Stück Code mal nach dem "Was aber nicht Wie" Prinzip. WAS soll der Code machen? Hm keine Ahnung. WIE macht er es? Er multipliziert x in Abhängigkeit von y. Okay das ist alles, was man wissen muss, um den Code umzuschreiben. Nennen wir die neue Funktion, welche obige Multiplikation ausführt als "toLLL". Damit sieht der Code so aus:

pure function toLLL(point) result(LLLPoint)
  implicit none
  type(Point2D), intent(in) :: point
  type(Point2D) LLLpoint
  
  LLLpoint%x = point%x * cos(point%y*DEGRAD)
  LLLpoint%y = point%y
end function

forall(i=1:size(data))
  data2(i) = toLLL(data(i))
end forall

Nun compiliert der Code ohne Fehler! Warum? Weil die vermeintliche Variable x nicht mehr temporär vorkommt. Lass es mich nochmal verdeutlichen. Im alten Code hatten wir so etwas:

xTemp = xalt*fac;
xneu = xTemp

Und im neuen Code haben wir so etwas

xneu= xalt*fac;

Es gibt also keine temporäre Zwischenvariablen die irgendwie von Hand verwaltet werden müsste. PLUS, der Code ist lesbarer geworden. PLUS der Code ist performanter, da er u.U. parallel ausgeführt werden kann. Ich will nicht behaupten, dass durch das Einführen einer neuen Funktion der Code mehr Richtung Funktionale Programmierung geht. Aber er ist definitiv besser. Das wäre es erst, wenn die Funktion toLLL als Funktionsparamter mit übergeben würde. Hier ist die Grenze von Fortran erreicht. Ohne Templates muss für jede generische Variante von Point2D eine neue Funktion toLLL geschrieben werden. Und das inlinen über mehrere Compilations Units hinweg ist auch nicht möglich.

In C++ gibt es natürlich die selbe Art von Problem, aber sie können leicht gelöst werden. Von Hand geschriebene Schleifen funktionieren natürlich immer. Aber, sobald Parallelisierung in's Spiel kommt. Egal auf welcher Art und Weise auch immer, ist es nicht mehr so einfach wie früher. Parallelisierung bedeutet auf viele Kleinigkeiten achten, auf die man vorher nicht achten musste. Aber wir leben ja nicht im Mittelalter. Wie haben Programme, welche die Checks für uns durchführen. Wir müssen sie nur benutzen.
Ein Beispiel dafür wäre std::transform Hier wird keine Schleife mehr explizit aufgeschrieben. Es wird nur noch gesagt WAS getan werden sollen. Das transformieren von data mit Funktion toLLL nach data2.
Hier das komplette C++ Programm

#include <vector>
#include <algorithm>
#include <cmath>
#include <iostream>

using namespace std;

constexpr const double DEGRAD = M_PI/180.0;

struct Point2D {
  double x, y;
};

Point2D toLLL(const Point2D& point) {
  return {point.x*cos(point.y*DEGRAD), point.y};
}
  

int main() {
  vector<Point2D> data(180), data2(180);
  for(size_t i=0; i < 180; ++i) {
    data[i].x = 180;
    data[i].y = -90.0+i;
  }
  
  transform(std::execution::par, begin(data), end(data), begin(data2), toLLL);
  
  for(const Point2D& p : data2) {
    cout << p.x << " " << p.y << "\n";
  }
}

Der ersten Parameter von transform bestimmt, ob der Datensatz parallel oder sequentiell abgearbeitet werden soll. Das Einführen der Funktion toLLL verhindert temporäre Variablen, und die Funktion wird als Parameter übergeben. Das generalisieren der Funktion toLLL auf einem Point Datentyp mit beliebiger Dimension und fundamentalen Type (int,float) ist nun leicht. Beispiel:

template<class Point>
Point toLLL(const Point& point) {
  Point LLLpoint = point;
  LLLpoint[0] *= cos(LLLpoint[1]*DEGRAD);
  return LLLpoint;
}

Die Idee ist von Stephen C. Dewhurst [1] die Abkürzung von mir

FVCRTTPPP - Fancy Variadic Curiously Recurring Template Template Parameter Pack Pattern

Praxisbeispiel:

using FVCRTTPPPScalar = Scalar<int, Ctr, Eq, Rel, Stream>; ! A countable, equal compareable, relational-able, streamable, scalar integer type.
FVCRTTPPPScalar x = 0;

[1] http://stevedewhurst.com/once_weakly/once-weakly20170328/once-weakly20170328.pdf

Der normale dynamic programing Zuker Algorithmus von Faltungen von Sequenzen hat eine Laufzeit von O(n^4), mit zuätzlichen Cache noch O(n^3). Mit der allgemeinen Four-Russians Speedup Methode für dynamic programing Algorithmen ist eine Laufzeit von O(n^3/log n) möglich. Für kleines n also absolut irrelevant. Für großes n halbiert, oder drittelt sogar die Laufzeit. Was dann aber immer noch unheimlich zu lange ist.
Aber was ist die Idee dahinter?
Als erstes wird festgestellt, dass es nicht möglich ist, die Four-Russians Technik auf die klassische Implementation des Zuker anzuwenden. Da die Auswertungsreihenfolge der Matrizen nicht passt. Also wird die Reihenfolge erst mal umgestellt. Erst kommt eine Schleife die nicht von der derzeitigen Spalte j in der Matrix abhängt. Und dann kommt eine 3fach Schleife die abhängt.
In der innersten Schleife findet nun die Optimierung statt. Die Schleife sieht so aus:

S(i,j)=max( S(i,j), S(i,k-1)+S(k,j) )

Für fixes i,j iteriert nur k von i+1 zu j-1. Die Idee ist nun, Indizes zu Gruppe der Größe q zusammen zu fassen. So dass die innerste Schleife nur noch n/q oft läuft.
Nun ja wenn sie meinen. Um das Maximum zu ermitteln, muss man so oder so jeden Wert in der Matrix in diesem Bereich anfassen. Wie die das verhindern wollen, das lese ich aus dem Paper nicht raus.
Zeitverschwendung.

[1] A Simple, Practical and Complete Time Algorithm for RNA Folding Using the Four-Russians Speedup - Yelena Frid and Dan Gusfield

Für die, die es immer noch nicht kapiert haben: Hier noch ein einfachereres Beispiel. Temporäre Objekte, alles klar? Speichert keine Referenzen, es gibt kein Grund das zu tun. Sonst schießt ich euch in den Fuß. In jede Zehe einzeln.

bool alive = true;

struct Foo {
  Foo() { }
  ~Foo() { alive = false; }
};

struct Bar {
  const Foo& dangling;
  
  Bar(const Foo& foo) : dangling(foo) { }
  
  void func() {
    if(alive) {
      cout << "dangling is okay\n";
    } else {
      cout << "dangling is dangling\n";
    }
  }
};

int main() {
  Bar bar{Foo{}};
  bar.func();
  return 0;
}