Re: Do loop beschleunigen

[Udo und Susanne Krause <su.krause@XXXXXXX.ch>]

Hallo Ulrich

> Das Problem ist, daß die Do-Schleife sehr langsam ist. Ein
> Ausweg ist natürlich ein cpp-Program mit Mathlink, das wird
> mittelfristig sowieso geschrieben werden. Aber für viele ähnliche
> Probleme würde mich der prinzipielle Weg interessieren, solche
> primitiven Schleifen zu beschleunigen.

<snip>

> 2.) Ist es schneller, eine Liste aufzubauen, in der man zu jedem
> Matrixelement den aus dem Radius berechneten Index für die radiale Liste
> berechnet und anschließend über diese Matrix als Index nutzt und damit
> Do-Schleifen ganz vermeiden? Aber mit welchem Befehl? Die Elemente der
> zweidimensionale Liste müßen ja auf die einer eindimensionale Liste
> abgebildet werden.

Das ist kein Dimensionsproblem, denn wenn Sie die den Abstand des Bildpunktes
haben,
koennen Sie die Matrix der Bildpunkte zu einer flachen Liste machen. Das
Problem ist,
dass man jeden Bildpunkt durchlaufen muss und dann entsprechend zählen und
umspeichern muss.

> 3.) Oder ist hier ein ganz anderer Weg einzuschlagen?
>
> Hier meine sehr cpp-artige Schleife:
> c1=500; c2=492;
> rad = Table[{i, 0, 0}, {i, 1, 1024}];     (* Tabelle fuer radiale
> Intensitaetsabhaengigkeit *)
> Do[
>          r = Sqrt[(i - c1)^2 + (j - c2)^2]; (* berechnet Dadius zum Zentrum
> (c1,c2) *)
>          buc = Round[r] + 1;                       (* daraus Index "buc" wie
> "bucket" berechnen *)
>          rad[[buc, 2]] += 1;                          (* Zaehlt, wieviele
> Pixel in diesen Topf fallen *)
>     rad[[buc, 3]] += bm[[i, j]],                (* addiert Pixel-Intensitaet
> zu diesem Topf *)
> {i, 1, Length[bm]}, {j, 1, Length[bm[[1]]]}]; (* Schleife ueber 1024x1024
> Pixel *)
> Die Schleife braucht auf einem schnelle Rechner ca. eine Minute.

Nimmt man das zum Ausgangspunkt, so hat man auf einem 1998-NT-Rechner (64 MB
Arbeitsspeicher):
In[1]:=
Remove[metrik]
metrik[i1_Integer, i2_Integer, c1_Integer, c2_Integer] :=
    Round[Sqrt[(i1 - c1)^2 + (i2 - c2)^2]] + 1;
In[5]:=
Remove[schwarzLoopBrute];
schwarzLoopBrute[seed_Integer, origX_Integer, origY_Integer, diagL_Integer] :=
Module[{c1 = origX(*= 500*), c2 = origY (*= 492*), bm, buc, rad, maxI},
  SeedRandom[seed];
  bm = Table[Random[Integer, {0, 100}], {i1, diagL}, {i2, diagL}];
  maxI = Max[metrik[1, 1, origX, origY], metrik[diagL, 1, origX, origY],
      metrik[1, diagL, origX, origY], metrik[diagL, diagL, origX, origY]];
  rad = Table[{i, 0, 0}, {i, 1, maxI}];
  Do[
    (* If[j == 1 && Mod[i, diagL/4] == 0, Print[" i: ", i]]; *)
    buc = Round[Sqrt[(i - c1)^2 + (j - c2)^2]] + 1;
    rad[[buc,2]] += 1;
    rad[[buc,3]] += bm[[i,j]],
    {i, 1, Length[bm]}, {j, 1, Length[bm[[1]]]}
  ];
  rad
] /; diagL > 0
In[7]:=
res0 = Timing[schwarzLoopBrute[20020712, 127, 129, 256]];
In[8]:=
res0[[1]]
Out[8]=
56.752 Second

Also fast 1 Minute für eine Aufgabe mit 1/16 der Datenmenge mit einem alten
Rechner.

Die üblichen Mittel, um schneller zu werden, sind funktionale Programmierung
und reine Funktionen (inlining). Wenn es sehr viele Daten sind, muss auch
Fold[] und Sort[] vermieden werden.

Eine Möglichkeit ist:
In[3]:=
Remove[schwarzNoLoop];
schwarzNoLoop[seed_Integer, origX_Integer, origY_Integer, rowL_Integer,
    colL_Integer] :=
  Module[{c1 = origX, c2 = origY, hL, o, maxI},
      SeedRandom[seed];
      maxI =
        Max[metrik[1, 1, origX, origY], metrik[rowL, 1, origX, origY],
          metrik[1, colL, origX, origY], metrik[rowL, colL, origX, origY]];
      Do[hL[i0] = List[i0], {i0, maxI}];
      o =
        MapThread[(hL[#1] = {hL[#1], #2}) &,
          Transpose[
            Flatten[Table[{Round[Sqrt[(i1 - c1)^2 + (i2 - c2)^2]] + 1,
                  Random[Integer, {0, 100}]}, {i1, colL}, {i2, rowL}], 1]]];
      Clear[o];
      ({First[#], Length[#] - 1, Plus @@  Drop[#, 1]}) &  /@
        Table[Flatten[hL[i0]], {i0, maxI}]
      ] /; rowL > 0 && colL > 0
In[10]:=
res1 = Timing[schwarzNoLoop[20020712, 127, 129, 256, 256]];
In[11]:=
res1[[1]]
Out[11]=
40.168 Second
In[12]:=
SameQ[res0[[2]], res1[[2]]]
Out[12]=
True

Das spart immerhin fast 1/3 der Zeit. Es wird eine Anzahl Hashlisten hL
definiert, die mit einer reinen Funktion gefüllt werden. Die
Auswerte-Do-Schleife ist das MapThread geworden. Die Auswertung ist auch
funktional.

Mit 1024 x 1024 Punkten konnte ich keine Zeitmessung machen, weil dann der
Rechner schon swappt hier.

Mit den besten Grüssen
Udo.