Compile a code evaluating a condition only if another condition is positive

Question

Consider the following three tables:

thzrange = 
  Flatten[Table[{-Pi + i*0.1, 40 + j*1}, {i, 0, 20, 1}, {j, 0, 30, 
     1}], {1, 2}];
phrange = RandomReal[{-Pi, Pi}, 10^4];
Exrange = Table[Ex, {Ex, 5, 100, 5}];

Using their elements, one may obtain values of some functions:

coord1[th_, z_, ph_] = z*Cos[th]*Sin[ph];
coord2[Ex_, th_, z_, ph_] = Ex*Exp[-z/50]*Cos[th]^2*Sin[ph]^5;

I would like to obtain the following table:

{Ex,th,z,acc1,acc2}

where acc1,acc2 are some acceptances depending on Ex,th,z (boolean conditions on correspondingly coord1,coord2) that are averaged over values of phrange. Moreover, acc2 may be non-zero (i.e., starts evaluating) only if acc1 is non-zero. The test conditions are

condition1[x1_] = Boole[-5 < x1 < 5];
condition2[x2_] = Boole[-15 < x2 < 0];

where x1, x2 are the values of coord1 and coord2.In reality, they may be much more complicated.

So this is my code: AccComp which computes the acceptances acc1, acc2 for kth combination of th,z and lth value of Ex:

AccComp = 
  Hold@Compile[{{thzrange, _Real, 2}, {phrange, _Real, 
           1}, {Exrange, _Real, 1}, {k, _Integer}, {l, _Integer}}, 
         Module[{x1val, x2val, thval, zval, Exval, phval, acc1v, acc1,
            acc2, acc2counter, xproj1val},
          acc1 = 0;
          acc2 = 0.;
          acc2counter = 0;
          Do[
           phval = Compile`GetElement[phrange, m];
           thval = Compile`GetElement[thzrange, k, 1];
           zval = Compile`GetElement[thzrange, k, 2];
           x1val = coord1[thval, zval, phval];
           acc1v = condition1[x1val];
           acc1 += acc1v;
           If[acc1v == 1,
            Exval = Exrange[[l]];
            acc2counter += 1;
            x2val = coord2[Exval, thval, zval, phval];
            acc2 += condition2[x2val];
            ], {m, 1, Length[phrange], 1}];
          {acc1/Length[phrange], 
           If[acc2counter != 0, acc2/acc2counter, 0]}], 
         CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
       DownValues@coord1 /. DownValues@coord2 /. 
     DownValues@condition1 /. DownValues@condition2 // ReleaseHold;

and the final table:

thzExtable = 
  Flatten[Table[{Exrange[[l]],thzrange[[k]][[1]], thzrange[[k]][[2]]},{l, 1, 
     Length[Exrange], 1}, {k, 1, Length[thzrange], 1}], {1, 2}];
tableaccvalsComp = 
 Hold@Compile[{{thzrange, _Real, 2}, {phrange, _Real, 
       1}, {Exrange, _Real, 1}}, 
     Table[AccComp[thzrange, phrange, Exrange, k, l],{l, 1, Length[Exrange], 1}, {k, 1, 
       Length[thzrange], 1}], 
     CompilationTarget -> "C", RuntimeOptions -> "Speed"] /. 
   DownValues@AccComp // ReleaseHold

It is not fast:

Join[thzExtable, 
   Flatten[tableaccvalsComp[thzrange, phrange, Exrange], {1, 2}], 
   2]; // AbsoluteTiming

{2.48617,Null}

Could you please tell me how to optimize it further?

Answer 1

You can remove the random sampling of $\phi$ completely by computing the limits of the acceptable $\phi$ values. For this, we can use that your expressions only depend on a factor $ain\phi$, which is monotonous in $\phi$ over $[-\pi/2,\pi/2]$ (values outside that interval don't matter, since $\sin\phi$ is simply mirrored at $\pm\pi/2$. The nested condition for x2 corresponds then simply to the intersection of the respective intervals.

Below the code for this "analytic approach". Note that I set phrange to include 10^5 values for better sampling of the original version of the code (this variable is not used in the new version):

thzrange = 
  Flatten[Table[{-Pi + i*0.1, 40 + j*1}, {i, 0, 20, 1}, {j, 0, 30, 
     1}], {1, 2}];
phrange = RandomReal[{-Pi, Pi}, 10^5];
Exrange = Table[Ex, {Ex, 5, 100, 20}];

thzExtable = 
  Flatten[Table[{Exrange[[l]], thzrange[[k]][[1]], 
     thzrange[[k]][[2]]}, {l, 1, Length[Exrange], 1}, {k, 1, 
     Length[thzrange], 1}], {1, 2}];

tableaccvals[thzrange_, Exrange_] :=
 Module[{thval, zval, x1fac, int1, x2fac, int2},
  Transpose@Table[
    thval = thzrange[[k, 1]];
    zval = thzrange[[k, 2]];
    x1fac = zval Cos[thval];
    int1 = If[x1fac == 0,
      {-\[Pi]/2., \[Pi]/2.},
      Sort@N@ArcSin@Clip[{-5., 5.}/zval/Cos[thval]]
      ];
    Table[
     x2fac = E^(-zval/50) Exrange[[l]] Cos[thval]^2;
     int2 = If[x2fac == 0,
       {-\[Pi]/2., \[Pi]/2.},
       Sort@N@ArcSin@Clip[Surd[{0., -15.}/x2fac, 5]]
       ];
     {-Subtract @@ int1/Pi
      , If[-Subtract @@ int1 != 0, (
       Min[int1[[2]], int2[[2]]] - 
        Max[int1[[1]], int2[[1]]])/-Subtract @@ int1, 0]},
     {l, Length[Exrange]}
     ],
    {k, Length[thzrange]}
    ]
  ]

cresult = 
   Join[thzExtable, 
    Flatten[tableaccvalsComp[thzrange, phrangeL, Exrange], {1, 2}], 
    2]; // AbsoluteTiming
(* {8.95623, Null} *)

aresult = 
   Join[thzExtable, Flatten[tableaccvals[thzrange, Exrange], {1, 2}], 
    2]; // AbsoluteTiming
(* {0.0897391, Null} *)

Max@Abs[aresult - cresult]
(* 0.00756831 *)

Note the 100x speedup, while simulatenously improving the accuracy of the result (the error will get smaller if you increase the number of ph samples, and thus runtime, of the original version)