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