0
$\begingroup$

I apologize for the very long question with a relatively complicated code. I am dealing with a real problem, and the question below appeared from my attempt to simplify this problem to ask about it here, while simultaneously not simplifying too much such that its scope will be lost. In addition, the reason for the problem is relatively simple: how Break[] compares (in particular, inside a compiled code).

What I need to do

So, consider a 2D table tab2 that has 4*n columns, where $n\geq 2$ is an arbitrary integer. The columns may be grouped by 4-columns,

1,2,3,4,...,4*(n-1)+1,4*(n-1)+2,4*(n-1)+3,4*n = (1,2,3,4),...(4*(n-1)+1,4*(n-1)+2,4*(n-1)+3,4*n)

The data from each of the 4-columns are used to calculate some boolean condition condition. The arguments of condition depend also on some external parameters contained in a table tab1. For the given row of tab2, let us define the acceptance: it is =1 if the product condition(4-column1)*condition(4-column2) = 1 for at least one pair of the 4-columns, or zero otherwise.

(The tables tab1, tab2, and the condition condition are defined below.)

I need to calculate the acceptance averaged over the rows of tab2 for the given row of tab1. So the output should be the list of acceptances with the dimensionality of tab1.

My attempts

My attempt to solve this task resulted in two codes: acceptanceNotOptimized and acceptanceOptimized.

(See these codes below)

The first one, acceptanceNotOptimized, loops over all pair combinations and the rows, then sums the boolean conditions.

It has two problems:

  1. Consider, some 4-column, say, (1,2,3,4). For n 4-columns, it enters m<n-1 4-column pair combinations: e.g.,

(1,2,3,4)+(5,6,7,8), (1,2,3,4)+(9,10,11,12),...

In the code, the functions projection1,projection2,condition are evaluated m times for the given 4-column. In the truly optimized code, this has to be done only once.

  1. The code loops over all pair combinations and does not care whether the acceptance is 1 for one of them. Because of this, even if for the given pair combination for the given row j the pair acceptance is 1, the code does not stop. It evaluates the acceptances for all possible pair combinations.

(Because of this, in the end, the code may return acceptances > 1, which I will then replace with 1 using acceptanceMaxRepl defined below)

acceptanceOptimized should partially deal with these problems. It differs from acceptanceNotOptimized literally by adding two Break[] inside looping over the first 4-column and the second 4-column. The goal of Break[] is to stop once at least for one pair of 4-columns the acceptance = 1. For large n, it works considerably faster than acceptanceNotOptimized.

However, it is still not the optimal one. For instance, in the marginal case, when there is no pair for which the acceptance = 1, it will compute all the boolean conditions for each 4-column m times, where m is the number of its occurrences in the pair combinations.

I have an idea for further optimization, but I got stuck on reproducing the results of the first code with the second code, which is crucial.

The question

The problem is that acceptanceOptimized and acceptanceNotOptimized give different acceptances (see the illustrations below). It may be related to the second Break[], but I do not see any problem with it. What may be the reason? Also, is it possible to further optimize the acceptanceOptimized?

Pre-definitions required to launch the code

(*Dimensionality of tab1*)
nvals = 10^6;
(*tab1*)
tab1Temp := 
  Join[RandomReal[{0, 1}, {nvals, 1}], RandomReal[{0, 1}, {nvals, 1}],
    RandomReal[{1, 5}, {nvals, 1}], 2];
