| 查看: 297 | 回复: 1 | |||
| 当前主题已经存档。 | |||
lovein金虫 (小有名气)
|
[交流]
【求助】麻烦帮忙看下一个“随机游走”的程序，为什么运行时老出错！！！
|
||
|
整个程序如下： (*----------------------------------------------------------------------*) BeginPackage["RW2D`"] RW2D::usage = "RW2D[maxstep, ktry, range, startarea, datadirectory] reads the two - dimensional digital image file (text file) and exports the mean - square displacement of random walkers as a text file, 'msdlist2d.out'. 'maxstep' is the final value of the time step of the random walker. 'ktry' is the number of average (i.e., number of walkers). 'range' is the lower and upper limit of the voxel intensity of the pore space. 'startarea' is the dimension of a square as the percentage to the system size. The random walker starts to diffuse from inside the square which is located at the center of the image. 'datadirectory' is the directory or folder name where a single two-dimensional digital image file (text file) is located. " ShowTrajectory::usage = "ShowTrajectory[kk] is for the graphic display of the random walk trajectory of the 'kk' - th trial. " Begin["`Private`"] (* Reading the two - dimensional data in 'datadirectory' *) readdata2d[range_List, datadirectory_] := Module[{dataname, datasize, x, m2d, rmin = range[[1]], rmax = \ range[[2]]}, dataname = First[FileNames["*", {datadirectory}]]; datasize = Round[Sqrt[Length[ReadList[ToString[dataname], Real]]]]; m2d = Partition[ Map[If[rmin (* random walk *) (* 'xylist' is the trajectory of 'maxstep' - step from 'o2d' *) tra2d[maxstep_, o2d_List, m2d_List] := Module[{n2d = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}}, now = o2d, xylist = {o2d}, outoffield = 0, xmax = ymax = Length[m2d], step, new, x, y}, ran2d := Random[Integer, {1, 4}]; For[i = 1, i (* random walk of 'ktry' - walkers *) (* 'tra2dktry' is list of the trajectory of 'ktry' - walkers *) rw2d[maxstep_, ktry_, startarea_, m2d_List] := Module[{xdatasize = ydatasize = Length[m2d], xylist, outoffield, tra2dktry = {}, failure = 0}, rano2d[lx_, ly_, sa_] := {Random[Integer, {(1 - sa/100)/2 lx, (1 + sa/100)/2lx}], Random[Integer, {(1 - sa/100)/2 ly, (1 + sa/100)/2ly}]}; o2d := (For[i = 1; {x0, y0} = rano2d[xdatasize, ydatasize, startarea], m2d[[x0, y0]] != 0, i++, {x0, y0} = rano2d[xdatasize, ydatasize, startarea]]; {x0, y0}); Do[{xylist, outoffield} = tra2d[maxstep, o2d, m2d]; AppendTo[tra2dktry, xylist]; AddTo[failure, outoffield], {i, 1, \ ktry}]; Return[{tra2dktry, failure}]; ]; (* calculation of the mean - square displacement *) msd2d[data_List] := Module[{ktry = Dimensions[data][[1]], maxstep = Dimensions[data][[2]] - \ 1, d1, d2, d3, dd, ii}, d1 = Transpose[data, {2, 1}]; d2 = Map[(# - First[d1]) &, d1]; d3 = Transpose[Map[(#^2) &, d2]]; dd = Apply[Plus, Transpose[Apply[Plus, d3]/ktry]]; ii = Range[0, maxstep]; msdlist2d = Transpose[{ii, dd}]; Return[msdlist2d]; ]; RW2D[maxstep_, ktry_, range_List, startarea_, datadirectory_] := Module[{failure, msdlist2d}, (* reading the two - dimensional data in 'datadirectory' *) m2d = readdata2d[range, datadirectory]; (* execution of the random walk *) {tra2dktry, failure} = rw2d[maxstep, ktry, startarea, m2d]; Print["Failure = ", failure, " / Success = ", ktry]; (* calculation of the mean - square displacement of 'ktry' walkers *) msdlist2d = msd2d[tra2dktry]; (* file output of the mean - square displacement *) (* The file name is 'msdlist2d.out' *) Export["msdlist2d.out", N[msdlist2d], "Table"]; (* graphic output of the mean - square displacement *) ListPlot[msdlist2d, PlotJoined -> True, AxesLabel -> {"time \n step", "mean-square \n displacement"}]; ]; (* graphic output of the trajectory *) (* Show the trajectory of a single walker as an example. *) (* RGBColor[] of 0, 1, 2 are colors of the pore, obstacle, and trajectory, respectively. You can change the color as you like. *) gtra2d[tra2d_List, m2d_List] := Module[{color0 = RGBColor[0.9, 0.9, 0.9], color1 = RGBColor[0.5, 0.5, 0.5], color2 = RGBColor[0.8, 0.8, 0.8], plmin = 1 - 2, plmax = Length[m2d] + 1}, tra2dshift = Map[(# - {0.5, 0.5}) &, tra2d]; tra2dm2d = (aa2 = m2d; Map[(aa2[[#[[1]], #[[2]]]] = 2) &, tra2d]; aa2); Show[ Graphics[ RasterArray[tra2dm2d /. {0 -> color0, 1 -> color1, 2 -> color2}]], Graphics[{Line[Map[Reverse, tra2dshift]], AbsolutePointSize[5.0], Point[Reverse[First[tra2dshift]]], Point[Reverse[Last[tra2dshift]]]}], PlotRange -> {{plmin, plmax}, {plmin, plmax}}, AspectRatio -> 1]; ]; (* Only for the 'kk' - th trial, the trajectory of the random walk is shown graphically *) ShowTrajectory[kk_] := gtra2d[tra2dktry[[kk]], m2d]; End[]; EndPackage[]; (*----------------------------------------------------------------------*) (* example of usage *) (*----------------------------------------------------------------------*) (* As noted below, first of all, you put the command '<< RW2D.m'. *) (* Then you can put the command set, 'RW2D[500, 200, {0, 34000}, 10, "DATAsingle"];' and 'ShowTrajectory[2];'. \ *) (* 'RW2D[maxstep, ktry, range, startarea, datadirectory]' is the command for the performance of the random walk and the calculation of mean - square displacement. *) (* The command 'ShowTrajectory[kk]' is for the graphic display of the trajectory of the 'kk' - th trial. *) (* In the example below, a trajectory of the second trial will be shown graphically. *) (* In[1] := << RW2D.m *) (* In[2] := RW2D[500, 200, {0, 34000}, 10, "DATAsingle"]; ShowTrajectory[2]; \ *) 不知道为啥，运行的时候老是出错，我又看不大懂这程序，请求帮忙 [ Last edited by lovein on 2009-10-26 at 13:35 ] |
回复此楼» 猜你喜欢
国社科又开始会评了，不知道这次命运如何
已经有13人回复
-
麻烦专家们看看评委们的意见（F口面上）
已经有13人回复
-
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有4人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有6人回复
-
为什么资助数各大高校都创新高，自己申请怎么就这么难
已经有15人回复
-
科研人应该花精力去思考如何解决问题，而不是去凝练问题
已经有14人回复
-
要骂人了，新模版改版就是要淡化问题凝练这种虚的东西，结果有个评委还在说凝练得不够
已经有17人回复
-
基金系统什么内容也没有
已经有10人回复
-
学科评审组评审是指会评吗？
已经有5人回复
-
面上函评意见出来了，像什么等级？
已经有17人回复
lovein
金虫 (小有名气)
- 应助: 1 (幼儿园)
- 金币: 1419.6
- 散金: 4
- 红花: 1
- 帖子: 254
- 在线: 38.3小时
- 虫号: 534931
- 注册: 2008-03-28
- 性别: MM
- 专业: 其他无机非金属材料
2楼2009-10-26 13:39:16
30