0
|
8534 Views
|
7 Replies
|
4 Total Likes
View groups...
Share

# Make code faster and in C-compatible manner?

Posted 8 years ago
 Hello, I need to make the following code very fast: B[a_, m_, t_] := Module[{n, p, q, seq, s, f, b}, n = 1 - m; p = m/a; q = m (1 - 1/a); seq[i_, j_] := Sequence @@ Flatten[Table[{i == K || j == t - K, s[i, j, K]}, {K, 1, Round[t/2]}], 1]; s[i_, j_, K_] := Sum[Binomial[i, i - k] Binomial[t - i, j - i + k] (p q)^ k ((n + p) (n + q))^(K - k), {k, 0, K}]; f[i_, j_] := If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)]; b[i_, j_] := Which[ j == t, p^(j - i) f[i, j] Binomial[t, i], i == t, q^(i - j) f[i, j] Binomial[t, j], i <= j, p^(j - i) f[i, j] Binomial[t, i] Which[Evaluate[seq[i, j]]], i > j, q^(i - j) f[i, j] Binomial[t, j] Which[Evaluate[seq[j, i]]]]; Return[Table[b[i, j], {i, 1, t}, {j, 1, t}]]];  It generates a t x t numerical matrix for an integer a and a real number m between 0 and 1. For low t it is fast enough but I need it for t=2000. I estimated that this would take years on a desktop computer... So I guess I need to compile the whole thing to C-code and possibly parallelize it. The problem is the seq function. It says that it cannot be compiled because the expression "Sequence" is not allowed. It generates a sequence of conditions of length t/2 and also uses the function s inside it. How can it be written in C-compatible manner? And is it possible to compile a function that uses another inside it? So I have a lot of questions. Unfortunately, my programming skills are only basic so any help would greatly be appreciated.
7 Replies
Sort By:
Posted 23 days ago
 I rewrote a bit. This still uses seq but not as a Sequence object. It has a final clause that is never actually needed. B2[a_, m_, t_] := Module[{n, p, q, seq, s, f, b}, n = 1 - m; p = m/a; q = m - p; s[i_, j_, K_] := Sum[Binomial[i, i - k] Binomial[t - i, j - i + k] (p q)^ k ((n + p) (n + q))^(K - k), {k, 0, K}]; f[i_, j_] := If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)]; seq[i_, j_] := Which[ i <= Round[t/2] && t - j <= Round[t/2], s[i, j, Min[i, t - j]] , i <= Round[t/2], s[i, j, i] , t - j <= Round[t/2], s[i, j, t - j] , True, 0]; b[i_, j_] := Which[ j == t, p^(j - i) f[i, j] Binomial[t, i] , i == t, q^(i - j) f[i, j] Binomial[t, j] , i <= j, p^(j - i) f[i, j] Binomial[t, i] seq[i, j] , i > j, q^(i - j) f[i, j] Binomial[t, j] seq[j, i]]; Table[b[i, j], {i, 1, t}, {j, 1, t}]] As you note, this can get into trouble with machine numbers. Timing[B2[3, .22, 400];] (* During evaluation of In[127]:= General::munfl: 3.75211*10^-299 5.02373*10^-10 is too small to represent as a normalized machine number; precision may be lost. During evaluation of In[127]:= General::munfl: 2.75154*10^-300 5.88718*10^-10 is too small to represent as a normalized machine number; precision may be lost. During evaluation of In[127]:= General::munfl: 2.0178*10^-301 6.89904*10^-10 is too small to represent as a normalized machine number; precision may be lost. During evaluation of In[127]:= General::stop: Further output of General::munfl will be suppressed during this calculation. Out[127]= {30.3564, Null} *) It handles higher precision without signs of distress though, and this is not hugely slower. Timing[b2big = B2[3, .2220, 400];] (* Out[129]= {51.2323, Null} *) 
Posted 23 days ago
 If you change Sum to: Total@Table** code is 20%** faster.
