Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

This is un updated version of a previous (now deleted) post of mine.

I need to marginalize a distribution. I was suggested to use the code in this thread. This works just fine, but the problem is that the construction of tab (see below) is very slow: I tried with a pTab[i] with 100 elements and it took 2 days on my dual core 2GHz laptop to build it, and my actual pTab[i] are 180000 elements long.

I guess you understand why I am worried.

Does anyone have any idea about how to speed it up?

Or a different code, maybe?

Here is the code I am using (only the slow part):

dat = (*Import data from dat.txt*)
mat = (*Import matrix from matrix.txt*)
nosys = Inverse[mat];
z = dat[[All, 1]];
mu = dat[[All, 2]];
H = ConstantArray[1, 580];
chi[{a_, b_, c_, 
   e_}] = (5/
      Log[10]*(Log[z] + 0.5*z*(3 - b) + 
        1/24*z^2*(21 + 9*b^2 - 2*b - 4*c)) - 
     mu).nosys.(5/
      Log[10]*(Log[z] + 0.5*z*(3 - b) + 
        1/24*z^2*(21 + 9*b^2 - 2*b - 4*c)) - 
     mu) - (H.nosys.(5/
         Log[10]*(Log[z] + 0.5*z*(3 - b) + 
           1/24*z^2*(21 + 9*b^2 - 2*b - 4*c)) - 
        mu))^2/(H.nosys.H) + ((a - 0.742)/0.024)^2 + ((e - 0.315)/
     0.017)^2;

ClearAll[pTab];

val = (*Import values needed for Ptab from ptab_values.txt*)

d = 4;(*number of parameters*)
Do[pTab[i] = val[[All, i + 1]], {i, d}];
tab = Module[{T, n}, 
   With[{F = chi}, 
    Replace[Hold[
       Array[Append[T, F[T]] &, Length[pTab[#]] & /@ Range[d]]] /. 
      T -> Table[Hold[pTab[n][[Slot[n]]]] /. n -> i, {i, d}], 
     Hold[x_] :> x, {0, Infinity}, Heads -> True]]];
Print["packing array"]
tab // Developer`PackedArrayQ;

The needed files are:

for ptab (these are the numbers needed for building pTab (skip the first column: that is needed for another work))

matrix (a 580 x 580 matrix)

for dat (only the first 2 columns are needed in this code)

Those are rapidshare links: I hope you allow it because the matrix is 580 X 580, and dat is 580 rows long: I can't write them in a post.

EDIT: I was thinking can Parallelize be of any help in this case?

share|improve this question
    
Do you think that Parallelize can be of any help? –  mattiav27 May 22 at 10:01
    
You could provide fake data like this, I suppose: nrec = 580; dat = RandomReal[1, {nrec, 2}]; mat = RandomReal[{-1, 1}, {nrec, nrec}]; etc. –  Michael E2 May 22 at 14:07
add comment

1 Answer 1

Perhaps this?

The main issue is that the chi function takes almost 3 seconds to evaluate once, which I did not investigate further. Compiling it speeds it up by factor of about 17000! Making the compiled function Listable and setting Parallelization -> True reduces the time of the whole program by another 30%.

Second, it's usually better to deal with arrays as arrays. The tab code below takes about half of the time as the code in the OP, but some of that is due to the listability and parallelization of chiCF0.

The whole program takes a little less than 0.03 seconds. It may still not be fast enough to switch from 580 elements to 180000, though.

ClearAll[pTab];
pTab = Transpose[val[[All, 2 ;; d + 1]]];

chiCF0 = Compile[{{arg, _Real, 1},
                  {z, _Real, 1}, {nosys, _Real, 2}, {H, _Real, 1}, {mu, _Real, 1}},
   Module[{a = arg[[1]], b = arg[[2]], c = arg[[3]], e = arg[[4]]},
    (5/Log[10]*(Log[z] + 0.5*z*(3 - b) + 
           1/24*z^2*(21 + 9*b^2 - 2*b - 4*c)) - 
        mu).nosys.(5/
          Log[10]*(Log[z] + 0.5*z*(3 - b) + 
           1/24*z^2*(21 + 9*b^2 - 2*b - 4*c)) - 
        mu) - (H.nosys.(5/
             Log[10]*(Log[z] + 0.5*z*(3 - b) + 
              1/24*z^2*(21 + 9*b^2 - 2*b - 4*c)) - 
           mu))^2/(H.nosys.H) + ((a - 0.742)/0.024)^2 + ((e - 0.315)/
        0.017)^2
    ],
   RuntimeAttributes -> {Listable}, Parallelization -> True
   ];
chiCF = chiCF0[#, z, nosys, H, mu] &;

tab = Module[{params = Tuples[pTab]},
  ArrayReshape[
   Transpose@Append[Transpose[params], chiCF[params]],
   Append[Table[d, {d}], d + 1]]
  ]
share|improve this answer
    
Thank you! I'll give it a try tonight. –  mattiav27 May 22 at 15:28
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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