(*Dimensionality of tab2*)
nvals2 = 10^3;
(*tab2. Each 4-th column is a "charge". The 4-columns may form a pair \
only if they have the opposite charges*)
tab2temp[n_] := Block[{}, charge := Exp[-I*Pi*RandomInteger[{0, 1}]];
  chargelist = Table[0, n];
  While[Length[DeleteDuplicates[chargelist]] == 1, 
   chargelist = Table[charge, n]];
  Join[##, 2] & @@ 
   Table[Join[RandomReal[{0., 0.3}, {nvals2, 1}], 
     RandomReal[{0., 0.2}, {nvals2, 1}], 
     RandomReal[{0.8, 0.99}, {nvals2, 1}], 
     Table[{chargelist[[i]]}, nvals2], 2], {i, 1, n, 1}]]
(*This code returns tab1,tab2*)
datapreptemp[ncols_] := Block[{}, tt = tab2temp[ncols];
  {tab1Temp, tt, chargelist, ncols}]
(*Functions of elements of tab1 (x0,y0,z0) and the first 3 columns of \
a 4-column of tab2 (vx,vy,vz)*)
projection1[x0_, z0_, vx_, vz_] = x0 + (10 - z0)*vx/vz;
projection2[y0_, z0_, vy_, vz_] = y0 + (10 - z0)*vy/vz;
(*The boolean condition where the functions above should be inserted*)
condition[projx_, projy_] = 
  Boole[-2 < projx < 2 && -1.5 < projy < 1.5];

The code - not optimized

acceptanceNotOptimized = 
  Hold@Compile[{{tab1, _Real, 1}, {tab2, _Real, 
           2}, {ncols, _Integer}}, 
         Module[{count, x0vals, y0vals, z0vals, indexprod2, vxval1, 
           vyval1, vzval1, vxval2, vyval2, vzval2, xprod1, yprod1, 
           xprod2, yprod2, acc, productindex, charge1, charge2, acc1, 
           acc2},
          count = 0.;
          (*Extracting elements of tab1*)
          x0vals = Compile`GetElement[tab1, 1];
          y0vals = Compile`GetElement[tab1, 2];
          z0vals = Compile`GetElement[tab1, 3];
          (*Looping over rows of tab2*)
          Do[
           acc = 0.;
           (*Looping over 4-columns of tab2*)
           Do[
            (*Extracting elements of tab2*)
            vxval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 1];
            vyval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 2];
            vzval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 3];
            (*Calculating the boolean arguments*)
            xprod1 = projection1[x0vals, z0vals, vxval1, vzval1];
            yprod1 = projection2[y0vals, z0vals, vyval1, vzval1];
            (*Condition for the first 4-column*)
            acc = condition[xprod1, yprod1];
            acc2 = 0.;
            (*If the condition is !=0, looping over 4-
            columns that may create a pair with the first 4-column*)
            If[acc != 0.,
             charge1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 4];
             Do[
              charge2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 4];
              (*Only if charge1 = -charge2, the given 4-columns form a pair*)
              If[charge1 == -charge2,
               vxval2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 1];
               vyval2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 2];
               vzval2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 3];
               xprod2 = projection1[x0vals, z0vals, vxval2, vzval2];
               yprod2 = projection2[y0vals, z0vals, vyval2, vzval2];
               (*Condition for the second 4-column - this is the 4-
               column pair acceptance*)
               acc2 = condition[xprod2, yprod2];
               count += acc2];
              , {f2, f1 + 1, ncols, 1}]];
            , {f1, 1, ncols - 1, 1}];
           , {j, 1, nvals2, 1}];
          {count/nvals2}], CompilationTarget -> "C", 
         RuntimeOptions -> "Speed", 
         RuntimeAttributes -> {Listable}] /. OwnValues@nvals2 /. 
      DownValues@projection1 /. DownValues@projection2 /. 
    DownValues@condition // ReleaseHold;
acceptanceMaxRepl = 
  Compile[{{tab, _Real, 2}}, 
   Table[{Min[tab[[i]][[1]], 1]}, {i, 1, Length[tab], 1}], 
   CompilationTarget -> "C", RuntimeOptions -> "Speed"];

The code - optimized

acceptanceOptimized = 
  Hold@Compile[{{tab1, _Real, 1}, {tab2, _Real, 
           2}, {ncols, _Integer}}, 
         Module[{count, x0vals, y0vals, z0vals, indexprod2, vxval1, 
           vyval1, vzval1, vxval2, vyval2, vzval2, xprod1, yprod1, 
           xprod2, yprod2, acc, productindex, charge1, charge2, acc1, 
           acc2},
          count = 0.;
          (*Extracting elements of tab1*)
          x0vals = Compile`GetElement[tab1, 1];
          y0vals = Compile`GetElement[tab1, 2];
          z0vals = Compile`GetElement[tab1, 3];
          (*Looping over rows of tab2*)
          Do[
           acc = 0.;
           (*Looping over 4-columns of tab2*)
           Do[
            (*Extracting elements of tab2*)
            vxval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 1];
            vyval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 2];
            vzval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 3];
            (*Calculating the boolean arguments*)
            xprod1 = projection1[x0vals, z0vals, vxval1, vzval1];
            yprod1 = projection2[y0vals, z0vals, vyval1, vzval1];
            (*Condition for the first 4-column*)
            acc = condition[xprod1, yprod1];
            acc2 = 0.;
            (*If the condition is !=0, looping over 4-
            columns that may create a pair with the first 4-column*)
            If[acc != 0.,
             charge1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 4];
             Do[
              charge2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 4];
              (*Only if charge1 = -charge2, the given 4-columns form a pair*)
              If[charge1 == -charge2,
               vxval2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 1];
               vyval2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 2];
               vzval2 = Compile`GetElement[tab2, j, 4*(f2 - 1) + 3];
               xprod2 = projection1[x0vals, z0vals, vxval2, vzval2];
               yprod2 = projection2[y0vals, z0vals, vyval2, vzval2];
               (*Condition for the second 4-column - this is the 4-
               column pair acceptance*)
               acc2 = condition[xprod2, yprod2];
               count += acc2];
              (*If for the given second 4-
              column the pair acceptance is 1, then break*)
              If[acc2 > 0., Break[]], {f2, f1 + 1, ncols, 1}]];
            (*If for the given first 4-
            column the pair acceptance is 1, then break*)
            If[acc2 > 0., Break[]], {f1, 1, ncols - 1, 1}];
           , {j, 1, nvals2, 1}];
          {count/nvals2}], CompilationTarget -> "C", 
         RuntimeOptions -> "Speed", 
         RuntimeAttributes -> {Listable}] /. OwnValues@nvals2 /. 
      DownValues@projection1 /. DownValues@projection2 /. 
    DownValues@condition // ReleaseHold;

How to launch the code

First, one should generate tab1 and tab2 for the given n; let us call the corresponding list dataprep. This is done using (here n = 4)

nval = 4;
dataprep = datapreptemp[nval];

Once dataprep is generated, the codes are launched using this command (where Optimized is "True" if one wants to use acceptanceOptimized, and something else otherwise):

acctab[dataprep_, Optimized_] := Module[{tab1, tab2, ncols},
  tab1 = dataprep[[1]];
  tab2 = dataprep[[2]];
  ncols = dataprep[[4]];
  tabtemp = 
   If[Optimized == "True", acceptanceOptimized[tab1, tab2, ncols], 
    acceptanceNotOptimized[tab1, tab2, ncols]]
  ]

The differences: illustration

Consider n = 4:

nval = 4;
dataprep = datapreptemp[nval];
data1 = acctab[dataprep, "True"]; // AbsoluteTiming
data2 = acctab[dataprep, "False"]; // AbsoluteTiming

data1, data2 are different:

Style[Row[{Histogram[data1 // Flatten, 100, ScalingFunctions -> {"Linear", "Log"}], Histogram[acceptanceMaxRepl[data2] // Flatten, 100, ScalingFunctions -> {"Linear", "Log"}]}], ImageSizeMultipliers -> {1, 1}]

enter image description here

The jump in unit value in data2 is caused by the fact that there are many rows for which more than one pair combinations have the unit boolean condition. In data1, it is also observed, although much less explicitly. The shape of the distribution for the lower values is similar.

It looks like that the problem sits in the second Break[]: if removing it from acceptanceOptimized, it shows the same distribution. I do not understand what the issue is with this Break[]. It should just break the cycle over the first 4-column once the cycle over the second 4-column was broken by the first Break[]. Otherwise, it should do nothing as acc2 = 0 in the beginning...

It is interesting that the predictions of the two codes agree if nvals2 = 1, but the differences appear for any other nvals. I tried to put the break condition in different places of the code, but it did not work...

Maybe it has something to do with the option RuntimeAttributes->{Listable}...

Improve this question
$\endgroup$
5
  • $\begingroup$ Do you know which code is wrong? $\endgroup$
    – Michael E2
    Jul 8 at 19:35
  • $\begingroup$ @MichaelE2 : I do not know, but naively I expect that acceptanceNotOptimized is correct. Break[] may in principle work incorrectly inside the code. $\endgroup$
    – John Taylor
    Jul 8 at 19:50
  • 1
    $\begingroup$ One thing that has helped me debug things is Trace. $\endgroup$
    – Peter Burbery
    Jul 10 at 21:55
  • 1
    $\begingroup$ I strongly suggest replacing Break with alternative conditional code and see what the results are. $\endgroup$
    – Somos
    Jul 11 at 21:03
  • $\begingroup$ @MichaelE2 : it turns out that the non-optimized code was wrong. So I am currently looking for a way to further optimize acceptanceOptimized. $\endgroup$
    – John Taylor
    Jul 12 at 13:12

1 Answer 1

Reset to default
0
$\begingroup$

It seems that I understand why the results of the codes are different and which one is correct.

Let us assume that we have, e.g., ncols=4 4-columns and 4 different 4-columns pairs. Assuming that tab2 has 1000 rows, acceptanceNotOptimized evaluates the acceptance (per row of tab1) as $$ \text{min}\left(1,\sum_{i = 1}^{1000}\frac{1}{1000}(\text{acc}_{\text{row}_{i}})\right), $$ where $\text{acc}_{\text{row}_{i}}\in (0,4)$. In reality, it must be $$ \tag 2 \sum_{i = 1}^{1000}\frac{1}{1000}\text{min}[\text{acc}_{\text{row}_{i}},1] $$ In acceptanceOptimized, the evaluation for each row stops once at least for one pair acc2 is 1, so $\text{acc}_{\text{row}_{i}}$ is always $0$ or $1$. This way, acceptanceNotOptimized always produced a higher acceptance which is wrong.

Regarding the optimization, I came up with the possibility of creating a table acccolumns having 1000 rows and ncols columns. In the beginning, it is filled with 1. However, once for the given row j, I sum over f2, I then replace the corresponding element acccolumns[[j]][[f2]] = 0 if acc2 = 0. Then, once the cycle over f1 reaches this f2, it is just skipped. The corresponding code acceptanceOptimized2 (see below the code) should work faster than acceptanceOptimized. In particular, if most of the columns return zero acc2, it should be much more efficient.

In reality, however, it just always evaluates slightly slower than it (see below). I am trying to understand why. Most likely, it is related to the way I implement this "optimization"

I would appreciate any help with the speedup.

The code acceptanceOptimized2

acceptanceOptimized2 = 
  Hold@Compile[{{tab1, _Real, 1}, {tab2, _Real, 
           2}, {ncols, _Integer}, {acccolumns, _Real, 2}}, 
         Module[{count, x0vals, y0vals, z0vals, indexprod2, vxval1, 
           vyval1, vzval1, vxval2, vyval2, vzval2, xprod1, yprod1, 
           xprod2, yprod2, acc, productindex, charge1, charge2, acc2, 
           acc3, acccols},
          count = 0.;
          (*Extracting elements of tab1*)
          x0vals = Compile`GetElement[tab1, 1];
          y0vals = Compile`GetElement[tab1, 2];
          z0vals = Compile`GetElement[tab1, 3];
          acccols = acccolumns;
          (*Looping over rows of tab2*)
          Do[
           acc = 0.;
           acc2 = 0.;
           (*Looping over 4-columns of tab2 (index f1)*)
           Do[
            (*If the acceptance for the second 4-
            column generated for the previous f1 is zero, 
            then skip this 4-column*)
            If[Compile`GetElement[acccols, j, f1] == 1.,
             (*Extracting elements of tab2*)
             vxval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 1];
             vyval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 2];
             vzval1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 3];
             (*Calculating the boolean arguments*)
             xprod1 = projection1[x0vals, z0vals, vxval1, vzval1];
             yprod1 = projection2[y0vals, z0vals, vyval1, vzval1];
             (*Condition for the first 4-column*)
             acc = condition[xprod1, yprod1];
             (*If the condition is !=0, looping over 4-
             columns that may create a pair with the first 4-
             column (index f2)*)
             If[acc != 0.,
              charge1 = Compile`GetElement[tab2, j, 4*(f1 - 1) + 4];
              Do[
               If[Compile`GetElement[acccols, j, f2] == 1.,
                
                charge2 = 
                 Compile`GetElement[tab2, j, 4*(f2 - 1) + 4];
                (*The pair is formed only if charge1 = -charge2*)
                If[charge1 == -charge2,
                 
                 vxval2 = 
                  Compile`GetElement[tab2, j, 4*(f2 - 1) + 1];
                 
                 vyval2 = 
                  Compile`GetElement[tab2, j, 4*(f2 - 1) + 2];
                 
                 vzval2 = 
                  Compile`GetElement[tab2, j, 4*(f2 - 1) + 3];
                 
                 xprod2 = 
                  projection1[x0vals, z0vals, vxval2, vzval2];
                 
                 yprod2 = 
                  projection2[y0vals, z0vals, vyval2, vzval2];
                 (*Condition for the second 4-column - this is the 4-
                 column pair acceptance*)
                 acc2 = condition[xprod2, yprod2];
                 (*If the Boolean condition for the second 4-
                 column is zero, 
                 then for the cycle over f1 will skip f1=f2. 
                 Otherwise, break the cycle*)
                 
                 If[acc2 == 0.,(*Compile`GetElement[acccols,f2]*)
                  acccols[[j, f2]] = 0., count += acc2; Break[]];
                 ]], {f2, f1 + 1, ncols, 1}];
              ];
             (*Once for the given first 4-
             column the pair acceptance is 1, then break the cycle*)
             If[acc2 == 1, Break[]]], {f1, 1, ncols - 1, 1}];
           , {j, 1, nvals2, 1}];
          {count/nvals2}], CompilationTarget -> "C", 
         RuntimeOptions -> "Speed", 
         RuntimeAttributes -> {Listable}] /. OwnValues@nvals2 /. 
      DownValues@projection1 /. DownValues@projection2 /. 
    DownValues@condition // ReleaseHold;
