0
$\begingroup$

I'm converting a code from MatLab to Mathematica. It takes MatLab less than a second to compute it while its more than ten minutes in Mathematica. I was hoping to improve the performance.

Here's a skeleton of the code

finout= ParallelTable[Quiet[Max[
             Table[myFunc[var1, var2, var3, var4, var5,var6], 
                  {var6, 1, 30}, {var5, 1, 60}, {var4, 1,2}, 
                  {var3, 1, 2}, {var2, 1, 2}]]], 
        {var1, 1, 200}]

myFunc is function taking six variables . finout gives the maximum value of myFunc for all values of var1.

I'm only working with real numbers. I've already wrapped myFunc with Compile. Is there anything more that could be done?

EDIT

Heres the complete code that I use. There are three parts; first I define some functions and variables, then I define keyrateFunc (corresponding to myFunc in the skeleton code above), in the final part I evaluate keyrateFunc over the input range.

Can I improve the performance of keyrateFunc?

(*Defined functions and parameter values*)

functiongamma[a_, b_, c_, d_] := Module[{aux1, aux2, gamma},
 
 aux1 = ((c + d)*(1 - b)*b)/(c*d*Log[2]);
 aux2 = Log[2, ((c + d)/(c*d*(1 - b)*b))*(21^2/a^2)];
 gamma = Sqrt[aux1*aux2]]

functionh[x_] := Module[{hx},
 
 If[x === Indeterminate, 0, If[
   Element[x, Reals] && 0 < x < 1,
   hx = -x*Log[2, x] - (1 - x) Log[2, 1 - x], 0
   ]]]


myMin[x_] := 
Module[{}, 
 If[x === Indeterminate, 0.5, 
  If[Element[x, Reals], Min[0.5, x], 0.5]]]


myMax[x_] := 
Module[{}, 
 If[x === Indeterminate, 0, 
  If[Element[x, Reals], Max[0, x], 0]]]

SKRMin = 10.;
RepRate = 500000000.;
mu3 = 2.000000000000000 10^-04;
Y0 = 6.007379741251949 10^-07;
eps = 1.000000000000000 10^-25;
betasec = 57.564627324851145;
ed = 0.01;
fEC = 1.16;
Nb = 10^11.;
RateMin = SKRMin/RepRate;
NdBpoints = 200; 
dBmin = 0; 
dBmax = 65.;
step = 0.01;
index = 0;
bitsDecoy = 2.;
bitsBasis = 2.;
stepBasis = 0.5^bitsBasis;
stepDecoy = 0.5^bitsDecoy;






(*The main function *)

