- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我正在尝试用 Python 实现一个简单的蒙特卡洛(我对此还很陌生)。来自 C 我可能走的是最错误的道路,因为我的代码对于我所要求的来说太慢了:对于 60 个 3d 粒子和周期性边界条件(PBC),我有一个潜在的硬球状(请参阅代码中的 V_pot(r)
) ),所以我定义了以下函数
import timeit
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from numpy import inf
#
L, kb, d, eps, DIM = 100, 1, 1, 1, 3
r_c, T = L/2, eps/(.5*kb)
beta = 1/(kb*T)
#
def dist(A, B):
d = A - B
d -= L*np.around(d/L)
return np.sqrt(np.sum(d**2))
#
def V_pot(r):
V = -eps*(d**6/r**6 - d**6/r_c**6)
if r > r_c:
V = 0
elif r < d:
V = inf
return V
#
def ener(config):
V_jk_val, j = 0, N
#
while (j > 0):
j -= 1
i = 0
while (i < j):
V_jk_val += V_pot(dist(config[j,:], config[i,:]))
i += 1
#
return V_jk_val
#
def acc(en_n, en_o):
d_en = en_n-en_o
if (d_en <= 0):
acc_val = 1
else:
acc_val = np.exp(-beta*(d_en))
return acc_val
#
然后,从配置开始(其中数组的每一行代表 3D 粒子的坐标)
config = np.array([[16.24155657, 57.41672173, 94.39565792],
[76.38121764, 55.88334066, 5.72255163],
[38.41393783, 58.09432145, 6.26448054],
[86.44286438, 61.37100899, 91.97737383],
[37.7315366 , 44.52697269, 23.86320444],
[ 0.59231801, 39.20183376, 89.63974115],
[38.00998141, 3.84363202, 52.74021401],
[99.53480756, 69.97688928, 21.43528924],
[49.62030291, 93.60889503, 15.73723259],
[54.49195524, 0.6431965 , 25.37401196],
[33.82527814, 25.37776021, 67.4320553 ],
[64.61952893, 46.8407798 , 4.93960443],
[60.47322732, 16.48140136, 33.26481306],
[19.71667792, 46.56999616, 35.61044526],
[ 5.33252557, 4.44393836, 60.55759256],
[44.95897856, 7.81728046, 10.26000715],
[86.5548395 , 49.74079452, 4.80480133],
[52.47965686, 42.831448 , 22.03890639],
[ 2.88752006, 59.84605062, 22.75760029],
[ 9.49231045, 42.08653603, 40.63380097],
[13.90093641, 74.40377984, 32.62917915],
[97.44839233, 90.47695772, 91.60794836],
[51.29501624, 27.03796277, 57.09525454],
[10.30180295, 21.977336 , 69.54173272],
[59.61327648, 14.29582325, 11.70942289],
[89.52722796, 26.87758644, 76.34934637],
[82.03736088, 78.5665713 , 23.23587395],
[79.77571695, 66.140968 , 53.6784269 ],
[82.86070472, 40.82189833, 51.48739072],
[99.05647523, 98.63386809, 6.33888993],
[31.02997123, 66.99709163, 95.88332332],
[97.71654767, 59.24793618, 5.20183793],
[ 6.79964473, 45.01258652, 48.69477807],
[93.34977049, 55.20537774, 82.35693526],
[17.35577815, 20.45936211, 29.27981422],
[55.51942207, 52.22875901, 3.6616131 ],
[61.45612224, 36.50170405, 62.89796773],
[23.55822368, 7.09069623, 37.38274914],
[39.57082799, 58.95457592, 48.0304924 ],
[93.94997617, 64.34383203, 77.63346308],
[17.47989107, 90.01113402, 81.00648645],
[86.79068539, 66.35768515, 56.64402907],
[98.71924121, 38.33749023, 73.4715132 ],
[ 0.42356139, 78.32172925, 15.19883322],
[77.75572529, 2.60088767, 56.4683935 ],
[49.76486142, 3.01800153, 93.48019286],
[42.54483899, 4.27174457, 4.38942325],
[66.75777178, 41.1220603 , 19.64484167],
[19.69520773, 41.09230171, 2.51986091],
[73.20493772, 73.16590392, 99.19174281],
[94.16756184, 72.77653334, 10.32128552],
[29.95281655, 27.58596604, 85.12791195],
[ 2.44803886, 32.82333962, 41.6654683 ],
[23.9665915 , 49.94906612, 37.42701059],
[30.40282934, 39.63854309, 47.16572743],
[56.04809276, 30.19705527, 29.15729635],
[ 2.50566522, 70.37965564, 16.78016719],
[28.39713572, 4.04948368, 27.72615789],
[26.11873563, 41.49557167, 14.38703697],
[81.91731981, 12.10514972, 12.03083427]])
我使用以下代码进行了 5000 个时间步的模拟
N = 60
TIME_MC = 5000
DELTA_LIST = [d]
#d/6, d/3, d, 2*d, 3*d
np.random.seed(19680801)
en_mc_delta = np.zeros((TIME_MC, len(DELTA_LIST)))
start = timeit.default_timer()
config_tmp = config
#
for iD, Delta in enumerate(DELTA_LIST):
t=0
while (t < TIME_MC):
for k in range(N):
RND = np.random.rand()
config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)
en_o, en_n = ener(config), ener(config_tmp)
ACC = acc(en_n, en_o)
if (RND < ACC):
config[k,:] = config_tmp[k,:]
en_o = en_n
en_mc_delta[t][iD] = en_o
t += 1
stop = timeit.default_timer()
print('Time: ', stop-start)
遵循 Metropolis 算法的规则来接受使用 config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)
提取的提议的移动。
我做了一些尝试来检查代码卡在哪里,我发现函数 ener
(也是因为函数 dist
)非常慢:它需要像 ~0.02s
这样的东西来计算配置的能量,这意味着周围的东西~6000s
运行完整的模拟(60 个粒子,5000 个建议的移动)。
其外部只是计算 Delta
不同值的结果。
使用 TIME_MC=60
运行此代码可以让您了解此代码 ( ~218s
) 有多慢,如果用 C 实现,只需几秒钟。我读过一些有关如何加速 Python 代码的其他问题,但我不能了解如何在此处执行此操作。
我现在几乎可以肯定问题出在函数 dist
中,因为仅计算两个 3D 向量之间的 PBC 距离就需要 ~0.0012s
左右的时间,当您计算 5000*60 次时,这会带来疯狂的长时间。
最佳答案
请注意,这是对原始问题的评论的部分答案。
下面是一个示例,说明“展开”numpy 函数在替换为更直接的距离计算时如何提高性能。请注意,这尚未被验证为等效,特别是在舍入方面。我认为这个原则仍然适用。
import random
import time
import numpy as np
L = 100
inv_L = 0.01
vec_length = 10
repetitions = 100000
def dist_np(A, B):
d = A - B
d -= L*np.around(d/L)
return np.sqrt(np.sum(d**2))
def dist_direct(A, B):
sum = 0
for i in range(0, len(A)):
diff = (A[0,i] - B[0,i])
diff -= L * int(diff * inv_L)
sum += diff * diff
return np.sqrt(sum)
vec1 = np.zeros((1,vec_length))
vec2 = np.zeros((1,vec_length))
for i in range(0, vec_length):
vec1[0,i] = random.random()
vec2[0,i] = random.random()
print("with numpy method:")
start = time.time()
for i in range(0, repetitions):
dist_np(vec1, vec2)
print("done in {}".format(time.time() - start))
print("with direct method:")
start = time.time()
for i in range(0, repetitions):
dist_direct(vec1, vec2)
print("done in {}".format(time.time() - start))
输出:
with numpy method:
done in 6.332799911499023
with direct method:
done in 1.0938000679016113
研究平均向量长度和重复次数,看看最佳位置在哪里。我预计当改变这些元参数时,性能增益并不是恒定的。
关于python - Python 中的 Monte Carlo 和 Metropolis 算法非常慢,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/55901584/
滑动窗口限流 滑动窗口限流是一种常用的限流算法,通过维护一个固定大小的窗口,在单位时间内允许通过的请求次数不超过设定的阈值。具体来说,滑动窗口限流算法通常包括以下几个步骤: 初始化:设置窗口
表达式求值:一个只有+,-,*,/的表达式,没有括号 一种神奇的做法:使用数组存储数字和运算符,先把优先级别高的乘法和除法计算出来,再计算加法和减法 int GetVal(string s){
【算法】前缀和 题目 先来看一道题目:(前缀和模板题) 已知一个数组A[],现在想要求出其中一些数字的和。 输入格式: 先是整数N,M,表示一共有N个数字,有M组询问 接下来有N个数,表示A[1]..
1.前序遍历 根-左-右的顺序遍历,可以使用递归 void preOrder(Node *u){ if(u==NULL)return; printf("%d ",u->val);
先看题目 物品不能分隔,必须全部取走或者留下,因此称为01背包 (只有不取和取两种状态) 看第一个样例 我们需要把4个物品装入一个容量为10的背包 我们可以简化问题,从小到大入手分析 weightva
我最近在一次采访中遇到了这个问题: 给出以下矩阵: [[ R R R R R R], [ R B B B R R], [ B R R R B B], [ R B R R R R]] 找出是否有任
我正在尝试通过 C++ 算法从我的 outlook 帐户发送一封电子邮件,该帐户已经打开并记录,但真的不知道从哪里开始(对于 outlook-c++ 集成),谷歌也没有帮我这么多。任何提示将不胜感激。
我发现自己像这样编写了一个手工制作的 while 循环: std::list foo; // In my case, map, but list is simpler auto currentPoin
我有用于检测正方形的 opencv 代码。现在我想在检测正方形后,代码运行另一个命令。 代码如下: #include "cv.h" #include "cxcore.h" #include "high
我正在尝试模拟一个 matlab 函数“imfill”来填充二进制图像(1 和 0 的二维矩阵)。 我想在矩阵中指定一个起点,并像 imfill 的 4 连接版本那样进行洪水填充。 这是否已经存在于
我正在阅读 Robert Sedgewick 的《C++ 算法》。 Basic recurrences section it was mentioned as 这种循环出现在循环输入以消除一个项目的递
我正在思考如何在我的日历中生成代表任务的数据结构(仅供我个人使用)。我有来自 DBMS 的按日期排序的任务记录,如下所示: 买牛奶(18.1.2013) 任务日期 (2013-01-15) 任务标签(
输入一个未排序的整数数组A[1..n]只有 O(d) :(d int) 计算每个元素在单次迭代中出现在列表中的次数。 map 是balanced Binary Search Tree基于确保 O(nl
我遇到了一个问题,但我仍然不知道如何解决。我想出了如何用蛮力的方式来做到这一点,但是当有成千上万的元素时它就不起作用了。 Problem: Say you are given the followin
我有一个列表列表。 L1= [[...][...][.......].......]如果我在展平列表后获取所有元素并从中提取唯一值,那么我会得到一个列表 L2。我有另一个列表 L3,它是 L2 的某个
我们得到二维矩阵数组(假设长度为 i 和宽度为 j)和整数 k我们必须找到包含这个或更大总和的最小矩形的大小F.e k=7 4 1 1 1 1 1 4 4 Anwser是2,因为4+4=8 >= 7,
我实行 3 类倒制,每周换类。顺序为早类 (m)、晚类 (n) 和下午类 (a)。我固定的订单,即它永远不会改变,即使那个星期不工作也是如此。 我创建了一个函数来获取 ISO 周数。当我给它一个日期时
假设我们有一个输入,它是一个元素列表: {a, b, c, d, e, f} 还有不同的集合,可能包含这些元素的任意组合,也可能包含不在输入列表中的其他元素: A:{e,f} B:{d,f,a} C:
我有一个子集算法,可以找到给定集合的所有子集。原始集合的问题在于它是一个不断增长的集合,如果向其中添加元素,我需要再次重新计算它的子集。 有没有一种方法可以优化子集算法,该算法可以从最后一个计算点重新
我有一个包含 100 万个符号及其预期频率的表格。 我想通过为每个符号分配一个唯一(且前缀唯一)的可变长度位串来压缩这些符号的序列,然后将它们连接在一起以表示序列。 我想分配这些位串,以使编码序列的预
我是一名优秀的程序员,十分优秀!