acctab2[dataprep_] := 
 Module[{tab1, tab2, ncols, acccolumns1, acccolumns},
  tab1 = dataprep[[1]];
  tab2 = dataprep[[2]];
  ncols = dataprep[[4]];
  acccolumns1 = Table[1., ncols];
  acccolumns = Table[acccolumns1, nvals2];
  acceptanceOptimized2[tab1, tab2, ncols, acccolumns]
  ]
<< CompiledFunctionTools`
CompilePrint@acceptanceOptimized2;

Comparison with acceptanceOptimized

Before launching, please change nvals2 to 100.

dataprep = datapreptemp[25];
chargelist
data1 = acctab[dataprep, "False"]; // AbsoluteTiming
data2 = acctab[dataprep, "True"]; // AbsoluteTiming
data3 = acctab2[dataprep]; // AbsoluteTiming
{{"\!\(\*SubscriptBox[\(data\), \(1\)]\) = \
\!\(\*SubscriptBox[\(data\), \(2\)]\)", 
   "\!\(\*SubscriptBox[\(data\), \(2\)]\) = \!\(\*SubscriptBox[\(data\
\), \(3\)]\)"}, {If[acceptanceMaxRepl[data1] == data2, "True", 
    "False"], If[data2 == data3, "True", "False"]}} // TableForm
Histogram[{data1 // Flatten}, 100, 
 ScalingFunctions -> {"Linear", "Log"}]
Style[Row[{Histogram[{acceptanceMaxRepl[data1] // Flatten, 
     data2 // Flatten, data3 // Flatten}, 100, 
    ScalingFunctions -> {"Linear", "Log"}]}], 
 ImageSizeMultipliers -> {1, 1}]

So the timings for data2 and data3 are 1.44 and 1.62 s.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.