keyrateFunc = Compile[{dB, qx, pmu1, pmu2, mu1, mu2},

 pmu3 = 1. - pmu1 - pmu2;
 etasys = 10.^(-dB/10.);
 
 nxmu1 = Floor[Nb*qx*qx*pmu1*(1. - (1. - Y0)*Exp[-etasys*mu1])]; 
 nxmu2 = Floor[Nb*qx*qx*pmu2*(1. - (1. - Y0)*Exp[-etasys*mu2])]; 
 nxmu3 = Floor[Nb*qx*qx*pmu3*(1. - (1. - Y0)*Exp[-etasys*mu3])];
 
 Emu1 = ((Y0/2. - ed)*Exp[-etasys*mu1] + 
     ed)/(1. - (1. - Y0)*Exp[-etasys*mu1]); 
 Emu2 = ((Y0/2. - ed)*Exp[-etasys*mu2] + 
     ed)/(1. - (1. - Y0)*Exp[-etasys*mu2]); 
 Emu3 = ((Y0/2. - ed)*Exp[-etasys*mu3] + 
     ed)/(1. - (1. - Y0)*Exp[-etasys*mu3]); 
 mxmu1 = Ceiling[nxmu1*Emu1]; mxmu2 = Ceiling[nxmu2*Emu2]; 
 mxmu3 = Ceiling[nxmu3*Emu3];
 
 nzmu1 = 
  Floor[Nb*(1. - qx)*(1. - qx)*
    pmu1*(1. - (1. - Y0)*Exp[-etasys*mu1])]; 
 nzmu2 = Floor[
   Nb*(1. - qx)*(1. - qx)*pmu2*(1. - (1. - Y0)*Exp[-etasys*mu2])]; 
 nzmu3 = Floor[
   Nb*(1. - qx)*(1. - qx)*pmu3*(1. - (1. - Y0)*Exp[-etasys*mu3])];
 
 Emu1 = ((Y0/2. - ed)*Exp[-etasys*mu1] + 
     ed)/(1. - (1. - Y0)*Exp[-etasys*mu1]);
 Emu2 = ((Y0/2. - ed)*Exp[-etasys*mu2] + 
     ed)/(1. - (1. - Y0)*Exp[-etasys*mu2]);
 Emu3 = ((Y0/2. - ed)*Exp[-etasys*mu3] + 
     ed)/(1. - (1. - Y0)*Exp[-etasys*mu3]);
 mzmu1 = Ceiling[nzmu1*Emu1];
 mzmu2 = Ceiling[nzmu2*Emu2];
 mzmu3 = Ceiling[nzmu3*Emu3];
 
 nx = nxmu1 + nxmu2 + nxmu3;
 mx = mxmu1 + mxmu2 + mxmu3;
 nz = nzmu1 + nzmu2 + nzmu3;
 mz = mzmu1 + mzmu2 + mzmu3;
 eobs = mx/nx;
 
 tao0 = pmu1*Exp[-mu1] + pmu2*Exp[-mu2] + pmu3*Exp[-mu3];
 nxmu3m = (Exp[mu3]/
     pmu3)*(nxmu3/(1. + (3.*betasec + 
          Sqrt[8.*betasec*nxmu3 + betasec^2.])/(2.*(nxmu3 - 
            betasec))));
 nxmu2p = (Exp[mu2]/
     pmu2)*(nxmu2/(1. - (-betasec + 
          Sqrt[8.*betasec*nxmu2 + 9.*betasec^2.])/(2.*(nxmu2 + 
            betasec))));
 
 auxSx0 = tao0*(mu2*nxmu3m - mu3*nxmu2p)/(mu2 - mu3);
 
 Sx0 = myMax[auxSx0];
 
 Sx0true = Nb*qx^2*tao0*Y0;
 tao1 = 
  pmu1*Exp[-mu1]*mu1 + pmu2*Exp[-mu2]*mu2 + pmu3*Exp[-mu3]*mu3;
 nxmu2m = (Exp[mu2]/
     pmu2)*(nxmu2/(1. + (3.*betasec + 
          Sqrt[8.*betasec*nxmu2 + betasec^2.])/(2.*(nxmu2 - 
            betasec))));
 nxmu3p = (Exp[mu3]/
     pmu3)*(nxmu3/(1. - (-betasec + 
          Sqrt[8.*betasec*nxmu3 + 9.*betasec^2.])/(2.*(nxmu3 + 
            betasec))));
 nxmu1p = (Exp[mu1]/
     pmu1)*(nxmu1/(1. - (-betasec + 
          Sqrt[8.*betasec*nxmu1 + 9.*betasec^2.])/(2.*(nxmu1 + 
            betasec))));
 
 auxSx1 = (tao1*
     mu1*(nxmu2m - 
       nxmu3p - ((mu2^2 - mu3^2)/(mu1^2))*(nxmu1p - (Sx0/
            tao0))))/(mu1*(mu2 - mu3) - mu2^2 + mu3^2);
 
 Sx1 = myMax[auxSx1];
 
 Sx1true = Nb*qx^2*tao1*(1. - (1. - Y0)*(1. - etasys));
 nzmu3m = (Exp[mu3]/
     pmu3)*(nzmu3/(1. + (3.*betasec + 
          Sqrt[8.*betasec*nzmu3 + betasec^2])/(2*(nzmu3 - betasec))));
 nzmu2p = (Exp[mu2]/
     pmu2)*(nzmu2/(1. - (-betasec + 
          Sqrt[8.*betasec*nzmu2 + 9.*betasec^2])/(2*(nzmu2 + 
            betasec))));
 
 auxSz0 = tao0*(mu2*nzmu3m - mu3*nzmu2p)/(mu2 - mu3);
 
 Sz0 = myMax[auxSz0];
 Sz0true = Nb*(1. - qx)^2*tao0*Y0;
 
 
 
 nzmu2m = (Exp[mu2]/
     pmu2)*(nzmu2/(1. + (3.*betasec + 
          Sqrt[8.*betasec*nzmu2 + betasec^2.])/(2.*(nzmu2 - 
            betasec))));
 nzmu3p = (Exp[mu3]/
     pmu3)*(nzmu3/(1. - (-betasec + 
          Sqrt[8.*betasec*nzmu3 + 9.*betasec^2.])/(2.*(nzmu3 + 
            betasec))));
 nzmu1p = (Exp[mu1]/
     pmu1)*(nzmu1/(1. - (-betasec + 
          Sqrt[8.*betasec*nzmu1 + 9.*betasec^2.])/(2.*(nzmu1 + 
            betasec))));
 
 auxSz1 = (tao1*
     mu1*(nzmu2m - 
       nzmu3p - ((mu2^2. - mu3^2)/(mu1^2.))*(nzmu1p - (Sz0/
            tao0))))/(mu1*(mu2 - mu3) - mu2^2. + mu3^2.);
 
 Sz1 = myMax[auxSz1];
 Sz1true = Nb*(1. - qx)^2.*tao1*(1. - (1. - Y0)*(1. - etasys));
 
 
 mzmu2p = (Exp[mu2]/
     pmu2)*(mzmu2/(1. - (-betasec + 
          Sqrt[8.*betasec*mzmu2 + 9.*betasec^2.])/(2.*(mzmu2 + 
            betasec))));
 mzmu3m = (Exp[mu3]/
     pmu3)*(mzmu3/(1. + (3.*betasec + 
          Sqrt[8.*betasec*mzmu3 + betasec^2.])/(2.*(mzmu3 - 
            betasec))));
 
 auxVz1 = tao1*(mzmu2p - mzmu3m)/(mu2 - mu3);
 Vz1 = myMax[auxVz1];
 
 nzmu1singleerror = 
  Ceiling[Nb*(1. - qx)*(1. - qx)*Exp[-mu1]*
    mu1*((Y0/2. - ed)*(1. - etasys) + ed)];
 nzmu2singleerror = 
  Ceiling[Nb*(1. - qx)*(1. - qx)*Exp[-mu2]*
    mu2*((Y0/2. - ed)*(1. - etasys) + ed)];
 nzmu3singleerror = 
  Ceiling[Nb*(1. - qx)*(1. - qx)*Exp[-mu3]*
    mu3*((Y0/2. - ed)*(1. - etasys) + ed)];
 
 nzsingleerror = 
  nzmu1singleerror + nzmu2singleerror + nzmu3singleerror;
 
 
 auxphix = Vz1/Sz1 + functiongamma[eps, Vz1/Sz1, Sz1, Sx1];
 phix = myMin[auxphix];
 
  lambdaEC = fEC*functionh[eobs]; 
 keylength = 
  Floor[Max[0, 
    Sx0 + Sx1 - Sx1*functionh[phix] - nx*lambdaEC - 
     Log[1./(2.*eps^4)]]]; keyrate = keylength/Nb, 
 RuntimeOptions -> "Speed"]





