Skip to main content

scoping - function fails to use local variable when it is called inside a Module


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



Answer



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}]
  ]

Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more