1
$\begingroup$

I have these two functions fun and microstep.Fun makes use of a Module construct within which I define the Array I need to store the values of magnetization for different temperatures (each case stored in a different row). microstep is the function that store the data at the correct position at each step of the Monte Carlo algorithm. The monte Carlo procedure doesn't matter really much now, what bothers me is that when I define the magnetization array inside fun, the function doesn't work properly:

   fun [numbofsets_, nsteps_] := Module [{confinit, magnetization, index},

  index = 
   MapIndexed[   { #2[[1]], #   } &, numbofsets   ];     (* {index,
  temp} tuple*) 
  confinit = RandomChoice[{-1, 1}, {10, 10}];             (*  
  initial random matrix *)
  magnetization = 
   ConstantArray[ 0, {Length@numbofsets, nsteps}];


  Table[   

   NestList[
     microstep[ ##[[1]], ##[[2]], ##[[3]], ##[[4]]       ]         &  \
  , {   index[[i, 1]], index[[i, 2]]   ,   confinit, 2  }          ,
      nsteps];

   , {i, 1, Length@numbofsets}];

  ]

and

microstep[tindex_, temp_, matrix_, mcindex_] :=
 Module[{ tempmatrix = matrix, dimx, dimy, x , y , it = 1/temp, down, 
   up, left, right, spinsum , randnum, bool = False    , J = 1       },

  (* generic Metropolis Alghoritm *)
  dimx = Dimensions[matrix][[1]];
  dimy = Dimensions[matrix][[2]];
  x = RandomInteger[{1, dimx}];
  y = RandomInteger[{1, dimy}];
  randnum = RandomReal[];
  spinsum = 
   Plus[Compile`GetElement[matrix, Mod[x + 1, dimx, 1], y], 
    Compile`GetElement[matrix, Mod[x - 1, dimx, 1], y], 
    Compile`GetElement[matrix, x, Mod[y - 1, dimy, 1]], 
    Compile`GetElement[matrix, x, Mod[y + 1, dimy, 1]]];

  If[2*J *spinsum*tempmatrix[[x, y]] < 0 \[Or] 
    randnum < E^(- it*2*J*tempmatrix[[x, y]]*spinsum)
   ,
   tempmatrix[[x, y]] = -Compile`GetElement[matrix, x, y]; bool = True
   ];
  (*  tricky part  starts here *)

  If[bool,
   magnetization[[tindex, mcindex]] = 
     Abs[(magnetization[[tindex, mcindex - 1]] + 
        2 *tempmatrix[[x, y]])] ;                                      
   , 
   magnetization[[tindex, mcindex]] = 
     magnetization[[tindex, mcindex - 1]];

   ];
  {tindex, temp, tempmatrix, mcindex + 1}

  ]

now if i run

fun [{2, 3, 4}, 10]

i get

"Part specification magnetization[[1,1]] is longer than depth of \ object"

Meanwhile If I declare the magnetization array outside the Module, the function works properly giving me the correctly stored values, but it forces me to use global variables :

magnetization = ConstantArray[0, {3, 11}];
fun [{2, 3, 4}, 10];
magnetization

{ {0, 2, 0, 2, 4, 6, 8, 6, 8, 6, 6}, {0, 2, 2, 0, 2, 2, 4, 4, 2, 0, 2}, {0, 2, 4, 6, 8, 10, 8, 6, 6, 4, 4} })

I think the problem rise up from the fact that module is a scoping construct but I thought that a function called inside it would see the local variable but it doesn't and I don't know how to solve the problem. In c-like languages pointers can be used, is there anything similar in mathematica? Also, as always, any suggestion is appreciated

$\endgroup$
  • $\begingroup$ Have you tried reducing the example so the answer could be more generic and to get more attention? $\endgroup$ – Kuba May 16 '18 at 7:15
  • $\begingroup$ i guess i can simplify the step function a little bit. $\endgroup$ – Alucard May 16 '18 at 7:33
  • 1
    $\begingroup$ This is one of the cases where you would need the "dynamical scoping" of Block instead of the "lexical scoping" of Module. If you replace Module by Block in your definition of fun, everything should work fine. $\endgroup$ – mmeent May 16 '18 at 11:01
  • $\begingroup$ @mmeent uhm I will try it as soon as I reach home $\endgroup$ – Alucard May 16 '18 at 11:52
2
$\begingroup$

Does this work out for you?

Here I added magnetization as additional argument and gave microstep the attibute HoldAll to allow for call by reference.

SetAttributes[microstep, HoldAll];
microstep[tindex_, temp_, matrix_, mcindex_, magnetization_] := 
 Module[{tempmatrix = matrix, dimx, dimy, x, y, it = 1/temp, down, up,
    left, right, spinsum, randnum, bool = False, 
   J = 1},(*generic Metropolis Alghoritm*)

  dimx = Dimensions[matrix][[1]];
  dimy = Dimensions[matrix][[2]];
  x = RandomInteger[{1, dimx}];
  y = RandomInteger[{1, dimy}];
  randnum = RandomReal[];
  spinsum = Plus[
    Compile`GetElement[matrix, Mod[x + 1, dimx, 1], y],
    Compile`GetElement[matrix, Mod[x - 1, dimx, 1], y],
    Compile`GetElement[matrix, x, Mod[y - 1, dimy, 1]],
    Compile`GetElement[matrix, x, Mod[y + 1, dimy, 1]]
    ];
  If[2*J*spinsum*tempmatrix[[x, y]] < 0 ∨ 
    randnum < E^(-it*2*J*tempmatrix[[x, y]]*spinsum), 
   tempmatrix[[x, y]] = -Compile`GetElement[matrix, x, y]; 
   bool = True];
  (*tricky part starts here*)
  If[bool, magnetization[[tindex, mcindex]] = Abs[(magnetization[[tindex, mcindex - 1]] + 2*tempmatrix[[x, y]])];, 
   magnetization[[tindex, mcindex]] = magnetization[[tindex, mcindex - 1]];
   ];
  {tindex, temp, tempmatrix, mcindex + 1}]

There was also a second issue within fun: Apparently, the array magnetization was set up a bit too short, so I prolonged it by 1. Removing also some ;, the function fun executes without error and returns a result. Checking wether the result is correct is up to you.

fun[numbofsets_, nsteps_] := 
 Module[{confinit, index}, 
  index = MapIndexed[{#2[[1]], #} &, numbofsets];(*{index,temp} tuple*)

  confinit = RandomChoice[{-1, 1}, {10, 10}];(*initial random matrix*)

  magnetization = ConstantArray[0, {Length@numbofsets, nsteps + 1}];
  Table[
   NestList[
    microstep[##[[1]], ##[[2]], ##[[3]], ##[[4]], magnetization] &,
    {index[[i, 1]], index[[i, 2]], confinit, 2}, nsteps
    ]
   ,
   {i, 1, Length@numbofsets}]
  ]
$\endgroup$
  • $\begingroup$ I will let you know soon Henrik, the windows spring update hit my notebook with the strength of ten thousand suns :-( $\endgroup$ – Alucard May 16 '18 at 8:50
  • $\begingroup$ Whoa... Good luck, then! $\endgroup$ – Henrik Schumacher May 16 '18 at 9:11
  • $\begingroup$ Though I still have to analyze the data, the array is properly filled. Thanks again henrik $\endgroup$ – Alucard May 16 '18 at 11:49
  • $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher May 16 '18 at 12:03
2
$\begingroup$

Just give your function microstep[tindex_, temp_, matrix_, mcindex_] an extra arg, so make it microstep[tindex_, temp_, matrix_, mcindex_, magnetization_] and feed it in as an argument.

By default microstep is looking for magnetization in the global scope, so you need to either make magnetization Global or make feed it in as an argument.

Additionally, if you want to rewrite values (e.g pass by reference instead of pass by value, which is what the above would do) you can pass in arguments in their Held form.

So your downvalue would now be microstep[tindex_, temp_, matrix_, mcindex_, Hold[magnetization_]], and obviously change fun to call microstep with the extra argument.

If this is not clear, I can update with more code examples when I get to my machine that has Mma on it.

Cheers!

Edit -- It seems that if you look @Heinriks answer, he implemented what I was thinking.

$\endgroup$
  • $\begingroup$ i tried with hold and reporting the magnetization as argument, fun keeps giving me back an empty array. the error is ` set {{},....} in the part assignment is not a symbol`. Btw i had completely forgotten about Hold, damn few months of not using mathematica and i have already forgot many useful functions $\endgroup$ – Alucard May 16 '18 at 4:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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