Posted 22 days ago
 Thank you Daniel and Mariusz for your suggestions. Meanwhile I found that the following code is the fastest version:  Bmatrix[a_, t_, m_] := Block[{n, p, q, B}, n = 1 - m; p = m/a; q = m (1 - 1/a); B = ConstantArray[0., {t, t}]; Do[B[[i, j]] = j! /i! p^(-i + j) (n + p)^(-j + t) (n + q)^ i Hypergeometric2F1Regularized[-i, j - t, 1 - i + j, ( p q)/((n + p) (n + q))], {i, 1, Floor[t/2]}, {j, i, t - i}]; Do[B[[i, j]] = j! /i! p^(-i + j) (n + p)^(-j + t) (n + q)^ i Hypergeometric2F1Regularized[-i, j - t, 1 - i + j, ( p q)/((n + p) (n + q))], {j, Ceiling[t/2], t - 1}, {i, t - j + 1, j}]; Do[B[[i, t]] = p^(t - i) (n + q)^i Binomial[t, i], {i, 1, t}]; Do[B[[i, j]] = (q/p)^(i - j) Binomial[t, i]/ Binomial[t, j] B[[j, i]], {i, 2, t}, {j, 1, i - 1}]; Return[B]]; It can be C-compiled but j! /i! and similar expressions produce too big numbers. To prevent this, I could use ((p^j) j! )/(p^i i!) = FactorialPower[p j, j - i, p] in the first Do loop because there j>=i and 0
Posted 22 days ago
 You might experiment with taking logs, PowerExpanding using LogGamma instead of Log[Factorial[...]], and exponentiating at the end. If there is intermediate coefficient growth followed by cancelation, the LogGamma usage should ameliorate the issue numerically.
Posted 1 month ago
 After 8 years, I am again struggling with my matrix B to get it still faster... The problem is that the code involves big numbers generated by Binomial[t, i]. For instance, Binomial[t, t/2]=2.04815162699310^600 for t=2000. So this is well above double floating of 2^1023=8.9884710^307. In the code, there are also factorials t!, Anyone has an idea how such numbers can be handled to prevent the compiled code to generate errors and reverse back to the slow uncompiled code?
Posted 1 month ago
 Maybe post your question here. They will definitely answer your question and you will be satisfiedI'm not a good programmer who knows tricks to speed up the code and I won't help you in your situation.Regards M.I.
Posted 8 years ago
 In case anyone is interested: the main brake in the above code is the seq function. It creates useless redundancy because the list of conditions it generates can be integrated directly into Do-loops as follows: B[a_, m_, t_] := Module[{n, p, q, s, f, b}, n = 1 - m; p = m/a; q = m (1 - 1/a); s[i_, j_, K_] := Sum[Binomial[i, i - k] Binomial[t - i, j - i + k] (p q)^ k ((n + p) (n + q))^(K - k), {k, 0, K}]; f[i_, j_] := If[(t - i - j) >= 0, (n + p)^(t - i - j), (n + q)^-(t - i - j)]; Do[b[i, j] = p^(j - i) f[i, j] Binomial[t, i]/Binomial[t, j] s[i, j, i], {i, 1, Round[t/2] + 1}, {j, i, t - i}]; Do[b[i, j] = p^(j - i) f[i, j] Binomial[t, i]/Binomial[t, j] s[i, j, t - j], {j, Round[t/2], t - 1}, {i, t - j + 1, j}]; Do[b[i, t] = p^(t - i) f[i, t] Binomial[t, i], {i, 1, t}]; Do[b[j, i] = (q/p)^(j - i) b[i, j] Binomial[t, j]/Binomial[t, i], {i, 1, t}, {j, i + 1, t}]; Return[Table[b[i, j], {i, 1, t}, {j, 1, t}]]]; `This way it is possible to reduce the time for t=2000 from thousands of years to hours, still non-compiled...