(*The higher the value, the faster the code runs;Ideally should 
be equal to 1*)
stepmul = 
 5;
keytab = 
 ParallelTable[
  Quiet[ Max[
    Table[keyrateFunc[dB, qx, pmu1, pmu2, mu1, mu2], {mu1, 0.4, 1., 
      stepmul step}, {mu2, mu3 + 0.001, Min[0.3, mu1 - mu3 - 0.001], 
      stepmul step}, {pmu1, 0.25, 1. - 2.*stepDecoy, 
      stepDecoy}, {pmu2, stepDecoy, Min[1. - pmu1 - stepDecoy, 0.5], 
      stepDecoy}, {qx, 0.5, 1. - stepBasis, stepBasis}]]], {dB, 
   dBmin, 65, (dBmax - dBmin)/(NdBpoints - 1.)}];
ListLogPlot[keytab]






Link takes you to the.nb file in my drive

$\endgroup$
9
  • 3
    $\begingroup$ It's difficult to provide the best answer without myFunc. However, it will certainly be significantly faster if you use Table instead of many Dos that constantly rewrite alist: ParallelTable[Quiet[Max[Table[myFunc @@ {var1, var2, var3, var4, var5, var6}, {var6, 1, 30}, {var5, 1, 60}, {var4, 1, 2}, {var3, 1, 2}, {var2, 1, 2}]]], {var1, 1, 200}] $\endgroup$
    – Domen
    Aug 17, 2023 at 16:32
  • $\begingroup$ Its faster...thanks. However still too slow. myFunc is not my code so I dont know how much I can share. Ill try to add a model soon. $\endgroup$
    – Dotman
    Aug 17, 2023 at 16:50
  • 1
    $\begingroup$ @Dotman Please show your matlab code too! $\endgroup$ Aug 17, 2023 at 17:18
  • 1
    $\begingroup$ NMaximize supports constrained integer programming to find the max within your variable bounds. It might be faster than brute-force testing every single point in the region to find the max for each value of var1. $\endgroup$
    – ydd
    Aug 17, 2023 at 17:27
  • 1
    $\begingroup$ @Dotman, I have looked into the code. Your compilation doesn't help much because the code is not written "correctly" (that is: compiled code has to call the main kernel evaluator many times, which slows down the code). Therefore, the main problem is not in your ParrallelTable code, but in your keyrateFunc. If you'd like help on that, please edit the question accordingly :) Also, you can just include the whole code directly into this question (it is not very long anyways). $\endgroup$
    – Domen
    Aug 18, 2023 at 8:19

1 Answer 1

2
$\begingroup$

The main reason your compiled code is slow is that it is not actually being compiled very much. You can check this with:

Needs["CompiledFunctionTools`"]
CompilePrint[keyrateFunc]

Wherever you see MainEvaluate, the code is being slowed down because it calls Mathematica kernel to get some result. To actually produce compiled code, you have to make the code such that it doesn't depend or interfere with global symbols outside of it. All numeric constants (SKRMin, RepRate, mu3 ...) have to be inserted into the code at the compilation time, and all intermediate results (pmu3, etasys, nxmu1) have to be localized to the function.

Therefore, do it like this (I've included all constants inside, but you should move the ones that you need globally outside of Compile):

keyrateFuncBetter = Compile[{dB, qx, pmu1, pmu2, mu1, mu2},
  Module[{SKRMin = 10.,RepRate = 500000000.,mu3 = 2.000000000000000 10^-04,
    Y0 = 6.007379741251949 10^-07,eps = 1.000000000000000 10^-25,
    betasec = 57.564627324851145,ed = 0.01,fEC = 1.16,Nb = 10^11.,
    RateMin = SKRMin/RepRate, NdBpoints = 200,dBmin = 0,
    dBmax = 65.,step = 0.01,
    index = 0,bitsDecoy = 2.,bitsBasis = 2.,
    stepBasis = 0.5^bitsBasis,stepDecoy = 0.5^bitsDecoy,
    pmu3, etasys, nxmu1, nxmu2, nxmu3, Emu1, Emu2, Emu3, mxmu1,
    mxmu2, mxmu3, nzmu1, nzmu2, nzmu3, mzmu1, mzmu2, mzmu3, nx, mx, 
    nz, mz, eobs, tao0, nxmu3m, nxmu2p, auxSx0, Sx0, Sx0true, tao1, 
    nxmu2m, nxmu3p, nxmu1p, auxSx1, Sx1, Sx1true, nzmu3m, nzmu2p, 
    auxSz0, Sz0, Sz0true, nzmu2m, nzmu3p, nzmu1p, auxSz1, Sz1, 
    Sz1true, mzmu2p, mzmu3m, auxVz1, Vz1, nzmu1singleerror, 
    nzmu2singleerror, nzmu3singleerror, nzsingleerror, auxphix, phix,
    lambdaEC, keylength, keyrate},

    pmu3 = 1. - pmu1 - pmu2;
    etasys = 10.^(-dB/10.);
    (* ... remaining of the code ... *)
    keyrate = keylength/Nb

 ], RuntimeOptions -> "Speed", CompilationTarget -> "C", 
 RuntimeAttributes -> Listable]

Note that I have also set it to compile to C-code, and made the function listable (see comment by @LukasLang).

If you do CompilePrint[keyrateFuncBetter], you will now see that pretty much all code is being properly compiled. However, there are still some external evaluates which call your other functions (functiongamma, functionh, myMin ...). Long story short, you have to separately compile them, too, and called compiled versions of them. Note that functionh, myMin and myMax will probably not compile well because you use very high-level Mathematica functions. You should probably rewrite them.

$\endgroup$
1
  • $\begingroup$ Thanks..This does increase the speed considerably. Using CompilationOptions -> {"InlineExternalDefinitions" -> True} also helped, althoughit introduces some errors in the final result $\endgroup$
    – Dotman
    Aug 18, 2023 at 12:37

Your Answer

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

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